Land Systems Division, Norwegian Defence Research Establishment, P.O. Box 25, NO-2027 Kjeller, NorwayReceived 26 March 2015; revised 8 May 2015; accepted 12 May 2015 Available online 17 June 2015

Simulation of natural fragmentation of rings cut from warheads John F. MOXNES*, Steinar BØRVE1

Land Systems Division, Norwegian Defence Research Establishment, P.O. Box 25, NO-2027 Kjeller, Norway
Received 26 March 2015; revised 8 May 2015; accepted 12 May 2015 Available online 17 June 2015

Abstract

Natural fragmentation of warheads that detonates causes the casing of the warhead to split into various sized fragments through shear or radial fractures depending on the toughness, density, and grain size of the material. The best known formula for the prediction of the size distribution is the Mott formulae, which is further examined by Grady and Kipp by investigating more carefully the statistical most random way of portioning a given area into a number of entities. We examine the fragmentation behavior of radially expanding steel rings cut from a 25 mm warhead by using an in house smooth particle hydrodynamic (SPH) simulation code called REGULUS. Experimental results were compared with numerical results applying varying particle size and stochastic fracture strain. The numerically obtained number of fragments was consistent with experimental results. Increasing expansion velocity of the rings increases the number of fragments. Statistical variation of the material parameters influences the fragment characteristics, especially for low expansion velocities. A least square regression fit to the cumulative number of fragments by applying a generalized Mott distribution shows that the shape parameter is around 4 for the rings, which is in contrast to the Mott distribution with a shape parameter of½. For initially polar distributed particles, we see signs of a bimodal cumulative fragment distribution. Adding statistical variation in material parameters of the fracture model causes the velocity numerical solutions to become less sensitive to changes in resolution for Cartesian distributed particles.

Copyright©2015, China Ordnance Society. Production and hosting by Elsevier B.V. All rights reserved.

Keywords:Warhead; Fragmentation; Simulation; Fracture models; Expanding ring; Smooth particle

E-mail address: john-f.moxnes@ffi.no (J.F. MOXNES).

Peer review under responsibility of China Ordnance Society.

1Tel.: +47 63807676.

http://dx.doi.org/10.1016/j.dt.2015.05.005

2214-9147/Copyright©2015, China Ordnance Society. Production and hosting by Elsevier B.V. All rights reserved.

1. Introduction

Most conventional weapons contain some type of explosive charge encased in a steel metallic container. When the explosive filling of a shell or a bomb detonates, the casing is subjected to an extreme high pressure from the gaseous products of the detonation. The casing ruptures and a number of fragments that vary greatly in size are produced. During natural fragmentation the spatial shape and velocity distributions are probabilistic because of the strength of the casing and the nature of the shock wave that gives the detonation. Rigorous hydrocode calculations can offer insight into the physics of fragmentation. However, no first principle calculations can currently be used for calculating the correct fragment size distribution. This deficiency may be due to the lack of viable fracture models during high rates of strain or microstructure variations that are not accounted for. However, hydrocode accuracy may also play a role due to the stochastic or chaotic nature of the fragmentation process.

The fracture behavior of steel rings, taken from a 25 mm steel warhead is studied. To reach high strain rates around 104/ s, an expanding ring test was performed. A streak camera was used to examine the radial ring velocity, and a water tank was used to collect the fragments (Moxnes et al., 2014 [31], Moxnes et al., 2015 [32]). The fracture strain in the standard Johnson-Cook (hereafter abbreviated J-C) fracture model (1985) [25] is deterministic. Current research in the literature on fracture/failure models focus on the dependency of fracture/failure strain on triaxiality (the ratio of the invariant I1to J2) or even the third invariant, strain rate influence on ductility,element size and the connection to adiabatic shear bands during high rates of strain. Goto et al. (2008) [17] and Grady and Hightower (1990) [13] investigated the fragmentation size of explosively driven rings and cylinders. The data were used to determine relevant coefficients for the J-C fracture model. 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 (1983 [24], 1985 [25]) strength model. The parameters of the J-C fracture model were found using the results from quasi-static tensile tests on three different sample geometries (Moxnes et al., 2014 [31]). However, variations in the micromorphology of the material may lead to variations in fracture strain that may be important in fragmentation studies.

