Fragmentation processes in impact of spheres
FFragmentation processes in impact of spheres
H. A. Carmona , , F. K. Wittel , F. Kun , and H. J. Herrmann , Centro de Ciˆencias e Tecnologia, Universidade Estadual do Cear´a, 60740-903 Fortaleza, Cear´a, Brazil Computational Physics IfB, HIF, ETH, H¨onggerberg, 8093 Z¨urich, Switzerland Department of Theoretical Physics, University of Debrecen, P. O. Box:5, H-4010 Debrecen, Hungary and Departamento de F´ısica, Universidade Federal do Cear´a, 60451-970 Fortaleza, Cear´a, Brazil
We study the brittle fragmentation of spheres by using a three-dimensional Discrete ElementModel. Large scale computer simulations are performed with a model that consists of agglomeratesof many particles, interconnected by beam-truss elements. We focus on the detailed development ofthe fragmentation process and study several fragmentation mechanisms. The evolution of meridionalcracks is studied in detail. These cracks are found to initiate in the inside of the specimen withquasi-periodic angular distribution. The fragments that are formed when these cracks penetratethe specimen surface give a broad peak in the fragment mass distribution for large fragments thatcan be fitted by a two-parameter Weibull distribution. This mechanism can only be observed in3D models or experiments. The results prove to be independent of the degree of disorder in themodel. Our results significantly improve the understanding of the fragmentation process for impactfracture since besides reproducing the experimental observations of fragment shapes, impact energydependence and mass distribution, we also have full access to the failure conditions and evolution.
PACS numbers: 46.50.+a, 62.20.Mk, 05.10.-a
I. INTRODUCTION
Comminution is a very important step in many indus-trial applications, for which one desires to reduce theenergy necessary to achieve a given size reduction andminimizing the amount of fine powder resulting fromthe fragmentation process. Therefore a large amountof research has already been carried out to predict theoutcome of fragmentation processes. Today the mecha-nisms involved in the initiation and propagation of singlecracks are fairly well understood, and statistical modelshave been applied to describe macroscopic fragmenta-tion [1, 2]. However, when it comes to complex frag-mentation processes with dynamic growth of many com-peting cracks in three-dimensional space (3D), much lessis understood. Today computers allow for 3D simula-tions with many thousands of particles and interactionforces that are more realistic than simple central poten-tials. These give a good refined insight of what is reallyhappening inside the system, and how the predicted out-come of the fragmentation process depends on the systemproperties.Experimental and numerical studies of the fragmenta-tion of single brittle spheres have been largely applied tounderstand the elementary processes that govern com-minution [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,18, 19, 20, 21, 22]. Experiments that were carried out inthe 60s analyzed the fragment mass and size distributions[3, 4, 5] with the striking result that the mass distribu-tion in the range of small fragments follows a power lawwith exponents that are universal with respect to mate-rial, or the way energy is imparted to the system. Laterit became clear that the exponents depend on the di-mensionality of the object. These results were confirmedby numerical simulations that were mainly based on Dis-crete Element Models (DEM) [18, 23, 24, 25, 26]. For large fragment masses, deviation from the power law dis-tribution could be modeled by introducing an exponentialcut-off, and by using a bi-linear or Weibull distribution[12, 13, 15, 27, 28, 29, 30]. Another important findingwas, that fragmentation is only obtained above a certainmaterial dependent energy input [5, 6, 31]. Numericalsimulation could show that a phase transition at a crit-ical energy exists, with the fragmentation regime above,and the fracture or damaged regime below the criticalpoint [18, 19, 21].The fragmentation process itself became experimen-tally accessible with the availability of high speed cam-eras, giving a clear picture on the formation of the frag-ments [5, 6, 7, 8, 9, 10, 11, 12, 13]. Below the criticalpoint, only slight damage can be observed, but the spec-imen mainly keeps its integrity. Above but close to thecritical point, the specimen breaks into a small number offragments of the shape of wedges, formed by meridionalfracture planes, and additional cone-shaped fragments atthe specimen-target contact point. Way above the criti-cal point, additional oblique fracture planes develop, thatfurther reduce the size of the wedge shaped fragments.Numerical simulations can recover some of these find-ings, but while two-dimensional simulations cannot re-produce the meridional fracture planes that are respon-sible for the large fragments [14, 16, 17, 18, 20, 21], three-dimensional simulations have been restricted to relativelysmall systems, and have not focused their attention onthe mechanisms that initiate and drive these meridionalfracture planes [15, 19]. Therefore, their formation andpropagation is still not clarified, although the resultingtwo to four spherical wedged-shaped fragments are ob-served for a variety of materials and impact conditions[5, 8, 11, 20]. Arbiter et al. [5] argued, based on the anal-ysis of high speed photographs, that fracture starts fromthe periphery of the contact disc between the specimen a r X i v : . [ c ond - m a t . s t a t - m ec h ] N ov and the plane, due to the circumferential tension inducedby a highly compressed cone driven into the specimen.However, their experiments did not allow access to thedamage developed inside the specimen during impact.Using transparent acrylic resin, Majzoub and Chaudhri[8] observed damage initiation at the border of the con-tact disc, but in their experiments plastic flow and ma-terial imperfections may have a dominant role.In this paper we present three-dimensional simulationsof brittle solid spheres under impact with a hard plate.With our simulations, the time evolution of the fragmen-tation process and stress fields involved are directly ac-cessible. We have focused our attention on the processesinvolved in the initiation and development of fracture,and how they lead to different regimes in the resultingfragment mass distributions. Our results can reproduceexperimental observations on fragment shapes, impactenergy dependence, and mass distributions, significantlyimproving our understanding of the fragmentation pro-cess in impact fracture. II. MODEL AND SIMULATION
Discrete Element Models (DEM) have been success-fully used since they were introduced by Cundall andStrack to study rock mechanics [32]. Applications rangefrom static, to impact and explosive loading, usingelementary particles of various shapes that are con-nected by different types of massless cohesive elements[14, 17, 18, 19, 20, 24, 33, 34, 35, 36, 37, 38]. In gen-eral, Newton’s equation governs the translational androtational motion of the elements, that concentrate thewhole mass. Forces and torques arise from element in-teractions, from the cohesive elements, volumetric forces,and of course from interaction with boundaries like walls.Throughout this work we use a three-dimensional (3D)implementation of DEM where the solid is representedby an assembly of spheres of two different sizes. Theyare connected via beam-truss elements that deform byelongation, shear, bending, and torsion. The total forceand moment acting on each element consists of thecontact forces resulting from sphere-sphere interactions, (cid:126)F c = (cid:126)F ov + (cid:126)F diss , the stretching and bending forces (cid:126)F b = (cid:126)F elo + (cid:126)Q and moments (cid:126)M b transmitted by thebeams attached.The contact force has a repulsive term due to elastic in-teraction between overlapping spherical elements, whichis given by the Hertz theory [39] as a function of the ma-terial Young’s modulus E p , the Poisson ratio ν p , and thedeformation ξ . The force on element j at a distance (cid:126)r ij relative to element i (see Fig. 1(a) ) is given by (cid:126)F ovj = 43 E p √ R eff (1 − ν ) ξ / ij ˆ r ij , (1)where the overlapping distance ξ ij = R i + R i − | (cid:126)r ij | describes the deformation of the spheres, 1 /R eff = 1 /R i + 1 /R j , and ˆ r ij = (cid:126)r ij / | (cid:126)r ij | . The additional termsof the contact force include damping and friction forcesand torques in the same way as described in Refs.[1, 18, 24, 40]. FIG. 1: (a) Representation of the overlap interaction betweentwo elements. (b) Typical deformation of a beam in the x-y plane, showing the resulting bending and shear forces andtorques. The z-axis is perpendicular to the image.
The 3D representation of beams used in this work is anextension of the two-dimensional case of Euler-Bernoullibeams described in Ref. [41]. In 3D, however, the totaldeformation of a beam is calculated by the superpositionof elongation, torsion, as well as bending and shearingin two different planes. The restoring force acting onelement j connected by a beam to element i due to theelongation of the beam is given by (cid:126)F eloj = − E b A b ε ˆ r ij , (2)where E b is the beam stiffness, ε = ( | (cid:126)r ij | − l ) /l , withthe initial length of the beam l and its cross section A b .The flexural forces and moments transmitted by abeam are calculated from the change in the orientationson each beam end, relative to the body-fixed coordinatesystem of the beam (ˆ e bx , ˆ e by , ˆ e bz ). Figure 1(b) shows a ty-pical deformation due to rotation of both beam ends rel-ative to the ˆ e bz axis, with ˆ e bx oriented in the direction ofˆ r ij . Given the angular orientations θ zi and θ zj , the cor-responding bending force (cid:126)Q z,bj and moment (cid:126)M z,bj for theelastic deformation of the beam are given by [41]: (cid:126)Q z,bj = 3 E b I ( θ zi + θ zj ) L ˆ e by , (3a) (cid:126)M z,bj = E b I ( θ zi − θ zj ) L ˆ e bz + (cid:16) (cid:126)Q z,bi × | (cid:126)r ij | ˆ e bx (cid:17) , (3b)where I is the beam moment of inertia. Correspondingequations are written for general rotations around ˆ e by , andthe forces and moments are added up. Additional torsionmoments are added to consider a relative rotation of theelements around ˆ e bx : (cid:126)M x,bj = − G b I tor ( θ xj − θ xi ) L ˆ e bx , (4)with G b and I tor representing the shear modulus and mo-ment of inertia of the beams along the beam axis, respec-tively. The bending forces and moments are transformedto the global coordinate system before they are added tothe contact, volume and walls forces.Beams can break in order to explicitly model damage,fracture, and failure of the solid. The imposed breakingrule takes into account breaking due to stretching andbending of a beam [21, 23, 24, 40, 42], which breaks if (cid:18) εε th (cid:19) + max ( | θ i | , | θ j | ) θ th ≥ , (5)where ε = ∆ l/l is the longitudinal strain, and θ i and θ j are the general rotation angles at the beam ends be-tween elements i and j , respectively. Here cos θ i = ˆ e ibx · ˆ e bx ,where (ˆ e ibx , ˆ e iby , ˆ e ibz ) define the i − particle’s orientation inthe beam body-fixed coordinate system, similar calcula-tion is performed to evaluate θ j . Equation (5) has theform of the von Mises yield criterion for metal plasticity[40, 43]. The first part of Eq. (5) refers to the breakingof the beam through stretching and the second throughbending, with ε th and θ th being the respective thresh-old values. The introduced threshold values are takenrandomly for each beam, according to the Weibull distri-butions: P ( ε th ) = kε o (cid:18) ε th ε o (cid:19) k − exp (cid:34) − (cid:18) ε th ε o (cid:19) k (cid:35) , (6a) P ( θ th ) = kθ o (cid:18) θ th θ o (cid:19) k − exp (cid:34) − (cid:18) θ th θ o (cid:19) k (cid:35) . (6b)Here k , ε o and θ o are parameters of the model, control-ling the width of the distributions and the average valuesfor ε th and θ th respectively. Low disorder is obtained byusing large k values, large disorder by small k . Disor-der is also introduced in the model by the different beamlengths in the discretization as described below.The time evolution of the system is followed by nu-merically solving the equations of motion for the transla-tion and rotation of all elements using a 6 th -order Gearpredictor-corrector algorithm, and the dynamics of the rotations of the elements is described using quaternions[41, 44]. The breaking rules are evaluated at each timestep. The beam breaking is irreversible, which meansthat broken beams are excluded from the force calcula-tions for all consecutive time steps. System formation and characterization
Special attention needs to be given to the discretiza-tion in order to prevent artifacts arising from the sys-tem topology, like anisotropic properties, leading to nonuniform propagation of elastic waves or preferred crackpaths. In our procedure we first start with 27000 spher-ical elements that we initially place on a cubic latticewith random velocities. The element diameters are oftwo different sizes, with D = 0 . D , that are randomlyassigned, leading to more or less equal fractions. Oncethe elements are placed, the system is left to evolve for50000 time steps, using periodic boundary conditions, ina volume that is about 8 times larger then the total vol-ume of the elements. This way we obtain truly randomand uniformly distributed positions.To compact the elements, a centripetal constant ac-celeration field, directed towards the center of the sim-ulation box, is imposed. Due to this field the elementsform a nearly spherical agglomerate at the center of thebox. The system is allowed to evolve until all particlevelocities are reduced to nearly zero due to dissipativeforces.With the elements compacted, the next stage is to con-nect them by beam-truss elements. This is achieved inour model through a Delaunay triangulation of their po-sitions. As a consequence, not only spherical elementsthat are initially in contact or nearly in contact witheach other are connected, but the resulting beam lat-tice is equivalent to a discretization of the material us-ing a dual Voronoi tessellation of the material domain[43, 45, 46]. After the bonds have been positioned, theirYoung’s moduli are slowly increased while the centripetalfield is reduced to zero. During this process the materialexpands to an equilibrium state, reducing the contactforces. The bond lengths and orientations are then resetso that no initial residual stresses are present in the beamlattice. The final solid fraction obtained is approximately0.65. We have compared impact simulations of specimenscompacted as described above with specimens using ran-dom packings of spheres as reported in Ref. [47], whichhave no preferential direction in the packing process suchas the one that could be imposed by the centripetal field.No significant difference was found in the simulation re-sults, indicating that possible radially aligned locked-inforce chains are not relevant.Once the system is formed, the specimen is shaped tothe desired geometry by removing particles and beamsthat are situated outside the chosen volume. The mi-croscopic properties, namely the elastic properties ofthe elements and bonds, as well as the bond breakingthresholds, are chosen to attain the desired macroscopicYoung’s modulus, Poisson’s ratio, as well as the tensileand compressive strength. Table 1 summarizes the in-put values used in the simulations presented in this pa-per. These were chosen to obtain macroscopic proper-ties close to the mechanical properties of polymers likePMMA, PA, and nylon at low temperatures. Figure 2displays the stress-strain curve measured by quasi-static,uni-axial tensile loading of a bar, as depicted in the inset.The microscopic and resulting macroscopic properties areresumed in Table 1, for a sample size (16 × × .
004 s − . Thestress-strain curve is basically linear until the strength isreached where rapid brittle fracture of the material takesplace. Oscillations in the broken specimen fractions canbe seen after the system is completely unloaded due toelastic waves. The Young’s modulus measured from theslope of the curve is 7 . ± . FIG. 2: Stress-strain curve for specimen under quasi-staticloading. The inset shows the load geometry. Abrupt brittlefracture behavior can be observed at about ε = 0 . In order to simulate the impact of a sphere on a fric-tionless hard plate, a spherical specimen with diame-ter D = 16 mm is constructed, and a fixed plane withYoung’s modulus 70 GPa is added to the simulation. Thespherical specimen has a total of approximately 22000 el-ements, with around 32 across the sample diameter. Thecontact interaction between the elements and the plateis identical to the element-element contact interaction,only with ξ = R i − r ip , where r ip is the distance be-tween the particle center and the plate. The specimenis placed close to the plate with an impact velocity v i , in the negative z-direction, assigned to all its composingelements. The computation continues until no additionalbonds are broken for at least 50 µs .For comparative reasons we calculate the evolution ofthe stress field using an explicit Finite Element (FE)analysis. The FE model is composed of axisymmetric,linear 4-node elements with macroscopic properties takenfrom the results of the DEM simulations (see Table 1).Along the central axis through the sphere and groundplate, symmetry boundary conditions are imposed, thebottom of the target plate is encastred and contact sur-faces for the sphere and plate are defined. Figure 3(a)shows a comparison between the impact simulation usingour DEM model and a Finite Element Model simulation.In Fig. 3(a), the DEM elements are colored according tothe amplitude of their accelerations to show the propaga-tion of a longitudinal shock wave that was initiated at thecontact point. The wave speed can be estimated to be ap-proximately 2200 ± FIG. 3: (a) Comparison of deformations and shock-wavepropagation obtained between DEM and FEM simulationsfor v i = 117 m/s. (b) Time evolution of the elastic potentialenergy stored in the system for the same velocity obtained byDEM (solid line) and FEM (dashed line) simulations. lTABLE I: Micro- and macroscopic material model properties.Typical model properties (DEM):Beams:stiffness E b /G b L d ε θ ◦ shape parameter k E p D ρ Hard plate:stiffness E w
70 GPaInteraction:friction coefficient µ γ n − friction coefficient γ t − System:time increment ∆ t N p N b D
16 mmMacroscopic properties (DEM):system stiffness E . ± . ν ρ system strength σ c
110 MPaComparison: DEM FEMlongitudinal 2210 ±
100 2270 ±
20 m/swave speedcontact time 31.4 31.4 µ s III. FRAGMENTATION MECHANISMS
In this section we explore the different fragmentationmechanisms in the order of occurrence and increasingenergy input. The first yield that arises in the materialis diffuse damage that occurs in the region above thecontact disc. It can be seen from Fig. 4(a), that thisdamage region is centered in the load axis, at a distanceapproximately D/ FIG. 4: Initial damage due to bi-axial stress state. (a) Ver-tical cut through the center of the sphere from the DEM sim-ulation showing broken bonds represented by dark color. (b)Stress fields calculated with FEM model. Left side are shearstresses in global coordinates from 0 to -400MPa (black towhite) while the right side shows circumferential stresses inlocal spherical coordinates with the center in the sphere centerranging from 0 to 130MPa (black to white).
As time evolves, meridional cracks start to appear.The origination of this type of cracks is explored in Fig. 5,where we plot in Fig. 5(a) the positions of the brokenbonds in two different projections, showing well definedmeridional crack planes that propagate towards the lat-eral and upper free surfaces of the specimen. In Fig. 5(b)we plot the angular distribution of the broken bonds fordifferent times. Here g ( θ ) is the probability of findingtwo broken bonds with an angular separation θ . Notethat their positions are projected into the plane perpen-dicular to the load axis. The evident peaks in g ( θ ) area clear indication that the cracks are meridional planesthat include the load axis. In this particular case, thecracks are separated by an average angle of about 60 de-grees, and they become evident 13 to 15 µ s after impact( v i = 120 m/s).In order to understand what governs the orientationand angular separation of these meridional cracks we per-formed many different realizations with different seeds ofthe random number generator and impact points. Forall cases the orientation of the cracks can change, butnot their average angular separation. We observe thatfor strong disorder (Eq. (6)), a larger amount of uncorre-lated damage occurs, but the average angular separationof the primary cracks does not change. This suggeststhat the formation of these cracks arises due to a com-bination of the existence of local disorder and the stressfield in the material, but does not depend on the degreeof disorder.As we can see from the FE calculations and from thedamage orientation correlation plot (Fig. 5) inside of theuniform biaxial tensile zone, no preferred crack orienta-tion is evidenced. Many microcracks weaken this ma-terial zone, decreasing the effective stiffness of the core.Around the weakened core the material is intact and un-der high circumferential stresses. It is in this ring shaped FIG. 5: (a) Colored dots display the positions of the brokenbonds according to the time of breaking. (b) Angular distri-bution function of broken bonds as a function of the angularseparation when their positions are projected into the planeperpendicular to the load axis. zone, that we observe to be the onset of the meridionalcracks when we trace them back. For increasing impactvelocity we observe a decrease in the angular separationof crack planes and thus more wedge-shaped fragments.Therefore this fragment formation mechanism can not beexplained by a quasi-static stress analysis. The observa-tion is in agreement with experimental findings and canbe explained by the basic ideas of Mott’s fragmentationtheory for expanding rings [48]. Due to the stress releasefront for circumferential stresses, once a meridional crackforms, stress is released in its neighborhood; the fracturedregions spread with a constant velocity and the probabil-ity for fracture in neighboring regions decreases. On theother hand in the unstressed regions, the strain still in-creases, and so does the fracture probability along with it.The average size of the wedge shaped fragments thereforeis determined by the relationship between the velocity ofthe stress release wave and the rate at which cracks nucle-ate. Thus the higher the strain rate, the higher the cracknucleation rate and the more fragments are formed. Wemeasured the strain rate at different positions inside thebi-axially loaded zone, finding a clear correlation betweenimpact velocity and strain rate. Even though we frag-ment a compact sphere and not a ring, when it comes tothe formation of meridional cracks, we observe that theyform in a ring shaped region and that Mott’s theory canqualitatively explain the decrease of angular separation between wedge shaped fragments with increasing impactvelocity.If enough energy is still available, some of the merid-ional plane cracks grow outwards and upwards, breakingthe sample into wedge shaped fragments that resemble orange slices .As the sphere continues moving towards the plate, aring of broken bonds forms at the border of the con-tact disc due to shear failure (compare Fig. 6(a)). Whenthe sphere begins to detach from the plate, the cone hasbeen formed by high shear stresses in the contact zone(see Fig. 4(b) left) by a ring crack that was able to growfrom the surface to the inside of the material under ap-proximately 45 ◦ (Fig. 6(b)). It detaches, leaving a smallnumber of cone shaped fragments that have a smallerrebound velocity than the rest of the fragments due todissipated elastic energy (Fig. 6(c)). FIG. 6: Vertical meridional cut of the sphere at differentstages during impact, showing the separation of lower frag-ments. (a) Formation of a ring of broken bonds due to shearfailure. (b) These broken bonds evolve into cracks that prop-agate inside the sample. (c) Finally these cracks lead to thedetachment of the lower fragments.
Oblique plane cracks may still break the large frag-ments further, if the initial energy given to the system ishigh enough. Therefore they are called secondary cracks.Figure 7(a) shows a vertical meridional cut of a samplewhere these cracks can be seen. The intact bonds arecolored according the final fragment they belong to.
FIG. 7: (a) DEM simulation at v i = 140 m/s exemplifyingthe secondary cracks. The bonds are colored according tothe final fragment they belong to. (b) 2D simulations usingpolygons as elementary particles [21]. These secondary cracks are very similar to the obliquecracks observed in 2D simulations [14, 21]. For compari-son we show in Fig.7(b) the crack patterns obtained froma 2D DEM simulation that uses polygons as elementaryparticles. In the 2D case, we observe a cone of numeroussingle element fragments and meridional cracks cannotbe observed.
IV. FRAGMENTATION REGIMES
The amount of energy necessary to fragment a materialis a parameter that is very important for practical appli-cations in comminution. In fragmentation experimentstwo distinct regimes for damage and fragmentation canbe identified depending on the impact energy: below acritical energy [4, 6, 19] damage takes place, while abovefragments are formed. Figure 8 compares the final crackpatterns after impact with different initial velocities. Theintact bounds are colored according to the final fragmentthey belong to, and gray dots display the positions of bro-ken bonds. The fragments have been reassembled to theirinitial positions to provide a clearer picture of the result-ing crack patterns. For the smaller impact velocities it ispossible to identify meridional cracks that reach the sam-ple surface above the contact point, but fragmentation isnot complete and one large piece remains (Figs. 8(a) and8(b)). We call these meridional cracks primary cracks ,since as one increases the initial energy given to the sys-tem, some of them are the first to reach the top freesurface of the sphere, fragmenting the material into afew large pieces, typically two or three fragments withwedge shapes (Fig. 8(c)). When we increase the initialenergy secondary oblique plane cracks break the orangeslice shaped fragments further (Fig. 8(d)). Additionalincrease in the impact velocity causes more secondarycracks and consequently the reduction of the fragmentsizes (Figs. 8(e) and 8(f)).The shape and number of large fragments resultingfrom the numerical model for smaller impact energies, aswell as the location and orientation of oblique secondarycracks for larger energies, are in agreement with experi-mental findings [11, 12, 20].We can identify that for velocities smaller then athreshold value, the sample is damaged by the impactbut not fragmented. This threshold velocity for frag-mentation has been found experimentally and numeri-cally [6, 18, 21]. In particular, it has been found from2D simulations that a continuous phase transition fromdamaged to fragmented outcome of impact fragmenta-tion can be tuned by varying the initial energy impartedto the system [18, 21].Following the analysis in references [18, 21] the finalstate of the system after impact is analyzed by observingthe evolution of the mass of the two largest fragments,as well as the average fragment size (shown in Fig. 9).The average mass M /M , with M k = (cid:80) N f i M ki − M kmax excludes the largest fragment. It can be observed that FIG. 8: Front view of the reconstructed spheres, showingthe final crack patterns at the surfaces for different initialvelocities. Gray dots are placed at the positions of brokenbeams while different colors are chosen for different fragments. below the threshold value v th = 115 m/s the largest frag-ment has almost the total mass of the system, while thesecond largest is orders of magnitude smaller. This be-haviour implies that for v i < v th the system does notfragment, it only gets damaged. For velocities larger then v th the mass of the largest fragment decreases rapidly.The second largest and average fragment masses increase,having their maximum at 117 . FIG. 9: The mass of first and second largest fragment and theaverage fragment mass as a function of the impact velocity.
V. RESULTING FRAGMENT MASSDISTRIBUTION
One of the first and still most important character-izations for fragmentation processes are fragment massdistributions. Experimental and numerical studies onfragmentation show that the mass distribution follow apower law in the range of small fragments, whose expo-nent depends on the fragmentation mechanisms, whilethe mass distribution for large fragments is usually rep-resented by an exponential cut-off of the power law. Thefragment mass distributions that are obtained from ourthree-dimensional simulations are given in Fig. 10(a) fordifferent impact velocities v i . Here F ( m ) represents theprobability density of finding a fragment with mass m between m and m + ∆ m . Where m is the fragmentmass normalized by the total mass of the sphere M tot .The values are averaged over 36 simulations, changingthe random breaking thresholds and randomly rotatingthe sample to obtain different impact points. For ve-locities below the critical velocity v th of our model, thefragment mass distribution shows a peak at low frag-ment masses, corresponding to some small fragments.The pronounced isolated peaks near the total mass ofthe system correspond to large damaged, but still unbro-ken system (see also Figs. 8(a) and (b)). Fragments atintermediate mass range are not found for small initialvelocities. At and above v th , the fragment mass distri-bution exhibits a power law dependence for intermediatemasses, F ( m ) ∼ m − τ , (dashed line in Fig. 10(a)) with τ = 1 . ± . Q , for the same velocities. Q is cal-culated by summing the mass of all the fragments smallerthen a given size s . The size of a fragment is estimated asthe diameter of a sphere with identical mass, the valuesare normalized by the sample diameter. By this repre-sentation the large fragments are better resolved. Wecan see that the shape of the size distribution for largefragments can be described by a two-parameter Weibulldistribution, Q ( s ) = 1 − exp[ − ( s/s c ) k s ] (dashed line inFig. 10(b), with s c = 0 .
75 and k s = 5 . k is in the breaking thresholds distributionsin Eq. (7)). Near the critical velocity v th we can iden-tify two main parts in the fragment mass distribution.For m up to around 1 /
40 (approximately 550 elements),the power law F ( m ) ∼ m − τ f ( m/ ¯ m o ) with the cutofffunction f ( m/ ¯ m o ) containing an exponential componentexp ( − m/ ¯ m o ) can be used like in Ref. [49]. However, forlarger m , F ( m ) can also be described by a two-parameterWeibull distribution F ( m ) ∼ (cid:18) k l ¯ m l (cid:19) (cid:18) m ¯ m l (cid:19) k l − exp (cid:34) − (cid:18) m ¯ m l (cid:19) k l (cid:35) . (7)The dashed line in Fig. 11 corresponds to a fit to thedata using ¯ m o = 0 . ± . m l = 0 . ± .
02 and k l =1 . ± .
1. The good quality of the fit allows for a betterestimation of the exponent of the power-law distributionin the small fragment mass range τ = 2 . ± . FIG. 10: (a) Fragment mass distribution for different initialvelocities. The straight line corresponds to a power-law withexponent -1.9 (b) Fragment size distribution weighted by massfor initial velocities with identical legend as above.
VI. CONCLUSIONS
We studied a brittle, disordered fragmenting solidsphere. We performed 3D DEM simulations with 3Dbeam-truss elements for the particle cohesion. Due tothis computationally more laborious approach as com-pared to previous works, we were able to obtain a clearerpicture of the fragmentation process, the evolution offragmentation mechanisms, and its consequences for thefragment mass distribution. To get a clearer insight intothe fracture initiation, we used continuum solutions forthe stress field, obtained by the Finite Element Method.We were able to show, that 2D simulations for frag-menting systems are not capable of capturing fragmenta-tion by meridional cracks, that are the primary crackingmechanism. We found that cracks form inside the samplein the region above a compressive cone long time before
FIG. 11: Fragment mass distribution at v = 122 . they are experimentally observable from the outside, if atall. They grow to form meridional fracture planes thatresult in a small number of large wedge shaped fragments,typically two to four. The increasing tensile radial andcircumferential stress in the ring-shaped region above thecontact plane gives rise to meridional cracks. The de-crease in the angular separation between these crackscould be explained by the Mott fragmentation model.Some of these cracks grow to form the meridional frac-ture planes that break the material in a small number oflarge fragments, and it is only then, that they becomevisible in experiments.The resulting mass distribution of the fragmentspresents a power law regime for small fragments and abroad peak for large fragments that can be fitted witha two-parameter Weibull distribution, in agreement withexperimental results [10, 13, 29, 30]. The fragment massdistribution is quite robust, independent on the macro-scopic material properties such as material strength anddisorder distribution. Only the large fragment range ofthe mass distribution happens to be energy dependent,due to additional fragmentation processes that arise asone increases the impact energy.Even though our results are valid for various materi-als with disorder, they are limited to the class of brittle,heterogeneous media. Extensions to ductile materials arein progress. Another class of interesting questions dealwith the problem of size effects, the influence of polydis-perse particles or the stiffness contrast of particles and0beam-elements. The ability of the model to reproducewell defined crack planes also opens up the possibilityto study other crack propagation problems in 3D. Fortechnological applications questions about the influenceof target shapes and the optimization potential to obtaindesired fragment size distributions or to reduce impactenergies are of broad interest. VII. ACKNOWLEDGMENTS
We thank the German Federation of Industrial Re-search Associations ”Otto von Guericke” e.V. (AiF) for financial support, under grant 14516N from the GermanFederal Ministry of Economics and Technology (BMWI).H. J. Herrmann thanks the Max Planck Prize. F. Kunwas supported by OTKA T049209. We thank Dr. JanBl¨omer and Prof. Jos´e Andrade Soares Jr. for helpfulldiscussions. [1] H. J. Herrmann and S. Roux, eds.,
Statistical Models forthe Fracture of Disordered Media (North-Holland, Ams-terdam, 1990).[2] J. A. ˚Astr¨om, Adv. Phys. , 247 (2006).[3] J. J. Gilvarry and B. H. Bergstrom, J. Appl. Phys. ,400 (1961).[4] J. J. Gilvarry and B. H. Bergstrom, J. Appl. Phys. ,3211 (1962).[5] N. Arbiter, C. C. Harris, and G. A. Stamboltzis, T. Soc.Min. Eng. , 118 (1969).[6] E. W. Andrews and K. S. Kim, Mech. Mater. , 161(1998).[7] J. Tomas, M. Schreier, T. Groger, and S. Ehlers, PowderTechnol. , 39 (1999).[8] R. Majzoub and M. M. Chaudhri, Philos. Mag. Lett. ,387 (2000).[9] K. T. Chau, X. X. Wei, R. H. C. Wong, and T. X. Yu,Mech. Mater. , 543 (2000).[10] A. D. Salman, C. A. Biggs, J. Fu, I. Angyal, M. Szabo,and M. J. Hounslow, Powder Technol. , 36 (2002).[11] S. Z. Wu, K. T. Chau, and T. X. Yu, Powder Technol. , 41 (2004).[12] W. Schubert, M. Khanal, and J. Tomas, Int. J. Miner.Process. , 41 (2005).[13] S. Antonyuk, M. Khanal, J. Tomas, S. Heinrich, andL. Morl, Chem. Eng. Process. , 838 (2006).[14] A. V. Potapov, M. A. Hopkins, and C. S. Campbell, Int.J. Mod. Phys. C. , 371 (1995).[15] A. V. Potapov and C. S. Campbell, Int. J. Mod. Phys.C. , 717 (1996).[16] A. V. Potapov and C. S. Campbell, Powder Technol. ,13 (1997).[17] C. Thornton, K. K. Yin, and M. J. Adams, J. Phys. D.Appl. Phys. , 424 (1996).[18] F. Kun and H. J. Herrmann, Phys. Rev. E. , 2623(1999).[19] C. Thornton, M. T. Ciomocos, and M. J. Adams, PowderTechnol. , 74 (1999).[20] M. Khanal, W. Schubert, and J. Tomas, Granul. Matter. , 177 (2004).[21] B. Behera, F. Kun, S. McNamara, and H. J. Herrmann,J. Phys-condens. Mat. , S2439 (2005).[22] H. J. Herrmann, F. K. Wittel, and F. Kun, Physica A , 59 (2006).[23] F. Kun and H. J. Herrmann, Int. J. Mod. Phys. C. , 837 (1996).[24] F. Kun and H. J. Herrmann, Comput. Method. Appl. M. , 3 (1996).[25] A. Diehl, H. A. Carmona, L. E. Araripe, J. S. Andrade,and G. A. Farias, Phys. Rev. E. , 4742 (2000).[26] J. A. ˚Astr¨om, B. L. Holian, and J. Timonen, Phys. Rev.Lett. , 3061 (2000).[27] L. Oddershede, P. Dimon, and J. Bohr, Phys. Rev. Lett. , 3107 (1993).[28] A. Meibom and I. Balslev, Phys. Rev. Lett. , 2492(1996).[29] C. S. Lu, R. Danzer, and F. D. Fischer, Phys. Rev. E. , 067102 (2002).[30] Y. S. Cheong, G. K. Reynolds, A. D. Salman, and M. J.Hounslow, Int. J. Miner. Process. , S227 (2004).[31] E. W. Andrews and K. S. Kim, Mech. Mater. , 689(1999).[32] P. A. Cundall and O. D. L. Strack, Geotechnique. , 47(1979).[33] B. K. Mishra and C. Thornton, Int. J. Miner. Process. , 225 (2001).[34] C. Thornton and L. F. Liu, Powder Technol. , 110(2004).[35] C. Thornton, M. T. Ciomocos, and M. J. Adams, PowderTechnol. , 258 (2004).[36] D. O. Potyondy and P. A. Cundall, Int. J. Rock. Mech.Min. , 1329 (2004).[37] G. A. D’Addetta and E. Ramm, Granul. Matter. , 159(2006).[38] H. A. Carmona, F. Kun, J. S. Andrade Jr, and H. J.Herrmann, Phys. Rev. E. , 046115 (2007).[39] L. D. Landau and E. M. Lifshitz, Theory of Elasticity ,vol. 7 of
Course of Theoretical Physics (Butterworth-Heinemann, London, 1986), 3rd ed.[40] H. J. Herrmann, A. Hansen, and S. Roux, Phys. Rev. B. , 637 (1989).[41] T. P¨oschel and T. Schwager, Computational Granu-lar Dynamics: Models and Algorithms (Springer-VerlagBerlin Heidelberg New York, 2005).[42] G. A. D’Addetta, F. Kun, E. Ramm, and H. J. Her-rmann, in
Continuous and Discontinuous Modelling ofCohesive-Frictional Materials , edited by P. Vermeer(Springer-Verlag, Berlin, 2001), vol. 568 of
Springer Lec-ture Notes in Physics , pp. 231–258.[43] G. Lilliu and J. G. M. Van Mier, Eng. Fract. Mech. ,
927 (2003).[44] D. C. Rapaport,
The Art of Molecular Dynamics Simu-lation (Cambridge University Press, Cambridge, 2004),2nd ed.[45] J. E. Bolander and N. Sukumar, Phys. Rev. B. , 094106(2005).[46] M. Yip, Z. Li, B. S. Liao, and J. E. Bolander, Int. J.Fracture. , 113 (2006).[47] R. M. Baram and H. J. Herrmann, Phys. Rev. Lett. ,224303 (2005). [48] N. F. Mott, Proceedings of the Royal Society of LondonA , 300 (1946).[49] F. K. Wittel, F. Kun, H. J. Herrmann, and B. H. Kroplin,Phys. Rev. E. , 016108 (2005).[50] D. L. Turcotte, J. Geophys. Res-solid. , 1921 (1986).[51] R. P. Linna, J. A. ˚Astr¨om, and J. Timonen, Phys. Rev.E. , 015601 (2005).[52] F. Wittel, F. Kun, H. J. Herrmann, and B. H. Kr¨oplin,Phys. Rev. Lett.93