# A simple artificial damping method for total Lagrangian smoothed particle
hydrodynamics

Chi Zhang<EMAIL_ADDRESS>Yujie Zhu<EMAIL_ADDRESS>Yongchuan Yu
<EMAIL_ADDRESS>Massoud Rezavand<EMAIL_ADDRESS>Xiangyu Hu
<EMAIL_ADDRESS>Department of Mechanical Engineering, Technical University
of Munich, 85748 Garching, Germany  Department of Aerospace and Geodesy,
Technical University of Munich, 82024 Taufkirchen, Germany

###### Abstract

In this paper, we present a simple artificial damping method to enhance the
robustness of total Lagrangian smoothed particle hydrodynamics (TL-SPH).
Specifically, an artificial damping stress based on the Kelvin-Voigt type
damper with a scaling factor imitating a von Neumann-Richtmyer type artificial
viscosity is introduced in the constitutive equation to alleviate the spurious
oscillation in the vicinity of the sharp spatial gradients. After validating
the robustness and accuracy of the present method with a set of benchmark
tests with very challenging cases, we demonstrate its potentials in the field
of bio-mechanics by simulating the deformation of complex stent structures.

###### keywords:

Total Lagrangian formulation, Smoothed particle hydrodynamics , Solid dynamics
, Kelvin-Voigt damper

††journal: J. X. X. X

##  1 Introduction

Numerical simulation of large-strain solid dynamics problems, where flexible
structures experience large deformations, plays a key role in a vast range of
engineering problems in the aerospace, automotive, manufacturing and
biomaterial industries. Besides the traditional mesh-based finite element
methods (FEM), smooth particle hydrodynamics (SPH), which is a mesh-free
method and originally designed for fluid dynamics [1, 2], has also been
adopted for such problems [3, 4], and received increasing attention in the
past decades. Generally, there are two types of SPH formulation for solid
dynamics: one is updated Lagrangian SPH (UL-SPH) formulation, in which the
current configuration is used as the reference [5, 6, 4, 7, 8, 9, 10, 11]; the
other is total Lagrangian SPH (TL-SPH) formulation [12], in which the initial
configuration is used as the reference. Compared with UL-SPH, TL-SPH shows
promising potential for solid dynamics due to its attractive advantages in
being free from tensile instability and ensuring $1$st-order consistency when
computing deformation gradient by introducing a kernel correction matrix.
Since its inception, TL-SPH has been applied for the problems of necking and
fracture in electromagnetically driven rings [13], thermomechanical
deformations [14], fluid-structure interaction (FSI) [15, 16, 17, 18] and bio-
mechanics [19, 20], among many others. It is known that without appropriate
stabilization the original TL-SPH exhibits spurious fluctuations in the
vicinity of sharp spatial gradients (as will also be shown in Section 4.1).
This deficiency may result in numerical instability and lead to wrongly
predicted deformation for problems involving large strain. To rectify this
deficiency, Lee et al. [21] proposed a Jameson-Schmidt-Turkel (JST) based
method, which shows good performance of eliminating spurious pressure
oscillations in nearly incompressible solids. In a recent work, Lee et al.
[22] further proposed a new stabilization method by introducing a
characteristic-based Riemann solver in conjunction with a linear
reconstruction procedure. This method also shows good performance in the
simulation of nearly and truly incompressible explicit fast solid dynamics
with large deformations.

The main objective of this paper is to present a simple and effective
artificial damping method for TL-SPH to enhance its numerical stability. In
particular, an artificial damping based on Kelvin–Voigt type damper is
introduced in the constitutive equation to alleviate the spurious oscillation
in the vicinity of the sharp spatial gradients. By imitating a von Neumann-
Richtmyer type artificial viscosity, a scaling factor is introduced to control
the damping force. Compared with the works of Refs. [21, 22], the present
method is very simple and easy to be implemented into the original TL-SPH
formulation in a straightforward way. A set of benchmark tests with very
challenging cases are investigated to validate the robustness and accuracy of
the present method. Furthermore, its versatility is demonstrated by modeling
the deformation of complex stent structures. Also, all the codes and data-sets
accompanying this work are available from the open-source library of SPHinXsys
[19] on GitHub at https://github.com/Xiangyu-Hu/SPHinXsys. The remainder of
this paper is arranged as follows: Section 2 briefly summarizes the governing
equations and original TL-SPH formulation. Then, TL-SPH-KV is detailed in
Section 3 and numerical validations and applications are presented and
discussed in Section 4. Finally, brief concluding remarks are given in Section
5.

