Proceedings of the ASME 2010 10th Biennial Conference on Engineering Systems Design and Analysis ESDA2010 July 12-14, 2010, Istanbul, Turkey The ASME 2010 10th Biennial Conference on Engineering Systems Design and Analysis ESDA 2010 12-14 July 2010, Istanbul, Turkey

ESDA2010-0 Put paper number here SIMULATION OF IMPACT AND FRAGMENTATION WITH THE MESHLESS METHODS Namık Kılıç, Atıl Erdik Otokar Otomotiv ve Savunma San. A.Ş, CAE Department Sakarya, Turkiye

Ass. Prof.Bülent Ekici Marmara University Mechanical Engineering Department İstanbul, Turkiye

ABSTRACT High velocity impact and penetration problems include large deformation, erosion, high strain rate dependent nonlinear material behavior and fragmentation. Therefore, meshless methods seem to be ideally suited for the modeling of penetration events as they allow unrestricted deformation and easy tracking of material interfaces and loading histories. In the first part of this study, a brief overview about meshless methods is given. Also the most important features of meshless methods with respect to mesh based approaches are compared. In the second part, numerical model is developed using one of the most frequently used meshless method, Smoothed Particle Hydrodynamics (SPH). 3D numerical simulations are performed on a high performance computer using MPP version of the explicit finite element code LS-Dyna. For reasonable behavior of material response under dynamic loading, Johnson Cook material models for armor steel target and 7,62 armor piercing projectile are derived using SHPB (Split Hopkinson Pressure Bar) test data. SPH computational investigation is compared with available experimental data such as penetration depth and impact crater diameter. For the future work, other potential meshless methods for ballistic impact problems are identified.

approaches are become quite popular in investigation of ballistic penetration. These are test based empirical formulations, analytical and numerical methods. Experimental investigation of ballistic penetration in development of armor and ammunition also requires quite complicated test and data analysis procedures, therefore new methods in experimental studies are still in progress. With the tests conducted in laboratory environment, numerous empirical formulations are derived and used in solution of problems [1]. Although, its difficulty in derivation, analytical models are useful from the view of their direct applicability. Since the qualitative laws determined by the use of these models can be considered as basis for the theoretical and experimental investigations, there are many recent studies on this field [2]. Multiple physical phenomena are involved during penetration such as fracture, failure, residual stresses and friction heating [3]. In reality since empirical and analytic approaches cannot capture all of these phenomena, numerical simulation has become a necessary tool for the study of ballistic penetration. Numerical methods and corresponding computing technologies have been recently evolved to the level where mentioned complex deformation and penetration pattern during ballistic impacts can be accurately predicted. The usage of numerical methods in the development of design alternatives not only introduces shorter development times, but also reduces the number of prototypes to minimize the number of the real scaled field tests required, and also help us to visualize the impact behavior. In a numerical model of a continuum, the material is discretised into finite sections over which the conservation and constitutive equations are solved. [4,5,6,12]. The way in which this spatial discretisation is performed leads to different numerical methods. A review of the literature on

INTRODUCTION In civil and military applications, ballistic penetration has become prime importance in development of armor solutions and continues to be a challenging research field for engineers. High velocity impact and penetration problems include large deformation, erosion and high strain rate dependent nonlinear material behavior. Due to complexity of the problem three

1

1

Copyright © 20xx by ASME

Copyright © 2010 by ASME

Although great success of finite element methods in many fields from academic to industrial application, there are some limitations especially arise in simulation of ballistic penetration. The meshless methods were born with the objective of eliminating below issues: [7]  In stress calculation, solutions are discontinuous and therefore less accurate  When dealing with large deformations considerably accuracy is lost because of element distortions  It is always difficult to simulate crack growth due to discontinuities do not coincide with the original nodes  It is difficult to simulate breakage of material into large number of fragments as finite element method essentially based on continuum mechanics, in which the elements formulated can not be broken. In FE solutions, the elements can be either eroded or stay as a whole piece. Serious error can occur because the nature of the problem is nonlinear and results are highly path dependent.  In remesh approach, the problem domain is remeshed at steps during the simulation process to prevent severe distortion. By the way the node lines allowed remaining coincident with the discontinuity boundaries. We can say that there are no reliable pre processing code is available for creating hexahedral meshes for 3D problems. The concept of meshless methods has great potential for solution of above mentioned problems of ballistic simulation.

