High velocity impact of metal sphere on thin metallic plate using smooth particle hydrodynamics (SPH) method

Abstract

The modeling of high velocity impact is an important topic in impact engineering. Due to various constraints, experimental data are extremely limited. Therefore, detailed numerical simulation can be considered as a desirable alternative. However, the physical processes involved in the impact are very sophisticated; hence a practical and complete reproduction of the phenomena involves complicated numerical models. In this paper, we present a smoothed particle hydrodynamics (SPH) method to model two-dimensional impact of metal sphere on thin metallic plate. The simulations are applied to different materials (Aluminum, Lead and Steel); however the target and projectile are formed of similar metals. A wide range of velocities (300, 1000, 2000, and 3100 m/s) are considered in this study. The goal is to study the most sensitive input parameters (impact velocity and plate thickness) on the longitudinal extension of the projectile, penetration depth and damage crater.

Keywords

smoothed particle hydrodynamics / high velocity impact / sensitivity analysis

Introduction

Impact in mechanical engineering problems can be defined as a high force or shock which is applied over a short time period. Modeling and simulating the collision of two or more bodies can help predict the consequences of the impact and protect the surrounding area against these effects.

Since experimental studies of the behavior of solid material or fluid are conducted with very expensive devices and equipment, computational modeling and numerical simulation have become the proper and accepted ways to investigate behavior of materials. Meshfree methods [1–3] are a powerful alternative to finite element methods to model large deformations, dynamic fracture and fragmentation.

SPH method is a fully Lagrangian, non-grid based computational technique that is typically used for continuum dynamics simulations. SPH was mainly used to simulate the motion of compressible flow. SPH has recently been extended to simulate incompressible fluid motions. In SPH modeling, particles of the fluid can be modeled as spheres, and the interacting forces between them can be calculated by applying dynamics and fluid dynamics formulas like gravity effect, inertia effect, viscosity effect, and pressure effect. Smoothed Particle Hydrodynamics method (SPH) was first introduced in the 1970s to model astrophysical applications [4,5], but recently, SPH can be found in a wide range of research areas. SPH method plays a key role in modeling of engineering problems including heat conduction [6], multi-phase flows [7], chemical explosions [8], as well as deformation and impact problems [9]. The SPH method can overcome the difficulties in large deformations, large inhomogeneities, and large discontinuities as well as in material interface tracking and therefore it is a good alternative to traditional grid-based numerical methods in simulating impact problems.

In the following, first a brief summary of the SPH method will be presented. The difficulty in modeling the primary steps, in particular the number of particles used for initial particles placement, is then described. The effects of initial velocity variations for three types of materials (aluminum, lead and steel) are then described. For the validation of our SPH code, our research is compared with previously published results. Finally a sensitivity analysis of the output variables is conducted with respect to the initial projectile velocity and plate thickness.

SPH methodology and numerical aspects

Governing equations

In SPH method, the fluid is represented by particles. In this technique the interpolation function A is defined as

(1)A_k(r)={\displaystyle \underset{\upsilon }{\int }A\left(r^{\prime }\right)}w\left(r-r^{\prime },h\right)\mathrm{d}{r}^{\prime },

where the integral is over the domain í and w(r - r′, h) signifies the smoothing kernel function. r and r′ are the position vectors at different points and h, the smoothing length, is the effective width of smoothing.

When using interpolation for a model, it is necessary to divide it into a set of small elements which are represented as mass elements. For instance, the properties of element j will be expressed in terms of m_j (mass of element j), p_j (density of element j), and r_j (position of element j). The value of parameter A at particle j is denoted by A_j.

Moreover, in numerical techniques the integral of interpolation function is generally approximated with a weight interpolating function. The integral can be estimated by summing up the mass elements. The summation interpolant can be written as

(2)A_s(r)={\displaystyle \sum _j{m}_j{\displaystyle \frac{A_j}{{\rho }_j}}w(r-r_j,h)}.

One of the most well-known kernels used in SPH is based on the cubic spline function.

(3)

where s={\displaystyle \frac{\left|r-r^{\prime }\right|}{h}}, \alpha is a normalized constant and v is the dimension of the problem. \alpha is equal to {\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}, {\displaystyle \frac{\mathrm{10}}{\mathrm{7}\mathrm{\pi }}}, \mathrm{1}/\mathrm{\pi } for one, two, and three dimension respectively.

The second derivative of the above kernel is continuous, and the dominant error term in the integral of interpolation function is o(h²). The compact support length of this kernel is 2h which means that there is no interaction between the particles in distance beyond 2h. Figure 1 shows how a kernel function acts on its compact support.