##  2 Preliminary

The kinematics of the finite-deformation solid dynamics can be characterized
by introducing a deformation map $\varphi$, which maps a material point
$\mathbf{\mathbf{r}^{0}}$ from the initial reference configuration
$\Omega^{0}\subset\mathbb{R}^{d}$ to the point
$\mathbf{r}=\mathbf{\varphi}\left(\mathbf{r}^{0},t\right)$ in the deformed
configuration $\Omega=\mathbf{\varphi}\left(\Omega^{0}\right)$. Here, we
denote the superscript $\left({\bullet}\right)^{0}$ as the quantities in the
initial reference configuration. Then, the deformation tensor $\mathbb{F}$ can
be defined as

$\mathbb{F}=\nabla^{0}{\varphi}=\frac{\partial\varphi}{\partial\mathbf{r}^{0}}=\frac{\partial\mathbf{r}}{\partial\mathbf{r}^{0}},$
(1)

where the derivative is evaluated with respect to the initial reference
configuration. Note that the deformation tensor $\mathbb{F}$ can also be
calculated from the point displacement $\mathbf{u}=\mathbf{r}-\mathbf{r}^{0}$
through

$\mathbb{F}=\nabla^{0}{\mathbf{u}}+\mathbb{I},$ (2)

where $\mathbb{I}$ represents the unit matrix.

###  2.1 Governing equation for solid dynamics

In the total Lagrangian framework, the momentum conservation equation can be
expressed as

$\rho^{0}\frac{\text{d}\mathbf{v}}{\text{d}t}=\nabla^{0}\cdot\mathbb{P}^{T}+\rho^{0}\mathbf{g},$
(3)

where $\rho^{0}$ is the density in the reference configuration, $\mathbf{v}$
the velocity and $\mathbb{P}$ the first Piola-Kirchhoff stress tensor. For an
ideal elastic or Kirchhoff material, $\mathbb{P}$ is given by

$\mathbb{P}=\mathbb{F}\mathbb{S},$ (4)

where $\mathbb{S}$ represents the second Piola-Kirchhoff stress which is
evaluated via the constitutive equation relating $\mathbb{F}$ with the Green-
Lagrangian strain tensor $\mathbb{E}$ defined as

$\mathbb{E}=\frac{1}{2}\left(\mathbb{F}^{T}\mathbb{F}-\mathbb{I}\right).$ (5)

In particular, when the material is linear elastic and isotropic, the
constitutive equation is simply given by

$\displaystyle\mathbb{S}$ $\displaystyle=$ $\displaystyle
K\mathop{\mathrm{tr}}\left(\mathbb{E}\right)\mathbb{I}+2G\left(\mathbb{E}-\frac{1}{3}\mathop{\mathrm{tr}}\left(\mathbb{E}\right)\mathbb{I}\right)$
(6) $\displaystyle=$
$\displaystyle\lambda\mathop{\mathrm{tr}}\left(\mathbb{E}\right)\mathbb{I}+2\mu\mathbb{E},$

where $\lambda$ and $\mu$ are Lam$\acute{e}$ parameters, $K=\lambda+(2\mu/3)$
the bulk modulus and $G=\mu$ the shear modulus. The relation between the two
modulus is given by

$E=2G\left(1+2\nu\right)=3K\left(1-2\nu\right),$ (7)

where $E$ denotes the Young’s modulus and $\nu$ the Poisson’s ratio. Note that
the sound speed of solid structure is defined as $c=\sqrt{K/\rho^{0}}$. In the
present work, a neo-Hookean material model defined by the strain-energy
density function

$\displaystyle W=\mu\mathop{\mathrm{tr}}\left(\mathbb{E}\right)-\mu\ln
J+\frac{\lambda}{2}(\ln J)^{2},$ (8)

is also applied for predicting the nonlinear stress-strain behavior of
materials undergoing large deformations. For neo-Hookean material, the second
Piola-Kirchhoff stress $\mathbb{S}$ can be derived as

