CC BY
82
Вычислительные технологии
Том 7, № 1, 2002
SMOOTH PARTICLE HYDRODYNAMICS:
SOME RESULTS *
S. K. BURUCHENKO Russian Federal Nuclear Center — Scientific Research Institute of Technical Physics, Snezhinsk e-mail: S.K.Buruchenko@vniitf.ru
Рассматривается бессеточный лагранжев метод сглаженных частиц в применении к многомерным задачам гидродинамики. Отсутствие сетки позволяет естественно моделировать произвольные закрученные и сдвиговые течения, отделение односвяз-ных зон и соединение многосвязных зон. Обсуждаются также некоторые двумерные гидродинамические расчеты, выполненные в декартовых координатах и иллюстрирующие достоинства и недостатки предлагаемого метода. Результаты моделирования высокоскоростного соударения достаточно хорошо согласуются с экспериментальными данными.
Introduction
Smooth Particle Hydrodynamics (SPH) is a pure Lagrangian method that employs no spatial mesh. It was proposed by Lucy [1], Gingold and Monaghan [2] in 1977 for astrophysical problems. In early 90s Libersky [3-5] and Benz [6-8] extended the method to the problems of continuum mechanics in which material strength is important. Besides astrophysics, the method is widely used for impact simulations (see, for example HVIS-1994, 1996, 1998) [9-15]. This is probably due to the fact that the method cannot treat arbitrary boundary terms but handles well the boundary terms typical of impact problems. Sometimes the SPH modeling is not sufficiently accurate as, for example, in the case of shock waves. The limitations of the method are discussed in [9].
The foundation of SPH is interpolation theory that allows defining a function by its values in a set of points. SPH basics can be found in [2, 16-21]. A set of moving points (particles) is considered. Mass, density, velocity, pressure and energy are known at these points. The moving particles interact and can change their neighbors. The particles are "smeared" in space by a spherically symmetric interpolation kernel with smoothing length h. The conservation laws of continuum mechanics, in the form of partial differential equations, are transformed into integral equations through the use of an interpolation function that gives the "kernel estimate" (mean value) of the field variables at a point.
For hydrodynamic equations, the transformation from continuous to discrete form is implemented for arbitrary function fj and its derivative by introducing a discrete analogue V/j
*The authors are responsible for possible misprints and the quality of translation.
© S.K. Buruchenko, 2002.
defined from the equations:
N
<f (x) mfiW(x - Xi,h), (1)
i= 1
N
< Vf (x) >= J] -X-- f XiVW(x - Xi, h). (2)
i= 1
The SPH equations of continuum mechanics are as follows (their derivation can be found, for example, in [22]:
the continuity equation
dp = Pi ^ m (V/ - j ; (3)
the momentum equation
IT = ^ j + j + %) ^; (4)
the energy equation
(5)
Here p is density, V is velocity, u is specific internal energy, P is pressure, Wij = W(xi — Xj, h) is the interpolation kernel; h is a measure of the width of the kernel or smoothing length:
dWij = W
Our version of the SPH-method uses the kernel proposed in [22]. According to [18], the artificial viscosity has the form:
n = i (-aCijßij + ßßij2)/pij, if (Vi - Vj)(xi - xj) < 0, (6) ij
= h(Vi - Vj)(x - xj) (7)
pij = (xi - xj)2 + eh2 (7)
and
Cij = (ci + Cj )/2, Pij = (Pi + Pj )/2, (8)
a = 1.0, ß = 0.5, e = 0.01, ci and Cj are speeds of sound at points i and j, respectively. We computed using the Mie — Gruneisen equation of state in the form:
0, otherwise
where
P = (kip + k2p2 + fcs^3) f 1 - ^ + rp0E, where p = P - 1. (9)
V 2 / po
The time step was:
At = 0.3- h
c + vmax
Vmax is the maximum particle velocity, h is smoothing length, is sound speed.
1. Calculations
All calculations were performed on the SGI R10 000 processor. The parameters of the equation of state (9) are given in Tab. 1. Two-dimensional calculations were performed with a constant smoothing length.
Table 1
State equations constants
Material K1, K2, K3, r P>
Kbar Kbar Kbar g/cm3
Copper 1407 2871 2335 2.04 8.93
Aluminum 765 1659 428 2.13 2.71
Zinc 662 1577 1242 2.38 11.346
1.1. The shock tube problem
The problem presented here has an exact solution. The initial conditions are by Sod [23]:
0 <x< 0.5, p = 1, P = 1, E = 2.5, V = 0; 0.5 <x< 1.0, p = 0.125, P = 0.1, E = 2.0, V = 0, where p is density, P is pressure, E is energy, V is velocity. The equation of state for ideal gas is P = Ep(Y — 1), where y = 1.4.
Fig. 1, a, b and c show velocity, pressure and density profiles at t = 0.23. The calculation used 400 particles. Initially the particles were uniformly distributed in space. Three runs were made: one with a constant smoothing length and two with a variable one. In the first case we used the algorithm of Benz [21] and that one of Steinmetz — Muller [24] in the second case. The shock wave was smeared within the range of 6-8 particles and the contact discontinuity within that of 14-16 particles.
1,2 - . . ........,........,........,........,........,........,........,........,........,........,........,........,.......,........,........,......., ,........,........,.....
/ 1
r / i 1 1
! 1 f t
Be rix h ..... Constant 1» --'itpimiii'tv.Miilli'f I-J 1 1
I ----I:ixa.ct solution
i_■
1 1 1 1 11 111 1 1 1 1 1 1 1 1
0.1 0:2 o.l 0.4 0.Í 0.G 0.7 OS 0.9 1.0
0.0
Ren/ li Constant ti Steinmetz-Muller li Exact solution
0.0
0.2
0.4
0.6
0.8
1.0
0 8 -
0 4
02
00
Benz li Constant li Steinmetz-Vlul lor 1i Exact solution
0.0 0.1 0,2 0.3 0 4 0 5 0.6 0.7 0.8 0.9 10
Fig. 1. Velocity (a), pressure (b), density (c) profile (t = 0.23).
b
c
The calculation with the constant smoothing length gave the worst result. We supposed the reason was that the masses of particles in the right and left parts of the region differed by the factor of 8. To prove this we performed calculations using the above algorithms for 44 particles in the right part and 356 particles in the left one that provided equal particle masses in the parts. Calculations show that in this case using a constant smoothing length gives a good result, which is close to that one obtained by the Steinmetz — Muller algorithm.
1.2. The development of Rayleigh — Taylor instability
Consider a plane system of two layers of heavy and light fluids. An acceleration g = 100 is pointed from the heavy fluid to the light one. The computational region L in Cartesian coordinates (XY) is:
L = (x,y) | 0 < x < 1; -2 < y < 2.
The contact surface Y = 0 divides the region into two zones: with heavy and light fluids. The contact surface is not perturbed. Instead we set the initial distribution of velocity vector u(u,v) from the condition Vu = 0. Below are formulas for the components of velocity vector, which satisfy this condition:
U = Uo sin(kx)(2H(y) - 1) exp(-k | y |), V = VOcos(kx) exp(-k | y |),
where
H (y)
0,
if y>0, if y<0,
k = 2n/A is the wave number, A is the wave length, UO and VO are the initial amplitudes of perturbation. Heavy material: density = 10, sound speed squared = 2000. Light material: density = 1, sound speed squared = 200. Equation of state: isothermal gas.
Fig. 2. The development of Rayleigh — Taylor instability p'/p = 10 : 1.
The results of calculations are shown in Fig. 2. One can see linear and non-linear phases in the growth of perturbation. The calculation used 50 x 200 particles. This seems insufficient for the detailed simulation with a constant smoothing length, however, the calculation captures the essential features of flow except the inward eddying of "mushroom cap" by the heavy fluid. We can expect that using a variable smoothing length and a greater number of particles will ensure more exact simulation.
1.3. Cumulative jet formation
A 7.07-mm-thick copper plate impacts on a rigid surface at 4 km/s. The velocity is normal to the plate. The angle between the plate and the rigid surface is 72°. The calculation used 1522 particles. The initial distance between particles was 0.3535 mm. The smoothing length was taken to be 0.707 mm. The calculation was run to 14 ^s (Fig. 3). It took 1425 cycles and 6 min on the SGI Power Challenge 10 000 computer. Fig. 4 shows how the maximum and minimum velocities of jet particles vary with time. In accord to the hydrodynamic theory of cumulation which basically proposes that jet matter is incompressible, the jet velocity is defined as V = U x ctg(a/2) = 5505 m/s. Our calculation gave the maximum jet velocity of 6580 m/s that differs by 19.5 % from this estimate.
Fig. 3. The generation of cumulative jet.
7000
6000
5000
4000
3000
2000
maximum v minimum ve ilOCIty locity
\ \ *
\ \ \ \
0.0 2.5 5.0 7 5 10.0 12.5 15.0 17.5
Fig. 4. Maximum and minimum velocities vs time.
20 o
1.4. Explosion in the center of sphere
A symmetric motion of gas resulted from an explosion in homogeneous matter is considered. In a 0.1-radius-sphere, the internal energy per unit mass is 10 ■ 106 (Fig. 5, a). In a spherical layer from 0.1 to 1 the energy is equal to zero. Boundary terms: the velocity component normal to the sphere of the unit radius is equal to zero. The matter is ideal gas of unit density and y = 1. We used 100 particles along radius. On the circle, particles were 0.01 distant from each other. So, distances between all particles were almost equal. Totally 15 916 particles were used and the calculation was run to t = 3.2 ■ 10-3. Fig. 5, b, 5, c and 5, d show the compression wave at t = 0.202 ♦ 10-3, 1.623 ♦ 10-3, 3.255 ♦ 10-3, respectively in the form of density distribution as it plott< re for
y
Fig. 5. The compression wave at: (a) t = 0; (b) t = 0.202 • 10 3; (c) t = 1.623 • 10-3; (d) t = 3.255 • 10-3.
Fig. 6. Density profile.
Fig. 7. Velocity profile.
y = 0 and x > 0. It is seen that the compression wave moves from the center of the 0.1-radius sphere. The wave is not spherical with strongly varying velocities in the region of expanding gas.
1.5. Copper disk impacting on aluminum plate
The problem is taken from [4]. A 3-g copper disk (11.18-mm diameter x 3.45-mm thick) normally impacts on a 2.87-mm thick aluminum plate at 5.55 km/s. The calculation used 16 272 particles: 96 x 30 = 2880 in the disk and 558 x 24 = 13 392 in the plate. We did not use symmetry, i e. the calculation was made for real geometry. The initial distance between particles was set to be 0.115 mm, the smoothing length was 0.23 mm. The calculation was run to 6.4ps. It took 3551 cycles and 3 h on the SGI Power Challenge 10 000 computer.
The purpose of the calculation was to compare the hydrodynamic solution with the solution obtained in [4] by a model that allows for strength and also with experiment [25]. Since the impact velocity is rather high and we compare the shape of debris cloud below the bumper plate, calculated results are expected to be close. Our result agrees well with that one from [4] and experiment.
There is a good experiment for this problem in [25]. It is a radiograph of the disk and the plate 6.4 ps after impact (Fig. 8, a). Both the matters can be clearly distinguished because the density of copper is about three times higher than that of aluminum. Fig. 8, b illustrates the calculation. The peculiar shape of the cloud consisting of aluminum and the remnants of the copper disk is nicely captured by the simulation. The experimental and calculated figures are scaled equally and one can use a ruler to compare various dimensions. Such measurements show superb agreement between experiment and calculation.
1.6. The impact of zinc cylinder on zinc plate
Ref. [26] provides experimental data on the impacts of differently shaped zinc projectiles on a zinc plate at different velocities. We chose test T4-1554 to verify our modeling of debris cloud.
ffigpHRR
mmm
Fig. 8. Impact of disk (Cu) on plate (Al), v = 5.55 km/s, t = 6.4 ^s: (a) — experiment; (b) — calculation [25].
1.7. Experimental data
Test Т4-1554 [26] uses commercially pure zinc for projectile and bumper plate. Projectile is a 3.98-mm-diam-cylinder 14.15 mm long. A light-gas gun accelerates it to 4.97 km/s. The bumper plate is 0.965 mm thick. The numerical parameters of debris cloud were defined at three times (Tab. 2).
Experimental data
Table 2
Time, ps Debris cloud length, mm Debris cloud front velocity, km/s
10.5 55.4 —
18.2 96.8 5.38
25.9 135.4 5.01
a
>
/
f-
Fig. 9. The impact of zinc rod on zinc slab (10.5 ps after impact): (a) — experiment; (b) — calculation [26].
Fig. 9, a and 10, a depict the experimental radiographs of debris cloud at t = 10.5 and 18.2 ^s, respectively after the impact of the rod on the plate. Fig. 9, b and 10, b show the results of calculation. The cloud consists of three parts: the main part is like a balloon, the remnant of the rod in the balloon and a cone below the balloon. The length of debris cloud given in Tab. 2 is measured from the plate side opposite to the impacting rod.
1.8. SPH calculations
In the experiment described in [26], there was a second plate used for determining the size of hole. We simulated only one plate because our goal was to compare the shape of debris clouds obtained in calculation and experiment.
Two calculations were made with different numbers of particles. The plate was 40 mm long. The first calculation totally used 9980 particles: 40 x 146 = 5840 in the projectile and 414 x 10 = 4140 in the plate. In the second calculations these figures were 19 554, 56 x 204 = 11 424 and 580 x 14 = 8120, respectively.
Fig. 10. The impact of zinc rod on zinc slab (18.2 ^s after impact): (a) — experiment [26]; (b) — calculation.
Experimental and calculated results are shown in Tab. 3 and 4.
Table 3
Experimental and ca culated data
Time, Experiment Calculation 9880 particles
Debris cloud Debris cloud Debris cloud Debris cloud
ßs length, front velocity, length, front velocity,
mm (%) km/s (%) mm (%) km/s (%)
61.63 (11)
10.5 55.4 56.54* (2.06 )
106.03 (9.53) 5.94(10.4)
18.2 96.8 5.38 97.24* (0.46) 5.28* (1.89)
151.03 (11.54) 5.94 (18.56)
25.9 135.4 5.01 139.63* (3.12) 5.505* (9.88)
*At the cone vertex.
Table 4
Experimental and calculated data
Time, Experiment Calculation 19 544 particles
Debris cloud Debris cloud Debris cloud Debris cloud
ßS length, front velocity, length, front velocity,
mm (%) km/s (%) mm (%) km/s (%)
10.5 55.4 59.2 (6.58) 5.73
56.48* (0.19)
18.2 96.8 5.38 103(6.4) 5.73 (6.5)
97.265* (0.48)
25.9 135.4 5.01 148.03 (9.3) 5.73 (14.3)
138.1* (1.99)
*At the cone vertex.
it "V
t
» ' « V
f
Fig. 11. The impact of zinc rod on zinc slab at 4.97 km/s (t = 10.5 ^s after impact): (a) — SPH calculation with 9800 particles; (b) — SPH calculation with 19554 particles.
Fig. 9 and 10 show the results of numerical simulation at 10.5 and 18.2 ^s. Fig. 11 shows the rod and the plate at 10.5 ^s simulated with the different numbers of particles.
If compare experiment (Fig. 9, a and 10, a) and calculation (Fig. 9, b and 10, b), it becomes clear that all the important details seen in the radiograph are reproduced by the calculation. See, for example, the cone part of debris cloud. In the SPH calculation, there is a small drop at the very front of the cone, which is absent in experiment. The cone is absent in the CTH calculation described in [26]. In experiment and the CTH calculation, the diameter of hole in the plate is equal to 10 mm while in the SPH calculation it is equal to 13.709 mm. This disagreement can be due to not accounting for material strength.
In general, the calculation reproduces all essential features of experiment.
Conclusion
A computer code implementing 2D Smooth Particle Hydrodynamics has been developed. Results for hypervelocity impact calculations agree well with experimental data. Several calculations demonstrate the key features of the technique.
Further development proposes the allowance for material strength and 3D simulation.
Acknowledgements to: Joe Monaghan, Monash University, Australia; Willy Benz, University Bern, Switzerland; Gordon R. Johnson, Alliant Techsystems Inc., MN, USA; Matthias Steinmetz, Universirty of Arizona, USA; Ewald Muller, Max-Planck-Institut fur Astrophysik, Germany; for sending their papers. Special thanks go to Willy Benz for valuable consultations and to Roland Speith, Institut fur Astronomie und Astrophysik, Universitat Tubingen, Germany for his code SPH3D publicated in Internet, which served as the base of the code developed by author; also to A. V. Polionov (RFNC) for his interested support to the work. Author extends thanks as well to Dmitry Mogilenskikh and Margarita Shishkina for help in presentation of the
results.
References
[1] LUCY L. B. A numerical approach to the testing of the fission hypothesis // The Astronomical J. 1977. Vol. 82, No. 12. P. 1013-1024.
[2] GlNGOLD R. A., MONAGHAN J.J. Smoothed particle hydrodynamics: theory and application to non-spherical stars // Mon. Not. R. Astr. Soc. 1977. Vol. 181. P. 375-389.
[3] Libersky L.D., PETSCHEK A. G. Smooth particle hydrodynamics with strength of materials: Advances in the Free Lagrange Method // Lecture Notes in Phys. 1990. Vol. 395. P. 248-257.
[4] Libersky L.D., Petschek A.G., Carney T. C. et al. High strain Lagrangian hydrodynamics. A tree-dimensional SPH code for dynamic material response //J. Comput. Phys. 1993. Vol. 109. P. 67-75.
[5] Petschek A. G., Libersky L. D. Cylindrical smoothed particle hydrodynamics // Ibid. P. 76-83.
[6] BENZ W., ASPHAUG E., RYAN E. V. Numerical simulations of catastrophic disruption: recent results // Planet. Space Sci. 1994. Vol. 42, No. 12. P. 1053-1066.
[7] BENZ W., ASPHAUG E. Impact Simulations with Fracture. Pt I: Method and Tests // Icarus. 1994. Vol. 107. P. 98-116.
[8] BENZ W., ASPHAUG E. Simulations of brittle solids using smooth particle hydrodynamics // Comput. Phys. Commun. 1995. Vol. 87. P. 253-265.
[9] FULK D. A., QuiNN D. W. Hybrid formulations of smoothed particle hydrodynamics // Int. J. Impact Eng. 1995. Vol. 17. P. 329-340.
[10] HAYHURST C. J., CLEGG R. A. Cylindrically symmetric SPH simulations of hypervelocity impacts on thin plates // Ibid. 1997. Vol. 20. P. 337-348.
[11] GROENENBOOM P. H. L. Numerical simulation of 2D and 3D hypervelocity impact using the SPH option in PAM-SHOCK // Ibid. P. 309-323.
[12] Hiermaier S., Konke D., Stilp A.J., THOMA K. Computational simulation of the hypervelocity impact of Al-spheres on thin plates of different materials // Ibid. P. 363374.
[13] Libersky L.D., Randles Ph.W., Carney T.C., Dickson D. L. Recent improvements in SPH modeling of hypervelocity impact // Ibid. P. 525-532.
[14] VlNGJEVIC R., CAMPBELL J. A penalty approach for contact in smoothed particle hydrodynamics // Ibid. 1999. Vol. 23. P. 945-956.
[15] FARAUD M., Destefanis R., Palmieri D., Marchetti M. SPH simulations of debris impacts using two different computer codes // Ibid. P. 249-260.
[16] MüNAGHAN J.J. Particle methods for hydrodynamics / / Comput. Phys. Rep. 1985. Vol. 3. P. 71-124.
[17] MüNAGHAN J.J. Smoothed particle hydrodynamics / / Ann. Rev. Astron. and Astrophys. 1992. Vol. 30. P. 543-574.
[18] MüNAGHAN J.J. An introduction to SPH // Comput. Phys. Communications. 1988. Vol. 48. P. 89-96.
[19] MüNAGHAN J.J, GlNGOLD R.A. Shock Simulation by the Particle Method SPH // J. Comput. Phys. 1983. Vol. 52. P. 374-389.
[20] MüNAGHAN J.J. Why particle methods work // SIAM J. Stat. Comput. 1982. Vol. 3, 4.
[21] BENZ W. Smoothed Particle Hydrodynamics: a Review // The numerical modeling of Nonlinear Stellar Pulsation / J. R. Buchler (ed.). N.Y.: Kluwer Acad. Publ., 1990. P. 269-288.
[22] MüNAGHAN J. J., LATTANZIü J. C. // Astron. and Astrophys. 1985. Vol. 149. P. 135.
[23] SüD G. A. A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws //J. Comput. Phys. 1978. Vol. 27. P. 1-31.
[24] STEINMETZ M., Müller E., On the capabilities and limits of smoothed particle hydrodynamics // Astron. and Astrophys. 1993. P. 268, 391.
[25] PlEKUTüWSKl A. J. // Int. J. Impact Eng. 1990. Vol. 10. P. 453.
[26] Müllin S.A., Littlefield D.L., Andersün C.E., Jr., Chhabildas L.C. An examination of velocity scaling // Ibid. 1995. Vol. 17. P. 571-581.
Received for publication May 3, 2001
