Experimental and numerical study of the fragmentation of expanding warhead casings by using different numerical codes and solution techniques

John F.MOXNES*,Anne K.PRYTZ,?yvind FR?YLAND1,Siri KLOKKEHAUG,Stian SKRIUDALEN1,Eva FRIIS,Jan A.TELAND1,Cato D?RUM1,Gard ?DEG?RDSTUEN

aLand Systems Division,Norwegian Defence Research Establishment,P.O.Box 25,NO-2027 Kjeller,Norway

bNammo Raufoss AS,P.O.Box 162,NO-2831 Raufoss,Norway

1.Introduction

When a warhead detonates,the resulting terminal effects in the target are dependent on the velocity and physical characteristics of the fragments.However,it is notable that,for blast-frag warheads,the total damage to the target is the results of blast and fragments acting together in a synergistic manner.

The initial fragment velocity is not very dependent of material properties,such as hardness,strength and ductility.It is well known that the same is true for changes in the geometrical aspects of casing design,such as 1)the use of shear-control grids on either the inner or outer surface of the casing;2)different control sizes in the dimensions of the diamond pattern grids;and 3)variations in the cross sectional profiles of the grid elements.

The main controlling parameter for initial fragment velocity from expanding warhead is the charge over mass relation of the warhead,as described by Gurney,based on the energy balance approach[1].For expanding thin walled casing,the strain rate is typically in the range of 104-105/s.The natural fragmentation process for a steel casing may start with the initiation of shear fracture at the outer or inner surface of the casing.The fracture then proceeds through the wall of the casing following the trajectories of maximum shear.However,the radial fractures could be initiated and dominate in the outer tensile region.

An elasto-plastic body that is homogeneously deformed will ultimately fail at a single random located imperfection.After failure is initiated,the static loads decrease so that the quasi-static stresses are no longer sufficient to trigger multiple fracture surfaces.However,when the same body is homogeneously deformed at high strain rates fragmentation can occur since initiation sites are isolated by wave propagation effects.A fracture that develops at one location can only influence the stress or strain at a neighbouring location after a finite delay time.This delayed interaction between initiation sites provides time for crack growth at neighbouring sites.Although the fractures resulting in fragmentation are randomly located,the fragmentation size distribution depends on the deformation,the fracture characteristics of the material and the interactions with neighbouring fractures.

In metals,the nucleation of micro voids and/or micro cracks can grow under tension.The fracture process in metals is now known to be caused by nucleation of micro voids or micro cracks,micro void and micro crack growth,and finally coalescence of voids or cracks to form macro cracks.Brittle materials form planar,circular micro cracks.Ductile metals show spherical micro voids[2-5].Under compression,micro shear bands can nucleate and grow into shear bands.

Changes through wide ranges of such properties of steel casing,such as hardness,strength and ductility,can produce the marked changes in the geometrical configuration of the fragments produced.To simulate fragmentation of,for instance,an expanding warhead,the accurate material strength and fracture models are necessary.Recent experiments on metals have shown that both pressure and the Lode[6]variable may be included in the flow curve for some materials.Bai and Wierzbicki[7]developed a model for metal plasticity and fracture with pressure and Lode dependency.The nucleation and fracture coalescence depends on pressure or triaxiality,that means on the proportion of the invariant I1to J2[8,9].It has been questioned whether triaxiality fully can describe isotropic ductility.It has been shown that the internal necking of ligament between voids that have grown significantly dominates at high triaxiality,while the internal shear localization of plastic strain ligament between voids that have experienced limited growth dominates at low stress triaxiality[10].

To quantify the influence of stress triaxiality on ductility,the different experiments on smoothed and notched bars are traditionally utilized[11].In general,the larger the triaxiality is,the smaller the strains are at fracture.This is in agreement with theoretical models for void growth[12,13].However,McClintock[12]and Johnson and Cook[14]found that the plastic strain to fracture was smaller in torsion(no triaxiality)compared to tension(larger triaxiality)for many materials.Recently,Bao and Wierzbicki[15]and Bao and Wierzbicki[16]compared the different models to cover the influence of triaxiality.They concluded that none of the models were able to capture the behaviour in the entire triaxiality range.For large triaxialities(above 0.4),void growth was the dominating fracture mode,while at low triaxialities the shear of voids dominated.The main conclusion was that there is a possible slope discontinuity in the fracture locus corresponding to the point of fracture transition[16].Bao and Wierzbicki[38]found a cut-off value of the stress triaxiality to 1/3 below which fracture never occurs.Teng and Wierzbicki[39]evaluated six fracture models in high velocity perforation.

A Lode parameter dependency has therefore been proposed.Wilkins et al.[17]noticed that the combination of hydrostatic tension and shear is related to the mechanism of ductile fracture.They concluded that the order of the applied loads,i.e.,hydrostatic load followed by shear load or vice versa,should be important in fracture modelling.