$\mathbb{S}=\frac{\partial W}{\partial\mathbb{E}}.$ (9)

###  2.2 TL-SPH formulation

In TL-SPH, the kernel correction or normalization technique [23, 24, 4] has
demonstrated its effects to improve the accuracy and consistency of SPH
method. The correction matrix $\mathbb{B}^{0}$ is introduced as [12]

$\mathbb{B}^{0}_{i}=\left(-\sum_{j}V^{0}_{j}\mathbf{r}_{ij}^{0}\otimes\nabla^{0}_{i}W_{ij}\right)^{-1},$
(10)

where

$\nabla^{0}_{i}W_{ij}=\frac{\partial
W\left(|\mathbf{r}^{0}_{ij}|,h\right)}{\partial|\mathbf{r}^{0}_{ij}|}\mathbf{e}^{0}_{ij},$
(11)

denoting the gradient of the kernel function evaluated at the initial
reference configuration. Note that the correction matrix is computed in the
initial configuration and, therefore, it is calculated only once before the
simulation. Then, the momentum conservation equation, Eq.(3), can be
discretized as

$\frac{\text{d}\mathbf{v}_{i}}{\text{d}t}=\frac{2}{m_{i}}\sum_{j}V^{0}_{i}V^{0}_{j}\tilde{\mathbb{P}}_{ij}\nabla^{0}_{i}W_{ij}+\mathbf{g},$
(12)

where the inter-particle averaged first Piola-Kirchhoff stress
$\tilde{\mathbb{P}}$ is defined as

$\tilde{\mathbb{P}}_{ij}=\frac{1}{2}\left(\mathbb{P}_{i}\mathbb{B}^{0}_{i}+\mathbb{P}_{j}\mathbb{B}^{0}_{j}\right).$
(13)

Here, the first Piola-Kirchhoff stress tensor is computed with the
constitutive law where the deformation tensor $\mathbb{F}$ is updated by the
change rate evaluated through

$\frac{\text{d}\mathbb{F}_{i}}{\text{d}t}=\left(-\sum_{j}V^{0}_{j}\mathbf{v}_{ij}\otimes\nabla^{0}_{i}W_{ij}\right)\mathbb{B}^{0}_{i}.$
(14)

##  3 Kelvin-Voigt type damping

An elastic solid undergoing large strains can be modeled with the mechanical
components of springs and dashpots. The former represents the restorative
force component and the later denotes the damping component. The Kelvin-Voigt
(KV) model can be represented by a purely viscous damper and purely elastic
spring connected in parallel, and the total stress is decomposed into two
parts

$\sigma_{total}=\sigma_{S}+\sigma_{D}.$ (15)

Here, $\sigma_{total}$ is the total stress, $\sigma_{S}$ the elastic stress
and $\sigma_{D}$ represent the damper stress as

$\sigma_{D}=\eta\frac{\text{d}\epsilon(t)}{\text{d}t}$ (16)

where $\text{d}\epsilon(t)/\text{d}t$ is the strain rate and $\eta$ the
physical viscosity. Applying the KV model to TL-SPH formulation, the second
Piola-Kirchhoff stress $\mathbb{S}$ can be rewritten as

$\mathbb{S}=\mathbb{S}_{S}+\mathbb{S}_{D},$ (17)

where $\mathbb{S}_{S}$ is given by the constitutive equation of Eq. (6) or Eq.
(8), and the damper $\mathbb{S}_{D}$ is defined as

$\mathbb{S}_{D}=\pi\frac{\text{d}\epsilon(t)}{\text{d}t}=\frac{\pi}{2}\left[\left(\frac{\text{d}\mathbb{F}}{\text{d}t}\right)^{T}\mathbb{F}+\mathbb{F}^{T}\left(\frac{\text{d}\mathbb{F}}{\text{d}t}\right)\right],$
(18)

where $\pi$ represents an artificial viscosity. As the main objective of
introducing the KV-type damper is to enhance the robustness of original TL-
SPH, we introduce a von Neumann-Richtmyer type scaling factor with the speed
of sound $c$ for Eq. (18) as

$\pi=\alpha\rho ch.$ (19)

