SPH method with tensor form of smoothing parameter is proposed for high speed impact modelling. Calculation of tensor smoothing 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 integration of motion equations scheme and procedure of restoring of particle consistency is used for calculation of spatial derivatives. Boundary conditions for free surface and contact surface are realized.
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 fragmentation 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 local 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 accuracy when variable smoothing length is used. This problem arises from the fact, 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 proposed too. 1.1
Basis of SPH approximation is equation
where h is smoothing parameter, which define a radius of fluency for points, – is a space coordinate, - smoothing function.
where - integration constant, ( ) - cubical b-spline
Integration constant defined as where V is 1D, 2D or 3D volume. Function (u) usially is cubical B -spline: ∼= ︁∫
( ) ( − , ℎ ) ( , ℎ ) =
(| | /ℎ ) ℎ 3 ︁∫
=
( ) ( ) = ⎨ ⎧ 0, ≥ 2;
0.25 · (2 − )3 , ∈ (
Corresponding to (5) particle approximation is written as:
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)
Approximation of test function (
−1; 1; = −1;
Matrix of approximations for this function at n-th point, found via uncorrected SPH approximation (6):
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): , ︀( − , ℎ )︀ ︃} space step in FEM/FDM. grangian of system: Smoothing length ℎ defines maximum interparticle distance and therefore defines particle count. From the other hand, particle count is related to interparticle distance and smoothing length- greater particle count allow to use smaller smoothing length. There is some analogy between smoothing length in SPH and For calculation of elastic plastic flow, approximation (11) is applied to La(7) (8) (9) (10) (11) (12) (13) where - is particle speed, - is internal energy (per mass). Internal energy change can be written in form: where is stress tensor and ˙ is strain rate tensor of particle k respectively. Lets define correction matrix as = ︁∑ ︂( 2 − ︂)
= : ˙
This implies: ities is linear: ︃) Since the relation between a small deformations rate tensor and nodal veloc ˙ = 1 (︂ + 2 ︂) , ︁( + )︁
Because of energy conservation (14) and (15) leads to the relation between the accelerations of nodes and stress field: ( , ℎ ) = ︁( · [ℎ ]−1)︁ 2
︃) |ℎ | ︃( √︃ ︁∑
˙ ℎ
= ℎ ,
For scalar smoothing parameter spatial derivatives of smoothing function have the form: ( , ℎ ) = ℎ 3 (| /ℎ |) / = ℎ 3
2 · ℎ −1 ⃒
⃒ ⃒ ⃒ =(| ·ℎ −1|) for tensor smoothing ones have similar form: ( , ℎ ) = |ℎ | · [ℎ ]−1 ⃒
⃒ 2 ⃒ ⃒ =(| ·ℎ −1|) = ︁∑ (︀
− : ˙ )︀ ︃( ︁∑ = ︁∑ = − , / treat a boundary. was from 11 to 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 operations, such as ”ghost particles” or special integration procedures for (5) to
Smoothing length h was selected to provide sufficient count of neighbors for each particle for correction matrix to be well defined. Usually neighbor count Time step is defined via Courant–Friedrichs–Lewy condition: < min [ℎ/ ], where C is the sound speed. At practice stability condition for heat transfer problem 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 :
Initial value ℎ = ℎ , equation of evolution of ℎ is: (14) (15) (16) (17) (18) (19) (20)
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, 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
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.
Acknowledgments. This work was supported by the Russian Science Foundation (RSF) No. 16-19- 10264.