Phenomenological continuum damage models have been developed to account for failure and fracture.To account for the order of the applied loads,the cumulative damage criterion has been applied[17].It is assumed that the failure occurs at a point of the material where a weighted measure of the accumulated plastic strain reaches a critical value.The weighing function depends on the triaxiality.An appropriate weighting function is still an active field of research[7,18].In the Johnson-Cook(J-C)model used herein,an uncoupled(passive)damage evolution formulation with no Lode dependency is adopted,which entails that there is no coupling between the stress-strain behaviour and the damage evolution until a failure occurs at the critical damage.

Shear localization due to adiabatic shear band that is much smaller than the element size can soften the material.The greater the shear strain rate is,the larger is a number of these shears bands generated,and hence there is a lower stress for a given strain.This unstable thermoplastic shear occurs locally in the shear bands when the local flow stress decreases with the increase in strain.This happens when the rate of thermal softening due to the internally generated heat exceeds the rate of isothermal work hardening.The shearing deformation could even be so intense as to cause melting of the material in the bands.The temperature rise in the adiabatic shear bands should not be confused with the bulk temperature rise of the metal on the element size undergoing deformation under adiabatic conditions[19].

The Lagrangian processor,in which the numerical grid distorts with the material,is computationally fast and gives good definition of material interfaces.However,the ability of the Lagrangian processors to simulate the explosive events can usually only be enhanced by use of an erosion algorithm.The erosion algorithm works by removing Lagrangian zones which have reached a user-specific strain,typically in the order of 100-150%.The Eulerian processor,which uses a fixed grid through which material flows,is computationally much more expensive then the Lagrangian process,but is often better when modelling larger deformations and fluid flow.See Refs.[34]and[33]for use of Eulerian CTH code,Refs.[35]and[36]for use of the Arbitrary-Lagrangian-Eulerian ALE3D/CALE codes,and Ref.[37]for semi-empirical-numerical methods.

The smooth particle hydrodynamics(SPH)method is a Lagrangian technique[20].It is a gridless technique so it does not suffer from the problem associated with the Lagrangian technique of grid tangling in large deformation problems.SPH is based on two main approximations of the continuum equations.First,an arbitrary scalar field variable is described by an integral over the space that is only approximate due to the use of a kernel approximation.A smoothed kernel is used in the integral instead of the exact Dirac delta function.Second,this integral is approximated by a discrete sum of a finite set of interpolation points(the particles).In AUTODYN and LS-DYNA,SPH nodes interact with Lagrangian surfaces,and this allows the user to model the regions which undergo small deformations using the Lagrangian processor while those regions experiencing large deformations(i.e.explosive)can be modelled using SPH.The most well-known problem with SPH is loss of stability due to particle splitting.

In IMPETUS Afea,a corpuscular Lagrangian method has also been implemented.The method does not start from the continuum equations,but postulates a number of corpuscles(particles)that interact by collisions.Simulations show that they in fact follow the Maxwell's kinetic molecular theory with a Maxwell-Boltzmann distributed speed.The particles may have the same size as the SPH particles,and each particle can thus represent 1015-1020molecules.Since most fluid constitutive models start from continuum mechanics,the continuum constitutive equations must match the constitutive interaction parameters of the corpuscles.Initially the volume,the density and the chosen number of particles give the mass of the particles.The initial internal energy of the continuous fluid is matched to the kinetic and vibration/spin energy of the corpuscles.The vector components of the velocity of the corpuscles are initially set to follow the Gaussian distribution.Thus the speed of the corpuscles is Maxwell-Boltzmann distributed.The balance between kinetic energy of the particles,Wt,,and the vibration/spin of the particle,Ws,is given by the ratio Ws/Wt=(5 3γ)/(3γ 3),where γ is the ratio of the specific heats.Thus for γ=5/3,corresponding to only translatory degrees of freedom,all the internal energy of the fluid defines the kinetic energy of the particles.However,it is notable that the corpuscular theory uses two basic parameters called b and γ.They are fitted to the equation of state of the fluid and can be associated(not equal)to the co-volume and the ratio of the specific heats[21,22].The parameter b thus gives the interaction distance of the particles.

In this paper,the fracture behaviour of the expanding steel casing of a 25 mm warhead was studied,and the applicability of the J-C strength and fracture model was evaluated by comparing the simulated fracture behaviour of an expanding warhead with experiments.However,first,a quasi-static strength model of the steel was established by using a smooth uniaxial tensile test to find the von Mises flow plastic function in a J-C strength model.The parameters of a J-C fracture/failure model were found using the results from a smoothed bar and two notched bars.The fracture strain as a function of the triaxiality was fitted to the two notched tensile experiments and the smoothed tensile test by inverse modelling by comparing the experimental results with simulation results of the bars.It is commonly seen that the initial or average triaxiality is used for construction of parametric values.However,using initial triaxial values on the axis of the bar or average values over the necked region does not necessarily give the correct parameters.See also Refs.[16],[23],[24],and[18]for the study on triaxiality variations during tension test.