The first intriguing stochastic approach to the fragmentation problem was to investigate the statistical most random way of portioning a given topology into a number of discrete entities. Mott and Linfoot (1943) [34] referred to an earlier work of Lineau (1936) [29] for the description of fragmentation of ammunition based on this purely stochastic principle. The best known formula for the prediction of the size distribution is the Mott cumulative fragment mass formulae (Mott and Linfoot 1943 [34]; see Appendix A for details), which is further developed by Grady and Kipp (1985) [15] by investigating more carefully the statistically most random way of portioning a given area into a number of discrete entities. The derivative of the Mott formulae is a two-parameter Weibull distribution which is infinite at zero mass. The literature gives possible modifications to account for maximum and minimum fragment size. See Cohen (1980) [5], Grady and Kipp (1985) [15], Strømsøe et al. (1987) [39], Grady (1990) [14] and Baker et al. (1992) [1] for further studies in geometric fragmentation statistics. Grady (1990) [14] developed a theory for the statistical fragmentation based on the Poisson process for masses or for areas. The Weibull distribution has been used for fracture statistics for brittle materials (Lu et al., 2002) [27]. A method for analyzing the mass distribution was developed where the cumulative fragment mass is plotted as a function of cumulative fragment number, i.e. beginning with the heaviest fragment (Held and Ku¨hl 1976 [20], Held 1990 [21], Held 1991 [22]).

Metals have a microstructure whose details may create variations in material strength and strain to fracture. Rather than explicitely model the microstructure one attempts to calculate some effects of material inhomogeneity by a physical based statistical description. In a later work, Mott (1947) [33] assumed that fractures occurred at random around the circumference of the ring of a casing at a frequency governed by a strain dependent hazard function. Mott used tensile test data on steels to investigate parameters of the hazard function and provided some analytical solutions. It should be emphasized that this second model by Mott bears no relationship to the earlier Mott distribution derived by Mott and Linfoot (1943) [34]. Moreover, computer codes can provide details on fragmentation behavior, and more easily use the stochastic fracture approach as developed by Mott (1947) [33]. However, the spatial scale of the microstructure is typically of the order of micrometers and is currently not readily accessible to computational tools and resources for system level calculations. To account for microstructure physics or even adiabatic shear banding at the sub-grid level, a statistical approach may indeed be useful. A current research area is whether statistical fracture in constitutive models predicts the size distribution of fragments better than a homogeneous fracture model (Haldorsen and Moxnes 1998a [18] and 1998b [19], Glansville et al., 2010 [11], Hopson et al., 2011 [23], Meyer and Brannon 2012 [30], Rakvåg et al., 2014 [38], Moxnes et al., 2015 [32]). Hopson et al. (2011) [23] concluded by using a Eulerian code that a statistically compensated J-C fracture model substantially improved the fragmentation mass distribution for an explosive loaded cylinder. The homogeneous solution produced larger fragments in comparison to the Weibull solutions and the test data. Meyer and Brannon (2012) [30] concluded that using a Eulerian code with inherent variability in the continuum mechanics simulations lead to more realistic predictions. The statistical J-C fracture model achieved better predictions of the intermediate-sized fragments. It was concluded that proper ways to incorporate sub scale physical effects in strength and fracture models remains a subject of research. Moxnes et al. (2014) [32] show by using smoothed particle hydrodynamics (SPH) that randomness increased the number of fragments.

Another research issue is whether numerical noise is useful without any stochastic fracture model (Diep et al., 2000 [7], Diep et al., 2004 [8], Prytz and Ødegårdstuen 2011 [36], Cullis et al., 2014 [6], Moxnes et al., 2014 [31], 2015 [32]). The interesting aspect of numerical noise is that it does not require any artificial seeding of fracture sites within the material as a part of the initial conditions of the problem. However, mesh sensitivity makes results of fracture models difficult to validate (Brannon et al., 2007 [3]). Glansville et al. (2010) [11] found that mesh sensitivity was significantly reduced in the explicit Autodyn Lagrangian code when applying volume scaling in a Weibull distribution. Meyer and Brannon (2012) [30] concluded that further studies were warranted to ensure mesh independence of the predictions and accuracy in a variety of applications. Moxnes et al. (2014) [31] show by using SPH that increasing the resolution (i.e. reducing the particle size) increased the number of fragments.

