=Paper= {{Paper |id=Vol-1839/MIT2016-p24 |storemode=property |title= Tensor smooth length for SPH modelling of high speed impact |pdfUrl=https://ceur-ws.org/Vol-1839/MIT2016-p24.pdf |volume=Vol-1839 |authors=Roman Cherepanov,Alexander Gerasimov }} == Tensor smooth length for SPH modelling of high speed impact== https://ceur-ws.org/Vol-1839/MIT2016-p24.pdf 
Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling

      Tensor Smooth Length for SPH Modelling of
                 High Speed Impact

                  Roman Cherepanov and Alexander Gerasimov

       Institute of Applied mathematics and mechanics, Tomsk State University
                         634050, Lenina av. 36, Tomsk, Russia
                    RCherepanov82@gmail.com,Ger@niipmm.tsu.ru
                               http://niipmm.tsu.ru



       Abstract. SPH method with tensor form of smoothing parameter is
       proposed for high speed impact modelling. Calculation of tensor smooth-
       ing parameter is based on deformation of local coordinate system and
       can be obtained from strain tensor evolution. Strain anisotropy in such
       an approach does not cause mixing of the particles while maintaining
       the uniform distribution of particles and the resulting solution is more
       smooth. Weak variational formulation is used to construct numerical in-
       tegration of motion equations scheme and procedure of restoring of par-
       ticle consistency is used for calculation of spatial derivatives. Boundary
       conditions for free surface and contact surface are realized.

       Keywords: smoothed particles, variable smoothing length, high speed
       impact, SPH


1     Introduction
Smoothed particle hydrodynamics was introduced as a numerical method for
study of the motion of compressible gas with self-gravity of arbitrary geometry
in three dimensions in astro-physics [1], [2]. Particles are used to represent a
sub-set of the fluid elements in the Lagrangian description of a fluid, and spatial
derivatives are calculated from analytical interpolation formulae, SPH is well
suited to problems in which large deformations can occur, such as the fragmen-
tation of self-gravitating gas clouds. The spatial resolution of SPH is determined
by particle density and the smoothing length, h. Originally h was taken to be
constant, and was a globally defined function of the average density of particles
in the system. More recently, however, it has become common for each particle
to have its own time dependent smoothing length which corresponds to the lo-
cal neighbors count. Full advantage is then taken of the Lagrangian nature of
SPH, and the dynamic range in spatial resolution of the method is increased.
Time-dependent smoothing lengths can, however, lead to problems with energy
conservation in certain situations [4]. The problem arises from the fact that the
use of variable smoothing lengths means that additional terms should appear
in the particle equations of motion. Other important problem is lack of accu-
racy when variable smoothing length is used. This problem arises from the fact,

271
               Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling

that SPH approcsimation originally has form of analytical interpolation, based
on integration, and there is some difference between analytical integration and
its approximation via summation over number of particles. Nonuniform particle
distribution lead to particle inconsistency problem, as a presence of boundary
and variable smoothing length lead to one too. In should be noted, that high
speed impact is accompanied by large deformations. When smoothing length is a
scalar, smoothing kernel should have spherical symmetry, and large anisotropic
deformations in this case can cause non-physical numerical fracture and particle
mixing. In this paper we propose conservative SPH with variable time-dependent
tensor smooth parameter, and free surface boundary condition algorithm is pro-
posed too.


1.1     SPH basics

Basis of SPH approximation is equation
                             ∫︁
                         𝑖 ∼
                       𝑓 = 𝑓 (π‘₯)π‘Š (π‘₯ βˆ’ π‘₯𝑖 , β„Ž)𝑑π‘₯                                      (1)

   where h is smoothing parameter, which define a radius of fluency for points,
π‘₯– is a space coordinate, π‘Š - smoothing function.

                                          𝐢
                                π‘Š (π‘Ÿ, β„Ž) =   πœ‘ (|π‘Ÿ| /β„Ž)                               (2)
                                          β„Ž3
      where 𝐢 - integration constant, πœ‘ (𝑒) - cubical b-spline
      Integration constant defined as
                                       ∫︁
                                  𝐢=      πœ‘(𝑒)𝑑𝑉                                      (3)
                                           𝑉

where V is 1D, 2D or 3D volume. Function πœ‘(u) usially is cubical B -spline:
                       ⎧
                       ⎨ 0, 𝑒 β‰₯ 2;
                                      3
                πœ‘ (𝑒) = 0.25 Β· (2
                               (οΈ€ βˆ’ 𝑒) 2,)︀𝑒 ∈ (1, 2) ;                     (4)
                       ⎩
                         0.75 Β· 1.0 βˆ’ 𝑒 (2 βˆ’ 𝑒) , 𝑒 ∈ [0, 1] ;

      Spatial derivatives are defined via:
                                      ∫︁
                         𝑖     πœ•π‘“ 𝑖 ∼
                        𝑓,𝛼 =       = 𝑓 (π‘₯)π‘Š,𝛼 (π‘₯ βˆ’ π‘₯𝑖 , β„Ž)𝑑π‘₯                         (5)
                              πœ•π‘₯𝛼
      Corresponding to (5) particle approximation is written as:

                        𝑖      πœ•π‘“ 𝑖 ∼ βˆ‘οΈ π‘˜                 π‘šπ‘˜
                       𝑓,𝛼 =        =   𝑓 π‘Š,𝛼 (π‘₯π‘˜ βˆ’ π‘₯𝑖 , β„Ž) π‘˜                         (6)
                               πœ•π‘₯𝛼                         𝜌
                                       π‘˜

   where π‘₯π‘˜ , 𝑓 π‘˜ , π‘šπ‘˜ , πœŒπ‘˜ – radius-vector, approximated function value, mass and
density at k-th point.

                                                                                      272
Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling

   Equation (6) have 𝐢 0 -consistency [3] near surfaces and boundaries or at
non-uniform particle distribution, and for restoring particle consistency special
approaches is needed. According [3], a test vector-function is introduced in form:

                                    π›₯ (π‘₯) = (1, π‘₯0 , π‘₯1 , π‘₯2 )                    (7)
     Approximation of test function (1) and its derivatives (5) can be found with
𝐢 0 -consistency, but due to the fact that this values are known, this approxima-
tion can be used to construct correction procedure and restore particle consis-
tency. Let’s define :
                                        {οΈ‚
                                           π‘₯𝛼 ; 𝛼 > βˆ’1;
                             π›₯(π‘₯)𝛼 =                                          (8)
                                           1; 𝛼 = βˆ’1;
    Matrix of approximations for this function at n-th point, found via uncor-
rected SPH approximation (6):
                               βˆ‘οΈ                                     π‘šπ‘˜
                 𝑇𝛽𝛼 (π‘₯𝑛 ) =        π›₯𝛼 (π‘₯π‘š βˆ’ π‘₯𝑛 ) π‘Š,𝛽 (π‘₯π‘š βˆ’ π‘₯𝑛 , β„Ž)      ;        (9)
                               π‘š
                                                                      πœŒπ‘˜
      Lets define correction matrix as

                  𝑛
                                  [οΈ€        ]οΈ€βˆ’1
                 𝐡𝛼𝛽 = 𝐡𝛼𝛽 (π‘₯𝑛 ) = 𝑇𝛼𝛽 (π‘₯𝑛 )     ; 𝛼, 𝛽 = βˆ’1, 0, 1, 2;          (10)
   Now, corrected approximation of function f(x) value or its derivatives can
be found, if we build full matrix of uncorrected approximations (6) and apply
correction matrix (10):
                                   {οΈƒ                          }οΈƒ
                      πœ•π‘“ 𝑖 ∼ 𝑖       βˆ‘οΈ π‘šπ‘˜        (οΈ€ π‘˜      )οΈ€
                𝑖                            π‘˜           𝑖
               𝑓,𝛼 =       = 𝑇𝛼𝛽 Β·         𝑓 π‘Š,𝛽 π‘₯ βˆ’ π‘₯ , β„Ž               (11)
                     πœ•π‘₯𝛼                πœŒπ‘˜
                                             π‘˜

Smoothing length β„Ž defines maximum interparticle distance and therefore de-
fines particle count. From the other hand, particle count is related to interpar-
ticle distance and smoothing length- greater particle count allow to use smaller
smoothing length. There is some analogy between smoothing length in SPH and
space step in FEM/FDM.
    For calculation of elastic plastic flow, approximation (11) is applied to La-
grangian of system:
                                 βˆ‘οΈ       (οΈ‚ π‘˜ π‘˜      )οΈ‚
                                        π‘˜   𝑣 𝑣     π‘˜
                            𝐿=       π‘š           βˆ’π‘’                          (12)
                                              2
                                        π‘˜
             π‘˜                          π‘˜
   where 𝑣 - is particle speed, 𝑒 - is internal energy (per mass). Internal energy
change can be written in form:

                                         π‘‘π‘’π‘˜
                                             = 𝜎 π‘˜ : πœ€Λ™π‘˜                        (13)
                                          𝑑𝑑
where 𝜎 π‘˜ is stress tensor and πœ€Λ™π‘˜ is strain rate tensor of particle k respectively.

273
             Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling

   This implies:
                            𝑑𝐿 βˆ‘οΈ π‘˜ (οΈ€ π‘˜ π‘˜         )οΈ€
                               = π‘š 𝑣 π‘Ž βˆ’ 𝜎 π‘˜ : πœ€Λ™π‘˜                                 (14)
                            𝑑𝑑
                                     π‘˜
    Since the relation between a small deformations rate tensor and nodal veloc-
ities is linear:
                        (οΈ‚         )οΈ‚ βˆ‘οΈ 𝑛          (︁                   )︁
                       1 πœ•π‘£π‘–   πœ•π‘£π‘—        π‘š
               πœ€Λ™π‘π‘–π‘— =       +       =       𝑛
                                                 𝑛𝑝
                                               π‘Š,π‘˜          𝑝
                                                       𝑣𝑖𝑛 π‘‡π‘—π‘˜        𝑝
                                                               + 𝑣𝑗𝑛 π‘‡π‘–π‘˜    (15)
                       2 πœ•π‘₯𝑗   πœ•π‘₯𝑖      𝑛
                                           𝜌

   Because of energy conservation (14) and (15) leads to the relation between
the accelerations of nodes and stress field:
                              (οΈƒ                  )οΈƒ
                                βˆ‘οΈ π‘šπ‘˜
                        𝑛                  π‘˜ π‘˜ π‘˜π‘›
                       π‘Žπ›Ύ = βˆ’            𝜎 𝑇 π‘Š       /πœŒπ‘›                  (16)
                                     πœŒπ‘˜ 𝛾𝑗 𝑗𝛽 ,𝛽
                                         π‘˜
    For impact problem we use a condition of zero normal stress at free surface
and equal normal stress and equal speeds on contact surface.
    In described method this boundary conditions does not need additional op-
erations, such as ”ghost particles” or special integration procedures for (5) to
treat a boundary.
    Smoothing length h was selected to provide sufficient count of neighbors for
each particle for correction matrix to be well defined. Usually neighbor count
was from 11 to 16.
    Time step is defined via Courant–Friedrichs–Lewy condition: π›₯𝑑 < min [β„Ž/𝐢],
where C is the sound speed. At practice stability condition for heat transfer prob-
lem is more weak while droplet radius exceed 1 mm.
    When smoothing length h is scalar parameter intensive deformations lead to
particle mixing and artificial fragmentation. To avoid this effects we introduced
a tensor smoothing parameter h𝑖𝑗 :
                                             (οΈƒβˆšοΈƒ                      )οΈƒ
                                    𝐢             βˆ‘οΈ (︁         βˆ’1
                                                                   )︁2
                  π‘Š (π‘Ÿ, β„Žπ‘–π‘— ) =            πœ‘            π‘Ÿπ‘— Β· [β„Ž]𝑖𝑗             (17)
                                 𝑑𝑒𝑑|β„Žπ‘–π‘— |         𝑖

   Initial value β„Žπ‘–π‘— = β„ŽπΌ, equation of evolution of β„Žπ‘–π‘— is:

                                         β„ŽΜ‡π‘–π‘— = β„Žπ‘–π‘˜ 𝑣𝑗,π‘˜                           (18)

   For scalar smoothing parameter spatial derivatives of smoothing function
have the form:
                                                        βƒ’
           πœ•π‘Š (π‘Ÿ, β„Ž)   𝐢                  𝐢 π‘Ÿ Β· β„Žβˆ’1 πœ•πœ‘ βƒ’βƒ’
                     = 3 πœ•πœ‘ (|π‘Ÿ/β„Ž|) /πœ•π‘Ÿ = 3                            (19)
               πœ•π‘Ÿ     β„Ž                  β„Ž     π‘Ÿ2   πœ•π‘’ ⃒𝑒=(|π‘ŸΒ·β„Žβˆ’1 |)
   for tensor smoothing ones have similar form:
                                                             βƒ’
                   πœ•π‘Š (π‘Ÿ, β„Žπ‘–π‘— )      𝐢      π‘Ÿ Β· [β„Žπ‘–π‘— ]βˆ’1 πœ•πœ‘ βƒ’βƒ’
                                =                                                  (20)
                      πœ•π‘Ÿ          𝑑𝑒𝑑|β„Žπ‘–π‘— |      π‘Ÿ2      πœ•π‘’ ⃒𝑒=(|π‘ŸΒ·β„Žβˆ’1 |)

                                                                                    274
Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling




        Fig. 1. Deformation of local coordinate system and its relation to β„ŽΜ‡π‘–π‘—



    3 Results
    A 38 πœ‡m radius stainless steel cylinder with length of 190 πœ‡m is impacting
with a velocity of 200 m/s on stainless steel substrate 64 πœ‡m thin. Parameters
of steel: density 7.87 g/cm3, E=200GPa, G=70 GPa, Ypl=1.2MPa. As it shown
on Fig. 2, when impact speed is relatively small, results are similar.




Fig. 2. Scalar smoothing (I) and tensor smoothing (II) comparison at t=0.55 πœ‡s, 𝑣=200
m/s.



    Fig. 2, when impact speed is relatively small, results are similar. Smoothing
tensor remains almost spherical in this condition and smoothing kernel with
tensor parameter differs slightly from standard kernel with scalar smoothing
length. More different results are obtained for impact speed 800 m/s (Fig. 3).
Tensile instability and particle mixing are not observed when tensor smoothing is
used, and shape of material ejection is tracked more accurately. At the same time,
usage of constant scalar smoothing length lead to particle mixing, numerical
fracture (as it seen on bottom part of Fig. 3-I) and intensive particle clustering.
Also, symmetry of results obtained with tensor smoothing is much better.


1.2   Conclusion

Introduction of tensor smoothing length in variational SPH-code is quite simple
and does not need significant calculation cost, and provide more accurate results
for high speed impact modelling.

275
              Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling




Fig. 3. Scalar smoothing (I) and tensor smoothing (II) comparison at t=0.03 πœ‡s, 𝑣=800
m/s.


Acknowledgments. This work was supported by the Russian Science Founda-
tion (RSF) No. 16-19- 10264.


References
1. Lucy, L.B.: A numerical approach to the testing of fusion hypothesis. Astronomical
   Journal 82, 1013-1024 (1977)
2. Gingold, R.A. and Monaghan, J.J.: Smoothed particle hydro-dynamics: theory and
   applications to non-spherical stars. Monthly Notices of the Royal Astronomical So-
   ciety 181, 375-389 (1977)
3. Liu, M.B. and Liu, G.R.: Restoring particle consistency in smoothed particle hy-
   drodynamics. Applied Numerical Mathematics 56(1), 19-36 (2006)
4. Richard P. Nelson, John C.B. Papaloizou:Variable Smoothing Lengths and Energy
   Conservation in Smoothed Particle Hydrodynamics. Mon.Not.Roy.Astron.Soc. 270
   (1994) 1




                                                                                     276