The numerical codes and different numerical solution techniques,such as Eulerian,Lagrangian and smooth particle hydrodynamics(SPH),and the corpuscular models recently implemented in LS-DYNA and IMPETUS Afea,were compared.

Section 2 gives the quasi-static strength and fracture model established by studying the tension of one smoothed and two notched bars.In Section 3,the fragmentation results of the expanding warhead are validated by comparing simulation results with flash X-ray results and water tank measurements of the fragmentation characteristic.Section 4 compares computer codes and solution techniques.We study the possible strain rate dependencies in the J-C strength and failure models.Section 5 presents the conclusions and discussion.

2.Strength and fracture/failure model

Uniaxial tensile test specimens and notched tensile specimens were extracted from a heat-treated steel material to establish a J-C strength and fracture model.The steel alloy composition is provided in Table 1.The steel is first casted,then rolled and heat-treated by quenching.Finally it is tempered.The hardness is 530 Vickers(5.6 GPa).

The tests were carried out at room temperature in a hydraulic test machine with a strain rate of approximately 5?104s1(quasi-static condition).Next,the numerical simulations of the mechanical tests were performed,assuming isotropic material properties-and the results are compared with the experimental data.The configuration of the steel bar is recorded by a mounted camera during the tension.The diameter was used to give the true strain(logarithmic strain calculated as Ln(A0/A),where Ln means the natural logarithm and A0is the initial area and A is the current area).In addition,the force and the engineering strain were recorded.The engineering strain was recoded using an extensometer.

Implicit FE simulations of the tensile tests were performed with the software Solid Works Simulation(www.solidswork.com)with a mesh of 1.1?105-1.4?105parabolic tetrahedral solid elements,depending on the specimen.Fig.1 shows the finite element(FE)-mesh and the geometry of the three bars that are used in the present study.Typical element length is 1.0 mm for the smooth specimen,0.2 for the notched K2 specimen,and 0.1 mm for the notched K08 specimen.The calculation time was around 24 h.The bars were given a displacement at the upper boundary.The lower boundary was fixed.The solution algorithm was quasi-static.The consistency of the numerical predictions was controlled by comparing with simulations using LS-DYNA.It was found that the two numerical codes gave insignificant differences.

Table 1Steel alloy composition.

Fig.1.Smooth specimen-Type G(left),K2(top right)and K08(bottom right).Colour plot is Von Mises stress at failure/fracture with unstressed state shown as grey.

Fig.2(a)-(c)show the comparison between the experimental results from the camera recoding with the simulations overlaid.The background image is from the experimental tests showing the last frame before failure of the specimens.The grey shadow corresponds to the initial configuration of the test specimen,and the red shadow show the configuration of the test specimen in the simulation at the same stress state as in the experiment.Fig.3 shows a photo of the fracture surface.A ductile fracture is observed.For the outer edge a shear fracture is found.

We inversely model the plastic flow curve(strength)to match the measured values of force and diameter for a smoothed bar.However,the well-known Bridgman-LeRoy correction was found to give a fast track to establish the Mises stress as a function of the equivalent plastic strain.Fig.4 shows the simulated and measured stresses as function of the strain.The Bridgeman corrected true stress gives the Mises stress which is plotted as a function of plastic strain defined by Ln(A0/A) εe,where εeis the elastic strain.The measured engineering stress vs engineering strain(measured by an extensometer)is in excellent agreement with the simulations.The measured and simulated true stress(defined as the force over the current area)is plotted as a function of the plastic strain.A small deviation is seen for the true stress.The reason for this may be that the camera recoding is not quite exact,or that the diameter is not simulated exactly.

The J-C(1985)strength model is

where Y1(εp)is a piecewise linear function of the plastic strain;c is the strain hardening parameter;Troomis the reference temperature set to 300 K;and Tmeltis the melting temperature set to 1800 K.The heat capacity is set to 500 J/kg.E is the Young's modulus,ν is the Poisson ratio,and ρ is the density.For the tensile tests we set that T=Troomand the strength hardening constant c is set to zero.

The well-known current J-C fracture/failure model does not account for Lode dependency.The Lode parameter μ =(2σ2σ1σ3)/(σ1σ3)= 1 for all three tensile tests since σ2σ3for the principal stresses.So any Lode angle parameter dependency would not be found anyway.The damage development is given as

Fig.2.The configuration of the test specimen according to the camera recording compared with experimental results.

Fig.3.The fracture surface of the smoothed bar.