The two-dimensional SPH expressions for equations of continuum mechanics, conservation of mass, momentum and energy, are given by

(4){\displaystyle \frac{\mathrm{d}{\rho }_i}{\mathrm{d}t}}=\rho {\displaystyle \sum _j{\displaystyle \frac{m_j}{{\rho }_j}}(V_i-V_j)}\cdot \nabla w_{ij},

(5){\displaystyle \frac{\mathrm{d}{\mathit{V}}_{\mathit{t}}}{\mathrm{d}\mathit{t}}}={\displaystyle \sum _j{m}_j}\left({\displaystyle \frac{{\sigma }_i}{p_i^{\mathrm{2}}}}+{\displaystyle \frac{{\sigma }_j}{p_j^{\mathrm{2}}}}-{\mathrm{\Pi }}_{ij}\right).\nabla w_{ij},

(6){\displaystyle \frac{\mathrm{d}{\mathit{e}}_{\mathit{i}}}{\mathrm{d}\mathit{t}}}=-{\displaystyle \sum _j{m}_j}(V_i-V_j)\left({\displaystyle \frac{{\sigma }_i}{{\mathit{\rho }}_i^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathrm{\Pi }}_{ij}\right).\nabla w_{ij}.

In above equations {\rho }_i, m_i, V_i, {\sigma }_i, and e_i are the density, mass, velocity, stress and specific internal energy of the ith particle respectively. {\mathrm{\Pi }}_{ij} denotes the artificial viscosity which is introduced to capture the shock wave. For more comprehensive details about effects caused by variation of artificial viscosity, one may refer to [10].

In order to reduce the numerical noise (numerical error) in simulation, some spatial filtering practices are often employed which improve the stability and convergence of the model. Generally there are various correction methods in the SPH method and XSPH is one of them. Particularly in the simulation of high speed flow or impact this technique must be considered. For instance, numerical noise in velocity regardless of direction may even cause unwanted penetration between particles resulting in the failure of the simulation. To clarify, consider closely spaced particles which have similar but not equal velocity. The particle with higher speed has priority and passes over the other particles; this situation results in particle penetration which is not allowed for the continuum material. The main purpose of using the XSPH correction technique is to correct and modify velocity to the proper value which is related to the velocities of the neighbor particles. In this case, particles can move smoothly [11]. The equations below are used to calculate the velocity.

(7){\displaystyle \frac{\mathrm{d}{\upsilon }_i}{\mathrm{d}t}}={\upsilon }_i+ϵ\mathrm{\Delta }{\upsilon }_i,

where,

(8)\mathrm{\Delta }{\upsilon }_i={\displaystyle \sum _j{\displaystyle \frac{m_j}{{\rho }_j}}({\upsilon }_j-{\upsilon }_i)}W(r_i-r_j,h).

In this study, the Von-Mises yield criterion is used. The following steps show the procedure of applying governing equations:

1) First, reduced tensions must be compared with the yield condition. If S_x^{\mathrm{2}}+S_y^{\mathrm{2}}+\mathrm{2}{S}_{xy}^{\mathrm{2}}\le {\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{(\sigma )}^{\mathrm{2}}, then the value of stress corresponds to the yield stress and particles show elastic behavior. But if S_x^{\mathrm{2}}+S_y^{\mathrm{2}}+\mathrm{2}{S}_{xy}^{\mathrm{2}}>{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{(\sigma )}^{\mathrm{2}}, then the material yields and it is in the plastic region.

2) If the material yields in the first stage, in this step all stress elements are multiplied by the following factor to return to the yield level

(9){\displaystyle \frac{\sqrt{\mathrm{2}/\mathrm{3}}\sigma }{\sqrt{S_x^{\mathrm{2}}+S_y^{\mathrm{2}}+\mathrm{2}{S}_{xy}^{\mathrm{2}}}}}.

3) Stresses can be calculated from the equation below

(10){\sigma }^{ij}=-P{\delta }^{ij}+S^{ij},

where {\sigma }^{ij} is the stress tensor, P is pressure and S^{ij} is the deviatoric stress tensor. In Eq. (10) {\delta }^{ij} is the Kronecker delta. It is noted that the mean stress σ⁰ is given by σ⁰ = -P. Assuming Hooke’s law with shear modulus, the evolution equation for the deviatoric stress S [12] is:

(11){\displaystyle \frac{\mathrm{d}{S}^{ij}}{\mathrm{d}\mathrm{t}}}=\mathrm{2}\mu \left({\dot{\in }}^{-ij}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\delta }^{ij}{\dot{\in }}^{-ij}\right)+{\mathrm{\Omega }}^{jk}{S}^{ik}+{\mathrm{\Omega }}^{ik}{S}^{jk},