Here, $\alpha$ is a suitable and constant parameter and $h$ the smoothing
length. Note that a similar parameter $\pi$ is also widely used in the
artificial viscosity for Eulerain shock-capturing schemes in modeling
compressible flow. We suggest $\alpha=0.5$ according to numerical experiments
and use it throughout this paper.

##  4 Numerical examples

In this section, we first study three benchmark tests, where the structures
may experience large deformation, to validate the robustness and performance
of the proposed method (denoted as TL-SPH-KV). We also compare numerical
results with those obtained by the original TL-SPH in which no damping stress
is applied (denoted as “TL-SPH”). Having the validation studies presented, we
then demonstrate the versatility of the method for applications in bio-
mechanical system, i.e. stent structures. In all the following examples, the
$5th$-order Wendland smoothing kernel function with an smoothing length of
$h=1.15dp$ is employed, where $dp$ represents the initial inter-particle
spacing. The position-based Verlet scheme, which is a two-step explicit
algorithm proposed in the work of Zhang et al. [18], is used for time
integration. The CFL number applied here is $CFL=0.6$, which is twice as that
used in the work of Lee et al . [22].

###  4.1 Wave propagation in a cable

Figure 1: Wave propagation in a cable: Initial configuration.

In the first benchmark test, we investigate a simple elastic wave propagation
in an elastic cable where an analytical solution is available for quantitative
comparison. This problem is also studied by Refs. [25, 26, 27, 21] by using
either updated or total Lagrangian SPH method. Following Refs. [25, 26, 27,
21], we consider a rod of dimensions $10\times 0.2\times 0.2~{}(m)$ with the
left end fixed and the right end free as shown in Figure 1. The wave
propagation is initialized by imposing a velocity $v=5~{}\text{m}/\text{s}$
along the length direction on the right quarter of the rod. We consider a
linear elastic material of density $\rho=8000~{}\text{kg}/\text{m}^{3}$,
Young’s modulus $E=200~{}\text{GPa}$ and Poisson’s ratio $\nu=0.0$.

Figure 2 shows the time histories of velocity and displacement in the length
direction at the right tip end of the rod and the comparison with analytical
solution. Similar to the reports in Ref. [27], TL-SPH exhibits excessive
oscillation, and similar overshoots in the velocity and displacement profiles,
respectively. As expected, TL-SPH-KV predicts the correct velocity and
displacement in both profiles. The convergence study of the present method is
also presented in Fig. 2 and shows clear convergence as the spatial resolution
increases. Note that, the comparison between the present results and those of
Ref. [21] is also reported in Fig.2. The main velocity and displacement
plateaus of these results are in good agreement except that small overshoots
are observed in those of Ref. [21].

Figure 2: Wave propagation in a cable: Velocity (upper panel) and displacement
(bottom panel) profiles at the right tip end. A linear elastic material with
density $\rho=8000~{}\text{kg}/\text{m}^{3}$, Young’s modulus
$E=200~{}\text{GPa}$ and Poisson’s ratio $\nu=0.0$ is applied.

###  4.2 Bending column

Figure 3: Bending column: Initial configuration. Note that point $S$ is
located at $(1,1,6)^{T}\text{m}$.

In this second benchmark test, we consider the bending deformation of a column
of span $\text{L}=6\text{m}$ and square cross section (height $h=1\text{m}$),
whose bottom end is clamped to the ground and its body is allowed to bend
freely by imposing an initial uniform velocity
$\mathbf{v}=(5\sqrt{3},5,0)^{T}\text{m}\cdot\text{s}^{-1}$ as shown in Fig. 3.
The neo-Hookean material with density $\rho=1.1\times
10^{3}~{}\text{kg}\cdot\text{m}^{-3}$, Young’s modulus $E=1.7\times
10^{7}~{}\text{Pa}$ and Poisson’s ratio $\nu=0.45$ is applied. This test is
also investigated by Aguirre et al.[28] and we set their results as a
reference for a rigorous comparison.

Figure 4 shows the deformed configuration colored with von Mises stress
contours for TL-SPH and TL-SPH-KV. Obviously, oscillations in the von Mieses
stress field is obtained by TL-SPH due to insufficient stabilization while
these oscillations are suppressed by TL-SPH-KV. Similar to the last test, TL-
SPH produces noisy oscillations in the velocity profile which are eliminated
by the present method as shown in Fig. 5. Note that, Fig. 5 also gives the
time history of the vertical displacement of point $S$ and its comparison with
that of Ref. [28]. Compared with the results reported in latter, good
agreements in the deformation are observed. Also note that a 2nd-order
convergence of the solution is achieved by the present method with increasing
spatial resolution, even though no linear reconstruction procedure is applied
as in Ref. [22]. Compared with the present method, the TL-SPH shows overshoots
in the displacement profile which is similar to fluctuations produced in the
velocity field.