where σ*is the so-called triaxiality(the negative value of pressure over the Mises stress).We assume that the bars start to fail at the centre of the neck on the axis.Inverse modelling is used to find the parameters D1,D2,D3such that damage equals 1 at failure for the three bars as seen in Fig.5.D4is set to zero.At failure the simulated cross section area,force and engineering stress equal the measured values.For the two notched bars we did not know the area or the engineering stress at failure since the camera was not good enough and the extensometer was not useful.So only the force at failure was used to calibrate the failure model for these two notches.Results are shown in Figs.5 and 6.

Fig.3 shows a sheared fractured region at the outer edge of the specimen.It is notable that we were not able to simulate this sheared fractured region in the tension test.The simulation of the stretched bar G shows a straight plane fracture surface when the J-C fracture model with our parameters is used.The reason for this could be related to the strength surface that may not be able to trigger such sheared region,or to the fracture model.

3.Validation

Fig.7(a)and(b)shows the different sections of the warhead that we study.The 25 mm warhead consists of a steel casing,tracer,fuze,incendiary(Zr)and explosive.The warhead also has a nose that is modelled as a steel part.

Figs.8 and 9 show the experimental results for the steel casing and base where fragments are collected in the water tank.Figs.8 and 9 shows two different shots.

All the simulations of the expanding warhead were performed by using the quarter symmetry.The material properties for the casing(shell body)are from Section 2 and we use quasi-static values for all parameters.Strain rate and temperature may influence the results.We have not measured any of the parameters that relate to strain rate and temperature effects.We prefer to keep the quasi-static parameters as our baseline parameters.Our baseline parameters are c=0,T=Troom,D1=0.069,D2=10.80,D3=4.8,D4=0.

Fig.4.Simulated and measured stresses vs strain.

Tracer and explosive are both modelled by the SPH solver.For the steel casing we always apply the Lagrangian solver.We also performed some simulations applying SPH in the steel casing but the results were not very good since Sections A,C and D showed instability that gave particle splitting.

Simulations are compared to experiments using flash X-ray measurements.We first scale the experimental pictures so that the simulated diameter of the nose of the steel casing fit the simulations.Due to an uncertainty of the triggering time for the experimental pictures,we chose simulation time by matching the outer radius of the fragments in zone B with the measurements.The velocity of these fragments should be mostly related to the charge over mass relation of the warhead as described by Gurney based on the energy balance approach[1].

Fig.10a,b and c show AUTODYN 3D simulations compared to flash X-ray measurements applying the baseline parameters.The simulations fit the measurements quite well,although discrepancies are observed in Sections A,C and D where no fractures are seen in the simulations but are evident when looking at the findings in Fig.8 and flash X-ray in Fig.10c.It should be noted that these fractures are hardly seen on the flash X-ray pictures.

We erode the elements at failure.However,due to grid tangling it was also applied global erosion criteria for strain equal to 1.We find it physically reasonable that the cell should not erode before failure.However,the visual inspection of the grid shows that some cells close to the corner section of the base erode before damage.

Fig.5.The J-C fracture/failure function and the plastic strain at the centre of the neck as a function of the triaxiallity at the centre of the neck for the three bars.D1=0.069,D2=10.80,D3=4.8.

Fig.6.The damage as a function of plastic strain for the three bars.

Fig.11 shows the cumulative number of fragments after collecting the fragments in the water tank.We plot the cumulative number N(m)as the number of fragments which is larger than m,where m is the mass of the fragments in grams.The mass of all fragments is 35 g in Region B.The total mass of the steel casing is 125.6 g.The smallest fragments in the experiment are not counted and correspond to a loss of 0.4%in Section B.We observe that the simulations give too small number of the medium fragments for IMPETUS Afea,although the overall agreement is fairly good.

The simulations in Fig.12 show that the damage of the base depends on the solver.Fig.12(1)shows the solution when applying the SPH for the tracer explosive and the incendiary.We observe that a part of the base is eroded and we find no sign of the fragmentation pattern observed in Fig.9.Fig.12(2)shows the result where SPH is used for HE and incendiary only.The base fractures along the lines of symmetry which we believe is a numerical artefact.If the Eulerian solver is use for the tracer,explosive and incendiary,the base is too much damaged(Fig.12(3)).The simulation in Fig.12(4)is similar to the simulation in Fig.12(1),but with SPH also for the fuze.We do not know the reason for this difference.In Fig.12(5),the Eulerian solver is use for the explosive and the incendiary.The damaged zone seems to be somewhat too large.Fig.12(6)shows the LS-DYNA result which is somewhat different since small fragments appear in the centre.Fig.12(7)shows the results when applying the IMPETUS Afea with baseline parameters.All parts are modelled using standard Lagrangian solver,except for the HE which is modelled using discrete particles.The base is too much eroded,mostly due to material failure.The result is very similar to the result obtained by AUTODYN.Fig.12(8)shows the results when applying the IMPETUS Afea with the node-splitting technique(to be discussed in the last section),which is an alternative to the erosion at material failure technique.As a result,we now do not see any erosion but only fragmentation.However,it still does not show the fragmentation pattern in Fig.9.Instead we see more or less a complete pulverization of the base.This is somewhat in contrast to the well description of the fragmentation in Section B as will be shown later.