Particle methods such as smoothed particle hydrodynamics (SPH) show tremendous potential for fragmentation simulations since they support both arbitrary large deformations and Lagrangian state variable tracking that avoids corruption by advection errors. SPH is a Lagrangian technique (Gingold and Monaghan 1977 [10], Lucy 1977 [28], Benz 1990 [2]) based on two main assumptions: First, an arbitrary scalar field variable can be estimated at any point in space by multiplying the variable by a suitable weight kernel and integrating over the entire simulation domain. The scale length of the weight kernel is referred to as the smoothing length, and in practice, the kernel has compact support so that the integral can be restricted to a relatively small volume. Secondly, the continuous integral is replaced by a discrete sum over a finite set of interpolation points (the particles). The gradient of the variablein question can be found by differentiating the interpolated estimate. The most well-known problems with SPH, in particular when simulating solid dynamics, is loss of stability due to tensile instability (Gray et al., 2001 [16]) and numerical fragmentation due to large particle spacing relative to the smoothing length (Liu et al., 2006 [26], Feldman and Bonet 2007 [9]). In the current work, we use an in house research simulation code called REGULUS. REGULUS is used with a state of the art handling of tensile instability (Gray et al., 2001 [16]). In addition REGULUS includes a technique to minimize numerical fragmentation referred to as regularized smoothed particle hydrodynamics (RSPH), which was originally developed to increase accuracy in shock wave modeling (Børve et al., 2005 [4]). However, in the current work, REGULUS has been used without regularization and the initial particle spacing is chosen so as to minimize the problem of numerical fragmentation. In future, the plan is to adapt the regularization technique to solid dynamics in order to reduce the problem of numerical fragmentation even further. In the current work, we let the fracture strain of SPH particles vary randomly according to a Weibull distribution, and examine the fragmentation behavior change as a function of the variance. To examine resolution dependency, the smoothing length (and particle size) is varied. The initial particle distribution is either Cartesian or polar structured. The discretization errors depend on the ratio of smoothing length to particle spacing and how the particle distribution is structured. In particular, the level of numerical noise will depend on how well the structuring of the particles matches the geometry of the problem being solved. In the current work, the initial particle distribution is either Cartesian or polar structured so that we can investigate to what extent the initial particle structuring affects the fragmentation behavior. The experimental part of the study is mostly reported in Moxnes et al. (2015) [32].

2. The experimental set up and geometrical data

Fig. 1 shows the set up. A brass tube with constant outside diameter and variable inside diameter is used to control the radial expansion velocity of the steel rings. The steel rings were manufactured from projectile bodies of in-service rounds. The test item is placed such that the expansion of the ring is perpendicular to the axis of a rotating mirror camera that is used to find the expansion velocity. The fragmentation studies were a duplication of the streak camera studies. However, in this case the fragments were collected in a water tank. To be able to repeat the actual velocity-time conditions, the tubes and rings were allowed to expand first in a thin plastic bag filled with air that was submerged underwater. Thus the expansions and break up occurred in air. The water barrel was then emptied and sieved, and the fragments collected with a magnet. More than 95% of the total mass was collected.

Fig. 1. The material locations and the geometrical set up for the expanding ring test.

The explosive is ignited at time zero at one end of the cylinder. The density of the explosive is 1.87 g/cm3. The total length of the cylinder with explosive is 10.2 cm. The length of the steel ring is 1 cm and the thickness is 0.33 cm. Two different shots (loadings) were studied numerically and experimentally by varying the thickness of the brass cylinder to achieve different expansions velocities of the steel rings. The steady state numerical velocities were found to be 190 m/s and 630 m/s and in good agreement with the measurements (Moxnes et al., 2015 [32]). The parameters of the two different loadings are seen in Table 1.

Uniaxial tensile test specimens and two 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 2. The steel is first casted, then rolled and heat-treated by quenching. Finally it is tempered. The hardness is 530 Vickers which corresponds to 5.2 GPa, or to 5.6 GPa when defined as force per projected contact area of the indenter. The tests were carried out at room temperature in a hydraulic test machine with a strain rate of approximately 5×10-4s-1(quasi-static condition). The numerical simulations of the mechanical tests were performed, assuming isotropic material properties. The results were compared with the experimental results.

The J-C (1983 [24], 1985 [25]) strength model is

In the original model, Y(εp) = A + B εpn(Johnson and Cook 1983 [24]), where εpis the plastic strain. In the current work, Y(εp) was set as a piecewise linear function of εp, as shown in Fig. 2.

.εpis the plastic strain rate and .ε*pis the nominal plastic strain rate of 1/s”. mtparameterizes the strength dependency of the temperature. Troomis the reference temperature set to 300 K and Tmeltis the melting temperature set to 1800 K. For the quasi-static tensile tests we set T = Troom. Other properties given for this steel is E = 210 GPa as the elastic modulus, ν= 0.33 as the Poisson ratio and ρ= 7850 kg/m3as the density. Strain rate parameter c and mtin equation (1) are set to zero for the quasi-static tests and as the baseline values in this article.

The J-C (1985 [25]) fracture model is, when not accounting for temperature dependency or strain rate dependency in the fracture strain, given as

Table 1Dimensions for brass cylinder and steel ring.

Table 2Steel alloy composition in percent.