impact simulations shows that the most research in this field has focused on the development and application of continuum hydrocodes using either Eulerian, Lagrangian or Arbitrary Lagrangian Eulerian (ALE) formulations. In the Lagrange solver, the numerical mesh moves and distorts with the physical material. This formulation is widely used because of its advantages, such as low level of computational cost [3], being able to track accurately and efficiently material interfaces and incorporate complex material models. However, it is very sensitive to distortions resulting small time step and possible loss of accuracy. The well known negative volume errors occur as a result of mesh tangling. Numerical codes can handle these problems with adaptive re-meshing algorithms or by eroding the highly distorted elements by artificial manipulation which can cause loss of accuracy. In Euler solver, the numerical mesh is fixed in space and physical material flows through the mesh [11]. This formulation is generally used to represent large deformations and flow problems. Free surfaces of material interfaces can move with fixed mesh, which also needs to model void mesh. In order to describe solid behavior, solid stress tensor and history of the material must be transported from cell to cell, therefore has less applicability in ballistic penetration. In the ALE solver, solver allows for a type of automatic rezoning which can be quite useful when Lagrangian mesh structure starts to distort [17]. It is particularly well suited for a variety of fluid structure interaction problems.

Meshless Backgrounds

Despite the fact that mesh based numerical methods are the primary computational methodology in engineering computational mechanics, they still have limited applications in many complex problems. Accurate and realistic simulation of ballistic penetration with mesh based methods requires complicated contact and erosion algorithms. When large deformation occurs, the accuracy of results significantly decreases. Additionally it is difficult to simulate structural behavior of material including fragmentation since mesh based formulation is based on connectivity of nodes. This is especially will be a problem in accurate simulation of ceramic materials which are highly used in composite armor structures. Mesh based formulation development continues to overcome numerical weakness experienced in ballistic simulation. Especially in Lagrange method adaptive remeshing is recommend to fix mesh tangling issues and successfully implemented by many researchers. On the other hand, high computational cost of adaptive remeshing algorithms and its limited applicability on three dimensional problems decrease the attractiveness.

Over the last decade a number of meshless methods have been proposed. They can be classified with the definition of shape functions or in accordance with the minimization way of approximation. The two most commonly used approximation theories in meshless methods are the moving least squares (MLS) and the Reproducing Kernel Method (RKM) [4]. The MLS approach uses a local least squares approximation around a fixed point. Reproducing Kernel Methods are in general a class of operators that reproduce the function itself through integration over the domain. From the continuous form of the RKM approximation a discretised form, the Reproducing Kernel Particle Method (RKPM), was derived. One of the first and the simplest meshless methods was the Smooth Particle Hydrodynamics (SPH) method which is similar to the RKPM concept. In the Meshfree Local Petrov-Galerkin (MLPG) method (Atluri and Zhu 1998) a local weak form is used instead of the global form of a standard Galerkin approach. The local weak form is formulated over all local sub-domains and the integration is performed over these domains. Thus this method is called “truly meshless”, since no background mesh is required. Several special cases have been derived for the MLPG concept depending on the choice of the test functions.

MESHLESS METHODS

2

2

Copyright © 20xx by ASME

Copyright © 2010 by ASME

Alternative “non-standard” approaches to construct the partition of unity are the Sibsonian and Non-Sibsonian which are used in the Natural Neighbor Galerkin Method (Sukumar et al. 1998) also called Natural Element Method (NEM) [9,10]. In this method a Voronoi diagram of the nodal domain is used to determine the neighbor nodes of an interpolation point.



u  x    u  y    x  y  dy

(1)



Where  x  is the Dirac delta function. One can define delta function as follows:

  x  y  dy  1