Fig.7.The different sections of the warhead.(a)Projectile section definitions,(b)Baseline setup for a quarter model of the projectile in AUTODYN.

We examined the base fragments obtained from the experiments more carefully and found the sign of adiabatic shear bands,as shown in Fig.13.It appears that the larger fragments of the base are caused by cleavage.A 45?shear fracture is seen.Voids are not clearly seen,and a smooth fracture surface appears with a clear evidence of brittle fracture and a highly localized shear zone.In general,it appears that the fragmentation process of the steel casing is due to a combined ductilebrittle fracture process since both dimples and cleavage are present.

Fig.8.The collected fragment sizes in a water tank.

We conclude that the simulations do not fit exactly to the experiments.Regions A,C and D should be more fragmented compared with Figs.8 and 9.The fragmentation of the base is not correctly simulated.

Fig.9.The masses of the base fragments.

Fig.10.Flash X-ray pictures compared to AUTODYN simulations using baseline parameters.(a)25μs,(b)50μs,(c)100μs.

Fig. 11. Fragment count in Section B. Theoretical function:N（m）=N0Exp（ （m/μ）n）n=0.5，N0=350.

4.Comparison between different numerical models,solution techniques and fracture/failure models

Different codes may give different solutions due to different numerical approximation techniques.We compare AUTODYN,LS-DYNA and IMPETUS Afea for our set of baseline parameters.We use the very same grid and the same number of particles.The tracer is modelled in Eulerian,Lagrange,and SPH.The explosive is modelled in Eulerian,SPH and corpuscles.

4.1.Comparison between LS-DYNA and AUTODYN

We initially had a problem with the contact algorithm of LS-DYNA.We have compared LS-DYNA with AUTODYN using SPH for the explosive and the incendiary.Fig.14 shows that many SPH particles leak out of the casing compared to the AUTODYN simulations.However,the overall agreement is quite good.It should be noted that the particles used in LSDYNA is 4 times of those in AUTODYN.

4.2.Comparison between IMPETUS and AUTODYN

IMPETUS has no Eulerian solver,hence the Lagrangian solver is used for the tracer and the fuze.We apply SPH for the explosive in AUTODYN and the corpuscles for the explosive in IMPETUS.The comparison can be seen in Fig.15,and we find that the solutions are similar.

4.3.Different AUTODYN solvers

Fig.16 shows the time series for different solvers in AUTODYN.We find that the results are similar although differences are observed.The Eulerian processor needs about 5 time longer computation time.

4.4.Applying different parameters in the J-C fracture mode in AUTODYN

4.4.1.Strength and strain rate

Fig.12.The base of the projectile at 15 ms.1:SPH for tracer/HE/Zr;2:SPH for HE/Zr;3:Euler for tracer/HE/Zr;4:SPH for tracer/HE/Zr/fuze;5:Euler for HE/Zr;6:LS-DYNA,SPH for tracer/HE/Zr;7:IMPETUS,corpuscles for HE;8:IMPETUS,corpuscles for HE and node split.Lagrange is used for the rest of the parts.

Fig.13.45?shear zones observed for the base fragment.

Strain rate dependency has been observed for steel materials.For strain rates up to 103-104/s,a linear logarithmic relationship has traditionally been used[14].It is here believed that the dislocation motion is controlled by thermal activation.Above the rate of 103-104/s,the strength of materials are usually significantly enhanced[14,25,26].(Split Hopkinson bare experiments usually are un-useful in this range).The reason for this may be a change of the microstructure rate controlling mechanisms at high strain rates.Enhanced strain rate dependency may be due to resistance to dislocation motion in the lattice itself by phonon viscosity[27].The stress is here found to be linearly proportional to the strain rate[28].However,the conflicting evidence has been obtained regarding the postulated change from thermally activated to viscous damping mechanisms at high strain rates.Quite many different constitutive models,which are partly microstructural,have been developed(see Rohr et al.[25]and references therein).To improve our fracture simulations,the strength should probably decrease and not increase.Shear localization due to adiabatic shear band can soften the material.The greater the shear strain rate is,the larger is the number of these shears bands generated,and hence a stress for a given strain is lower.In general,the shearing deformation could even be so intense as to cause the melting of the material in the bands.

4.4.2.Ductility and strain rate

It has been reported that the ductility increase with the strain rate[14].However,high alloyed steel that achieves its strength by precipitation hardening process can behave differently with increased strength and decreased ductility with strain rate.High strength austenitic steel can also show significant strengthening effect with increasing strain rate,while the ductility decreases.Here again,a decrease in ductility at the higher strain rates can be explained by shear localization due to adiabatic shear bands[29,40,41].Although we find clear indications of shear bands in the base of the projectile,they are not clearly visible at any other places.In the base,45?shear zones were observed for the base fragments(Fig.13).