where σ* is the triaxiality (negative value of pressure/Mises stress ratio).εfis the fracture strain and D is the damage variable. When D≥1 the strength of the material is set to zero. The experimental results for our steel give D1= 0.069, D2= 10.8, and D3= 4.8 (Moxnes et al., 2014 [31], 2015 [32]).

We assume that fracture of the ring is dependent of the subscale microstructure. The fracture strain of the ensemble of elements/particles making up the ring is assumed to be Weibull distributed. Meyer and Brannon (2012) apply randomness to D1+ D2according to a Weibull distribution. Here we set D1and D3as fixed, while D2is set stochastic to account for subscale microstructure. We set D2= D2/D20, where D20= 10.8 is the average value of D2equal to the experimentally found value. The distribution ρ（D2）is a Weibull distribution, to read

It is readily found that the expectation E（D2）and variance Var（D2）are

where Γ() is the Gamma function. To ensure that the expectation E（D2）equals one, we set the scale parameter a = E（D2）/Γ（1 + 1/M）= 1/Γ（1 + 1/M）for varying values of the Weibull modulus M. When M approaches infinity the solution is deterministic since the variance becomes zero.

Fig. 2. Piecewise linear yield curve as a function of plastic strain εp.

3. Results

The rings were used in a numerical study. To compare with earlier result in Moxnes et al. (2015) [32] we keep the initial values the same and neglect the acceleration phase. To keep high numerical resolution in the rings we do note simulate the complete experimental set up. Only the steady state expansion velocities of the rings were input to most of the simulations. For the low and high velocity shots, the expansion velocity was 190 m/s and 630 m/s, respectively. The smoothing length determines the accuracy of the kernel estimation assuming a continuous description. The particle spacing relative to the smoothing length determines the level of error in going from a continuous to a discrete representation. Typically, the ratio of smoothing length to particle spacing is initially chosen to be in the range 1.0-1.5. We use an initial particle spacing of 0.01 cm as baseline and the smoothing length h in all simulations is chosen to be 1.5 times the initial particle spacing. All simulations are performed in 3D with a symmetry plan normal to the axial axis and through the middle of the ring. The numerical results obtained with the REGULUS code are compared with the experimental data. Figs. 3 and 4 show the ring after 30 μs when using initially Cartesian distributed particles. Figs. 5 and 6 show the corresponding results for initially polar distributed particles. The color coding indicates the mass of individual fragments through a function found as the logarithm of fragment mass normalized by a reference mass of 1 mg.

The number of fragments is much lower for the low velocity shot than for the high velocity shot. The reduced number of fragments can be explained. The fragmentation process starts with the initiation of shear or tensile fractures at some random points. After fractures are initiated, loads decrease so stresses are not sufficient to trigger multiple fracture surfaces. However, when the same ring is deformed at high rate of strain, fragmentation number increases since a fracture that develops at one location can only influence the stress and strain at a neighboring location after a finite delay time. This delayed interaction between initiation sites provides time for crack growth at neighboring sites.

In Figs. 3-6, the statistical variation in fracture strain increases in the panels from left to right. All simulations include a certain level of numerical noise which depends on resolution and initial structuring of the particles. Cartesian distributed particles results in more irregular angular particle position compared to polar distributed particles. If we compare the leftmost panel of Figs. 3 and 5, we see that the number of fragments in the low velocity case is about 4 times larger when using Cartesian instead of polar distributed particles. When the statistical variation in fracture strain is increased, the results using the polar distributed particles start to resemble the corresponding Cartesian results. For M = 1, simulations with both types of particle distributions agree on a total fragment count of about 12-15. In the high velocity case, the difference between the Cartesian and polar results is less pronounced even when M equals infinity. This can, at least partly, be due to increased numerical fragmentation.

Fig. 3. The fragmentation pattern at 30 μs with the ring expansion velocity of 190 m/s using Cartesian distributed particles with h = 0.015 cm. The color indicates the logarithmic fragment mass normalized by a reference mass of 1 mg.

Fig. 4. The fragmentation pattern at 30 μs with the ring expansion velocity of 630 m/s using Cartesian distributed particles with h = 0.015 cm. The color indicates the logarithmic fragment mass normalized by a reference mass of 1 mg.

Fig. 5. The fragmentation pattern at 30 μs with the ring expansion velocity of 190 m/s using polar distributed particles with h = 0.015 cm. The color indicates the logarithmic fragment mass normalized by a reference mass of 1 mg.

Fig. 6. The fragmentation pattern at 30 μs with the ring expansion velocity of 630 m/s using polar distributed particles with h = 0.015 cm. The color indicates the logarithmic fragment mass normalized by a reference mass of 1 mg.

