A simple artificial damping method for total Lagrangian smoothed particle hydrodynamics
Chi Zhang, Yujie Zhu, Yongchuan Yu, Massoud Rezavand, Xiangyu Hu
AA simple artificial damping method for totalLagrangian smoothed particle hydrodynamics
Chi Zhang a , Yujie Zhu a , Yongchuan Yu b , Massoud Rezavand a , Xiangyu Hu a, ∗ a Department of Mechanical Engineering, Technical University of Munich, 85748 Garching,Germany b Department of Aerospace and Geodesy, Technical University of Munich, 82024Taufkirchen, Germany
Abstract
In this paper, we present a simple artificial damping method to enhance therobustness of total Lagrangian smoothed particle hydrodynamics (TL-SPH).Specifically, an artificial damping stress based on the Kelvin-Voigt type damperwith a scaling factor imitating a von Neumann-Richtmyer type artificial viscosityis introduced in the constitutive equation to alleviate the spurious oscillationin the vicinity of the sharp spatial gradients. After validating the robustnessand accuracy of the present method with a set of benchmark tests with verychallenging cases, we demonstrate its potentials in the field of bio-mechanics bysimulating the deformation of complex stent structures.
Keywords:
Total Lagrangian formulation, Smoothed particle hydrodynamics,Solid dynamics, Kelvin-Voigt damper
1. Introduction
Numerical simulation of large-strain solid dynamics problems, where flexi-ble structures experience large deformations, plays a key role in a vast rangeof engineering problems in the aerospace, automotive, manufacturing and bio-material industries. Besides the traditional mesh-based finite element methods ∗ Corresponding author.
Email addresses: [email protected] (Chi Zhang ), [email protected] (Yujie Zhu), [email protected] (Yongchuan Yu), [email protected] (Massoud Rezavand), [email protected] (Xiangyu Hu )
Preprint submitted to J. X. X. X February 10, 2021 a r X i v : . [ c s . C E ] F e b FEM), smooth particle hydrodynamics (SPH), which is a mesh-free methodand originally designed for fluid dynamics [1, 2], has also been adopted for suchproblems [3, 4], and received increasing attention in the past decades. Gener-ally, there are two types of SPH formulation for solid dynamics: one is updatedLagrangian SPH (UL-SPH) formulation, in which the current configuration isused 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 ref-erence. Compared with UL-SPH, TL-SPH shows promising potential for soliddynamics due to its attractive advantages in being free from tensile instabilityand ensuring 1st-order consistency when computing deformation gradient byintroducing a kernel correction matrix. Since its inception, TL-SPH has beenapplied for the problems of necking and fracture in electromagnetically drivenrings [13], thermomechanical deformations [14], fluid-structure interaction (FSI)[15, 16, 17, 18] and bio-mechanics [19, 20], among many others. It is known thatwithout appropriate stabilization the original TL-SPH exhibits spurious fluctu-ations in the vicinity of sharp spatial gradients (as will also be shown in Section4.1). This deficiency may result in numerical instability and lead to wronglypredicted deformation for problems involving large strain. To rectify this defi-ciency, Lee et al. [21] proposed a Jameson-Schmidt-Turkel (JST) based method,which shows good performance of eliminating spurious pressure oscillations innearly incompressible solids. In a recent work, Lee et al. [22] further proposed anew stabilization method by introducing a characteristic-based Riemann solverin conjunction with a linear reconstruction procedure. This method also showsgood performance in the simulation of nearly and truly incompressible explicitfast solid dynamics with large deformations.The main objective of this paper is to present a simple and effective arti-ficial damping method for TL-SPH to enhance its numerical stability. In par-ticular, an artificial damping based on Kelvin–Voigt type damper is introducedin the constitutive equation to alleviate the spurious oscillation in the vicin-ity of the sharp spatial gradients. By imitating a von Neumann-Richtmyertype artificial viscosity, a scaling factor is introduced to control the damping2orce. Compared with the works of Refs. [21, 22], the present method is verysimple and easy to be implemented into the original TL-SPH formulation ina straightforward way. A set of benchmark tests with very challenging casesare investigated to validate the robustness and accuracy of the present method.Furthermore, its versatility is demonstrated by modeling the deformation ofcomplex stent structures. Also, all the codes and data-sets accompanying thiswork are available from the open-source library of SPHinXsys [19] on GitHub at https://github.com/Xiangyu-Hu/SPHinXsys . The remainder of this paper isarranged as follows: Section 2 briefly summarizes the governing equations andoriginal TL-SPH formulation. Then, TL-SPH-KV is detailed in Section 3 andnumerical validations and applications are presented and discussed in Section4. Finally, brief concluding remarks are given in Section 5.
2. Preliminary
The kinematics of the finite-deformation solid dynamics can be characterizedby introducing a deformation map ϕ , which maps a material point r from theinitial reference configuration Ω ⊂ R d to the point r = ϕ (cid:0) r , t (cid:1) in the deformedconfiguration Ω = ϕ (cid:0) Ω (cid:1) . Here, we denote the superscript ( • ) as the quantitiesin the initial reference configuration. Then, the deformation tensor F can bedefined as F = ∇ ϕ = ∂ϕ∂ r = ∂ r ∂ r , (1)where the derivative is evaluated with respect to the initial reference configura-tion. Note that the deformation tensor F can also be calculated from the pointdisplacement u = r − r through F = ∇ u + I , (2)where I represents the unit matrix. 3 .1. Governing equation for solid dynamics In the total Lagrangian framework, the momentum conservation equationcan be expressed as ρ d v d t = ∇ · P T + ρ g , (3)where ρ is the density in the reference configuration, v the velocity and P thefirst Piola-Kirchhoff stress tensor. For an ideal elastic or Kirchhoff material, P is given by P = FS , (4)where S represents the second Piola-Kirchhoff stress which is evaluated via theconstitutive equation relating F with the Green-Lagrangian strain tensor E de-fined as E = 12 (cid:0) F T F − I (cid:1) . (5)In particular, when the material is linear elastic and isotropic, the constitutiveequation is simply given by S = K tr ( E ) I + 2 G (cid:18) E −
13 tr ( E ) I (cid:19) = λ tr ( E ) I + 2 µ E , (6)where λ and µ are Lam´ e parameters, K = λ + (2 µ/
3) the bulk modulus and G = µ the shear modulus. The relation between the two modulus is given by E = 2 G (1 + 2 ν ) = 3 K (1 − ν ) , (7)where E denotes the Young’s modulus and ν the Poisson’s ratio. Note that thesound speed of solid structure is defined as c = (cid:112) K/ρ . In the present work, aneo-Hookean material model defined by the strain-energy density function W = µ tr ( E ) − µ ln J + λ J ) , (8)is also applied for predicting the nonlinear stress-strain behavior of materialsundergoing large deformations. For neo-Hookean material, the second Piola-Kirchhoff stress S can be derived as S = ∂W∂ E . (9)4 .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 SPHmethod. The correction matrix B is introduced as [12] B i = − (cid:88) j V j r ij ⊗ ∇ i W ij − , (10)where ∇ i W ij = ∂W (cid:0) | r ij | , h (cid:1) ∂ | r ij | e ij , (11)denoting the gradient of the kernel function evaluated at the initial referenceconfiguration. Note that the correction matrix is computed in the initial config-uration and, therefore, it is calculated only once before the simulation. Then,the momentum conservation equation, Eq.(3), can be discretized asd v i d t = 2 m i (cid:88) j V i V j ˜ P ij ∇ i W ij + g , (12)where the inter-particle averaged first Piola-Kirchhoff stress ˜ P is defined as˜ P ij = 12 (cid:0) P i B i + P j B j (cid:1) . (13)Here, the first Piola-Kirchhoff stress tensor is computed with the constitutivelaw where the deformation tensor F is updated by the change rate evaluatedthrough d F i d t = − (cid:88) j V j v ij ⊗ ∇ i W ij B i . (14)
3. Kelvin-Voigt type damping
An elastic solid undergoing large strains can be modeled with the mechanicalcomponents of springs and dashpots. The former represents the restorative forcecomponent and the later denotes the damping component. The Kelvin-Voigt(KV) model can be represented by a purely viscous damper and purely elasticspring connected in parallel, and the total stress is decomposed into two parts σ total = σ S + σ D . (15)5ere, σ total is the total stress, σ S the elastic stress and σ D represent the damperstress as σ D = η d (cid:15) ( t )d t (16)where d (cid:15) ( t ) / d t is the strain rate and η the physical viscosity. Applying theKV model to TL-SPH formulation, the second Piola-Kirchhoff stress S can berewritten as S = S S + S D , (17)where S S is given by the constitutive equation of Eq. (6) or Eq. (8), and thedamper S D is defined as S D = π d (cid:15) ( t )d t = π (cid:34)(cid:18) d F d t (cid:19) T F + F T (cid:18) d F d t (cid:19)(cid:35) , (18)where π represents an artificial viscosity. As the main objective of introducingthe KV-type damper is to enhance the robustness of original TL-SPH, we in-troduce a von Neumann-Richtmyer type scaling factor with the speed of sound c for Eq. (18) as π = αρch. (19)Here, α is a suitable and constant parameter and h the smoothing length. Notethat a similar parameter π is also widely used in the artificial viscosity forEulerain shock-capturing schemes in modeling compressible flow. We suggest α = 0 .
4. Numerical examples
In this section, we first study three benchmark tests, where the structuresmay experience large deformation, to validate the robustness and performanceof the proposed method (denoted as TL-SPH-KV). We also compare numeri-cal results with those obtained by the original TL-SPH in which no dampingstress is applied (denoted as “TL-SPH”). Having the validation studies pre-sented, we then demonstrate the versatility of the method for applications inbio-mechanical system, i.e. stent structures. In all the following examples, the6 .75L 0.25L v Figure 1: Wave propagation in a cable: Initial configuration. th -order Wendland smoothing kernel function with an smoothing length of h = 1 . dp is employed, where dp represents the initial inter-particle spacing.The position-based Verlet scheme, which is a two-step explicit algorithm pro-posed in the work of Zhang et al. [18], is used for time integration. The CFLnumber applied here is CF L = 0 .
6, which is twice as that used in the work ofLee et al . [22].
In the first benchmark test, we investigate a simple elastic wave propagationin an elastic cable where an analytical solution is available for quantitativecomparison. This problem is also studied by Refs. [25, 26, 27, 21] by using eitherupdated or total Lagrangian SPH method. Following Refs. [25, 26, 27, 21], weconsider a rod of dimensions 10 × . × . m ) with the left end fixed and theright end free as shown in Figure 1. The wave propagation is initialized byimposing a velocity v = 5 m / s along the length direction on the right quarterof the rod. We consider a linear elastic material of density ρ = 8000 kg / m ,Young’s modulus E = 200 GPa and Poisson’s ratio ν = 0 . ime(s) V e l o c i t y ( m / s ) AnalyticalJST SPH (h/dp = 3)TL SPH KV (h/dp = 3)TL SPH (h/dp = 5)TL SPH KV (h/dp = 5)TL SPH KV (h/dp = 10)TL SPH KV (h/dp = 20)
Time(s) D i s p l a ce m e n t( m ) AnalyticalJST SPH (h/dp = 3)TL SPH KV (h/dp = 3)TL SPH (h/dp = 5)TL SPH KV (h/dp = 5)TL SPH KV (h/dp = 10)TL SPH KV (h/dp = 20)
Figure 2: Wave propagation in a cable: Velocity (upper panel) and displacement (bottompanel) profiles at the right tip end. A linear elastic material with density ρ = 8000 kg / m ,Young’s modulus E = 200 GPa and Poisson’s ratio ν = 0 . comparison between the present results and those of Ref. [21] is also reportedin Fig.2. The main velocity and displacement plateaus of these results are ingood agreement except that small overshoots are observed in those of Ref. [21]. In this second benchmark test, we consider the bending deformation of acolumn of span L = 6m and square cross section (height h = 1m), whosebottom end is clamped to the ground and its body is allowed to bend freely byimposing an initial uniform velocity v = (5 √ , , T m · s − as shown in Fig.8 v xx zz yyss Figure 3: Bending column: Initial configuration. Note that point S is located at (1 , , T m.
3. The neo-Hookean material with density ρ = 1 . × kg · m − , Young’smodulus E = 1 . × Pa and Poisson’s ratio ν = 0 .
45 is applied. This test isalso investigated by Aguirre et al.[28] and we set their results as a reference fora rigorous comparison.Figure 4 shows the deformed configuration colored with von Mises stresscontours for TL-SPH and TL-SPH-KV. Obviously, oscillations in the von Miesesstress field is obtained by TL-SPH due to insufficient stabilization while theseoscillations are suppressed by TL-SPH-KV. Similar to the last test, TL-SPHproduces noisy oscillations in the velocity profile which are eliminated by thepresent method as shown in Fig. 5. Note that, Fig. 5 also gives the timehistory of the vertical displacement of point S and its comparison with that ofRef. [28]. Compared with the results reported in latter, good agreements in thedeformation are observed. Also note that a 2nd-order convergence of the solutionis achieved by the present method with increasing spatial resolution, even thoughno linear reconstruction procedure is applied as in Ref. [22]. Compared withthe present method, the TL-SPH shows overshoots in the displacement profile9 igure 4: Bending column: Deformed configuration at different time stances for TL-SPH-KVand TL-SPH. The neo-Hookean material with density ρ = 1 . × kg · m − , Young’s modulus E = 1 . × Pa and Poisson’s ratio ν = 0 .
45 is applied. which is similar to fluctuations produced in the velocity field.10 ime(s) V e l o c it y ( z a x i s )( m / s ) TL SPH (h/dp = 6)TL SPH KV (h/dp = 6)
Time(s) D i s p l ace m e n t ( x a x i s )( m ) Aguirre et al. (2014)TL SPH (h/dp = 6)TL SPH KV (h/dp = 6)TL SPH KV (h/dp = 12)TL SPH KV (h/dp = 24)
Figure 5: Bending column: Velocity (upper panel) and displacement (bottom panel) profiles atthe right tip end. The displacement profile is compared with data in Ref. [28] and convergencestudy is also presented. The neo-Hookean material with density ρ = 1 . × kg · m − , Young’smodulus E = 1 . × Pa and Poisson’s ratio ν = 0 .
45 is applied.
Following Refs.[21, 22, 29], the example of bending column reported in Sec-tion 4.2 can be extended to a more challenging test to assess the robustnessof TL-SPH-KV for predicting the extremely highly nonlinear deformations. Asshown in Fig. 6, the twisting column is also clamped on its bottom face andthe body is initialized with a sinusoidal rotational velocity field relative to theorigin given by v (cid:0) r (cid:1) = ω × r , ω = (cid:104) . , Ω sin (cid:16) πy L (cid:17) , . (cid:105) T (rad · s − ) . (20)11 xz Figure 6: Twisting column: Initial configuration.
The column material is modeled as nearly incompressible by using a neo-Hookeanconstitutive model with density ρ = 1100 Kg / m , Young’s modulus E =0 .
017 GPa and Poisson’s ratio ν = 0 . = 105rad · s − obtained by TL-SPH and TL-SPH-KV. Clearly, TL-SPH pro-duces non-physical stress fluctuations due to insufficient numerical stabilizationand then fails to capture the correct deformation pattern. On the contrary,TL-SPH-KV alleviates these discrepancies and predicts more accurate deforma-tion 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 ofRef. [22] (see their Fig. 23) reports about slightly more twist may be due tothe slight different neo-Hookean material models. It is worth noting that TL-SPH-KV preserves the axial rotation very well without introducing out-of-axischaracteristics, as shown in Fig. 8 monitoring the time history of the horizon-12al velocity components at the point r = [0 . , . , . T . As observed, TL-SPHgenerates much larger out-of-axis fluctuations. Figure 9 shows the convergencestudy with particle refinement. Both the deformation and von Mises stressresolution obtained exhibit good convergence property for TL-SPH-KV.To demonstrate the robustness of the present TL-SPH-KV, we consider morechallenging tests by increasing the initial rotational velocity to Ω = 200 rad / sand Ω = 300 rad / s (with Poisson’s ratio ν = 0 . In the last example, complex flexible structures are investigated to demon-strate the robustness and versatility of TL-SPH-KV. As shown in Fig. 11, twodifferent stent structures, viz. Palmaz-Schatz shaped (PS-shaped) and Cyphershaped (C-shaped), are considered. Note that the present stent structures arerealistic cardiovascular stent and widely used in biomedical applications. Tostudy the deformation pattern of the stent structures, we apply a velocity fieldwith magnitude of v = 5 . · s − at the top and bottom of the structure, whichis modeled as neo-Hookean material with density ρ = 1100 Kg / m , Young’smodulus E = 0 .
017 GPa and Poisson’s ratio ν = 0 .
45 .Figure 12 shows the overall deformation of the PS-shaped stent structureat time t = 0 .
02s with von Mises stress contour. It can be observed that theregions of high stress are concentrated at the four corners of the cells ratherthan in the middle of the struts or the bridging strut itself. This is due to thefact that the struts pull apart from each other to form a rhomboid shape ofcells during the deformation. Figure 13 presents the deformed configuration ofthe C-shaped stent structure at time t = 0 . igure 7: Twisting column: Comparison of deformed configuration plotted with von Misesstress at serial time instance using the TL-SPH (upper panel) and the TL-SPH-KV (bottompanel). Results obtained with a initial sinusoidal rotational velocity Ω =
105 rad / s. A neo-Hookean material with density ρ = 1100 Kg / m , Young’s modulus E = 0 .
017 GPa andPoisson’s ratio ν = 0 . Different with the PS-shaped stent, the regions of high stress in C-shaped stentare concentrated at the curved regions of the struts. Another notable differenceis the obvious expansion of the free-end strut in the C-shaped stent becauseof the antisymmetric constraint link. Generally, it is remarkable that for bothstructures the von Mises stress field is reasonably captured. To the best knowl-14 igure 8: Twisting column: Time history of the velocity at the middle point of the free endusing TL-SPH and TL-SPH-KV. Results obtained with a initial sinusoidal rotational velocityΩ = 105 rad / s. A neo-Hookean material with density ρ = 1100 Kg / m , Young’s modulus E = 0 .
017 GPa and Poisson’s ratio ν = 0 . edge of the authors, this is first time that a SPH-based method is successfullyextended to the simulation of realistic cardiovascular stent and this will openup interesting possibilities for modeling bio-mechanical applications, where thisconsideration is very relevant.
5. Concluding remarks
In this paper, we present a simple and robust artificial damping methodfor stabilize the TL-SPH simulation of the solid mechanics problems involvinglarge deformations. The proposed stabilization strategy is based on a Kelvin-Voigt type damper in the constitution equation and a scaling factor imitating a15 igure 9: Twisting column: A sequence of particle refinement analysis using TL-SPH-KV.Results obtained with a initial sinusoidal rotational velocity Ω = 105 rad / s. A neo-Hookeanmaterial with density ρ = 1100 Kg / m , Young’s modulus E = 0 .
017 GPa and Poisson’s ratio ν = 0 . von Neumann-Richtmyer type artificial viscosity. A set of numerical examplestogether with very challenging cases have been investigated, and it is shownthat the present method is free of the non-physical fluctuations suffered by theoriginal TL-SPH. Finally, its versatility is also demonstrated by the simulationof complex stent structures, which is a stepping stone to possible applicationsin the field of bio-mechanics. 16 igure 10: Twisting column: Deformed configuration plotted with von Mises stress for largeinitial rotational velocity obtained by using the present TL-SPH-KV. (a) - (c) particle refine-ment analysis for Ω = 200 rad / s and (d) for Ω = 300 rad / s. A neo-Hookean material withdensity ρ = 1100 Kg / m , Young’s modulus E = 0 .
017 GPa is applied. Note that Poisson’sratio is set as ν = 0 . ν = 0 . CRediT authorship contribution statementChi Zhang:
Investigation, Conceptualization, Methodology, Visualization,Validation, Formal analysis, Writing - original draft, Writing - review & edit-ing;
Yujie Zhu:
Investigation, Writing - review & editing;
Yongchuan Yu:
Investigation;
Massoud Rezavand:
Investigation, Writing - review & edit-ing;
Xiangyu Hu:
Supervision, Conceptualization, Methodology, Investiga-tion, Writing - review & editing. 17 a) PS-shaped stent (b) C-shaped stent
Figure 11: Computer-aided Design (CAD) geometries of stent structures: Palmaz-Schatzshaped (PS-shaped) stent and Cypher shaped (C-shaped) stent. The corresponding CADfiles in STL format can be downloaded from our code repository or GrabCAD.
Declaration of competing interest
The authors declare that they have no known competing financial interests orpersonal relationships that could have appeared to influence the work reportedin this paper.
6. Acknowledgement
The authors would like to express their gratitude to Deutsche Forschungs-gemeinschaft (DFG) for their sponsorship of this research under grant numbersDFG HU1527/10-1 and HU1527/12-1.18 igure 12: PS-shaped stent structure: Deformed configuration with von Mises stress contourat t = 0 . s with present TL-SPH-KV. Here, v = 5 . / s is applied to impose the initialcondition and the neo-Hookean material with density ρ = 1100 Kg / m , Young’s modulus E = 0 .
017 GPa and Poisson’s ratio ν = 0 .
45 is used. igure 13: C-shaped stent structure: Deformed configuration with von Mises stress contourat t = 0 . s with present TL-SPH-KV. Here, v = 5 . / s is applied to impose the initialcondition and the neo-Hookean material with density ρ = 1100 Kg / m , Young’s modulus E = 0 .
017 GPa and Poisson’s ratio ν = 0 .
45 is used. eferenceseferences