Fig.17 compares the simulation results and show the flash X-ray pictures.We fit the D4 parameter in the J-C fracture model to our warhead fragmentation experiments.D4 is chosen to be negative.This means that the strain at failure/fracture decrease with strain rate.Increasing the failure/fracture strain with strain rate(i.e.positive D4)worsen the simulations results.We find much better agreement with experimental results in Region A.However,still Sections C and D do not fit.We believe that this result is applicable for all codes and all solvers.

We next study the effect of strain rate on strength by the strain rate factor shown in Fig.18.That means that we multiply the quasi-static J-C strength model with this function of strain rate to account for the linear logarithmic relationship in strength for strain rates up to 103-104/s and the significant enhancement in strength above the strain rate of 103-104/s[14,25,26].

We find that the simulation results are influenced as seen in Fig.19.However,we could not fit any of the simulations results to the experiments since we never found a time where the length and the diameter nose of the projectile and the expansion of Section B were in agreement.We conclude that the strain rate enhancement factor in Fig.18 is not useful for our steel material.

Fig.14.Comparison between LS-DYNA and AUTODYN.

Fig.15.Comparison between IMPETUS and AUTODYN.

We also studied the influence of temperature in the J-C strength model but the difference was rather insignificant(not shown in the figures).

Finally we applied the IMPETUS Afea node-splitting algorithm(Fig.20).Instead of eroding cells that fails,the nodes can split,resulting in a sort of crack propagation.These cracks are constrained by the mesh,or cell,size.The erosion by geometric strain was set to 2.Due to the node-splitting,the mass of the casing is more or less preserved.Only about 0.2%is eroded due to geometric strain(which is basically used to avoid the numerical problems).The fracturing of the casing is thus solely due to the node-splitting algorithm.The nodesplitting algorithm implemented in IMPETUS Afea shows very promising results,although the fragment count does not fit exactly(Fig.11).

5.Conclusions and discussion

We studied the fracture behaviour of the steel casing of 25 mm warhead experimentally and numerically.The parameters of a J-C strength and fracture model were found using the results from tensile tests of the smoothed bar and two notched bars.

The simulation of the fracture behaviour was very promising but differences with experiments were observed.

Increase in strength with strain rate according to a suggested relationship that also accounts for increasing strength due to phonon viscosity for strain rates above 103-104/s worsens the simulations.

In an attempt to improve the simulation,we assumed that strain rate decreases the ductility.This improved the simulations.

Fig.16.The casing with damage plot and mesh at 10,20 and 100μs.1:SPH for tracer/HE/Zr;2:SPH for HE/Zr;3:Euler for tracer/HE/Zr;4:SPH for tracer/HE/Zr/fuze;5:Euler for HE/Zr.Lagrange is used for the rest of the parts.

We compared the numerical codes and the different numerical solution techniques,such as Eulerian,Lagrangian,smooth particle hydrodynamic(SPH),the corpuscular technique and the node-splitting technique as implemented in numerical code IMPETUS Afea.For the same solution technique and material model,we find that the codes gave similar results.The SPH technique or corpuscular technique is superior to the Eulerian technique and the Lagrangian technique(with erosion)when we apply it to the materials that have fluid like behaviour,such as the explosive or the tracer.The Eulerian technique simply gives too much computational time,and both the Lagrangian and the Eulerian technique also seem to give less agreement with the measurements.We achieved the best simulations and the most efficient computer time by applying SHP technique or the corpuscular technique in the tracer and the explosive,while keeping the Lagrangian solver for the steel casing.

The node-splitting algorithm implemented in IMPETUS Afea was very promising.It allows gases to go through the Lagrangian grid of the steel casing without erosion.Of the different technique tested,the node-splitting algorithm gave the best results in predicting the fracturing of the casing.We mention that LS-DYNA has implemented a corpuscular method for airbag deployment simulations,but the code also has a more general particle method,discrete element method(DEM),that appears to be well suited for simulation of multiphysics problems including interaction of bonded and loose particles.The applicability of the implemented DEM in LSDYNA for simulation of fragmentation of expanding warheads will be investigated in a further work.

Simulation of the stretched smoothed bar showed a straight plane fracture surface when applying the J-C fracture model with our parameters.However,the tensile test showed a sheared fractured region at the outer edge of the specimen.The reason for the discrepancies probably relates to the strength or to the fracture model itself.We do not know whether this weakness of the mathematical model could explain why the fragmentation pattern of the base and Regions A,B and C did not fit exactly to the experiments.A Lode dependency of the fracture model could also be an explanation[15,16].

Fig.17.Flash X-ray pictures compared to AUTODYN simulations using modified J-C parameters.(a)25μs,D4= 0.077,(b)50μs,D4= 0.077,(c)100μs,D4= 0.077.

