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