Fig. 7. The accumulated number of fragments for v = 630 m/s as a function of the fragment mass in milligrams for M = Infinity, M = 4, M = 1 and Cartesian coordinates.

Fig. 8. The accumulated number of fragments for v = 630 m/s as a function of the fragment mass in milligrams for M = Infinity, M = 4, M = 1 and polar coordinates.

The cumulative number of fragments as a function of the fragment mass is shown by the solid lines in Figs. 7 and 8 for the high velocity shot for the Cartesian and polar distributed particles, respectively. The Weibull modulus M influences the fragmentation characteristics in this case for both Cartesian and polar distributed particles, and the variability with the smoothing length h is significant. The smoothing length h = 0.08 cm (yellow curves) is clearly too large. The results for N(m) are reasonably well converged when h is smaller than 0.04 cm (orange curves). A plateau in N(m) is seen for masses less than around 300 mg. In addition to a distinct plateau, most of the numerical solutions for M = 1 and M = 4 exhibit an noticeable increase in N(m) for m < 50 mg. Mott and Linfoot (1943) [34] use a shape parameter of n = 0.5, which gives no plateau. However, the Mott (1947) [33] solution shows a plateau (see also Fig. 5 in Grady and Olsen 2003 [12]). We calculated the least square regression fit to the different numerical SPH solutions using a generalized Mott distribution (Appendix A). To achieve a plateau, the shape parameter n of the fragment distribution must be larger or equal to one. Using a single, generalized Mott distribution (see Appendix A), we usually achieve a value of n of around 4. But in some cases, in particular for some of the polar solutions plotted in Fig. 8, the sharp low-mass increase in fragment count forced the estimated n to be around 1 or less. In these cases, however, the overall fit to the numerically obtained distributions was typically poor.

We fit the numerical data to the bimodal distribution, N(m) = N1(m) + N2(m) where N1(m) and N2(m) are twogeneralized Mott distributions (Vogler et al., 2003 [40]). To simplify the calculation, we assume that fragments with masses greater than 50 mg belong to distribution 1, while the remaining fragments belong to distribution 2. In most simulations, a clear majority of fragments have masses greater than 50 mg and thus N1(m)»N2(m). Even so, the introduction of a second Mott distribution makes it possible to find a better fit to the distribution of fragments with masses larger than 50 mg. The dotted lines in Figs. 7 and 8 show the least square regression fit to the different numerical SPH solutions using the bimodal Mott distribution. In the Cartesian case, the lowmass Mott distribution is visible for M = 1. There are also signs of an extra low-mass plateau in the case of M = 4 and M = Infinity. In the polar case in Fig. 8, the low-mass Mott distribution is narrow but more visible than in the Cartesian case.

For each bimodal distribution, N(m) we achieve a set of 4 Mott parameters, (n1,n2,μ1,μ2), which defines the fitted distribution. Based on these parameters, we define the massaveraged bimodal Mott parameters for the bimodal distribution as n = (M1n1+ M2n2)/MTOTand μ= (M1μ1+ M2μ2)/ MTOT, where M1is the total mass of fragments with masses larger than 50 mg and M2is the total mass of the remaining fragments. MTOT= M1+ M2. Fig. 9 shows the mass-averaged bimodal Mott parameters n and μ for the high velocity shot as function of h in the Cartesian case. The black, green, and red curves correspond to the M = Infinity, M = 4, and M = 1 distributions, respectively. We see that n is stable around 4 for all solutions, while μ increases from around 0.5 g to around 0.8 g when h is sufficiently large. There is little or no evidence of n and μ being dependent on Weibull modulus M. Fig. 10 reveals somewhat larger fluctuations in n and μ in the polar case for larger h, which seems to be a resolution dependent variation. We observe a small difference in μ between the M = Infinity solution and the other two solution when h is small. The difference between the deterministic and the statistical solutions in μ only when using polar distributed particles indicates that the numerical noise is lower when using polar than when using Cartesian distributed particles.

Fig. 10. The shape parameter n and the scale parameter μ as a function of smoothing length h for polar distributed particles and v = 630 m/s.

Fig. 9. The shape parameter n and the scale parameter μ as a function of smoothing length h for Cartesian distributed particles and v = 630 m/s.

To check the applicability of volume scaling on mesh dependency of cumulative mass, we apply the volume scaling according to the relation of Meyer and Brannon (2012) [30] on an M = 1 simulation with polar distributed particles. The required reference volume was set to 10-5cm3. Fig. 11 shows that the scaling is not viable. The accumulated number of fragments is depended on the particle size.

