- Email: [email protected]
Development of modified SPH approach for modeling of high-velocity impact
International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
Contents lists available at SciVerse ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Development of modified SPH approach for modeling of high-velocity impact K.A. Alhussan a, V.A. Babenko b, I.M. Kozlov b, A.S. Smetannikov b,⇑ a b
Space Research Institute, King Abdulaziz City for Science and Technology, Riyadh, Saudi Arabia Heat and Mass Transfer Institute, Brovki Str., 15, 220072 Minsk, Belarus
a r t i c l e
i n f o
Article history: Received 29 February 2012 Accepted 24 May 2012 Available online 13 July 2012 Keywords: High-velocity impact Shock waves Hydrodynamics Numerical methods
a b s t r a c t The problem of high-velocity impact was modeled on the basis of two numerical methods, the smoothed particle hydrodynamics (SPH) and the Eulerian total variation diminishing (TVD) finite-difference scheme. The modification of SPH method was elaborated. In this modification, the summation approach and the continuous approach are equivalent, and the problem of boundary particle deficiency is successfully resolved. The comparison of modified SPH and Eulerian finite-difference approaches shows a good agreement between them. The calculations results can be useful for estimation of efficiency of anti-meteorite protection shields, giving information on the shape and size of formed apertures and jets. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction The study of dynamics of the processes, which take place at high-velocity impact (HVI), is of interest for many problems of space physics and astrophysics: the study of meteoric craters, the origin of planetary atmospheres, and the consequences of large space object fall on the Earth etc., [1]. The long stay of space vehicles at an orbit results in practically inevitable influence of meteorites on their surface. Thus, the probability of collision with small particles (micrometeorites) with mass less than 0.01 g is especially great. At high velocities (10 km/s) the initial pressure near to impact place makes about megabar, and the specific energy exceeds the sublimation heat. At such action the material melts and evaporates. The time of action of high pressure appears to be rather short (owing to the small size of a meteorite), the shock wave quickly damps. At moving away from the impact place the specific energy of meteorite decreases, to become possible only melting. Then it falls below the melting energy, and at this stage there can be only fracturing of material. This process is described by creation and growth of microcracks that is determined by strength characteristics of material. We shall consider an initial stage of impact, when strength effects are insignificant. The experimental data on HVI concern to the velocities not exceeding 10 km/s. Thereof, in the range of high collision velocities characteristic of anti-meteorite protection the theoretical research of impact dynamics (in particular numerical modeling) is very important. ⇑ Corresponding author. E-mail address: [email protected] (A.S. Smetannikov). 0017-9310/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2012.06.002
The majority of traditional approaches for numerical modeling in hydrodynamics use a computational mesh for calculation. The Eulerian approaches study the gas motion using space mesh, which is not changed in time (Eulerian mesh). The Lagrangian approaches use a mesh moving together with gas. These mesh based techniques solve the fluid dynamic equations by calculating interaction with adjacent cells. The Eulerian description is a spatial description, whereas the Lagrangian description is a material description. The fundamental difference of these two descriptions is that the Lagrangian description employs the total time derivative as the combination of local derivative and convective derivative. Methods of particles refer to a class of Lagrangian methods. The media in these methods is represented by a set of particles, and movement of these particles is investigated. Among such particle methods the method of smoothed particle hydrodynamics (SPH) holds a special place because it is mesh-free. In the SPH method the equation of gas dynamics are solved on the basis of interaction of each particle with its nearest neighbors. Due to Lagrangian description, the problem can be reduced to solution of ordinary differential equations, each of them being formulated for the velocity, density, position, and energy of a particle. Many of studies of HVI in the past were performed with Eulerian approach. Nevertheless, this class of numerical models has unavoidable difficulties as applied to HVI investigation. In this problem there are solid objects moving with high velocities and having, at least initially, well-defined boundaries. We are interested mainly on transformation of these boundaries with time. But these boundaries can slide one over another with very high velocities, and in frame of Eulerian approaches it is very difficult getting answer on the real position of the material
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
interfaces. Besides, this sliding creates very serious problems with the numerical stability. These difficulties can be resolved and overcome to some extent by using complicated Eulerian methods such as TVD schemes with markers describing positions of the interfaces. But at the later stages of impact when the topology of created cracks can be very complex, it is difficult to apply this technique. That is why Lagrangian approaches, and especially SPH method have received wide acceptance in HVI studies. SPH is a new effective technique for impact simulation, which combines harmonically the benefits of Lagrangian methods and particle approximation. The books by Zukas [2,3], the papers by Libersky with co-authors [4,5], by Randles et al. [6,7], and by Johnson et al. [8,9] have made outstanding contributions in the application of SPH to impact problems. Many papers were made afterwards in order to combine SPH with finite element methods (FEM). Different contact algorithms have been used with SPH when combining the SPH method with the FEM [10], but the only purpose of these algorithms was to allow SPH nodes to interact with FEM mesh. An improved hybrid particle-finite element method has been developed for hypervelocity impact simulation in [11]. The SPH method possesses some advantages and disadvantages over traditional methods of numerical prediction. The main privileges of SPH consist in that it is remarkably robustness in sense of numerical stability and relatively simplicity of programming. The main disadvantage is probably insufficient accuracy if the number of particles in modeling is not great enough. Just the number of particles is a key problem in SHP. An amount of computational work in reported SPH programs grows with the number of particles N as the factor between N2 and N log N. That is why even dealing with optimized computer problem, user should distribute the particles very economically. Actually, the number of particles is usually bounded by the value of several hundred of thousands or a few million. So, at actual computations the value of smoothing length h can be large enough as the total amount of particles is bounded by computational time. Even for this restricted number of particles, it is important keeping balances of mass, momentum and energy exactly for the system as whole and for its component parts, as it actually became a standard in the theory of finite-difference schemes. In SPH, however, writing fully conservative numerical algorithms is rather problematic by many reasons. It will be explained later after formulation of main equations in SPH. Here, it worth noting only that the density approximation is especially important in this connection, since the density basically determines the particle distribution and the smoothing length evolution. The error arising in continuity equation is transferred to the momentum and energy equations thus worsened the modeling quality. Another important problem in SPH is a so called boundary particle deficiency. Near the boundaries, the support domain for a particle with smoothing length h can intersect the whole domain of the problem, thus reducing smoothed density of the particle and leading to erroneous results. The authors of [12], for instance, showed that SPH introduces spurious pressure forces acting on particles in regions where there are steep density gradients. Due to reduced density the pressure at the particle can increase substantially thus creating a source of numerical instability. Boundary particle deficiency can contribute also into the error in continuity equation. In this paper, a modification of SPH methods is proposed and investigated in application to HVI problem. Comparison is made with the results of application of the TVD difference scheme with markers. Suggested in this paper an improvement of SPH enables to keep the mass balance exactly and to eliminate the boundary particle deficiency problem.
6341
2. Governing equations For the Eulerian description of axial symmetrical flow arising at impact, the system of gas dynamic equations in two-dimensional representation reads:
@u @u @u @P @v @v @v @P þv þu ¼ ; q þv þu ¼ ; @t @r @z @z @r @t @r @z 1 dq 1 @ðr n v Þ @u de 1 @ðr n v Þ @u þ ; q ¼ P n þ ; P ¼ Pðe; qÞ; ¼ n r @r @z r @r @z q dt dt
q
ð1Þ
where t is the time, q the density, P the pressure, e the specific internal energy, v, u are the radial and axial components of velocity vector, n is the index of symmetry (n = 0 for plane and n = 1 for cylindrical geometry, correspondingly). The equation of state P = P(e, q) closes this system of equations. The complexity of numerical calculation of high-velocity impact is caused by prominent features of arising flow. At the initial stage of impact the flow is sharply spatially non-uniform, and there are areas of materials sliding. Prediction of such flows in Lagrange representation is extremely complicated. For later times, the other type of complexity is connected to considerable growth of spatial scale of flow. The system of gasdynamic Equation (1) should be approximated by some difference scheme for numerical solution. For this, we use two difference schemes. The first one is a method of total variation diminishing (TVD) in Eulerian variables [13,14], the second one the smoothed particle hydrodynamics (SPH) method in Lagrangian variables [15,16]. In the first method, the spatial mesh is used, the second method is meshfree. More detailed description of these methods will be made in the next section. The calculations are carried out for impact of a striker (meteorite) on plane aluminum shields (with thickness D) with speed u0, directed normally to the shield. The meteorite is simulated by a sphere with diameter D. As a material of striker, aluminum is taken also. For closing the system of gas dynamic equations it is necessary to determine the equation of state for materials, i.e. the dependence of pressure P = P(e, q) from energy and density. In calculations we used the Mie–Gruneisen equation of state in form proposed by Tillotson [17]. For material density higher than the normal density of a solid state q P q0 for all specific energy e (and in small region with q < q0 and e < e1) we assumed the relation
" P ¼ Pc ¼ a þ
#
b 1 þ e=ðe0 d2 Þ
qe þ Al þ Bl2 :
ð2Þ
With density lower than that a solid state q 6 q0 and e > e2 we used the equation P ¼ P H ¼ aqe þ
" ! 2 # 1 1 exp a 1 þ Al exp b 1 : d d 1 þ e=ðe0 d Þ b qe
2
ð3Þ
In these expressions, d = q/q0, l = d 1, and e1, e2 are the energies of incipient and full vaporization. For density smaller than that for solid state q 6 q0 and energy e1 < e < e2 we used next interpolation by energy between these two relations
P ¼ ½Pc ðe2 eÞ þ PH ðe e1 Þ=ðe2 e1 Þ:
ð4Þ
In calculations, the following parameters of aluminum were taken: a = 0.5, b = 1.63, e0 = 5 MJ/kg, A = 75 GPa, B = 65 GPa, q0 = 2790 kg/m3, a = 5, b = 5, e1 = 3 MJ/kg, e2 = 15 MJ/kg. It should be noted that in the unloading region (at densities somewhat smaller than that for solid state q0) at low energies the Tillotson equations of state gives very large negative pressures. In fact, when the
6342
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
negative pressure exceeds a magnitude of the order of 0.1 GPa, there occurs a breaking (or splitting off) the material, and the pressure goes to zero. In this connection, when it turned out that the pressure was negative, we assumed that p = 0. A similar procedure in this region was employed in [18]. 3. Technique of calculation In this section we will describe shortly the techniques used in the present work for calculation of flow hydrodynamics. 3.1. TVD method The technique of numerical solution of the system of Euler equations in application to the problem of the frontal impact of a spherical striker on a plane shield of finite thickness in vacuum is considered. A characteristic feature of the method is that the Eulerian approach for solving the gas dynamics equations is combined here with the description of the motion of contact boundaries between various objects. The impermeable contact boundaries are described with the aid of marker particles moving together with the substance. The application of markers imparts to the Euler scheme some features typical of special importance in modeling collision of solid bodies with distinct boundaries. Moreover, information about contact boundaries makes it possible to easily identify various objects in space and, consequently, to appropriately describe their thermodynamic and strength properties. The setting of boundary conditions at interfaces between different media allows one to describe the interaction of the objects that possess the properties of both gases and solids. A numerical scheme is realized on rectangular spatial grids and can be used in modeling two-dimensional problems possessing an axial or plane symmetry. As a difference scheme, use is made of the scheme of the second order of accuracy in space obtained by the method of finite volumes. The second order accuracy in time is provided by ‘‘predictor–corrector’’-type scheme. The scheme does not employ an artificial viscosity, with the stability of solution being ensured by application of the numerical TVD dissipation that does not alter the second order of approximation. To calculate the TVD correction of fluxes, an approach based on the use of the Riemann invariants is applied making it possible to be got rid of time-consuming matrix operations. For each characteristic solution the correcting fluxes are calculated with the aid of the formula for the scalar upwind TVD scheme [19,20]. In order to describe contact boundaries with the aid of markers in the two-dimensional case, the notion of a marker object is introduced. The marker object is understood to be a certain region of space entirely occupied by a medium of any one species and bounded by one or several marker lines. The marker line is a continuous curve formed on a plane by segments of a straight line that connect a set of markers ordered in space. A straight-line segment that connects two successive markers is called a marker segment. Each marker line unites its own sequence of markers and has a certain orientation starting with the first marker and ending with the last marker of the sequence. The marker is a point on a plane assigned by a pair of coordinates and having a velocity vector. The velocity of a marker is determined by interpolation of velocity components from neighboring nodes of the grid. In the ‘‘predictor– corrector’’-type of the scheme the velocity field for markers is calculated at the ‘‘predictor’’ stage. In addition to markers, configuration and index arrays are used for describing a marker object. The former serves for determining which specific computational cell belongs to that or other marker object. The latter establishes the connection between the marker object and the material medium that fills this object. For each cell crossed by one or several marker
lines the mean direction of these lines inside a cell is calculated and is represented in the form of a unit vector n. In a two-dimensional case, when pulse exchange through the contact is calculated, it is possible to take into account the true calculation of the marker line relative to the coordinate lines of the rectangular grid. For this purpose, a unit vector n is assigned to each cell crossed by the marker line. Moreover, taking into account the direction n of the marker line, one can model the slipping between the media. The numerical description of the contact boundary is composed of two parts: 1. Determination of the position of the contact boundary in space. 2. Determination of the fluxes of gas-dynamical quantities through the contact boundary. The contact boundary corresponding to the marker line passes only along the boundaries of the cells, i.e. each cell entirely belongs to only one marker object or background. Whether a concrete cell, crossed by one or several marker lines, belongs to that or other marker object or background is determined by comparing the areas occupied by different objects inside a cell. A cell is ascribed to the marker object (or background) which has the largest area. The position of the contact boundary changes discretely accurate to the magnitude of a spatial cell. The status of a substance in the cell through which the marker line passes, also changes discretely. If at a certain moment the status of a contact cell changes, an analysis of its neighbors is carried out, on the basis of which the given cell is assigned some parameters average over the neighboring cells. 3.2. SPH method SPH is a numerical method that describes a fluid flow in continuous media as a motion of the collective of interacting particles. For each particle i its characteristics such as mass mj, position ri, velocity vi, density qi and energy density ei are determined. The modification of SPH method in application to HVI, which is developed in the given paper, concerns mainly to the SPH approximation of mass conservation law. In literature, various SPH approximations for density equation are described (see [16]). The first approach is the summation density (5), which directly applies the SPH approximations to the density itself. For a given particle i, the density with summation density approach is written in the form of
qi ¼
N X mj W ij ;
ð5Þ
j¼1
where N is the number of particles in the support domain of particle i, and mj is the mass associated with particle j. Wij is the smoothing function of particle i evaluated at particle j with the smoothing length h.
W ij ¼ W jxi xj j; h ¼ Wðrij ; hÞ; where Wij = W(rij, h) is the kernel function with smoothing length h = (hi + hj)/2, where hi denotes the smoothing length associated with the ith particle and rij = jxi xjj is the distance between particles i and j. As a main kernel function we use Gauss form. The Gauss kernel is actually the ideal function for the purpose, but it takes computational time approximately half as much again as the cubic spline kernel. The summation density approach (5) is, at a first glance, the most natural. It conserves the mass exactly since the integration of the density over the entire problem domain is exactly the total mass of all the particles. Nevertheless, in so doing, the accuracy of density approximation is lost near the boundaries, where the approximation becomes poor due to the mentioned above particle
6343
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
deficiency. To improve the situation a modified form of summation density approach in form of relation (6) was proposed by many authors. It consists in normalization of the RHS of Eq. (5) with the SPH summation of the smoothing function over the neighboring particles:
PN
mj W ij
qi ¼ Pj¼1 N mj j¼1 qj
W ij
:
ð6Þ
qsi ¼
Here, the correction of (5) is made, which takes into account a non equality of the denominator to unity. The possibility of such correction at improving the SPH computation results is not very impressive. It represents rather a tool for masking improper results than a real remedy. Another approach to particle approximation for density is the continuity density, which approximates the density according to the continuity equation using the concepts of SPH approximations plus some transformations. Different transformation on the RHS of the continuity equation may lead to different forms of density approximation equations, see [16]. The most frequently used form of continuity density equation and for our opinion the best one is obtained as N Dqi X @W ij mj v bij : ¼ Dt @xbi j¼1
between them will save during evolution. So, the important thing is to satisfy Eq. (5) in initial data. It can be done only by adjusting masses in (5) as being free parameters in this equation. In the given paper, an iteration algorithm is proposed for seeking the masses satisfied to Eq. (5). It is built by the following way. As an initial approach, the condition m0i ¼ 1 is taken, and then the iteration process starts:
ð7Þ
This equation shows clearly that the time rate of density change of a particle is closely related to the relative velocities v bij ¼ v bi v bj between ith particle and all the other particles in the support domain. Note, that when differentiating (5) over time, we obtain (7), if the smoothing lengths for all particles are not dependent on time. In this case, (5) and (7) are equivalent and respond to the same conservation law in algebraic and differential forms, correspondingly. They do not contradict to each other, if in initial distribution the density field satisfies to (5) exactly. If these conditions are not satisfied, the summation density approach and the continuity approach are actually different. As a whole, summation density approach is more popular in SPH applications. However, in applications devoted to high velocity impact, the continuity density approach is more attractive. In high velocity hydrodynamics the state variable is just density, so a continuity density approach is rather natural for this kind of problems. Besides, it is more efficient computationally. A disadvantage for the density summation approach is that it needs more computational effort since the density must be evaluated before other parameters can be evaluated, and it is required to calculate the smoothing function itself. The continuity density approach does not need to calculate density before other parameters, therefore can save considerable computational efforts. Many authors use the both approaches simultaneously. For instance, if continuity density was chosen as the main approach, every 50th or every 30th time step summation density is applied to deduce computational errors, see [16,22]. In our modeling, we also used the both approaches in equation for density. Ordinary differential equations were solved with predictor–corrector numerical scheme. In so doing, continuity approach (7) is applied at the predictor step, and the summation approach (5) at the corrector step. Density calculated with (5) depends on the positional relationship of the particles, so in such organized predictor–corrector scheme we need applying a time consuming search of neighbor particles only once, at the corrector stage. As has been mentioned above, expression (5) is simply the integral of (7), so at fulfillment of (5) in initial conditions, this relation
X msj W ij ; j
msþ1 ¼ msi i
qi : qsi
At the high velocity impact problem one can choose a set of some number of interacting objects, strikers and screens. Iteration process runs separately for every object were the particles are residing. Iterations go up to achievement of criterion qqsi 1 < e i
for all particles of the object, where could be taken e = 109. According to this iteration process the particles try occupying the equal volume, which is determined as the ratio of mass and density of the particle. Thus, particles take their equilibrium volumes initially prior to the beginning of time steps. As a result, during further time evolution the more uniform distribution of particles is attained. The necessity to arrange the particles at the objects as more uniform as possible was specified even in the early papers devoted to SPH method, [15]. By the construction, the given iteration method produces only the positive masses of particles. But it is necessary to take into account also, that particles can accept very small values of mass that it is undesirable. Therefore, particles with a mass smaller than some threshold value are deleted. The given inspection is constructed as follows. Obtained on s+1 iterations the masses of partimsþ1 cles are checked on requirement of their smallness msi < d, where i max msimax is the maximum value of a particle mass at the given object at the previous iteration. Value of threshold for rejection of a particle with small mass is chosen to be equal d = 104. The sense of this rejection is to avoid negative masses of the particles. If we solve linear algebraic system (5) with any exact method of linear algebra to find the masses of particles, we can obtain negative masses as a result. It indicates the bad particle adjustment inside the object. Just such bad positioned particles are excluded in this iteration algorithm. The given algorithm is built on the supposition on the constancy h over the object. For a constant smoothing length h, it is also possible to build algorithm of adding the particles in places of their greatest rarefaction, but in the given algorithm this possibility is not implemented. If, as usually done, to assign the masses of particles, and to set densities, the instability can further develop. The force rendered on a particle depends on the pressure, and it in turn depends on the density. If the computed field of density is close to equilibrium, pressures in the object are also near to equilibrium, that does not give to non-physical oscillations developing. One more advantage of the proposed algorithm consists in that the Monaghan velocity correction become unnecessary. The stability of particle movement can be achieved without this correction, which is rather artificial. Other SPH equation are written in the standard form:
! N X pj Dv ai p @W ij ; ¼ mj 2i þ 2 þ Pij Dt @xai q q i j j¼1 ! N pj Dei 1 X pi @W ij mj 2 þ 2 þ Pij v bij ; ¼ 2 j¼1 Dt qi qj @xbi Dxai ¼ v ai : Dt
ð8Þ ð9Þ ð10Þ
In Eqs. (8), (9) an artificial viscosity Pij is used to obtaining numerical stability in the regions of shock waves. It reads, see [16]:
6344
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
Fig. 1. Density distribution for several time moments in the first variant of calculation.
Pij ¼
/ij ¼
8 < aP cij /ij þbP /2ij : 0;
q ij
hij v ij xij jxij j2 þ u2
v ij ¼ v i v j ;
;
;
v ij xij < 0; v ij xij P 0;
cij ¼
1 ðci þ cj Þ; 2
ð11Þ
1 2
q ij ¼ ðqi þ qj Þ;
xij ¼ xi xj :
In 11, the constants aP, bP are the coefficients for linear and quadratic viscosity, c is the speed of sound, v is the particle velocity vector, the factor u = 0.1hij is inserted to prevent numerical instability when two particles are approaching each other. 4. Results We shall consider next two variant of calculations. In both the material of a shield and striker is aluminum. Spherical striker has the diameter D = 0.95 cm, the initial velocity u0 = 6.7 km/s, the thickness of plane shield D = 0.155 cm in the first variant. Initial
density of shield and striker is equal to q0. Geometry of collision in this variant is cylindrical. Geometry of impact is plane in the second variant, i.e. we study an impact of long rod with circular crosssection (its axes parallel to a shield plane) with velocity directed normally to a plane of shield. Diameter of rod, its initial velocity and thickness of shield are next: D = 1 cm, u0 = 6.18 km/s, D = 0.4 cm. A configuration of a first variant was taken realized in an experiment [21] and may be used as a test problem for checking the numerical model and calculation technique. Let’s discuss the results of calculations. The dynamics of flow in the first variant of calculation can trace on Fig. 1, where the field of density for several moments of time is shown. The Eulerian TVDmethod (in cylindrical r-z coordinates) in first variant of calculation is used. Calculation was performed on the uniform difference grid with constant cell sizes Dz = Dr = 0.015 cm. At the initial moment all the material objects (the striker and shields) are in equilibrium with the surrounding medium. Since the latter is under vacuum, the objects are in a state with a zero pressure. In this way the absence of motion inside each object is ensured. The shields are at rest relative to the laboratory coordinate system,
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
6345
Fig. 2. Material distribution (striker and shield) for several time moments in the first variant of calculation.
but at the same time the striker moves as a single whole toward them. In the initial time moment a striker there comes in contact with the shield – their boundary is located at z = 8.155 cm. In a place of impact the pressure quickly grows and reaches maximum in 0.3 mcs and maximal density reaches 4 g/cm3. There are two shock waves – one go deep into target, another propagates on striker’s body towards to its movement. To the moment 1 mcs the maximum of density makes 3.5 g/cm3. Striker goes deep into a target, the crater is formed and the material of striker spreads on a surface of a crater. On edge of a formed crater the jet going upwards is created. The shock wave, going on striker’s body, is increased quickly a specific energy of material. The second shock wave going downwards into material of shield gradually gets the half-spherical form. From this motion creates more one shock wave moves in a radial direction on a material of shield. After this striker begin expanded and its density is decreased. In a plate the
aperture is formed. In the moment t = 2 mcs the aperture radius makes 0.8 cm. The size of an aperture some time increases, and then its growth stops. Between edge of this aperture and striker’s body the thin shell (bubble) extending in radial and axial directions in process of moving of the striker is formed. In the subsequent the flow develops, the striker moves in a direction of impact and the sizes of a bubble grows. In the moment t = 7.2 mcs the forward edge of the striker has reached z = 3.7 cm, and the radial size of a bubble makes 1.8 cm. On the edge of generated aperture long whiskers from lateral side of shield were formed. Dynamics of material movement is given on Fig. 2, where the material of striker and shield are shown at the same time moments as in Fig. 1. It is necessary to notice that according to this drawing the thin shell of a bubble consists from a shield material. Comparison of result of calculation and experimental data [21] at the
6346
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
Fig. 3. Comparison of experimental radiograph [21] (a) and a calculation image (b) at the moment t = 7.2 mcs (first variant of calculation).
Fig. 4. Density distribution for several time moments in the second variant of calculation.
moment of time t = 7.2 mcs is presented on Fig. 3. As shows this comparison there is quite satisfactory agreement of computational
and experimental parameters of the striker, a bubble and whiskers (by position, form and sizes).
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
6347
Fig. 5. Comparison of calculation with SPH-method (left half) and TVD-method (right half) at time t = 5 mcs in the second variant of calculation (a – traditional SPH and b, c – modified SPH, a, b – density, c – material).
In the second variant of calculation the geometry of problem is plane. In this case two methods of calculation are used – Eulerian TVD and Lagrangian SPH. In SPH calculation was taken 80000 particles in shield and 7858 – in striker. The dynamics of flow in the second variant of calculation can trace on Fig. 4, where the field of density for several moment of time is shown (SPH-method of calculation). The picture of flow in this variant of calculation is qualitatively similar to the first variant, in this connection we will not discuss it in details. A comparison of calculation with SPHmethod and TVD-method is shown on Fig. 5 at time t = 5 mcs. This comparison shows the satisfactory agreement of the solutions received on the basis of given methods of calculation (modified SPH). At the same time, it is necessary to notice that in a traditional SPH-method [16] near to a leading edge of a formed bubble (in an
impact direction) on contact boundary between a substance of striker and shield the small crack is formed (Fig. 5a, left half). Such crack is not present in Eulerian method of calculation. It is connected, possibly, with some defect in a SPH-method. We will notice that in calculation of a problem about decay of arbitrary discontinuity (shock tube problem) there are some fluctuations of pressure near to contact boundary (practically in all modifications of SPHmethod discussed in the literature). 5. Conclusion The results of calculations have allowed one to receive detailed spatial-temporary picture of arising flow, to study creation, propagation and subsequent attenuation of shock waves in high-velocity
6348
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
impact problem. The calculation of two-dimensional gas dynamic problems on a basis of difference schemes in Eulerian and Lagrangian coordinates shows satisfactory quality of the solutions received on this technique. The numerical modeling is practically unique method for study of high-velocity impact with taking in account real properties of materials. The results of calculations allow estimate efficiency of anti-meteorite protection, to determine the shape and size of formed craters and apertures. Besides, the numerical modeling allows find a phase composition of destruction products, to study dynamics of formed jets etc. References [1] [2] [3] [4]
[5]
[6]
[7]
[8]
F.L. Whipple, Meteorites and space travel, Astron. J. 52 (1947) 137–148. J.A. Zukas, Impact Dynamics, John Wiley & Sons, New York, 1982. J.A. Zukas, High Velocity Impact, John Wiley & Sons, New York, 1990. L.D. Libersky, A.G. Petscheck, T.C. Carney, J.R. Hipp, F.A. Allahdadi, High strain Lagrangian hydrodynamics-a three-dimensional SPH code for dynamic material response, J. Comput. Phys. 109 (1993) 67–75. L.D. Libersky, R W. Randles, T.C. Carney, SPH calculations of fragmentation, in: Proceedings of the 3rd U.S. Congress on Computational Mechanics, Dallas, TX, 1995. P.W. Randles, L.D. Libersky, Smoothed particle hydrodynamics some recent improvements and applications, Comput. Methods Appl. Mech. Eng. 138 (1996) 375–408. P.W. Randles, T.C. Carney, L.D. Libersky, J.R. Renick, A.G. Petschek, Calculation of oblique impact and fracture of tungsten cubes using smoothed particle hydrodynamics, Int. J. Impact Eng. 17 (1995) 661–672. G.R. Johnson, E.H. Petersen, R.A. Stryk, Incorporation of an SPH option into the EPIC code for a wide range of high velocity impact computations, Int. J. Impact Eng. 14 (1993) 385–394.
[9] G.R. Johnson, R.A. Stryk, S.R. Beissel, SPH for high velocity impact computations, Comput. Methods Appl. Mech. Eng. 139 (1996) 347–373. [10] G.R. Johnson, Linking of Lagrangian particle methods to standard finite element methods for high velocity impact computations, Nucl. Eng. Design 150 (1994) 265–274. [11] E.P. Fahrenthold, B.A. Horban, An improved hybrid particle-element method for hypervelocity impact simulation, Int. J. Impact Eng. 26 (2001) 169–178. [12] O. Agertz, B. Moore, J. Stadel, et al., Fundamental Differences Between SPH and Grid Methods 1 (1996) 15.. [13] P.L. Roe, Approximate Riemann solvers, parameter vectors, and difference schemes, J. Comput. Phys. 43 (1981) 357–372. [14] H.C. Yee, Linearized form of implicit TVD schemes for multidimensional Euler and Navier–Stokes equations, Int. J. Comput. Math. Appl. 12A (4/5) (1986) 413–432. [15] J.J. Monaghan, Smoothed particle hydrodynamics, Ann. Rev. Astron. Astrophys. 30 (1992) 543–574. [16] G.R. Liu, M.B. Liu, Smoothed Particle Hydrodynamics: A Meshfree Particle Method, World Scientific Pub., Singapore, 2003. 449 p. [17] J.H. Tillotson Metallic equation of state for hypervelocity impact, General Atomic Report GA-3216, 1962, 142 p. [18] J.E. Dienes, J.M. Walsh, in: R. Kinslow (Ed.), High-Velocity Impact Phenomena, Academic Press, New York–London, 1970, p. 456. [19] J. Moon, H.C. Yee, Numerical simulation by TVD schemes of complex shock reflections from airfoils at high angle of attack, AIAA Paper, No. 350, 1987, pp. 1–17. [20] H.C. Yee, Linearized form of implicit TVD schemes for multidimensional Euler and Navier–Stokes equations, Int. J. Comput. Math. Appl. 12A (4/5) (1986) 413–432. [21] A.J. Piekutowski, Characteristics of debris clouds produced by hypervelocity impact of aluminum spheres with thin aluminum plates, Int. J. Impact Eng. 14 (1993) 573–586. [22] M.G. Gesteira, B.D. Rogers, R.A. Dalrymple, et al., User Guide for the SPHysics Code, Johns Hopkins University, 2009.
Contents lists available at SciVerse ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Development of modified SPH approach for modeling of high-velocity impact K.A. Alhussan a, V.A. Babenko b, I.M. Kozlov b, A.S. Smetannikov b,⇑ a b
Space Research Institute, King Abdulaziz City for Science and Technology, Riyadh, Saudi Arabia Heat and Mass Transfer Institute, Brovki Str., 15, 220072 Minsk, Belarus
a r t i c l e
i n f o
Article history: Received 29 February 2012 Accepted 24 May 2012 Available online 13 July 2012 Keywords: High-velocity impact Shock waves Hydrodynamics Numerical methods
a b s t r a c t The problem of high-velocity impact was modeled on the basis of two numerical methods, the smoothed particle hydrodynamics (SPH) and the Eulerian total variation diminishing (TVD) finite-difference scheme. The modification of SPH method was elaborated. In this modification, the summation approach and the continuous approach are equivalent, and the problem of boundary particle deficiency is successfully resolved. The comparison of modified SPH and Eulerian finite-difference approaches shows a good agreement between them. The calculations results can be useful for estimation of efficiency of anti-meteorite protection shields, giving information on the shape and size of formed apertures and jets. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction The study of dynamics of the processes, which take place at high-velocity impact (HVI), is of interest for many problems of space physics and astrophysics: the study of meteoric craters, the origin of planetary atmospheres, and the consequences of large space object fall on the Earth etc., [1]. The long stay of space vehicles at an orbit results in practically inevitable influence of meteorites on their surface. Thus, the probability of collision with small particles (micrometeorites) with mass less than 0.01 g is especially great. At high velocities (10 km/s) the initial pressure near to impact place makes about megabar, and the specific energy exceeds the sublimation heat. At such action the material melts and evaporates. The time of action of high pressure appears to be rather short (owing to the small size of a meteorite), the shock wave quickly damps. At moving away from the impact place the specific energy of meteorite decreases, to become possible only melting. Then it falls below the melting energy, and at this stage there can be only fracturing of material. This process is described by creation and growth of microcracks that is determined by strength characteristics of material. We shall consider an initial stage of impact, when strength effects are insignificant. The experimental data on HVI concern to the velocities not exceeding 10 km/s. Thereof, in the range of high collision velocities characteristic of anti-meteorite protection the theoretical research of impact dynamics (in particular numerical modeling) is very important. ⇑ Corresponding author. E-mail address: [email protected] (A.S. Smetannikov). 0017-9310/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2012.06.002
The majority of traditional approaches for numerical modeling in hydrodynamics use a computational mesh for calculation. The Eulerian approaches study the gas motion using space mesh, which is not changed in time (Eulerian mesh). The Lagrangian approaches use a mesh moving together with gas. These mesh based techniques solve the fluid dynamic equations by calculating interaction with adjacent cells. The Eulerian description is a spatial description, whereas the Lagrangian description is a material description. The fundamental difference of these two descriptions is that the Lagrangian description employs the total time derivative as the combination of local derivative and convective derivative. Methods of particles refer to a class of Lagrangian methods. The media in these methods is represented by a set of particles, and movement of these particles is investigated. Among such particle methods the method of smoothed particle hydrodynamics (SPH) holds a special place because it is mesh-free. In the SPH method the equation of gas dynamics are solved on the basis of interaction of each particle with its nearest neighbors. Due to Lagrangian description, the problem can be reduced to solution of ordinary differential equations, each of them being formulated for the velocity, density, position, and energy of a particle. Many of studies of HVI in the past were performed with Eulerian approach. Nevertheless, this class of numerical models has unavoidable difficulties as applied to HVI investigation. In this problem there are solid objects moving with high velocities and having, at least initially, well-defined boundaries. We are interested mainly on transformation of these boundaries with time. But these boundaries can slide one over another with very high velocities, and in frame of Eulerian approaches it is very difficult getting answer on the real position of the material
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
interfaces. Besides, this sliding creates very serious problems with the numerical stability. These difficulties can be resolved and overcome to some extent by using complicated Eulerian methods such as TVD schemes with markers describing positions of the interfaces. But at the later stages of impact when the topology of created cracks can be very complex, it is difficult to apply this technique. That is why Lagrangian approaches, and especially SPH method have received wide acceptance in HVI studies. SPH is a new effective technique for impact simulation, which combines harmonically the benefits of Lagrangian methods and particle approximation. The books by Zukas [2,3], the papers by Libersky with co-authors [4,5], by Randles et al. [6,7], and by Johnson et al. [8,9] have made outstanding contributions in the application of SPH to impact problems. Many papers were made afterwards in order to combine SPH with finite element methods (FEM). Different contact algorithms have been used with SPH when combining the SPH method with the FEM [10], but the only purpose of these algorithms was to allow SPH nodes to interact with FEM mesh. An improved hybrid particle-finite element method has been developed for hypervelocity impact simulation in [11]. The SPH method possesses some advantages and disadvantages over traditional methods of numerical prediction. The main privileges of SPH consist in that it is remarkably robustness in sense of numerical stability and relatively simplicity of programming. The main disadvantage is probably insufficient accuracy if the number of particles in modeling is not great enough. Just the number of particles is a key problem in SHP. An amount of computational work in reported SPH programs grows with the number of particles N as the factor between N2 and N log N. That is why even dealing with optimized computer problem, user should distribute the particles very economically. Actually, the number of particles is usually bounded by the value of several hundred of thousands or a few million. So, at actual computations the value of smoothing length h can be large enough as the total amount of particles is bounded by computational time. Even for this restricted number of particles, it is important keeping balances of mass, momentum and energy exactly for the system as whole and for its component parts, as it actually became a standard in the theory of finite-difference schemes. In SPH, however, writing fully conservative numerical algorithms is rather problematic by many reasons. It will be explained later after formulation of main equations in SPH. Here, it worth noting only that the density approximation is especially important in this connection, since the density basically determines the particle distribution and the smoothing length evolution. The error arising in continuity equation is transferred to the momentum and energy equations thus worsened the modeling quality. Another important problem in SPH is a so called boundary particle deficiency. Near the boundaries, the support domain for a particle with smoothing length h can intersect the whole domain of the problem, thus reducing smoothed density of the particle and leading to erroneous results. The authors of [12], for instance, showed that SPH introduces spurious pressure forces acting on particles in regions where there are steep density gradients. Due to reduced density the pressure at the particle can increase substantially thus creating a source of numerical instability. Boundary particle deficiency can contribute also into the error in continuity equation. In this paper, a modification of SPH methods is proposed and investigated in application to HVI problem. Comparison is made with the results of application of the TVD difference scheme with markers. Suggested in this paper an improvement of SPH enables to keep the mass balance exactly and to eliminate the boundary particle deficiency problem.
6341
2. Governing equations For the Eulerian description of axial symmetrical flow arising at impact, the system of gas dynamic equations in two-dimensional representation reads:
@u @u @u @P @v @v @v @P þv þu ¼ ; q þv þu ¼ ; @t @r @z @z @r @t @r @z 1 dq 1 @ðr n v Þ @u de 1 @ðr n v Þ @u þ ; q ¼ P n þ ; P ¼ Pðe; qÞ; ¼ n r @r @z r @r @z q dt dt
q
ð1Þ
where t is the time, q the density, P the pressure, e the specific internal energy, v, u are the radial and axial components of velocity vector, n is the index of symmetry (n = 0 for plane and n = 1 for cylindrical geometry, correspondingly). The equation of state P = P(e, q) closes this system of equations. The complexity of numerical calculation of high-velocity impact is caused by prominent features of arising flow. At the initial stage of impact the flow is sharply spatially non-uniform, and there are areas of materials sliding. Prediction of such flows in Lagrange representation is extremely complicated. For later times, the other type of complexity is connected to considerable growth of spatial scale of flow. The system of gasdynamic Equation (1) should be approximated by some difference scheme for numerical solution. For this, we use two difference schemes. The first one is a method of total variation diminishing (TVD) in Eulerian variables [13,14], the second one the smoothed particle hydrodynamics (SPH) method in Lagrangian variables [15,16]. In the first method, the spatial mesh is used, the second method is meshfree. More detailed description of these methods will be made in the next section. The calculations are carried out for impact of a striker (meteorite) on plane aluminum shields (with thickness D) with speed u0, directed normally to the shield. The meteorite is simulated by a sphere with diameter D. As a material of striker, aluminum is taken also. For closing the system of gas dynamic equations it is necessary to determine the equation of state for materials, i.e. the dependence of pressure P = P(e, q) from energy and density. In calculations we used the Mie–Gruneisen equation of state in form proposed by Tillotson [17]. For material density higher than the normal density of a solid state q P q0 for all specific energy e (and in small region with q < q0 and e < e1) we assumed the relation
" P ¼ Pc ¼ a þ
#
b 1 þ e=ðe0 d2 Þ
qe þ Al þ Bl2 :
ð2Þ
With density lower than that a solid state q 6 q0 and e > e2 we used the equation P ¼ P H ¼ aqe þ
" ! 2 # 1 1 exp a 1 þ Al exp b 1 : d d 1 þ e=ðe0 d Þ b qe
2
ð3Þ
In these expressions, d = q/q0, l = d 1, and e1, e2 are the energies of incipient and full vaporization. For density smaller than that for solid state q 6 q0 and energy e1 < e < e2 we used next interpolation by energy between these two relations
P ¼ ½Pc ðe2 eÞ þ PH ðe e1 Þ=ðe2 e1 Þ:
ð4Þ
In calculations, the following parameters of aluminum were taken: a = 0.5, b = 1.63, e0 = 5 MJ/kg, A = 75 GPa, B = 65 GPa, q0 = 2790 kg/m3, a = 5, b = 5, e1 = 3 MJ/kg, e2 = 15 MJ/kg. It should be noted that in the unloading region (at densities somewhat smaller than that for solid state q0) at low energies the Tillotson equations of state gives very large negative pressures. In fact, when the
6342
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
negative pressure exceeds a magnitude of the order of 0.1 GPa, there occurs a breaking (or splitting off) the material, and the pressure goes to zero. In this connection, when it turned out that the pressure was negative, we assumed that p = 0. A similar procedure in this region was employed in [18]. 3. Technique of calculation In this section we will describe shortly the techniques used in the present work for calculation of flow hydrodynamics. 3.1. TVD method The technique of numerical solution of the system of Euler equations in application to the problem of the frontal impact of a spherical striker on a plane shield of finite thickness in vacuum is considered. A characteristic feature of the method is that the Eulerian approach for solving the gas dynamics equations is combined here with the description of the motion of contact boundaries between various objects. The impermeable contact boundaries are described with the aid of marker particles moving together with the substance. The application of markers imparts to the Euler scheme some features typical of special importance in modeling collision of solid bodies with distinct boundaries. Moreover, information about contact boundaries makes it possible to easily identify various objects in space and, consequently, to appropriately describe their thermodynamic and strength properties. The setting of boundary conditions at interfaces between different media allows one to describe the interaction of the objects that possess the properties of both gases and solids. A numerical scheme is realized on rectangular spatial grids and can be used in modeling two-dimensional problems possessing an axial or plane symmetry. As a difference scheme, use is made of the scheme of the second order of accuracy in space obtained by the method of finite volumes. The second order accuracy in time is provided by ‘‘predictor–corrector’’-type scheme. The scheme does not employ an artificial viscosity, with the stability of solution being ensured by application of the numerical TVD dissipation that does not alter the second order of approximation. To calculate the TVD correction of fluxes, an approach based on the use of the Riemann invariants is applied making it possible to be got rid of time-consuming matrix operations. For each characteristic solution the correcting fluxes are calculated with the aid of the formula for the scalar upwind TVD scheme [19,20]. In order to describe contact boundaries with the aid of markers in the two-dimensional case, the notion of a marker object is introduced. The marker object is understood to be a certain region of space entirely occupied by a medium of any one species and bounded by one or several marker lines. The marker line is a continuous curve formed on a plane by segments of a straight line that connect a set of markers ordered in space. A straight-line segment that connects two successive markers is called a marker segment. Each marker line unites its own sequence of markers and has a certain orientation starting with the first marker and ending with the last marker of the sequence. The marker is a point on a plane assigned by a pair of coordinates and having a velocity vector. The velocity of a marker is determined by interpolation of velocity components from neighboring nodes of the grid. In the ‘‘predictor– corrector’’-type of the scheme the velocity field for markers is calculated at the ‘‘predictor’’ stage. In addition to markers, configuration and index arrays are used for describing a marker object. The former serves for determining which specific computational cell belongs to that or other marker object. The latter establishes the connection between the marker object and the material medium that fills this object. For each cell crossed by one or several marker
lines the mean direction of these lines inside a cell is calculated and is represented in the form of a unit vector n. In a two-dimensional case, when pulse exchange through the contact is calculated, it is possible to take into account the true calculation of the marker line relative to the coordinate lines of the rectangular grid. For this purpose, a unit vector n is assigned to each cell crossed by the marker line. Moreover, taking into account the direction n of the marker line, one can model the slipping between the media. The numerical description of the contact boundary is composed of two parts: 1. Determination of the position of the contact boundary in space. 2. Determination of the fluxes of gas-dynamical quantities through the contact boundary. The contact boundary corresponding to the marker line passes only along the boundaries of the cells, i.e. each cell entirely belongs to only one marker object or background. Whether a concrete cell, crossed by one or several marker lines, belongs to that or other marker object or background is determined by comparing the areas occupied by different objects inside a cell. A cell is ascribed to the marker object (or background) which has the largest area. The position of the contact boundary changes discretely accurate to the magnitude of a spatial cell. The status of a substance in the cell through which the marker line passes, also changes discretely. If at a certain moment the status of a contact cell changes, an analysis of its neighbors is carried out, on the basis of which the given cell is assigned some parameters average over the neighboring cells. 3.2. SPH method SPH is a numerical method that describes a fluid flow in continuous media as a motion of the collective of interacting particles. For each particle i its characteristics such as mass mj, position ri, velocity vi, density qi and energy density ei are determined. The modification of SPH method in application to HVI, which is developed in the given paper, concerns mainly to the SPH approximation of mass conservation law. In literature, various SPH approximations for density equation are described (see [16]). The first approach is the summation density (5), which directly applies the SPH approximations to the density itself. For a given particle i, the density with summation density approach is written in the form of
qi ¼
N X mj W ij ;
ð5Þ
j¼1
where N is the number of particles in the support domain of particle i, and mj is the mass associated with particle j. Wij is the smoothing function of particle i evaluated at particle j with the smoothing length h.
W ij ¼ W jxi xj j; h ¼ Wðrij ; hÞ; where Wij = W(rij, h) is the kernel function with smoothing length h = (hi + hj)/2, where hi denotes the smoothing length associated with the ith particle and rij = jxi xjj is the distance between particles i and j. As a main kernel function we use Gauss form. The Gauss kernel is actually the ideal function for the purpose, but it takes computational time approximately half as much again as the cubic spline kernel. The summation density approach (5) is, at a first glance, the most natural. It conserves the mass exactly since the integration of the density over the entire problem domain is exactly the total mass of all the particles. Nevertheless, in so doing, the accuracy of density approximation is lost near the boundaries, where the approximation becomes poor due to the mentioned above particle
6343
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
deficiency. To improve the situation a modified form of summation density approach in form of relation (6) was proposed by many authors. It consists in normalization of the RHS of Eq. (5) with the SPH summation of the smoothing function over the neighboring particles:
PN
mj W ij
qi ¼ Pj¼1 N mj j¼1 qj
W ij
:
ð6Þ
qsi ¼
Here, the correction of (5) is made, which takes into account a non equality of the denominator to unity. The possibility of such correction at improving the SPH computation results is not very impressive. It represents rather a tool for masking improper results than a real remedy. Another approach to particle approximation for density is the continuity density, which approximates the density according to the continuity equation using the concepts of SPH approximations plus some transformations. Different transformation on the RHS of the continuity equation may lead to different forms of density approximation equations, see [16]. The most frequently used form of continuity density equation and for our opinion the best one is obtained as N Dqi X @W ij mj v bij : ¼ Dt @xbi j¼1
between them will save during evolution. So, the important thing is to satisfy Eq. (5) in initial data. It can be done only by adjusting masses in (5) as being free parameters in this equation. In the given paper, an iteration algorithm is proposed for seeking the masses satisfied to Eq. (5). It is built by the following way. As an initial approach, the condition m0i ¼ 1 is taken, and then the iteration process starts:
ð7Þ
This equation shows clearly that the time rate of density change of a particle is closely related to the relative velocities v bij ¼ v bi v bj between ith particle and all the other particles in the support domain. Note, that when differentiating (5) over time, we obtain (7), if the smoothing lengths for all particles are not dependent on time. In this case, (5) and (7) are equivalent and respond to the same conservation law in algebraic and differential forms, correspondingly. They do not contradict to each other, if in initial distribution the density field satisfies to (5) exactly. If these conditions are not satisfied, the summation density approach and the continuity approach are actually different. As a whole, summation density approach is more popular in SPH applications. However, in applications devoted to high velocity impact, the continuity density approach is more attractive. In high velocity hydrodynamics the state variable is just density, so a continuity density approach is rather natural for this kind of problems. Besides, it is more efficient computationally. A disadvantage for the density summation approach is that it needs more computational effort since the density must be evaluated before other parameters can be evaluated, and it is required to calculate the smoothing function itself. The continuity density approach does not need to calculate density before other parameters, therefore can save considerable computational efforts. Many authors use the both approaches simultaneously. For instance, if continuity density was chosen as the main approach, every 50th or every 30th time step summation density is applied to deduce computational errors, see [16,22]. In our modeling, we also used the both approaches in equation for density. Ordinary differential equations were solved with predictor–corrector numerical scheme. In so doing, continuity approach (7) is applied at the predictor step, and the summation approach (5) at the corrector step. Density calculated with (5) depends on the positional relationship of the particles, so in such organized predictor–corrector scheme we need applying a time consuming search of neighbor particles only once, at the corrector stage. As has been mentioned above, expression (5) is simply the integral of (7), so at fulfillment of (5) in initial conditions, this relation
X msj W ij ; j
msþ1 ¼ msi i
qi : qsi
At the high velocity impact problem one can choose a set of some number of interacting objects, strikers and screens. Iteration process runs separately for every object were the particles are residing. Iterations go up to achievement of criterion qqsi 1 < e i
for all particles of the object, where could be taken e = 109. According to this iteration process the particles try occupying the equal volume, which is determined as the ratio of mass and density of the particle. Thus, particles take their equilibrium volumes initially prior to the beginning of time steps. As a result, during further time evolution the more uniform distribution of particles is attained. The necessity to arrange the particles at the objects as more uniform as possible was specified even in the early papers devoted to SPH method, [15]. By the construction, the given iteration method produces only the positive masses of particles. But it is necessary to take into account also, that particles can accept very small values of mass that it is undesirable. Therefore, particles with a mass smaller than some threshold value are deleted. The given inspection is constructed as follows. Obtained on s+1 iterations the masses of partimsþ1 cles are checked on requirement of their smallness msi < d, where i max msimax is the maximum value of a particle mass at the given object at the previous iteration. Value of threshold for rejection of a particle with small mass is chosen to be equal d = 104. The sense of this rejection is to avoid negative masses of the particles. If we solve linear algebraic system (5) with any exact method of linear algebra to find the masses of particles, we can obtain negative masses as a result. It indicates the bad particle adjustment inside the object. Just such bad positioned particles are excluded in this iteration algorithm. The given algorithm is built on the supposition on the constancy h over the object. For a constant smoothing length h, it is also possible to build algorithm of adding the particles in places of their greatest rarefaction, but in the given algorithm this possibility is not implemented. If, as usually done, to assign the masses of particles, and to set densities, the instability can further develop. The force rendered on a particle depends on the pressure, and it in turn depends on the density. If the computed field of density is close to equilibrium, pressures in the object are also near to equilibrium, that does not give to non-physical oscillations developing. One more advantage of the proposed algorithm consists in that the Monaghan velocity correction become unnecessary. The stability of particle movement can be achieved without this correction, which is rather artificial. Other SPH equation are written in the standard form:
! N X pj Dv ai p @W ij ; ¼ mj 2i þ 2 þ Pij Dt @xai q q i j j¼1 ! N pj Dei 1 X pi @W ij mj 2 þ 2 þ Pij v bij ; ¼ 2 j¼1 Dt qi qj @xbi Dxai ¼ v ai : Dt
ð8Þ ð9Þ ð10Þ
In Eqs. (8), (9) an artificial viscosity Pij is used to obtaining numerical stability in the regions of shock waves. It reads, see [16]:
6344
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
Fig. 1. Density distribution for several time moments in the first variant of calculation.
Pij ¼
/ij ¼
8 < aP cij /ij þbP /2ij : 0;
q ij
hij v ij xij jxij j2 þ u2
v ij ¼ v i v j ;
;
;
v ij xij < 0; v ij xij P 0;
cij ¼
1 ðci þ cj Þ; 2
ð11Þ
1 2
q ij ¼ ðqi þ qj Þ;
xij ¼ xi xj :
In 11, the constants aP, bP are the coefficients for linear and quadratic viscosity, c is the speed of sound, v is the particle velocity vector, the factor u = 0.1hij is inserted to prevent numerical instability when two particles are approaching each other. 4. Results We shall consider next two variant of calculations. In both the material of a shield and striker is aluminum. Spherical striker has the diameter D = 0.95 cm, the initial velocity u0 = 6.7 km/s, the thickness of plane shield D = 0.155 cm in the first variant. Initial
density of shield and striker is equal to q0. Geometry of collision in this variant is cylindrical. Geometry of impact is plane in the second variant, i.e. we study an impact of long rod with circular crosssection (its axes parallel to a shield plane) with velocity directed normally to a plane of shield. Diameter of rod, its initial velocity and thickness of shield are next: D = 1 cm, u0 = 6.18 km/s, D = 0.4 cm. A configuration of a first variant was taken realized in an experiment [21] and may be used as a test problem for checking the numerical model and calculation technique. Let’s discuss the results of calculations. The dynamics of flow in the first variant of calculation can trace on Fig. 1, where the field of density for several moments of time is shown. The Eulerian TVDmethod (in cylindrical r-z coordinates) in first variant of calculation is used. Calculation was performed on the uniform difference grid with constant cell sizes Dz = Dr = 0.015 cm. At the initial moment all the material objects (the striker and shields) are in equilibrium with the surrounding medium. Since the latter is under vacuum, the objects are in a state with a zero pressure. In this way the absence of motion inside each object is ensured. The shields are at rest relative to the laboratory coordinate system,
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
6345
Fig. 2. Material distribution (striker and shield) for several time moments in the first variant of calculation.
but at the same time the striker moves as a single whole toward them. In the initial time moment a striker there comes in contact with the shield – their boundary is located at z = 8.155 cm. In a place of impact the pressure quickly grows and reaches maximum in 0.3 mcs and maximal density reaches 4 g/cm3. There are two shock waves – one go deep into target, another propagates on striker’s body towards to its movement. To the moment 1 mcs the maximum of density makes 3.5 g/cm3. Striker goes deep into a target, the crater is formed and the material of striker spreads on a surface of a crater. On edge of a formed crater the jet going upwards is created. The shock wave, going on striker’s body, is increased quickly a specific energy of material. The second shock wave going downwards into material of shield gradually gets the half-spherical form. From this motion creates more one shock wave moves in a radial direction on a material of shield. After this striker begin expanded and its density is decreased. In a plate the
aperture is formed. In the moment t = 2 mcs the aperture radius makes 0.8 cm. The size of an aperture some time increases, and then its growth stops. Between edge of this aperture and striker’s body the thin shell (bubble) extending in radial and axial directions in process of moving of the striker is formed. In the subsequent the flow develops, the striker moves in a direction of impact and the sizes of a bubble grows. In the moment t = 7.2 mcs the forward edge of the striker has reached z = 3.7 cm, and the radial size of a bubble makes 1.8 cm. On the edge of generated aperture long whiskers from lateral side of shield were formed. Dynamics of material movement is given on Fig. 2, where the material of striker and shield are shown at the same time moments as in Fig. 1. It is necessary to notice that according to this drawing the thin shell of a bubble consists from a shield material. Comparison of result of calculation and experimental data [21] at the
6346
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
Fig. 3. Comparison of experimental radiograph [21] (a) and a calculation image (b) at the moment t = 7.2 mcs (first variant of calculation).
Fig. 4. Density distribution for several time moments in the second variant of calculation.
moment of time t = 7.2 mcs is presented on Fig. 3. As shows this comparison there is quite satisfactory agreement of computational
and experimental parameters of the striker, a bubble and whiskers (by position, form and sizes).
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
6347
Fig. 5. Comparison of calculation with SPH-method (left half) and TVD-method (right half) at time t = 5 mcs in the second variant of calculation (a – traditional SPH and b, c – modified SPH, a, b – density, c – material).
In the second variant of calculation the geometry of problem is plane. In this case two methods of calculation are used – Eulerian TVD and Lagrangian SPH. In SPH calculation was taken 80000 particles in shield and 7858 – in striker. The dynamics of flow in the second variant of calculation can trace on Fig. 4, where the field of density for several moment of time is shown (SPH-method of calculation). The picture of flow in this variant of calculation is qualitatively similar to the first variant, in this connection we will not discuss it in details. A comparison of calculation with SPHmethod and TVD-method is shown on Fig. 5 at time t = 5 mcs. This comparison shows the satisfactory agreement of the solutions received on the basis of given methods of calculation (modified SPH). At the same time, it is necessary to notice that in a traditional SPH-method [16] near to a leading edge of a formed bubble (in an
impact direction) on contact boundary between a substance of striker and shield the small crack is formed (Fig. 5a, left half). Such crack is not present in Eulerian method of calculation. It is connected, possibly, with some defect in a SPH-method. We will notice that in calculation of a problem about decay of arbitrary discontinuity (shock tube problem) there are some fluctuations of pressure near to contact boundary (practically in all modifications of SPHmethod discussed in the literature). 5. Conclusion The results of calculations have allowed one to receive detailed spatial-temporary picture of arising flow, to study creation, propagation and subsequent attenuation of shock waves in high-velocity
6348
K.A. Alhussan et al. / International Journal of Heat and Mass Transfer 55 (2012) 6340–6348
impact problem. The calculation of two-dimensional gas dynamic problems on a basis of difference schemes in Eulerian and Lagrangian coordinates shows satisfactory quality of the solutions received on this technique. The numerical modeling is practically unique method for study of high-velocity impact with taking in account real properties of materials. The results of calculations allow estimate efficiency of anti-meteorite protection, to determine the shape and size of formed craters and apertures. Besides, the numerical modeling allows find a phase composition of destruction products, to study dynamics of formed jets etc. References [1] [2] [3] [4]
[5]
[6]
[7]
[8]
F.L. Whipple, Meteorites and space travel, Astron. J. 52 (1947) 137–148. J.A. Zukas, Impact Dynamics, John Wiley & Sons, New York, 1982. J.A. Zukas, High Velocity Impact, John Wiley & Sons, New York, 1990. L.D. Libersky, A.G. Petscheck, T.C. Carney, J.R. Hipp, F.A. Allahdadi, High strain Lagrangian hydrodynamics-a three-dimensional SPH code for dynamic material response, J. Comput. Phys. 109 (1993) 67–75. L.D. Libersky, R W. Randles, T.C. Carney, SPH calculations of fragmentation, in: Proceedings of the 3rd U.S. Congress on Computational Mechanics, Dallas, TX, 1995. P.W. Randles, L.D. Libersky, Smoothed particle hydrodynamics some recent improvements and applications, Comput. Methods Appl. Mech. Eng. 138 (1996) 375–408. P.W. Randles, T.C. Carney, L.D. Libersky, J.R. Renick, A.G. Petschek, Calculation of oblique impact and fracture of tungsten cubes using smoothed particle hydrodynamics, Int. J. Impact Eng. 17 (1995) 661–672. G.R. Johnson, E.H. Petersen, R.A. Stryk, Incorporation of an SPH option into the EPIC code for a wide range of high velocity impact computations, Int. J. Impact Eng. 14 (1993) 385–394.
[9] G.R. Johnson, R.A. Stryk, S.R. Beissel, SPH for high velocity impact computations, Comput. Methods Appl. Mech. Eng. 139 (1996) 347–373. [10] G.R. Johnson, Linking of Lagrangian particle methods to standard finite element methods for high velocity impact computations, Nucl. Eng. Design 150 (1994) 265–274. [11] E.P. Fahrenthold, B.A. Horban, An improved hybrid particle-element method for hypervelocity impact simulation, Int. J. Impact Eng. 26 (2001) 169–178. [12] O. Agertz, B. Moore, J. Stadel, et al., Fundamental Differences Between SPH and Grid Methods 1 (1996) 15.