Figure 4: Bending column: Deformed configuration at different time stances for
TL-SPH-KV and TL-SPH. The neo-Hookean material with density $\rho=1.1\times
10^{3}~{}\text{kg}\cdot\text{m}^{-3}$, Young’s modulus $E=1.7\times
10^{7}~{}\text{Pa}$ and Poisson’s ratio $\nu=0.45$ is applied.

Figure 5: Bending column: Velocity (upper panel) and displacement (bottom
panel) profiles at the right tip end. The displacement profile is compared
with data in Ref. [28] and convergence study is also presented. The neo-
Hookean material with density $\rho=1.1\times
10^{3}~{}\text{kg}\cdot\text{m}^{-3}$, Young’s modulus $E=1.7\times
10^{7}~{}\text{Pa}$ and Poisson’s ratio $\nu=0.45$ is applied.

###  4.3 Twisting column

Figure 6: Twisting column: Initial configuration.

Following Refs.[21, 22, 29], the example of bending column reported in Section
4.2 can be extended to a more challenging test to assess the robustness of TL-
SPH-KV for predicting the extremely highly nonlinear deformations. As shown in
Fig. 6, the twisting column is also clamped on its bottom face and the body is
initialized with a sinusoidal rotational velocity field relative to the origin
given by

$\mathbf{v}\left(\mathbf{r}^{0}\right)=\omega\times\mathbf{r}^{0},\omega=\left[0.0,\Omega_{0}\sin\left(\frac{\pi
y}{2L}\right),0.0\right]^{T}(\text{rad}\cdot\text{s}^{-1}).$ (20)

The column material is modeled as nearly incompressible by using a neo-Hookean
constitutive model with density $\rho_{0}=1100~{}\text{Kg}/\text{m}^{3}$,
Young’s modulus $E=0.017~{}\text{GPa}$ and Poisson’s ratio $\nu=0.4995$.

Figure 7 shows the deformed configurations colored by von Mises stress with
$\Omega_{0}=105\text{rad}\cdot\text{s}^{-1}$ obtained by TL-SPH and TL-SPH-KV.
Clearly, TL-SPH produces non-physical stress fluctuations due to insufficient
numerical stabilization and then fails to capture the correct deformation
pattern. On the contrary, TL-SPH-KV alleviates these discrepancies and
predicts more accurate deformation patterns as those in the literature [21,
22, 29] (see Figs 23 in Ref. [22]), demonstrating the robustness of the
proposed method. Note that the results of Ref. [22] (see their Fig. 23)
reports about slightly more twist may be due to the slight different neo-
Hookean material models. It is worth noting that TL-SPH-KV preserves the axial
rotation very well without introducing out-of-axis characteristics, as shown
in Fig. 8 monitoring the time history of the horizontal velocity components at
the point $\mathbf{r}=\left[0.0,6.0,0.0\right]^{T}$. As observed, TL-SPH
generates much larger out-of-axis fluctuations. Figure 9 shows the convergence
study with particle refinement. Both the deformation and von Mises stress
resolution obtained exhibit good convergence property for TL-SPH-KV.

To demonstrate the robustness of the present TL-SPH-KV, we consider more
challenging tests by increasing the initial rotational velocity to
$\Omega_{0}=200~{}\text{rad}/\text{s}$ and
$\Omega_{0}=300~{}\text{rad}/\text{s}$ (with Poisson’s ratio $\nu=0.49995$)
which induces extreme deformation and leads to high requirement for
robustness. It is observed from Fig. 10 (a) - (c) that the extremely large
deformation is well captured and, numerical convergence is achieved with the
increase of spatial resolution. Also, with the even higher initial angular
velocity, one more twist is obtained as shown in Fig. 10 (d), which gives
high-resolution results on the deformed configuration.