where,

(12){\dot{\in }}^{-ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial {\upsilon }^i}{\partial x^j}}+{\displaystyle \frac{\partial {\upsilon }^j}{\partial x^i}}\right),{\mathrm{\Omega }}^{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial {\upsilon }^i}{\partial x^j}}+{\displaystyle \frac{\partial {\upsilon }^j}{\partial x^i}}\right).

Modeling and results

In this section we discuss about the modeling of our case study. The codes used in this study are written in FORTRAN V.90 and the results are simulated in TECPLOT V.360 software package. The simulations are done for different velocities and materials, however, for each case the projectile and the target plate are made of the same metals. As it is shown in Fig. 2, particles are initially located in a regular grid with the initial inter-particle distance of 0.01 cm. This involves a spherical projectile with a 5 mm radius impacting on 2 mm thick target plate. The simulation is started at the moment of impact; therefore particles of the projectile are assigned the projectile velocity which are 300, 1000, 2000 and 3100 m/s.

The total number of particles used in this study is 17985 and the ratio of h/d is equal to 1.5 (where h is the smoothing length and d is the initial particles distance). This simulation is done in 2D planar geometry carried out for the total duration of 8 µs after impact.

At the beginning of modeling, due to the insufficient number of particles (less than 10000 particles), no jet of particles is observed coming out of the wall where the projectile is impacting. A rupture also occurred in the wall. These issues are clearly seen in Fig. 3.

Three quantities are tracked during simulation. Results for crater diameter, final longitudinal diameter of projectile and projectile penetration length are presented in Tables 1, 2, and 3 relating to the AL-AL, Lead-Lead, and Steel-Steel impact respectively. Note that the steel alloy used in this project is Steel-SAE-1030.

According to the results obtained from the above models, it can be seen that the projectile penetration increases by increasing the velocity. The rate of change of the crater diameter, for all three types of materials, decreases with each step of velocity changes. According to the above tables, for all materials with velocity of 3100 m/s, the targets are entirely penetrated and a debris cloud is formed to the back of the target entailing material of both projectile and target. Comparing the above tables, the obtained results from (AL-AL, Lead-Lead and Steel-Steel) simulations are quite similar to each other.

Some of the above cases are simulated using the TECPLT 360 software package. Figures 4–7 depict impact of aluminum projectile hitting aluminum target plate at velocity of 3100 m/s. It can be observed that the particles’ continuity is maintained over the time.

Figure 8 displays a general view of the above figures.

In the following, Figs. 9-12 illustrate velocity variation in X direction related to the AL-AL impact with the initial velocity of 3100 m/s.

The upper half of the velocity variation in X direction for Lead-Lead impact (at 8 µs after impact) for the velocities of 300, 1000, 2000, 3100 m/s are shown in Figs. 13-16 respectively.

Validating data with published literature

In this section, our results for a specific case (AL-AL impact at the velocity of 3100 m/s) are validated by comparison with a previous paper published by Mehra and Chaturvedi [13]. They have studied high velocity impact of aluminum sphere on thin aluminum plate as a smooth particle hydrodynamics study. They used four different shock capturing patterns which are used in SPH, and then applied the equations in the high velocity impact and hyper velocity impact of projectile on thin aluminum plate which are around 3 and 6 km/s respectively. These two simulations are applied for the smooth geometrically plane. An elastic––perfectly plastic constitutive prototype has been used. Hence, it is possible to compare the simulation with simulation done by Howell and Ball [14]. They compared the results from the simulation with the data which had been previously examined and published. These patterns are different in the manner of treating of artificial viscosity (AV). These schemes are included regular SPH artificial viscosity, Monaghan and Morris’s artificial viscosity (MON), Balsara’s artificial viscosity (BAL), Riemann based contact algorithm (CON) of Parshikov. Mehra and Chaturvedi [13] found that in the impact problems with moderately high velocity, CON achieved best result overall. Figure 6 shows the arrangement of particles used in SAV, BAL, MON, and CON methods.

The main differences between their modeling and what is presented here are number of particles, h/d and how particles are initially distributed. As shown in Fig. 17, particles are initially located on a regular mesh. The mesh is rectangular except that the circularity of the projectile is surrounded by two rings of equidistant particles along the circumference [13]. In Mehra and Chaturvedi, the initial distance between particles is 0.01 cm. Wall thickness and radius of the projectile are the same as in the model discussed in this paper. They have a total of 17,850 particles in the simulation including 7820 particles in the projectile. The results are given for h/d = 1.4 for all AVs, with the exception of CON, where h/d = 1.7 is used.