Fig.18.Strain rate enhancement factor on strength.

Fig.19.Flash X-ray pictures compared to AUTODYN simulations using modified J-C parameters.(a)25μs,strength increased by strain rate.(b)50μs,strength increased by strain rate.

We did find shear failure in the base.The ratio of the tensile outer region(y)to the thickness of the steel casing(h)was given by Taylor to be y/h=Y/p,where Y is the yield stress of the casing and p is the inner pressure[31].For Taylor's hypothesis,the compressive tangential stress decreases from maximum at the inner surface to zero at the failure surface.Thus the radial fractures of the inner surface occur when the expansion of the casing is such that p=Y(corresponding to zero hoop stress at the inner surface).The tensile fractures will not propagate into the inner compression zone.However,here shear fractures approximately at 45?from the radial direction will prevail.Unstable adiabatic shear indeed transfers the entire burden of strain to a finite number of these shear planes(adiabatic shear bands)in the deformation stage.The inability to simulate shear band could explain why the simulations did not exactly fit to the experiments.In an attempt to improve the simulation,we assumed that strain rate decreases the ductility.It is notable that the decrease in failure strain with strain rate may not be realistic.Decreasing the failure strain with strain rate is likely a numerical method that accounts for shear bands that applies on a scale below the grid size.Due to restriction on computational time the grid was simply too coarse to resolve the shear bands by direct simulation.

Fig.20.Flash X-ray pictures compared to AUTODYN simulations when applying IMPETUS and node split.(a)25 μs,(b)50μs.

Meyer and Brannon[30]used a Weibull distribution to generate statistical failure for predicting the size distribution of fragments better than a homogeneous failure model.Applying this approach to our warhead may enhance the fragmentation of the steel casing.However,since we have performed no quasi-static tensile experiments to establish such statistical distribution and we did not use this approach.

Finally we mention Bridgman's studies related to the Second World War(WW II)and the post WW II eras that are summarized in Ref.[32].Bridgman demonstrated that,when high hydrostatic pressure is superimposed on the steel tensile specimen,the ductility of the specimen may be increased by several hundred folds as compared to a similar specimen pulled to fracture at the atmospheric pressure.It is notable that the stress in a steel tensile specimen simultaneously subjected to a high hydrostatic pressure is also essentially realized in the wall of an expanding steel casing,which is the most notable near the inner surface of the casing.The Bridgman effect may provide a natural check seal in the zone of ductile behaviour in the inner portion of the cylindrical wall where pressure is positive.Even though the outer portion of the steel casing may be fractured in tension,the inner ductile portion may still expand intactly and thus contain the detonation products.This Bridgeman effect may explain why the simulated expansion velocity and the fracture is too low in Sections A,C and D.

In general,we believe that our research gives a better understanding of the warhead case behaviour.Velocity-fragment size contribution can be factored better into the synergistic aspect of the total damage mechanisms of a warhead.

[1]Jones GE,Kennedy JE,Bertholf LD.Ballistics calculations of R.W.Gurney.Am J Phys 1980;48(4):264-9.

[2]Tvergaard V.Influence of voids on shear band instability under plain strain conditions,Int.J Fract 1981;17:389-407.

[3]Tvergaard V.On localization in ductile materials containing spherical voids.Int J Fract 1982;18:237-52.

[4]Tvergaard V.Material failure by void growth.Adv Appl Mech 1990;27:83-151.

[5]Tvergaard V.Shear deformation of voids with contact modeled by internal pressure.Int J Mech Sci 2008;50:1459-65.

[6]Lode W.Versuche uber den Einfluss der mittleren Haupspannung auf das Fliessen der Metalle Eisen Kupfer und Nickel.Z Phys 1926;36:913-39.

[7]Bai Y,Wierzbicki T.A new model of metal plasticity and fracture with pressure and lode dependency.Int J Plast 2010;24:1071-96.

[8]Van Stone RH,Cox TB,Low JR,Psioda JA.Micro structural aspects of fracture by dimpled rupture.Int Met Rev 1985;30:157-79.

[9]Garrison Jr WM,Moody NR.Ductile fracture.J Phys Chem Solids 1987;48:1035-74.

[10]Parodoen T,Brechet T.Influence of microstructure driven strain localization on the ductile fracture of metallic alloys.Philos Mag2004;84:269-97.

[11]Hancock JW,Mackenzie AC.On the mechanisms of ductile failure in high-strength steels subjected to multi-axial stress states.J Mech Phys Solids 1976;24:147-69.

[12]McClintock FA.A criterion for ductile fracture growth of holes.J Appl Mech 1968;35:363.

[13]Rice JR,Tracey DM.On the ductile enlargement of voids in triaxial stress fields.J Mech Phys Solids 1969;17:201-17.