Figure 7: Twisting column: Comparison of deformed configuration plotted with
von Mises stress at serial time instance using the TL-SPH (upper panel) and
the TL-SPH-KV (bottom panel). Results obtained with a initial sinusoidal
rotational velocity $\Omega_{=}105~{}\text{rad}/\text{s}$. A neo-Hookean
material with density $\rho_{0}=1100~{}\text{Kg}/\text{m}^{3}$, Young’s
modulus $E=0.017~{}\text{GPa}$ and Poisson’s ratio $\nu=0.4995$ is applied.

Figure 8: Twisting column: Time history of the velocity at the middle point of
the free end using TL-SPH and TL-SPH-KV. Results obtained with a initial
sinusoidal rotational velocity $\Omega=105~{}\text{rad}/\text{s}$. A neo-
Hookean material with density $\rho_{0}=1100~{}\text{Kg}/\text{m}^{3}$,
Young’s modulus $E=0.017~{}\text{GPa}$ and Poisson’s ratio $\nu=0.4995$ is
applied. Figure 9: Twisting column: A sequence of particle refinement analysis
using TL-SPH-KV. Results obtained with a initial sinusoidal rotational
velocity $\Omega=105~{}\text{rad}/\text{s}$. A neo-Hookean material with
density $\rho_{0}=1100~{}\text{Kg}/\text{m}^{3}$, Young’s modulus
$E=0.017~{}\text{GPa}$ and Poisson’s ratio $\nu=0.4995$ is applied. Figure 10:
Twisting column: Deformed configuration plotted with von Mises stress for
large initial rotational velocity obtained by using the present TL-SPH-KV. (a)
- (c) particle refinement analysis for $\Omega_{0}=200~{}\text{rad}/\text{s}$
and (d) for $\Omega_{0}=300~{}\text{rad}/\text{s}$. A neo-Hookean material
with density $\rho_{0}=1100~{}\text{Kg}/\text{m}^{3}$, Young’s modulus
$E=0.017~{}\text{GPa}$ is applied. Note that Poisson’s ratio is set as
$\nu=0.4995$ for (a) - (c) and $\nu=0.49995$ for (d).

###  4.4 Stent structures

(a) PS-shaped stent

(b) C-shaped stent

Figure 11: Computer-aided Design (CAD) geometries of stent structures: Palmaz-
Schatz shaped (PS-shaped) stent and Cypher shaped (C-shaped) stent. The
corresponding CAD files in STL format can be downloaded from our code
repository or GrabCAD.

In the last example, complex flexible structures are investigated to
demonstrate the robustness and versatility of TL-SPH-KV. As shown in Fig. 11,
two different stent structures, viz. Palmaz-Schatz shaped (PS-shaped) and
Cypher shaped (C-shaped), are considered. Note that the present stent
structures are realistic cardiovascular stent and widely used in biomedical
applications. To study the deformation pattern of the stent structures, we
apply a velocity field with magnitude of $v_{0}=5.0\text{m}\cdot\text{s}^{-1}$
at the top and bottom of the structure, which is modeled as neo-Hookean
material with density $\rho_{0}=1100~{}\text{Kg}/\text{m}^{3}$, Young’s
modulus $E=0.017~{}\text{GPa}$ and Poisson’s ratio $\nu=0.45$ .

Figure 12 shows the overall deformation of the PS-shaped stent structure at
time $t=0.02\text{s}$ with von Mises stress contour. It can be observed that
the regions of high stress are concentrated at the four corners of the cells
rather than in the middle of the struts or the bridging strut itself. This is
due to the fact that the struts pull apart from each other to form a rhomboid
shape of cells during the deformation. Figure 13 presents the deformed
configuration of the C-shaped stent structure at time $t=0.025\text{s}$ with
von Mises stress contour. Different with the PS-shaped stent, the regions of
high stress in C-shaped stent are concentrated at the curved regions of the
struts. Another notable difference is the obvious expansion of the free-end
strut in the C-shaped stent because of the antisymmetric constraint link.
Generally, it is remarkable that for both structures the von Mises stress
field is reasonably captured. To the best knowledge of the authors, this is
first time that a SPH-based method is successfully extended to the simulation
of realistic cardiovascular stent and this will open up interesting
possibilities for modeling bio-mechanical applications, where this
consideration is very relevant.