Smoothed Particle Hydrodynamics (SPH) – Theory

(2)



Note that this integral representation is exact but it is difficult to use in numerical analysis. u  x  is approximated by the following integral form of representation

Smoothed Particle Hydrodynamics is a gridless Lagrangian hydrodynamics computational technique, which possess individual material properties and move according to the governing conservation equations. SPH was first developed by Lucy as Monte-Carlo approach to solve hydrodynamic time evolution equations. Subsequently Gingold and Monaghan, reformulated the derivation in terms of an interpolation theory. According to interpolation derivation each SPH particle represents a mathematical interpolation point at which the fluid properties are known. The complete solution is obtained at all points in space by application of an interpolation function. It has some special advantages over the traditional mesh-based numerical methods. The most significant is the adaptive nature of the SPH method, which is achieved at the very early stage of the field variable (i.e. density,velocity, energy) approximation that is performed at each time step based on a current local set of arbitrarily distributed particles. Because of the adaptive nature of the SPH approximation, the formulation of the SPH is not affected by the arbitrariness of the particle distribution. Therefore, it can handle problems with extremely large deformations very well. Unlike the meshfree nodes in other meshfree methods, which are only used as interpolation points, the SPH particles also carry material properties, functioning as both approximation points and material components. These particles are capable of moving in space, carry all computed information, and thus form the computational frame for solving the partial differential equations describing the conservation laws. The numerical solution procedure of the SPH formulation consists of the following steps [7]: (i) generation of the meshless numerical model, (ii) integral representation (kernel approximation), (iii) particle approximation, (iv) adaptation and (v) dynamic analysis.

u h  x    u  y  W x  y, h  dy

(3)



where u  x  represents the approximation of function u  x  , is a kernel or the smoothing function and h is the smoothing length determining the influence domain of the smoothing function. The particle approximation of the function (Equation 2) can be defined by below equation for moving particles xi t  i  1 N  where xi t  is the location of particle i: h

N

u h xi    w j u xi  W xi  x j , h 

(4)

j i

Where

wj 

mj

j

is the weight of the particle. The

formulation for conservation equations and derivation for weak forms are available in Reference 4. During the entire computational simulation it is important to know which particles interact with each other and which particles are within the influence domain. The search of neighboring particles in a model with N particles requests to perform N-1 distance calculations for only one particle, which consequently increases computational time. This can be overcome by the bucket sort method, which is based on splitting the analyzed domain into several boxes of a given size (Fig.1). This method reduces the number of distance calculations therefore reduces the computational time and increases its efficiency [4].

First point to be considered in explaining theoretical fundamentals of SPH is the interpolation theory; the second point will be the conservation laws of continuum in the form of partial differential equations. The SPH method uses integral representation of a function. Consider a function u  x  at any point of by

x  x, y, z  . Its integral representation can be given

Fig. 1. Bucket sort and neighbor search [8].

3

3

Copyright © 20xx by ASME

Copyright © 2010 by ASME

Despite all described advantages, the SPH method has still to cope with some numerical difficulties, like particle inconsistency, inaccuracy at domain boundaries and instabilities at tensile stress state [7]. However, the SPH method accuracy and stability also depend on the particle number (particle density) within the influence domain and the time step of time integration scheme. In application, the principal potential advantage of the SPH method is that there is no need for connectivity between particles with a conventional mesh, hence avoiding element distortion problems at large deformations. In comparison with the Eulerian description, the SPH method offers higher efficiency in terms of modeling domain, since only the material domains have to be discretised and not also the areas through which the material might move during the simulation.

The total mass of AP bullet is about velocity of the ammunition is 854 m/s 4569. The simulation will use the hardened steel core projectile and discretization for steel armor target.

10.04 g. The muzzle as specified in Stanag Lagrange method for SPH and Lagrange

NUMERICAL SIMULATIONS Numerical Modelling With the recent developments, SPH method became a reliable tool providing accurate and stable results for large deformation problems and successfully implemented in many commercial hydrocodes. One of FE codes which also includes SPH is LSDyna. In this study, due to high computational cost of SPH, an alternative numerical solution technique of coupling FEM/SPH is used. Since SPH uses a Lagrangian formulation, a coupling between SPH and Lagrangian mesh is possible with defining contact interfaces. For impact problem, the large deformation section of the target is modeled with SPH and where the small deformation occurs is modeled with solid finite elements. Such coupling should be carefully arranged to simulate wave propagation.

Fig. 3. Discretised Model The circular Lagrange layers of the armor are divided into two regions. In radial direction mesh is gradually coarsening from inner to outer. In SPH section of armor, an equispaced discretization applied. In simulations, gradually coarsening SPH discretization also applied but the results are not improving therefore omitted.

Geometry The geometric model consist of an armor steel target and Stanag 4569 level 3 7.62x54R AP ammunition. For the AP bullet, the ogival nose hardened steel core inserted in a brass sabot, before the jacket is clamped onto it. In front of the core, a cap of lead-antimony is placed. The purpose of this cap is to stabilize the projectile during flight and in the initial stage of penetration.

Fig.4. SPH Discretised Target Model Material Model The Johnson-Cook model is used for target and projectile as it is well suited to the model metals that are subject to high strain rate loading. The yield stress at non zero strain rate depends on the strain hardening, strain rate hardening and temperature softening such that [16]

Fig. 2. Geometric cross-section of AP ammunition

4

4

Copyright © 20xx by ASME

Copyright © 2010 by ASME







 Y  A  B np 1  C log  *p 1  THm



(5)

Where p is the equivalent plastic strain,

 *p is the

dimensionless plastic strain rate, TH is homologous temperature { TH=(T-Troom)/(Tmelting-Troom) }. The five material constants are A, B, C, n and m. The expression in the first set * of brackets gives the stress as a function of  p =1 and TH =0.

The expression in the second and third sets of brackets represent the affects of strain rate and thermal softening respectively. Fracture in the JC material model is based on a cumulative damage law:

D

 f

(6)

in which





  D1  D2 exp( D3 *)1  D4 ln  1  D5TH f

* p



Fig.5. Tensile test data for 0,1 s-1 strain rate (7)

For determining J-C material strength model for bullet, Dynamic Test & Modeling Laboratory at Izmir Institute of Technology (IYTE) performed experiments. SHPB is a mechanical test instrument used to characterize the dynamic response of materials at high strain rates (Figure 6).

where * is the ratio of hydrostatic stress per effective stress. The constants D are experimental constants.

*  

H  eq ( x   y   z ) / 3 2 x

2 y

2 z

       x y   y z   z x  3( xy2   yz2   zx2 )

Failure of the elements assumed to occur when D=1. The failure strain and thus the accumulation of damage is a function of mean stress, strain rate and temperature. The failed eleemnts are removed from the FE model. The JC material model was used in conjuctio with the Mie-Gruneisen equation of state model. The armour material used in ballistic test is alloyed, liquidquenched and tempered high-strength special steel for civil use. The measured hardness values and mechanical properties suit the Armox 500 specified by Skoglund and coworkers. The JC parameters A = 1470 MPa, B = 702 MPa, n = 0.199, C = 0.00549 and finally m = 0.811 are used in simulations [18]. The failure parameters D1=0,068, D2=5,382 D3=-2,554 are also taken from Skoglund and coworkers and compared with tensile test results(Fig 5). In simulations it is assumed that triaxiality has a larger influence on the strain to fracture than the strain rate in the investigated loading range, therefore D4=0. Thermal softening affects are also ignored D5=0.

Fig.6. Split-Hopkinson Pressure Bar Aparatus A Gas gun launches a striker bar that impacted on the end of the incident bar. A stress wave is generated that travels down the bar and is recorded sequentially by the first and second strain gages mounted longitudinally on the bar. The stress wave then passes through the specimen and the specimen is compressed. Part of the stress wave is reflected in the form of a tensile pulse and the reminder is transmitted to the transmitter bar and recorded by the third strain gage mounted on the transmitter bar. The three readings are used to determine the time dependent stress state of the specimen. From the time dependent strain state data, a stress vs. strain plot can be obtained.

5

5

Copyright © 20xx by ASME

Copyright © 2010 by ASME

Fig.7. A Schematic representation of the SHPB Voltage-stress conversion formula gives stress values of dynamic test: 2   2C 0V .10 6 x x1000 / L0 (8) (Gain).( f )(1   )( Ec) Where : Stress  C0 : Wave speed V : Voltage values Gain : Gains of amplicator f : Factor of strain gage : Poisson ratio  Ec : Stress bridge voltage L0 : Initial lengths

Fig.9. AP Bullet Tensile Test Material strength model for bullet has compared with Buchar and co-workers. [14] The material model is close to Buchar results with alterated yield strength.

Equation of State The relation between pressure, volume and internal energy of a material is defined as the equation of state. In this numerical modeling, Mie-Gruneisen and linear polynomial equations of state are used.

 0 C 2  1  1   2   2  2  P 2 1  S1  1  S 2  2 /   1  S 3  3 / 1    (9)     E0



Fig.8. Notched Dynamic Test Specimen



Where E 0 is the internal energy per unit volume, C is the intercept of the shock and particle velocity curve, S1 , S 2 , S 3

The quasi-static stress-strain response was measured for the steel core material. Tensile specimens were machined, instrumented with strain gauges and strained until failure. The test results are plotted in Figure 8. It was found that the core material fails approximately 4% strain at a stress of 2.1 GPa. However, there is variation in strain to failure and associated ultimate stress.

are the cofficients of the slope of the shock and particle velocity curve,  is the Gruneisen coefficient.  is the volume correction factor and

    0  1 is the compression

factor. The Mie Gruneisen EOS data for the target plate is taken from Reference 15 as 4340 material. For the AP bullet, the linear polinomial EOS is used.

PK 

(10)

Where K is the bulk modulus for steel. Simulation Model

6

6

Copyright © 20xx by ASME

Copyright © 2010 by ASME

Using the geometry detailed above a full model was generated for the analysis. The target core is modeled as SPH with particle size of 0,2 and surrounding Lagrange regions are modeled with elements size of 0,5 and 1 mm. The bullet is modeled with 3744 Lagrange solid elements and target is modeled with 144.000 Lagrange and 786.000 SPH elements. The contact between SPH and Lagrange elements is developed using CONTACT_TIED_NODES_TO_SURFACE. The contact during the penetration of bullet into target is achieved using ERODING_NODES_TO_SURFACE. In SPH simulations, the artificial bulk viscosity ignored. All simulations including Lagrangian and SPH discretisation models have been performed using a 20 CPUs high performance compute cluster.

Fig.11. Steel Core After Shot Test

RESULTS

After increasing dynamic yield stress, the results become more realistic and predict the depth of penetration with reasonable approximation. (Figure 12). In test the depth of crater varies between 12,3-12,9 mm and the approximate diameter of crater varies between 16,1-17,5 mm. With Lagrangian simulation penetration is overestimated by 25% and crater diameter is underestimated by 20%. For 20 mm target, the back face protrusion occurred but height is ignorable. In simulations also the height of protrusion is overestimated.

In order to see potential problems of Lagrange discretisation, simulations performed with modeling target and projectile with Lagrange elements. The initial results were disappointing for the hardened steel core AP ammunition. The bullet appeared to be eroded prematurely and the penetration depth results are underestimated for 20 mm target (Figure 10). Contrary to the simulations, in tests, the hardened steel core is not eroded and do not have significant deformation on tip. (Figure 11) In order to improve predictions dynamic yield strength of bullet increased by changing strain hardening coefficient.

Fig.12. Lagrange Discretisation, Penetration on 20 mm Target In SPH discretisation, since the tip of steel core did not deformed, projectile is modeled with rigid elements. During shots, residual speed is measured using 42.000 fps high speed camera with 256x176 resolution.

Fig. 10. Simulation Results for Lagrange Discretisation

Fig.13. 20 mm Armor Steel After Shot Test

7

7

Copyright © 20xx by ASME

Copyright © 2010 by ASME

[8] Vesenjak N., Ren Z., “Application Aspects of the Meshless SPH Method”, Journal of the Serbian Society for Computational Mechanics, Vol. 1, 2007, pp. 74-86 [9] Sukumar N., Moran B., Belytschko T., “The Natural Element Methods in Solid Mechanics”, International Journal of Numerical Methods in Engineering, Vol. 43, 1998, pp. 839887 [10] Sukumar N., Cueto E.., Calvo B., “Overview and Recent Advances in Natural Neighbour Galerkin Methods”, International Journal of Numerical Methods in Engineering, Vol. 10, 2003, pp. 307-384 [11] Buyuk M., Kan S.C.D.,Bedevi N.E., Durmus A.,Ulku S., “Moving Betond the Finite Elements, a Comparison Between the FEM and Meshless Methıds for a Ballistic Impact Sİmulation”, 8 th International LS-DYNA Users Conference [12] Century Dynamics Inc., “Sph User Manual & Tutorial”, Revision 4.3, 2005 [13] Belytschko T., Lu Y., Gu L., “Elemet Free Galerkin Methods”, International Journal of Numerical Methods in Engineering, Vol. 37, 1994, pp. 229-256 [14] Buchar J., Voldrich J., Rolc S., Lisy J., “Ballistic Performance of Dual Hardness Armor”, 20th International Symposium on Ballistics Orlando, 23-27 Sep 2002 [15] Li J., Li X.J., Zhao Z., Ou Y.X., Jiang D.A., “Simulation on Projectile with High Rotating Speed Peentration into the Moving Vehicular Door”, Theoretical and Applied Fracture Mechanics Vol 47, 2007, pp113-119 [16] Johnson G.R., Cook W.H., “A Constitutive Model and Data for Metals Subjected to Large Strains, High Strain Rates and High Temperatures”, [17] Donea J., Huerta A., Phontot J., Ferran A.R., “Arbitrary Lagrangian Eulerian Methods”, Encyclopedia of Computational Mechanics Vol1, 2004, John Wiley&Sons [18] Skoglund P., Nilsson M., Tjernberg A., “Fracture Modelling of a High Perfromance Armour Steel”, J.of Physic IV France Vol 134, 2006, pp197-202

Fig.14. SPH Discretisation, Penetration on 20 mm Target CONCLUSION/DISCUSSION This paper has illustrated that some of the difficulties involved in modeling of complex fragmentation behavior of hardened steel core armor piercing projectile impacting on armor steel having hardness of 500 BHN. In general, Lagrange simulation issues are solved with finer meshing. For SPH, even with finer mesh the crater diameter and fragmentation of armor would not have necessarily predicted more accurately.

ACKNOWLEDGMENTS The financial and technical support of this work from Otokar Otomotiv ve Savunma Sanayi A.S. is gratefully acknowledged.

REFERENCES [1] Bangash M.Y.H., “Shock, Impact and Explosion Structural Analysis and Design”, Springer 2009 [2] Ben-Dor G., Dubinsky A., Elperin T., “Ballistic Impact: Recent Advances in Analytical Modelling of Plate Penetration Dynamics”, ASME Applied Mechanics Reviews Nov 2005 Vol.58 pp 355-371 [3] Zukas J. A., Nicholas T., Swift H. F., Greszczuk L. B., Curan D. R.; “Impact Dynamics”, J. Wiley, New York, 1982 [4] Livermore Software Technology Corporation (LSTC), “LSDyna Keyword User’s Manual”, May 2007, Version 971 [5] Abaqus Theory Manual, Version 6.7 [6] MSC.DYTRAN Theory Manual Version 2005, MSC.Software Corporation [7] Liu G. R., “Mesh Free Methods”, CRC Press, Florida 2003

8

8

Copyright © 20xx by ASME

Copyright © 2010 by ASME