The low velocity shot of 190 m/s shows a small number of fragments and a good statistical examination is difficult. Fig. 12 shows the cumulative number of fragments for h = 0.015 cm for Cartesian and polar distributed particles, respectively. The top panel shows the deterministic case with M = Infinity. The Cartesian case gives a total fragment count of 20 when M = Infinity. The polar case shows a small number of fragments of the same mass, which is not identifiable for larger h values. Applying statistically varying material properties polar distributed particles show identifiable fragments. As the Weibull modulus decreases, the total number of fragments in the Cartesian case decreases. As can be seen in the bottom panel of Fig. 12, the Cartesian and polar solutions largely agree when M = 1 with a total number of fragments of about 10-15. This is in good agreement with experimental data. The agreement between Cartesian and polar distributed particles only for M～1 is consistent with the assumption that the numerical noise is considerable larger in the Cartesian solutions than in the polar solutions.

Fig. 11. The accumulated number of fragments for v = 630 m/s as a function of the fragment mass in milligrams for M = 1 and polar coordinates. Volume scaling is used according to Meyer and Brannon (2012).

Figs. 13 and 14 show a summary of the simulated and experimental fragment count. For v = 190 m/s and the largerparticle size (h = 0.04 cm), the fragment count shows only a weak dependence on the Weibull modulus M for the Cartesian distributed particles. However, for the polar distributed particles the number of fragments increases with decreasing M. For the smaller particle size (h = 0.015 cm) the number of fragments is much the same, but decreases somewhat with decreasing M for the Cartesian distributed particles. For the polar distributed particles the fragments count once again increases with decreasing M. For v = 630 m/s, the larger particle size gives only about half the expected fragment count for the Cartesian distributed particles. Decreasing the Weibull modulus reduces the fragment count somewhat. Using polar distributed particles with h = 0.04 cm, a larger number of small fragments is achieved. The number of fragments increases with decreasing M. Using the smaller sized particles gives overall good agreement with experimental data both for the initially Cartesian and polar structured particles. The total number of fragments increases while the number of fragments larger than 100 mg decreases as Weibull modulus M decreases in the polar case. In conclusion, the Cartesian solutions fit better to the experimental data when the Weibull modulus M is assumed to be infinite, while the polar solutions fit better when M is assumed to be 4 or less. According to the simulation of the complete expanding ring experiment including the brass and the explosive by using the Impetus Afea and the Autodyn code, the duration of the main acceleration phase is in the order of microseconds (Moxnes et al., 2015) [32]. However, few reverberations are needed to get to the final velocity. It is observed in Fig. 14 that the acceleration phase of one or two micro seconds marginally influences the fragmentation pattern for the high velocity. For the low velocity the results are the same. A closer examination shows that no fragments are developed during this phase. Future research may simulate the complete set up with high resolution in the steel rings. This may increase agreement with experimental results.

Fig. 12. The accumulated number of fragments for v = 190 m/s as a function of the fragment mass in milligrams for h = 0.015 cm, M = Infinity, M = 4 and M = 1.

Fig. 13. The total number of fragments for v = 190 m/s.

Fig. 14. The total number of fragments for v = 630 m/s.

To further investigate to what extent fragmentation depends on resolution, we compare the cumulative fragment distributions for the high velocity case. First, the bimodal regression fit to the cumulative fragment solutions with h = 0.015 cm is taken as the reference. Next, we calculate the average difference, dNerr,between a given REGULUS distribution N(m) as plotted in Figs. 7 and 8, and the corresponding reference distribution. We do this by summing the squared difference between the two distributions for each discrete fragment mass value, and dividing the sum by the number of summands. The results are plotted in Fig. 15 as functions of h. The Cartesian case dNerincreases strongly when h is larger than around 0.04 cm, which can be considered as a threshold. For h < 0.04 cm, the difference between the numerical solutions and the reference solution drops more slowly as h is reduced. For the polar case a corresponding threshold value seems to be around 0.03 cm. With polar distributed particles, dNerappears to grow more quickly with h than with Cartesian distributed particles. The curves for the different values of Weibull modulus M show roughly the same increase in dNerrwith h.

4. Conclusions and discussion

The fracture behavior of the radially expanding steel rings made of a casing of 25 mm warhead was studied experimentally and numerically by using the SPH method. The parameters of a J-C strength and J-C fracture model were established using the results from tensile tests of the smoothed bar and two notched bars. The simulated expansion velocity of the velocity of rings matches a streak camera measurement.

Fig. 15. Average difference between a cumulative fragment distribution and the corresponding reference distribution as function of h.