Figure 12: PS-shaped stent structure: Deformed configuration with von Mises
stress contour at $t=0.02~{}s$ with present TL-SPH-KV. Here,
$\mathbf{v}_{0}=5.0~{}\text{m}/\text{s}$ is applied to impose the initial
condition and the neo-Hookean material with density
$\rho_{0}=1100~{}\text{Kg}/\text{m}^{3}$, Young’s modulus
$E=0.017~{}\text{GPa}$ and Poisson’s ratio $\nu=0.45$ is used.

Figure 13: C-shaped stent structure: Deformed configuration with von Mises
stress contour at $t=0.025~{}s$ with present TL-SPH-KV. Here,
$\mathbf{v}_{0}=5.0~{}\text{m}/\text{s}$ is applied to impose the initial
condition and the neo-Hookean material with density
$\rho_{0}=1100~{}\text{Kg}/\text{m}^{3}$, Young’s modulus
$E=0.017~{}\text{GPa}$ and Poisson’s ratio $\nu=0.45$ is used.

##  5 Concluding remarks

In this paper, we present a simple and robust artificial damping method for
stabilize the TL-SPH simulation of the solid mechanics problems involving
large deformations. The proposed stabilization strategy is based on a Kelvin-
Voigt type damper in the constitution equation and a scaling factor imitating
a von Neumann-Richtmyer type artificial viscosity. A set of numerical examples
together with very challenging cases have been investigated, and it is shown
that the present method is free of the non-physical fluctuations suffered by
the original TL-SPH. Finally, its versatility is also demonstrated by the
simulation of complex stent structures, which is a stepping stone to possible
applications in the field of bio-mechanics.

## CRediT authorship contribution statement

Chi Zhang: Investigation, Conceptualization, Methodology, Visualization,
Validation, Formal analysis, Writing - original draft, Writing - review &
editing; Yujie Zhu: Investigation, Writing - review & editing; Yongchuan Yu:
Investigation; Massoud Rezavand: Investigation, Writing - review & editing;
Xiangyu Hu: Supervision, Conceptualization, Methodology, Investigation,
Writing - review & editing.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

##  6 Acknowledgement

The authors would like to express their gratitude to Deutsche
Forschungsgemeinschaft (DFG) for their sponsorship of this research under
grant numbers DFG HU1527/10-1 and HU1527/12-1.