Upper half of the configuration obtained by AL-AL impact at a velocity of 3100 m/s, at the time, 8 µs after impact, are shown in Figs. 18 and 19 for different AVs and the method used in this paper respectively.

Comparing the position of particles in different methods, CON method and the SPH method used in this study have the least trouble with numerical fracture and the tendency to form clumps at moderately high velocity impact.

Results of three mentioned parameters produced by AL-AL impact at the velocity of 3100 m/s, for this specific case, are shown in Table 4.

As shown in Table 4, the SPH technique used in this project generally produces slightly larger crater diameters than those reported by Mehra and Chaturvedi [13] but slightly smaller final longitudinal diameter of projectile. All SPH simulations are in agreement among themselves on geometrical measurement of penetration length, crater diameter and final longitudinal diameter of projectile.

In summary, the SPH code used in this paper is validated for a specific case which is impact of AL-AL at the velocity of 3100 m/s by the work done by Mehra and Chaturvedi [13]. The results demonstrate that there is a positive correlation between the particles cohesion and the accuracy of the methods which means that CON method and the SPH method used in this study can produce the best results.

Sensitivity analysis

In this section, we will examine the sensitivity of the model with respect to small changes in the input variables. Sensitivity analysis (SA) is becoming an important component of model validation and is part of the more general subject of uncertainty quantification. All input data are subject to uncertainty of varying degrees, either due to assumptions made or due to measurement errors. To assess the validity of a model, it is important to quantify the effect of these uncertainties on the output variables. As a side benefit, SA requires several test runs of the model which can be used for validation and for discovering software defects.

There are several types of methods that can be employed for SA such as local methods, variance based methods, metamodeling or Monte Carlo filtering. Many of these methods require a large number (500-1000) of model runs which is prohibitively expensive for this problem since each run requires more than 10 hours of CPU time on a standard desktop computer. In the following, we will limit ourselves to estimating the local sensitivity for two input variables (initial velocity and plate thickness) and three outputs (crater diameter, final longitudinal diameter of the projectile and projectile penetration).

First we calculate the sensitivity index [15], which, for an input parameter P, and output parameter Q is defined by the equation:

(13)S{I}_{PQ}=\left(\left|{\displaystyle \frac{Q_{\mathrm{0}}-Q_L}{P_{\mathrm{0}}-P_L}}\right|+\left|{\displaystyle \frac{Q_{\mathrm{0}}-Q_H}{Q_{\mathrm{0}}-Q_H}}\right|\right)/\mathrm{2},

where index 0 represents the base value, and the index L and H represent an increase and decrease from the base value of the parameter P, or the corresponding outputs Q. In all of the calculations, we considered a 5% increase and decrease from the base value of the two input parameters (initial velocity and plate thickness). Note that SI_PQ is an approximation of the partial derivative \partial Q/\partial P at P₀, which measures the absolute change in the output and is dependent on the units used. A measure of sensitivity that is dimensionless and measures relative change is the elasticity index, defined as

(14)E{I}_{PQ}=S{I}_{PQ}\cdot (P_{\mathrm{0}}/Q_{\mathrm{0}}).

This index can be interpreted as an estimate of the proportional change in the output Q caused by a 1% change in the input P at P₀. We evaluate EI_PQ at the four base initial velocity values and a base plate thickness of 2 mm, and show the results in Table 5. We observe that the three outputs are more sensitive to the initial velocity than the plate thickness, especially in the case of projectile penetration. The elasticity index is generally less than 1, except for the elasticity index of the velocity with respect to penetration which is slightly greater than 1.

Next, we will consider a first-order sensitivity analysis that takes into account the variance of the input parameters. In particular, given a probability distribution of the input parameters, we would like to determine the variance of the output and how much each input contributes to the output variance. The percentage contribution of each input parameter P_i can be written as [16]

(15)PC[P_i,Q]={\displaystyle \frac{{(\partial Q/\partial P_i)}^{\mathrm{2}}\mathrm{V}\mathrm{a}\mathrm{r}[P_i]}{{\displaystyle {\sum }_j^{}{(\partial Q/\partial P_i)}^{\mathrm{2}}/\mathrm{V}\mathrm{a}\mathrm{r}[P_i]}}}\cdot \mathrm{100.}

Here the partial derivatives \partial Q/\partial P_i can be approximated from the sensitivity indices SI_PQ. In the following, we will assume that the range between the high and low values of the input parameters correspond to a 5-95 percentile range for a normally distributed variable (i.e., a 90% probability interval or 3.3 standard deviations). Then the variances of the input parameters can be calculated from the equation

