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.

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.

𝑊 (𝑟, ℎ) = 𝐶 ℎ 3 𝜑 (|𝑟| /ℎ)(2)

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; 0.25 • (2 − 𝑢) 3 , 𝑢 ∈ (1, 2) ; 0.75 • (︀ 1.0 − 𝑢 2 )︀ (2 − 𝑢) , 𝑢 ∈ [0, 1] ;(4)

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.

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 approximation can be used to construct correction procedure and restore particle consistency. Let's define :

𝛥(𝑥) 𝛼 = {︂ 𝑥 𝛼 ; 𝛼 > −1; 1; 𝛼 = −1;(8)

Matrix of approximations for this function at n-th point, found via uncorrected 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 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 space step in FEM/FDM. For calculation of elastic plastic flow, approximation (11) is applied to Lagrangian of system:

𝐿 = ∑︁ 𝑘 𝑚 𝑘 (︂ 𝑣 𝑘 𝑣 𝑘 2 − 𝑢 𝑘 )︂(12)

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.

Mathematical and Information Technologies, MIT-2016 -modeling

This implies:

𝑑𝐿 𝑑𝑡 = ∑︁ 𝑘 𝑚 𝑘 (︀ 𝑣 𝑘 𝑎 𝑘 − 𝜎 𝑘 : ε𝑘 )︀ (14)

Since the relation between a small deformations rate tensor and nodal velocities is linear:

ε𝑝 𝑖𝑗 = 1 2 (︂ 𝜕𝑣 𝑖 𝜕𝑥 𝑗 + 𝜕𝑣 𝑗 𝜕𝑥 𝑖 )︂ = ∑︁ 𝑛 𝑚 𝑛 𝜌 𝑛 𝑊 𝑛𝑝 ,𝑘 (︁ 𝑣 𝑛 𝑖 𝑇 𝑝 𝑗𝑘 + 𝑣 𝑛 𝑗 𝑇 𝑝 𝑖𝑘 )︁(15)

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 operations, 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 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 𝑖𝑗 :

𝑊 (𝑟, ℎ 𝑖𝑗 ) = 𝐶 𝑑𝑒𝑡|ℎ 𝑖𝑗 | 𝜑 (︃ √︃ ∑︁ 𝑖 (︁ 𝑟 𝑗 • [ℎ] −1 𝑖𝑗 )︁ 2 )︃(17)

Initial value ℎ 𝑖𝑗 = ℎ𝐼, equation of evolution of ℎ 𝑖𝑗 is:

ḣ𝑖𝑗 = ℎ 𝑖𝑘 𝑣 𝑗,𝑘(18)

For scalar smoothing parameter spatial derivatives of smoothing function have the form:

𝜕𝑊 (𝑟, ℎ) 𝜕𝑟 = 𝐶 ℎ 3 𝜕𝜑 (|𝑟/ℎ|) /𝜕𝑟 = 𝐶 ℎ 3 𝑟 • ℎ −1 𝑟 2 𝜕𝜑 𝜕𝑢 ⃒ ⃒ ⃒ ⃒ 𝑢=(|𝑟•ℎ −1 |)(19)

for tensor smoothing ones have similar form: 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 𝜕𝜑 𝜕𝑢 ⃒ ⃒ ⃒ ⃒ 𝑢=(|𝑟•ℎ −1 |)(20)

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.