## References

  * [1] L. B. Lucy, A numerical approach to the testing of the fission hypothesis, The Astronomical Journal 82 (1977) 1013–1024.
  * [2] R. A. Gingold, J. J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical stars, Mon. Not. R. Astron. Soc. 181 (3) (1977) 375–389.
  * [3] L. D. Libersky, A. G. Petschek, Smooth particle hydrodynamics with strength of materials, in: Advances in the free-Lagrange method including contributions on adaptive gridding and the smooth particle hydrodynamics method, Springer, 1991, pp. 248–257.
  * [4] P. Randles, L. Libersky, Smoothed particle hydrodynamics: some recent improvements and applications, Comput. Methods Appl. Mech. Eng. 139 (1-4) (1996) 375–408.
  * [5] J. Swegle, D. Hicks, S. Attaway, Smoothed particle hydrodynamics stability analysis, Journal of computational physics 116 (1) (1995) 123–134.
  * [6] G. R. Johnson, S. R. Beissel, Normalized smoothing functions for SPH impact computations, International Journal for Numerical Methods in Engineering 39 (16) (1996) 2725–2741.
  * [7] G. A. Dilts, Moving-least-squares-particle hydrodynamics—i. consistency and stability, International Journal for Numerical Methods in Engineering 44 (8) (1999) 1115–1155.
  * [8] J. J. Monaghan, Simulating free surface flows with SPH, J. Comput. Phys. 110 (2) (1994) 399–406.
  * [9] J. P. Gray, J. J. Monaghan, R. Swift, Sph elastic dynamics, Computer methods in applied mechanics and engineering 190 (49-50) (2001) 6641–6662.
  * [10] C. Dyka, P. Randles, R. Ingel, Stress points for tension instability in SPH, International Journal for Numerical Methods in Engineering 40 (13) (1997) 2325–2341.
  * [11] C. Zhang, X. Y. Hu, N. A. Adams, A generalized transport-velocity formulation for smoothed particle hydrodynamics, Journal of Computational Physics 337 (2017) 216–232.
  * [12] R. Vignjevic, J. R. Reveles, J. Campbell, SPH in a total lagrangian formalism, CMC-Tech Science Press- 4 (3) (2006) 181.
  * [13] T. De Vuyst, R. Vignjevic, Total lagrangian SPH modelling of necking and fracture in electromagnetically driven rings, International Journal of Fracture 180 (1) (2013) 53–70.
  * [14] K. Ba, A. Gakwaya, Thermomechanical total lagrangian SPH formulation for solid mechanics in large deformation problems, Computer Methods in Applied Mechanics and Engineering 342 (2018) 458–473.
  * [15] C. Antoci, M. Gallati, S. Sibilla, Numerical simulation of fluid–structure interaction by SPH, Computers & Structures 85 (11-14) (2007) 879–890.
  * [16] A. Khayyer, H. Gotoh, H. Falahaty, Y. Shimizu, An enhanced ISPH-SPH coupled method for simulation of incompressible fluid–elastic structure interactions, Computer Physics Communications 232 (2018) 139–164.
  * [17] M. Liu, Z. Zhang, Smoothed particle hydrodynamics (SPH) for modeling fluid-structure interactions, Science China Physics, Mechanics & Astronomy 62 (8) (2019) 984701.
  * [18] C. Zhang, M. Rezavand, X. Hu, A multi-resolution SPH method for fluid-structure interactions, Journal of Computational Physics (2020) 110028.
  * [19] C. Zhang, M. Rezavand, Y. Zhu, Y. Yu, D. Wu, W. Zhang, S. Zhang, J. Wang, X. Hu, SPHinXsys: An open-source meshless, multi-resolution and multi-physics library, Software Impacts (2020) 100033.
  * [20] C. Zhang, J. Wang, M. Rezavand, D. Wu, X. Hu, An integrative smoothed particle hydrodynamics framework for modeling cardiac function, arXiv preprint arXiv:2009.03759 (2020).
  * [21] C. H. Lee, A. J. Gil, G. Greto, S. Kulasegaram, J. Bonet, A new jameson–schmidt–turkel smooth particle hydrodynamics algorithm for large strain explicit fast dynamics, Computer Methods in Applied Mechanics and Engineering 311 (2016) 71–111.
  * [22] C. H. Lee, A. J. Gil, A. Ghavamian, J. Bonet, A total lagrangian upwind smooth particle hydrodynamics algorithm for large strain explicit solid dynamics, Computer Methods in Applied Mechanics and Engineering 344 (2019) 209–250.
  * [23] R. Vignjevic, J. Campbell, L. Libersky, A treatment of zero-energy modes in the smoothed particle hydrodynamics method, Computer methods in Applied mechanics and Engineering 184 (1) (2000) 67–85.
  * [24] J. Bonet, S. Kulasegaram, A simplified approach to enhance the performance of smooth particle hydrodynamics methods, Applied Mathematics and Computation 126 (2-3) (2002) 133–155.
  * [25] C. Dyka, R. Ingel, An approach for tension instability in smoothed particle hydrodynamics (SPH), Computers & structures 57 (4) (1995) 573–580.
  * [26] Y. Vidal, J. Bonet, A. Huerta, Stabilized updated lagrangian corrected SPH for explicit dynamic problems, International journal for numerical methods in engineering 69 (13) (2007) 2687–2710.
  * [27] J. Bonet, S. Kulasegaram, Remarks on tension instability of eulerian and lagrangian corrected smooth particle hydrodynamics (CSPH) methods, International Journal for Numerical Methods in Engineering 52 (11) (2001) 1203–1220.
  * [28] M. Aguirre, A. J. Gil, J. Bonet, A. A. Carreño, A vertex centred finite volume jameson–schmidt–turkel (JST) algorithm for a mixed conservation formulation in solid dynamics, Journal of Computational Physics 259 (2014) 672–699.
  * [29] A. J. Gil, C. H. Lee, J. Bonet, M. Aguirre, A stabilised petrov–galerkin formulation for linear tetrahedral elements in compressible, nearly incompressible and truly incompressible fast dynamics, Computer Methods in Applied Mechanics and Engineering 276 (2014) 659–690.