A minimum resolution was required to achieve the simulated number of fragments consistent with the experimental results. Increasing expansion velocity of the rings increases the number of fragments. Added randomness of the material parameters influences the fragment characteristics, especially for low expansion velocity shots. Increasing particle size in SPH decreased the number of medium and large fragment size.

Mass distribution dependency on the SPH particle size is larger when using polar distributed particles than when using Cartesian distributed particles. Least square regression fits to the cumulative number of fragments by applying a generalized Mott distribution show a shape parameter around 4. For the polar distributed particles we find signs of a bimodal cumulative fragment distribution.

Meyer and Brannon (2012) [30] used a Weibull distribution to generate statistical fracture that predicts the size distribution of fragments better than a homogeneous fracture model. We performed no quasi-static tensile experiments to establish a statistical distribution and the true Weibull modulus is uncertain for the used steel material. To explain the difference between the results obtained with Cartesian distributed particles and polar distributed particles, it can be suggested that a Cartesian representation gives more numerical noise than the corresponding polar representation. This numerical noise stimulates randomness and fragmentation. We identify a minimum resolution required for simulating fragmentation characteristics with acceptable accuracy. This critical resolution corresponds to a particle size small enough to resolving the fragments. We hypothesize that decreasing particle sizes in RSPH, that increase computer time, can to some extent be avoided by adding randomness in material parameters of the fracture model. This may be more like a general principle to be examined in future research for Cartesian and polar distributed particles. However, use of numerical noise is not a controlled method that may mislead is some cases.

The need for a stochastic material model to mimic subscale physics may not be strong for our steel material. Numerical noise in the Cartesian case appeared to be useful. This is in agreement with earlier research from Diep et al. (2000 [7], 2004 [8]), Prytz and Ødegårdstuen (2011) [36], Cullis et al. (2014) [6]), and Moxnes et al. (2014 [31], 2015 [32]). A numerical computed chaotic trajectory diverges exponentially from then true trajectory in phase space with the same initial condition. There exists an errorless trajectory (no computational error) with slightly different initial condition that stays near (shadows) the numerical computed one. Thus a computational solution (with error), with no variation in the initial conditions may mimic the true solution with variations in the initial conditions or the material parameter (Ott 1993 [35]). In future research it will be useful to compare the mass distribution with the best numerical results and analytical models.

Acknowledgment

The authors appreciate the comments from Chief Scientist Ove Dullum at Norwegian Defence Research Establishment, which have improved this paper.

Appendix A. The generalized Mott distribution

The generalized Mott distribution for the cumulative number of fragments as a function of the mass of fragments is

where NTOTis the total number of fragments. This gives the particle density ρn(m) as the two parameter Weibull distribution

The total mass MTOTbecomes

where Γ is the gamma function. Thus

Some explicit examples are

The average fragment mass is defined by mdef= MTOT/NTOT. This gives the relation between average mass and scale factor and some example as

The masses distribution ρmand the cumulative mass Mc(m) are

where ΓIis the upper incomplete gamma function. Say as an alternative that the cumulative mass is given byThis gives the mass density ρmand particle density ρnas

The cumulative number of fragments becomes

References

[1] Baker L, Giancola AJ, Allahdadi F. Fracture and spall ejecta mass distribution: lognormal and multi fractal distributions. J Appl Phys 1992;72(7):2724-31.

[2] Benz W. In: Buchler JR, editor. In numerical modelling of nonlinear Stellar pulsation: problems and prospects. Dordrecht: Kluwer Academic; 1990. 1990.

[3] Brannon RM, Wells JM, Strack OE. Validating theories for brittle damage. Metall Mater Trans A 2007. http://dx.doi.org/10.11007/s11661-007-9310-7. online.

[4] Børve S, Omang M, Trulsen J. Regularized smoothed particle hydrodynamics with improved multi-resolution handling. J Comput Phys 2005;208:345-67.

[5] Cohen EA. New formulas for predicting the size distribution of warhead fragments. Math Model 1980;2:19-32.

[6] Cullis IE, Dunsmore P, Harrison A, Lewtas I, Townsley R. Numerical simulation of the natural fragmentation of explosive loaded thick walled cylinders. Def Technol 2014;10:198-2010. http://dx.doi.org/10.1016/ j.dt.2014.06.003.

[7] Diep QB, Moxnes JF,Ødegårdstuen G. A fragmentation model in 3D based on stochastic noise. In: Autodyn user group meeting, England; 2000.

[8] Diep QB, Moxnes JF, Nevstad G. Fragmentation of projectiles and steel rings using numerical 3D simulations. In: 21th Int. symp. Ballistics, 19-23 April, Adelaine, Australia; 2004.