(16)\mathrm{V}\mathrm{a}\mathrm{r}[P]={\left({\displaystyle \frac{P_H-P_L}{\mathrm{3.3}}}\right)}^{\mathrm{2}}.

In Table 6, we summarize the percentage contribution of the two input variables for the four pairs of base values, corresponding to the four different initial velocities and the wall thickness of 2 mm.

We note again that more than 50% of the contribution comes from the variance in the initial velocity for all 3 input variables, and the contribution is particularly significant for the variance in projectile penetration (close to 95%). Therefore we can conclude that outputs are more sensitive to changes in initial velocity than wall thickness.

Conclusions

SPH is a particle based method for modeling fluid flow and solid deformation and it has been used successfully to model impact problems. For these impact problems, SPH has the advantages of being able to follow very high deformations (further than what is possible with finite element method and finite volume methods) and to keep track of the specific history of each part of the metal and potentially direct prediction of many types defects occurring during impact.

In this paper, an implementation of the SPH method to deal with two-dimensional impact has been presented. This work is mainly focused on recording three parameters which are crater diameter, final longitudinal diameter of projectile and projectile penetration; then this method is compared with four different types of SPH techniques. Elastic-perfectly plastic constitutive equation is used in order to compare with free Lagrangian simulation of [13].

The numerical simulation showed a good illustration of the basic effects such as penetration of projectile and formation of the crater.

A sensitivity analysis of the result showed that the elasticity index for the input/outpt pairs examined at various data points is generally less than or close to one, which indicates the outputs are not overly sensitive with respect to the inputs. It was also shown that most of the variance in the outputs is due to the variance in the initial velocity rather than the plate thickness.

References

[1] Rabczuk T, Belytschko T. Adaptivity for structured meshfree particle methods in 2D and 3D. International Journal for Numerical Methods in Engineering, 2005, 63(11): 1559–1582

[2] Rabczuk T, Areias P M A, Belytschko T. A simplified mesh-free method for shear bands with cohesive surfaces. International Journal for Numerical Methods in Engineering, 2007, 69(5): 993–1021

[3] Rabczuk T, Eibl J, Stempniewski L. Numerical analysis of high speed concrete fragmentation using a meshfree Lagrangian method. Engineering Fracture Mechanics, 2004, 71(4-6): 547–556

[4] Gingold R A, Monaghan J J. Smoothed particle hydrodynamics: Theory and application to non-spherical stars. Monthly Notices of the Royal Astronomical Society, 1977, 181: 375–389

[5] Lucy L B. A numerical approach to the testing of the fission hypothesis. Astronomical Journal, 1977, 82: 1013–1024

[6] Cleary P W, Monaghan J J. Conduction modelling using smoothed particle hydrodynamics. Journal of Computational Physics, 1999, 148(1): 227–264

[7] Hu X Y, Adams N A. A multi-phase SPH method for macroscopic and mesoscopic flows. Journal of Computational Physics, 2005 (in press)

[8] Liu M B, Liu G R, Zong Z, Computer simulation of high explosive explosion using smoothed particle hydrodynamics methodology. Computers & Fluids, 2003, 32(3): 305–322

[9] Johnson G R, Petersen E H, Stryk R A. Incorporation of an SPH option into the EPIC code for a wide range of high velocity impact computations. International Journal of Impact Engineering, 1993, 14(1–4): 385–394

[10] Greenberg M D. Advanced Engineering Mathematics. New Jersey: Prentice Hall, 1998

[11] Zou S. Coastal sediment transport simulation by smoothed particle hydrodynamics. Dissertation for the Doctoral Degree. Baltimore: John Hopkins University, 2007

[12] Gray J P, Monaghan J J, Swift R P. SPH elastic dynamics. Computer Methods in Applied Mechanics and Engineering, 2001, 190(49–50): 6641–6662

[13] Mehra V. Chaturvedi S. High velocity impact of metal sphere on thin metallic plates: A comparative smooth particle hydrodynamics study. Journal of Computational Physics, 2006, 212: 318–337

[14] Howell B P, Ball G J. A free-lagrange augmented Godunov method for the simulation of elasticplastic solids. Journal of Computational Physics, 2002, 175(1): 128–167

[15] McCuen R. The role of sensitivity analysis in hydrologic modeling. Journal of Hydrology (Amsterdam), 1973, 18(1): 37–53

[16] Benjamin J R, Cornell C A. Probability, Statistics, and Decision for Civil Engineers. New York: McGraw-Hill Book Company, 1970

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg