Int. J. Impact Engng, Vol. 20, pp. 143-152, 1997 ©1997 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0734-743X/97 $17.00+0.00
DEVELOPMENT OF LAGRANGIAN HYDROCODE M O D E L L I N G FOR DEBRIS IMPACT DAMAGE PREDICTION J. C A M P B E L L and R. VIGNJEVIC College of Aeronautics,CranfieldUniversity,Cranfield,Beds, MK43 0AL. UK.
Summary--Work on improving the Lawrence Livermore DYNA3D code to allow it to accurately model hypervelocityimpact is presented. The unmodifiedversion of DYNA3D was tested, results showedthat improvementswere required in three main areas: materialmodelling, the element erosion criterion, the ability to model the debris cloud formed by a thin plate impact. The SESAME Equation of State has been implemented to improve the material modelling. To improve the element erosion criterion two different criterion have been implemented and tested. The first is based on total element deformationand shows improved results for semi-infiniteand thin plate targets. The second is intendedto delete an elementwhen it becomes numericallyinaccurate. Initial results show the limitationof the Lagrangian finite element approach and that further improvementis required.
~TRODUCTION All spacecraft in orbit are exposed to impacts with meteoroids and debris. Some of these particles are large enough that an impact could cause significant damage to the spacecraft, or one of its sub-systems. Designers must be aware that impacts will occur on the spacecraft, and act to minimize the risk that an impact will damage a spacecraft system. The probability of an impact with a particle of a certain size and velocity can be predicted. The potential damage due to particle impact must also be found. A suitable hydrocode would be a useful complement to laboratory experiments for estimating impact damage, and investigating the effect of design changes. This paper covers work on developing tools to allow reliable and accurate modelling of hypervelocity impact on spacecraft structures using a Lagrangian hydrocode. The work has been based on an existing computer code, DYNA3D . Calculations were performed to evaluate the capability of DYNA3D to accurately model hypervelocity impact. Based on this experience improvements have been made to the code. The SESAME Equation of State was added. Two alternative element erosion criteria have been developed, implemented and tested in DYNA3D.
E V A L U A T I O N OF DYNA3D This work uses the computer code DYNA3D, developed by Lawrence Livermore National Laboratory. DYNA3D is an explicit non-linear finite element code, based on the Lagrangian formulation, that is widely used for modelling lower velocity impacts such as car crashes, bird strikes etc. It contains features that make it suitable for modelling hypervelocity impact, such as: bulk viscosity for shock wave modelling, an element erosion algorithm and suitable material models.
J. CAMPBELLand R. VIGNJEVIC
The first analyses used the unmodified code. This allowed experience to be gained in modelling hypervelocity impact problems and allowed the capabilities of DYNA3D to be assessed. The case modelled was normal impact of a spherical aluminium projectile on a thin aluminium plate at 8 km/s. Element erosion, using an effective plastic strain erosion criterion, was used to allow modelling of the perforation of the plate. The Steinberg-Guinan constitutive model  was used as this was the only suitable model that allowed element erosion. This material strength model was combined with a Gruneisen Equation of State. None of the calculations produced an accurate result for the final hole diameter in the plate, which in the calculations was always close to the initial projectile diameter. In addition element erosion was interfering with the formation of the impact generated shockwaves, leading to a maximum pressure behind the shockwave only 50% of the expected value. Furthermore the analyses showed that the code was not capable of successfully modelling impact on a multi-plate structure as much of the material that would form the debris cloud was eroded, Figure 1.
Fig. 1. Impact of aluminium sphere on aluminium plate at 8 km/s, 0.75 ~s after impact. Effective plastic strain erosion criterion. The major cause of the errors was the element erosion criterion, based on the effective plastic strain, which allowed elements to be eroded too easily early in the calculation. This caused an unrealistic void to be created between the projectile and plate, which interfered with the formation of the impact shockwaves. This lead to unrealistic wave propagation through the projectile and target. The problem was not affected by changing the magnitude of the erosion criterion. The cause of the problem was found to be with the implementation of the Steinberg-Guinan constitutive model. This model allows the material in an element to melt, when this occurs the shear modulus of the element is set to zero. For an element the effective plastic strain is updated each time-step using: --n+l
(s *-cr, )
ep = e p - t (]'3G+lxl0-'°'
When G, the shear modulus, is set to zero the effective plastic strain becomes very large. As a consequence the element was immediately eroded. From these initial calculations three main areas were identified where improvements to the code were necessary, these are: • material modelling • the erosion algorithm • the capability to model the material forming the debris cloud
Lagrangian hydrocode modelling for debris impact I M P R O V E M E N T OF THE M A T E R I A L MODEL
The Gruneisen Equation of State gives poor results for multiply shocked material and material in the vapour region , and so could cause inaccuracies especially when multiple plate impacts come to be modelled. There was no suitable alternative present in DYNA3D. As a consequence the SESAME Equation of State library  was implemented into DYNA3D. The advantage of SESAME library is that it contains data for a wide range of materials, and can account for material melting and vaporisation. The implementation of the SESAME Equation of State was validated by calculation of the Hugoniot curve for aluminium from DYNA results. Analyses were performed of flyer-plate impacts over a velocity range of 2-16 krn/s. The material used was SESAME 3717. The results from these analyses were compared with experimental results for aluminium [5,6], and the Hugoniot curve calculated directly from SESAME 3717, Figure 2. Good agreement was observed between the DYNA results and the curve calculated directly from the SESAME data, showing that SESAME is implemented correctly.
LASL data for AI2024  Data from Mitchell at a l  SESAME
0,25 Specific Volume crrP/g
025 Specific Vobme crrPIg
Fig. 2. Comparison of Hugoniot curve for SESAME 3717 with experimental data and DYNA results
IMPROVEMENT OF EROSION ALGORITHM
The erosion algorithm in DYNA3D required improvement, as element erosion was interfering with the formation of impact generated shockwaves. The initial analyses showed that this was due to the criterion used to determine when to erode an element, rather than with the contact algorithm. Two approaches were taken: the first was to implement an alternative element strain criterion; the second was to investigate whether there were any other suitable criteria. This section describes the new criteria, and presents initial results for them. Element Strain Criterion
A new criterion based on overall element deformation, calculated from the rate of deformation tensor, was added to DYNA. The erosion algorithm was also extended to work with the Johnson-Cook constitutive model . The new criterion was tested by modelling impacts with both thin and semi-infinite targets. For the semi-infinite target the ease modelled was of a 4.5mm diameter, 4.5ram long aluminium cylinder impacting normally at 5km/s into an aluminium target. The results could be
J. CAMPBELLand R. VIGNJEVIC
compared to published experimental and hydrocode results . Calculations were performed using different values of the erosion criterion as well as both constitutive models. Results from these calculations are shown in Table 1. The depth is the distance from the base of the crater to the original surface, the diameter is measured in the plane of the original surface. Table 1. Results from simulations of impact on a semi-infinite target
Valuefor erosion Experiment  3D Lagrangian calculation  Plastic strain erosion criterion New criterion (Steinberg-Guinan model) New criterion (Steinberg-Guinan model) New criterion (Johnson-Cook model) New criterion ~Johnson-Cook model) • Values not known
* 1.5 2.0 2.5 2.0 2.5
Crater Crater Depth (mm) Diameter(mm) 8.4 8.3 6.4 6.7 7.4 8.2 8.7
* 15 11.2 13.6 13.2 16.5 17.5
With the effective plastic strain erosion criterion the final crater is far too small. With the new criterion the results from the calculation using the Steinberg-Guinan model also show a too small crater depth. This constitutive model is not valid at the low strain rates that occur late in the calculation, so affecting the final crater size. The result from the calculation using the JohnsonCook model, for a criterion of 2.0, (Figure 3), shows good agreement with experiment in crater depth, but the final diameter is too large. It can be seen that the final results are dependent on the magnitude of the erosion criterion. The results for the new criterion are an improvement over those for the plastic strain criterion.
i1 i i i i i i i i ! ! ! ! ! ' l I H ~
24gs Fig. 3. Impact of 4.5mm diameter, 4.5ram long aluminium cylinder into aluminium target at 5km/s. Deformation erosion criterion.
Lagrangian hydrocode modelling for debris impact
For the thin plate target the case modelled was of a 4mm diameter aluminium sphere impacting a 0.8mm thick aluminium plate at 8km/s. Again an experimental result was available for comparison . Results for these calculations are shown in Table 2. Table 2. Results from simulation of impacts on a thin plate
Value for erosion
(ram) Experiment  Deformation criterion (Steinberg-Guinan model) New criterion (Johnson-Cook model) New criterion (Johnson-Cook model)
7.8 7.2 7.5 7.6
2.5 2.0 2.5
Again the calculations using the Johnson-Cook constitutive model produce better results. With the result for a criterion of 2.5 (Figure 4) giving the closest result. This is due to the Steinberg-Guinan model not being valid for the lower strain rates that occur later in the calculation. These simulations can only be used to predict the final hole diameter in the plate. By the end of the calculation most of the material that would form the debris cloud has eroded.
A 3~ts Fig.4. Impact of 4mm diameter aluminium sphere into 0.8mm thick aluminium plate at 8 km/s. Deformation erosion criterion.
Alternative Criterion Element erosion algorithms were developed and refined in the 1980's  to allow Lagrangian hydrocodes to model penetration of a projectile into a target. This type of problem cannot be modelled by a Lagrangian code without element erosion, or user intervention in the form of remeshing. This is due to large material deformation leading to serious element distortion which forces the calculation to be terminated.
J. CAMPBELLand R. VIGNJEVIC
Element erosion does not model a physical effect, it is a convenient method for avoiding the numerical problems that would otherwise occur due to large deformations. Erosion assumes that the material represented by an eroded element will no longer contribute to the physics of the problem. Experience has shown that reasonable results can be gained using this method . Highly distorted elements give rise to three main problems: • The elements can invert, leading to the calculation of a negative volume for the element. DYNA automatically terminates the calculation if this occurs. • The distortion gives rise to errors in the evaluation of the constitutive equation, especially when the element has become large and is undergoing rapid deformation. • The timestep is calculated from the smallest element dimension, distorted elements usually reduce the timestep, leading to calculations that take an impractical amount of CPU time Any useful erosion criterion will erode an element before the element inverts, and before the element becomes significantly inaccurate. It is reasonable to assume that an element about to invert will have become inaccurate. So the principal aim of an erosion criterion is to detect inaccurate elements. The solid element used in DYNA3D is the 8-node iso-parametric brick. This is a tri-linear element, integrated using 1 point Gaussian quadrature. The element assumes a linear variation in solution variables along any element axis, an assumption that requires the analyst to use a fine mesh in order to resolve rapid changes in solution variables. An improvement would be the use of higher order elements or integration rules, but this would lead to an impractical amount of computer time. The integration is performed using the element axes, in which the element is a perfect cube, the results have to be mapped from the element axes to the global axes. As 1 point integration is used the determinant of the Jacobian matrix need only be evaluated at one point. So the assumption that the Jacobian is constant over the element has been made. This assumption holds as long as the element is a parallelepiped. If there is a significant variation in the value of the Jacobian determinant across the element, the element can no longer be considered accurate. An erosion criterion has been implemented into DYNA3D which detects this condition. The code evaluates the determinant of the Jacobian matrix at eight points in the element, the eight points are equivalent to the Gauss points of a 2x2x2 integration rule. For each element the ratio between the maximum and minimum values of the determinant is calculated and compared with a ratio specified by the user. If the variation exceeds the set ratio the element is eroded. As this condition does not detect a skewed element, it has been linked with a skew criterion. If the variation of the determinant is small then the element is also checked for skew. To save CPU time only one comer angle is calculated for each element face. The user can specify a minimum skew angle. This criterion is being tested using the same cases as for the strain criterion. Only preliminary results are available at the time of writing. For the semi-infinite target case three analyses have been made. All use the Johnson-Cook constitutive model. Only the ratio between the Jacobian determinants has been varied. The skew angle for erosion is 0.5 radians in all calculations. The results are shown in Table 3. Table 3. Results for simulation for impact on a semi-infinite target
Experiment  3D Lagrangian Calculation  Jacobian Criterion Jacobian Criterion Jacobian Criterion * Value not known
1.5 2.0 2.5
8.4 8.3 7.1 7.6 7.7
* 15.0 9.7 11.6 12.8
Lagrangian hydrocode modelling for debris impact
In all cases the crater depth is less than in the experiment . The crater diameter changes significantly with the magnitude of the criterion. Results for the case with maximum ratio 2.0 are shown in Figure 5. It is clear that analyses with higher erosion values are required to properly assess the new criterion.
16gs Fig. 5. Impact of 4.5mm diameter, 4.5mm long aluminium cylinder into aluminium target at 5km/s. Jacobian erosion criterion.
As for the semi-infinite target case only the maximum ratio has been varied. The skew angle for erosion is 0.5. Results are shown in Table 4. The expansion angle is calculated from the outer node of the surviving projectile material, it is the angle between the nodes velocity vector and the centreline. The experimental value for the expansion angle is the main spray angle behind the plate. Table 4. Results for simulation of impact on a thin plate
Maximum Ratio Experiment  Jacobian Criterion Jacobian Criterion Jacobian Criterion
1.5 2.0 2.5
Hole Diameter (mm)
Expansion Angle (degrees)
7.8 6.3 7.1 7.1
29 18 19 23
It can be seen, Fig.6, that fewer elements have been eroded than with the strain criterion. Material forming the back splash can be clearly seen and a significant amount of the projectile has survived to form the debris cloud.
J. CAMPBELLand R. VIGNJEVIC
Fig. 6. Impact of 4mm diameter aluminium sphere into 0.8mm thick aluminium plate at 8 km/s. Jacobian erosion criterion.
The effect that this criterion has on the solution time was checked by analysing a model containing 2419 nodes and 1812 solid elements on a SPARC20 workstation. For these analyses the impact velocity was low, so no elements were eroded. The values for the CPU time are shown in Table 5. Table 5. CPU timings for new criteria
No Erosion Deformation criterion Jacobian Criterion
Number of timesteps
CPU time (seconds)
707 707 707
3036 3048 4468
The new Jacobian criterion has a significant effect on the CPU time, in this case it required 47% more CPU time to solve. The increase is due to the extra computation required for every element each timestep. This new criterion requires further improvement if it is to be useful or practical. The preliminary results are worse than those for the deformation criterion, and require significantly more CPU time to solve. For the case of a semi-infinite impact there is currently no advantage in using the Jacobian criterion over a strain criterion. For the case of a thin plate impact fewer elements are eroded during the impact resulting in more information about the debris cloud and back splash.
Lagrangian hydrocode modelling for debris impact
MODELLING THE DEBRIS CLOUD If a hydrocode is to be a useful tool for modelling impact on spacecraft structures, it should be capable of modelling the formation and propagation of a debris cloud formed by impact on a thin plate. This allows the code to be used to evaluate shield designs, and also to evaluate potential damage to the inside of a satellite from an impact which penetrates the outer wall. DYNA3D is not capable of modelling impact on a thin plate, the propagation of the debris cloud and the impact of the debris cloud on a subsequent plate as one continuous problem. Element erosion must be used to model the perforation of the plate and as a much of the material which would form the debris cloud is removed from the calculation. Results from any subsequent impact of the debris cloud could not be taken to be accurate. Three main options are available that would add the capability to model impact on a multiplate structure. These are: add an Eulerian capability, add an Arbitrary Lagrangian Eulerian (ALE) capability or add a gridless Lagrangian capability, such as Smoothed Particle Hydrodynamics (SPH). Adding an Eulerian or multi-material ALE capability would be an extremely complex approach. Modelling the impact on the plate using an Eulerian mesh would remove the need to delete material from the calculation. For a single material ALE formulation our experience with DYNA2D, which has this capability, shows that this is not, on its own, a sufficient solution . Implementing a gridless Lagrangian capability appears the most promising approach. It allows the modelling of large deformations without the problem of grid tangling. It has been shown that SPH can be successfully added to a Lagrangian code . Adding a gridless Lagrangian option to DYNA3D will give the code the capability to model impact on multi-plate structures.
CONCLUSION The aim of the work is to develop a Lagrangian hydrocode capable of modelling hypervelocity impact on spacecraft structures. The Lawrence Livermore DYNA3D code is being used. The main conclusions at this stage of the work are: • DYNA3D requires improvement to allow it to model hypervelocity impact accurately. Specific areas are constitutive models, element erosion, the capability to model large deformations. • The SESAME Equation of State has been implemented into DYNA to improve the material modelling • An element erosion criterion based on total element deformation has been implemented. Results show a better agreement with experiment than when using the effective plastic strain criterion. • An erosion criterion intended to erode an element when it is numerically inaccurate has been implemented, and is being tested. Preliminary results for a thin plate impact show that fewer elements are eroded than for the deformation criterion and more information is available on the debris cloud and backsplash.
REFERENCES R.G. Whirley and B.E. Engelmann, DYNA3D A Non-linear, Explicit, Three-dimensional Finite Element Code for Solid and Structural Mechanics - User Manual. UCRL-MA-107254,Lawrence Livermore National Laboratory, USA. November (1993).
J. CAMPBELL and R. VIGNJEVIC
D.J. Steinberg, S.G. Cochran and M.W. Guinan, A constitutive model for metals applicable at high-strain rate. J. Appl. Phys, 51(3) (1980). 3. J.R. Assay and G.I. Kerley, The response of materials to dynamic loading, Int. J. Impact Engng., 5, 69-99, (1987). 4. S.P. Lyon and J.D. Johnson (eds.), SESAME: The Los Alamos National Laboratory equation of state database. LA-UR-92-3407, Los Alamos National Laboratory, USA. (1992). 5. S.P. Marsh (ed.), LASL Shock Hugoniot Data, University of California Press, (1980). 6. A.C. Mitchell, W.J. Nellis, J.A. Moriarty, R.A. Heinle, N.C. Holmes, R.E. Tipton and G.W. Repp, Equation of State ofAl, Cu, Mo, and Pb at shock pressures up to 2.4 TPa (24 Mbar), J. Appl. Phys., 69(5), 2981-2986, (1991). 7. G.R. Johnson and W.H. Cook, A constitutive model for metals subjected to large strains, high strain rates and high temperatures, Proc. 7th Int. Syrup. Ballistics, The Hague, The Netherlands. (1983). 8. C.J. Hayhurst, H.J. Ranson, D.J. Gardner and N.K. Birnbaum, Modelling ofmicroparticle hypervelocity oblique impacts on thick targets, Int. J. Impact Engng., 17, 375-386, (1995). 9. J.D. Frey, F. Janicot, X. Garaud, P. Groenenboom and M. Lambert, The validation of hydrocodes for orbital debris impact simulation, Int. J. Impact Engng., 14, 255-265, (1993). 10. G.R. Johnson and R.A. Stryk, Eroding interface and improved tetrahedral element algorithms for high-velocity impact computations in three dimensions. Int. J. Impact Engng., 5, 411-421, (1987). 11. J. Campbell and R. Vignjevic, Numerical Simulation of Hypervelocity Impact on Spacecraft, DYNA3D User Group Conference, September 28-29, Nottingham, UK (1994). 12. G.R. Johnson, E.H. Petersen and R.A. Stryk, Incorporation of an SPH option into the EPIC code for a wide range of high velocity impact computations. Int. J. Impact Engng., 14, 385-394, (1993).