[9] Feldman J, Bonet J. Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems. Int J Numer Methods Eng 2007;72:295-324.

[10] Gingold R, Monaghan J. Smoothed particle hydrodynamics: theory and applications to non-spherical stars. Mon Not R Astron 1977;181:375-89.

[11] Glansville JP, Barnes I, Fairlie GE. Explicit dynamics simulation of the natural fragmentation of cased explosives. In: The first international conference of protective structures, manchester; 2010.

[12] Grady DE, Olsen ML. A statistics and energy based theory of dynamic fragmentation. Int J Impact Eng 2003;29:293-306.

[13] Grady DE, Hightower MM. Natural fragmentation of exploding cylinders. In: Int. conf. on the materials effects of shock-wave and high-strain rate phenomena; 1990. p. 713-21.

[14] Grady DE. Particle size statistics in dynamic fragmentation. J Appl Phys 1990;68(12):6099-105.

[15] Grady DE, Kipp ME. Geometric statistics and dynamic fragmentation. J Appl Phys 1985;58(3):1210-22.

[16] Gray JP, Monaghan JJ, Swift RP. SPH elastic dynamics. Comput Methods Appl Mech Eng 2001;190:6641-62.

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

[18] Haldorsen Tom, Moxnes John F. Fragmentation of a steel ring by use of C-4. FFI/Report-98/03004. Forsvarets Forskningsinstitutt; 1998 [Approved for public release].

[19] Haldorsen Tom, Moxnes John F. A fragmentation model for the multipurpose round. FFI/Report-98/04940. 1998.

[20] Held M, Ku¨hl P. Consideration to the mass distribution of fragments by natural-fragmentation in combination with preformed fragments. Propellants Explos 1976;1:20-3.

[21] Held MH. Fragment mass distribution of HE projectiles. Propellants Explos Pyrotech 1990;15:254-60.

[22] Held M. Fragment mass distribution of secondary fragments. Propellants Explos Pyrotech 1991;16:21-6.

[23] 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.

[24] Johnson GR, Cook WH. A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures. In: Proc. 7th Int. Symp. Ballistics 21; 1983. p. 451-7.

[25] 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.

[26] Liu MB, Liu GR, Lam KY. Adaptive smoothed particle hydrodynamics for high strain hydrodynamics with material strength. Shock Waves 2006;15:21-9.

[27] Lu C, Danzer R, Fischer FD. Fracture statistics of brittle materials. Phys Rev E 2002;65:067102.

[28] Lucy LB. A numerical approach to the testing of the fission hypothesis. Astron J 1977;82:1013-24.

[29] Lineau CC. Random fracture of a brittle solid. J Frankl Inst 1936;221:485-94. 674-686, 769-787.

[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] Moxnes JF, Prytz AK, Frøyland Ø, Klokkehaug S, Skriudalen S, Friis E, TelandJA,DørumC,ØdegårdstuenG.Experimentalandnumericalstudyof the fragmentation of expanding warhead casings by using different numerical codes and solution techniques. Def Technol 2014;10(2):161-76.

[32] Moxnes JF, Prytz AK, Frøyland Ø, Skriudalen S, Børve S, Ødegårdstuen G. Strain rate dependency and fragmentation pattern of expanding warheads. Def Technol 2015;11(1):1-9.

[33] Mott NF. Fragmentation of shell cases. Proc Roy Soc A 1947;189:300-8.

[34] Mott NF, Linfoot EH. A theory of fragmentation. Ministry of Supply, A.C.; 1943. p. 3348.

[35] Ott E. Shadowing theorem, chaos in dynamical systems. New York: Cambridge University Press; 1993. p. 18-9 [Referenced on Wolfram｜Alpha: Shadowing Theorem].

[36] Prytz AK, Odegardstuen G. Fragmentation of 155 mm artillery grenade, simulations and experiments. In: 25th Int. J. Ballistics, September 12-22, Miami, Flodida, USA; 2011.

[38] Rakvåg KG, Børvik T, Hopperstad OS. A numerical study on the deformation and fracture modes of steel projectiles during Taylor impact test. Int J Solids Struct 2014:808-21.

[39] Stømsøe E, Ingebrigtsen KO. A modification of the Mott formula for prediction of the fragment size of distribution. Propellants Explos Pyrotech 1987;12:175-8.

[40] Vogler TJ, Thornhill TF, Reinhart WR, Chhabildas LC, Grady DE, Wilson LT, Hurricane OA, Sunwoo A. Fragmentation of materials in expanding tube experiments. Int J Impact Eng 2003;29:735-46.

* Corresponding author. Tel.: +47 63 807514, +47 63 807000; fax: +47 63 807509.