[14]Johnson GR,Cook WH.Fracture characteristics of three metals subjected to various strains,strain rates,temperatures and pressures.Eng Fract Mech 1985;21(1):31-48.

[15]Bao Y,Wierzbicki T.A comparative study on various ductile crack formation criteria.J Eng Mater Tech 2004;126:314-24.

[16]Bao Y,Wierzbicki T.On fracture locus in the equivalent strain and stress triaxiality space.Int J Mech Sci 2004;46:81-98.

[17]Wilkins ML,Streit RD,Reaugh JE.Cumulative-strain damage model of ductile fracture:simulation and prediction of engineering fracture tests.UCRL-53058,R3763;October 1980.

[18]Gruben G,Hopperstad OS,B?rvik T.Evaluation of uncoupled ductile fracture criteria for the dual-phase steel docol 600DL.Int J Mech Sci 1985;62:133-46.

[19]Manion SA,Stock TAC.Measurement of strain in adiabatic shear bands.J Australas Inst Met 1969;14(3):190-1.

[20]Benz W.In:Buchler JR,editor.Numerical modelling of nonlinear stellar pulsation:problems and prospects.Dordrecht:Kluwer Academic;1990.

[21]Olovsson L,Hanssen AG,B?rvik T,Langseth M.A particle-based approach to close-range blastloadings.EurJMech A Solids 2010;29(1):1-6.

[22]B?rvik T,Olovsson L,Hanssen AG,Dharmasena KP,Hansson H,Wadley HNG.A discrete particle approach to simulate the combined effect of blast and sand impact loading of steel plates.J Mech Phys Solids 2011;59:940-58.

[23]Wierzbicki T,Bao Y,Lee YW,Bai Y.Calibration and evaluation of seven fracture models.Int J Mech Sci 2005;47:719-43.

[24]Gruben G,Fagerholt E,Hopperstad OS,B?rvik T.Fracture characteristics of a cold-rolled dual-phase steel.EurJMech A Solids 2011;30:204-18.

[25]Rule WK,Jones SE.A revised form for the Johnson-cook strength model.Int J Impact Eng 1997;21(8):609-24.

[26]Rohr I,Nahme H,Thoma K.Material characterization and constitutive modelling of ductile high strength steel for a wide range of strain rates.Int J Impact Eng 2005;31(4):401-33.

[27]Regazzoni G,Kocks UF,Follansbee PS.Dislocation kinetics at high strain rates.Acta Metall 1987;35(12):2865-75.

[28]Ferguson WG,Kumar A,Dorn JE.Dislocation damping in aluminum at high strain rates.J Appl Phys 1967;38:1863.

[29]Johnson GR,Hoegfeldt JM,Lindhom US,Nagy A.Response of various metals to large torsional strains over a large range of strain rates-part 2:less ductile metals.J Eng Mater Technol 1983;105(49):48-53.

[30]Meyer HW,Brannon RM.A model for statistical variation of fracture properties in a continuum mechanics code.Int J Impact Eng 2012;42:48-58.

[31]Taylor GI.Scientific papers of G.I.Taylor,vol.III.Cambridge:Cambridge University Press;1963.No.44.

[32]Bridgman PW.Studies of plastic flow and fracture.Cambridge,MA:Harvard University Press;1964.

[33]Hopson MV,Scott CM,Patel R.Computational comparison of homogeneous and statistical description of AerMet100 steel subjected to high strain rate loading.Int J Impact Eng 2011;38:451-5.

[34]Vogler TJ,Thornhill TF,Reinhart WD,Chhabildas LC,Grady DE,Wilson LT,et al.Fragmentation of materials in expanding tube experiments.Int J Impact Eng 2003;29:735-46.

[35]Wang P.Modeling material responses by arbitrary Lagrangian Eulerian formulation and adaptive mesh refinement method.J Comp Phys 2010;229:1573-99.

[36]Goto DM,Becker R,Orzechowski TJ,Springer HK,Sunwoo AJ,Syn CK.Investigation of the fracture and fragmentation of explosive driven rings and cylinders.Int J Impact Eng 2008;35:1547-56.

[37]Gold VM,Baker EL.A model for fracture of explosive driven metal shells.Eng Fract Mech 2008;75:275-89.

[38]Bao Y,Wierzbicki T.On the cut-off value of negative triaxiality for fracture.Eng Fract Mech 2005;72:1049-69.

[39]Teng X,Wierzbicki T.Evaluation of six fracture models in high velocity perforation.Eng Fract Mech 2006;73:1653-78.

[40]Batra RC,Ko KI.Analysis of shear bands in dynamic axisymmetric compression of a thermoviscoplastic cylinder.Int J Sci 1993;31:529-47.

[41]Bonnet-Lebouvier A-S,Molinary A,Lipinski P.Analysis of the dynamic propagation of adiabatic shear bands. Int J Solids Struct 2002;39:4249-69.