Computer simulation of fatigue under diametrical compression H. A. Carmona1,2, F. Kun3, J. S. Andrade Jr.4, and H. J. Herrmann2 1Centro de Ciˆencias e Tecnologia, Universidade Estadual do Cear´a, 60740-903 Fortaleza, Cear´a, Brazil 2IfB, HIF, E18, ETH, H¨onggerberg, 8093 Zu¨rich, Switzerland 3Department of Theoretical Physics, University of Debrecen, P. O. Box:5, H-4010 Debrecen, Hungary and 4Departamento de F´ısica, Universidade Federal do Cear´a, 60451-970 Fortaleza, Cear´a, Brazil 7 0 We study the fatigue fracture of disordered materials by means of computer simulations of a 0 discrete element model. We extend a two-dimensional fracture model to capture the microscopic 2 mechanisms relevant for fatigue, and we simulate the diametric compression of a disc shape speci- n men under a constant external force. The model allows to follow the development of the fracture a process on the macro- and micro-level varying the relative influence of the mechanisms of damage J accumulation over the load history and healing of microcracks. As a specific example we consider 5 recentexperimentalresultsonthefatiguefractureofasphalt. Ournumericalsimulationsshowthat 2 forintermediateappliedloadsthelifetimeofthespecimenpresentsapowerlawbehavior. Underthe effect of healing, more prominent for small loads compared to the tensile strength of the material, ] thelifetime ofthesampleincreases andafatigue limit emerges below which nomacroscopic failure n occurs. The numerical results are in a good qualitative agreement with theexperimental findings. n - s PACSnumbers: 46.50.+a,62.20.Mk,05.10.-a i d . t I. INTRODUCTION obtain a theoretical understanding of the fatigue perfor- a mance of asphalt specimens [16, 17, 18]. Besides dam- m age accumulation, experimental research has shown the Understanding the fatigue damage process of disor- - importance of the healing process to the lifetime of as- d dered materials is obviously important in many areas. n The time scales related to the failure process, the rel- phalt mixtures [9, 10, 11, 13, 19]. The healing mecha- o nism is relatedto the recoveryof microcracksdue to the evant damage mechanisms, and scaling laws associated c viscoelastic nature of the binder material in polymeric with this phenomenon have been widely studied both [ mixtures, resulting in an extended lifetime. Theoreti- experimentally and theoretically [1, 2, 3, 4]. Universal 2 laws have been found in the fracture process related to cal approaches turned out to have difficulties in captur- v ing all the mechanisms involved in the fatigue process, both system sizes and geometry [5, 6], and in particular 0 such as damage growth, relaxationdue to viscoelasticity damage accumulationhas been shown to play an impor- 8 and healing [10, 16, 17, 20, 21]. Although, these mod- tantroleintheformationanddevelopmentofcrackswith 0 els agreewell with the results of the experiments, taking 0 fractal structure [7]. intoaccountthestochasticnatureofthefractureprocess 1 This phenomenon is also important in practical ap- and the cumulative effect of the load history, they fail 6 plications, in particular when dealing with asphalt mix- 0 to describe explicitly the processes of damage, fracture tures subject to traffic loading. Fatigue cracking is one / and failure of the solid. A more complete understand- t of the main causes for asphalt layer failure in pavement a ing ofthe failure processdue to damage accumulationin structures, among moisture damage and thermal crack- m theBraziliantestconfigurationmaybenefitfromdetailed ing. Most of the knowledge about fatigue failure of as- - computer simulations confronted with the experiments, d phalt mixtures relies, however,on experimental observa- as already was successfully performed in many fracture n tions [8, 9, 10, 11, 12, 13]. Materialmodeling to improve processes [22, 23, 24, 25, 26, 27, 28, 29, 30]. o structural design is becoming an important tool to over- c come this mode of material distress. In this paper we enhance a 2D discrete element model : v Recently, fatigue life tests have been carried out to of the fracture of disordered materials in order to cap- i study the performance of asphalt mixtures [14], by mea- ture the relevant microscopic mechanisms of fatigue. In X suring the accumulation of deformation with time for themodelconvexpolygonssymbolizinggrainsofthema- r a cylindrical discs under diametrical compression applied terial are coupled by breakable beams, which can suffer periodically with a constant amplitude. Experiments re- stretching and bending deformation [1, 2, 22, 24, 25, 31, vealedthatthefatigueprocesshasthreedifferentregimes 32, 33, 34, 35, 36]. The beams of the model can repre- depending on the external load amplitude σ: when σ sentthepolymericbinderofasphaltbetweenaggregates. falls close to the tensile strength of the material a rapid Similarly to former studies, breaking of a beam can be failure occurs, while for low load values a so-called fa- caused by stretching and bending captured by the von tigue limit emerges below which the material does not Mises-type breaking criterion [22, 24, 25, 35]. Addition- suffer macroscopic breaking. In the intermediate load ally,weintroduceanageingmechanism,i.e.intactbeams regime the lifetime of the specimen t has a power law undergoadamageaccumulationprocesswhichcanagain f dependence on the load, also called the Basquin law of give rise to breaking. The accumulation process intro- fatigue [15]. Different models have been developed to duces memory over the loading history of the system. 2 TABLE I:Parameter values used in thesimulations. Parameter Value Numberof elements 5070 Density ρ 5 g/cm3 Bulk Youngmodulus Y 1×1010 dyn/cm2 Beam Youngmodulus E 5×1010 dyn/cm2 Time step δt 1×10−6 s Diameter of thedisk D 20 cm Typical size of a single element 0.5 cm Width of theload platen b 2.5 cm Memory factor f0 10−500 Range of memory τ 500−∞ plex multi-axis stress field regime present in the geome- tries used experimentally. This model also resolves the localtorquesandshearforceswhichnaturallyarisewhen dealing with materials with complex micro-structures, where grains or aggregates are embedded in a matrix FIG. 1: (Color online) Geometrical layout of the model. A material. Beams have been particularly used in the lit- disk shaped sample is subjected to diametrical compression. erature to prevent an unphysical failure of the system Thediskiscomposedofconvexpolygonsconnectedbyelastic undershear[35],whichhappenswhenusingcentralforce beams (magnified image). elements. A detailed description of the beam model is givenelsewhere[25,31,37]. Briefly,thecentersofmassof neighboring elements are connected by beams which can Damage recovery leading to healing is implemented in suffer stretching and bending deformation, and hence, the model by limiting the range of memory which con- exert an elastic force and torque on the polygons when tributes to the amount of damage of fibers. We study in stretched or bent [1, 2, 22, 24, 31, 32, 33, 34, 35, 36]. detailsthetimeevolutionofthesystemanddemonstrate In the model beams can break according to specified that the model provides a good quantitative description rules in order to explicitly model damage, fracture and oftheexperimentalfindings,wheredamageaccumulation failure of the solid [25, 37]. The imposed breaking rule and healing proved to be essential. must take into account both the immediate breaking of the beams by stretching and bending, as well as the ac- cumulation of damage and the healing mechanism. In II. MODEL order to achieve these conditions for each beam, at time t, we evaluate the quantities p(t) and q(t) defined by Our model is an extension of a realistic discrete ele- 2 ε max(|θ |,|θ |) ment model of disordered materials which has been suc- p(t)= + i j , (1a) cessfully appliedto study variousaspectsoffracture and (cid:18)εth(cid:19) θth fragmentation phenomena [22, 24, 25, 35]. The material t −(t−t′) ′ ′ is assumed to be composed of a large number of meso- q(t)=p(t)+f0 e τ p(t)dt, (1b) Z scopic elements that interact elastically with each other. 0 In our model the mesoscopic elements are convex poly- where ε = ∆l/l is the longitudinal deformation of the 0 gons randomly generated using a Voronoi construction beam, and θ and θ are the rotation angles at the ends i j [37]. of beams between sites i and j, respectively. Equation The mesoscopic elements are elastic and unbreakable (1a)hastheformofthevonMisesplasticitycriterionde- bodies, all with the same physical properties. Defor- scribingthemechanicalstrengthofbeamswithrespectto mations are represented in the model by the possible stretchingandbendingdeformation. ThefirstpartofEq. overlapof the elements whenpressedagainsteachother. (1a)referstothebreakingofthebeamthroughstretching Between overlapping polygons we introduce a repulsive and the second through bending, with εth and θth being force which is proportional to the element’s bulk Young the threshold values for elongation and bending, respec- modulusY andtheoverlappingarea,andhasadirection tively. Equation (1b) refers to the long term memory, perpendicular to the contact line [25, 37]. and accounts for the accumulation of damage and the Cohesion between the elements is introduced in the healing mechanism. The rate of damage accumulationis model by the inclusion of beams between neighboring assumed to have the form ∆c(t) = f0p(t), i.e. the dam- polygons, see Fig. 1. The use of beams, as opposed to age accumulated until time t is obtained as an integral t ′ ′ simple springs, takes into account tension and bending over the loading history of elements c(t)=f p(t)dt. 0 0 deformations,whichareespeciallyimportantinthecom- In our model healing is captured such that mRicrocracks 3 nucleated in beams can recover, a process which limits (a) (b) the total amount of damage. We describe this effect by introducingafiniterangeofmemoryoverwhichtheload experienced by beams has a contribution to the dam- age. For polymeric materials the dynamics of long chain molecules typically lead to an exponential form of the healing term. In Eq. (1b) the parameter f is a memory 0 factor and τ is the time range of the memory overwhich theloadinghistoryofthespecimencontributestotheac- cumulationofdamage[7],andthereforeaparameterthat controls the healing mechanism. The breaking condition Eqs. (1a) and (1b) are evaluated at each iteration time step and those beams for which q(t) ≥ 1 are considered (c) (d) to be broken, i.e., are removed from the simulation. For simplicity, all the beams have the same thresh- old values ε and θ , and disorder is introduced in the th th model solely through the mesh generation [22, 23, 24, 25, 37, 38]. The global material properties can be tuned by adjusting geometrical and microscopic physical pa- rameters of the model. The randomness of the Voronoi tessellation and the average size of the mesoscopic ele- ments, which define the average length and width of the beams are geometrical parameters that, together with theYoungmodulusofpolygonsY andthatofthe beams E determine the macroscopic response of the model ma- FIG. 2: (Color online) Fracture of a disk-shaped solid sub- terial [22, 23, 24, 25, 37, 39]. The time evolution of the jectedtoaconstantdiametricload. Snapshotsoftheevolving systemiscalculatedbynumericallysolvingtheequations damage. Simulations were carried out with the parameter ofmotionforthe translationandrotationofallelements values: σ/σc = 0.8, f0 = 100, and τ = 1300. For further using a sixth-order predictor-corrector algorithm. The parameters see Table I. breaking rule is evaluated at each iteration step. The breaking of beams is irreversible, which means that it is excludedfromtheforcecalculationsforallfollowingtime III. RESULTS AND DISCUSSION steps. Table I summarizes the parameter values used in the simulations. A. Fracture development Inthenumericalexperimentsthesampleisloadedwith a constant external stress σ. In order to avoid undesir- Since usually the fracturing of materials is highly de- able large elastic waves,the load applied to the opposite pendentonthestructuralenvironmentinwhichitisstud- platens is slowly increased from zero. When the load ied, one needs to reproduce the experimental geometry reaches the value σ, the numerical sample is allowed to asaccuratelyaspossible. StartingfromtheVoronoicon- equilibrate during thousands of time steps including a struction of a square, one cuts a disk-shaped specimen smalldampingbetweenthe grainsandfrictionaccording by removing those polygons outside the disk region and to Coulomb’s friction law [25, 37]. Only after equilibra- reshaping those at the border. The loading conditions tion the breaking rules are applied to the beams. With for the numerical specimen are shown in Fig. 1. This this slow loading and long equilibration time, the vibra- specimen is loaded in diametral compression by non- tions of the sample are drastically reduced compared to deformable platens, which are modeled by introducing the case when a constantloadis appliedinstantaneously two new polygons with the same physical properties of to the platens. The parameter values used in the simu- the loaded material. The load platens are cut with the lations are summarized in Table I. samecurvature ofthe disk. The ratiobetweenthe width First we carried out computer simulations in order to of the load platen b and the diameter of the disk D is determine the quasi-staticstrengthσ ofthe model solid c b/D = 0.125, as in the experiments of Ref. [14]. In the at the given parameters. Then the damage and fracture simulationspresentedinthispapereachdiskspecimenis ofthediscwasstudiedvaryingtheexternalloadσ below composed of 5070 mesoscopic elements. σ . Inordertounderstandthe failureprocesswestudied c 4 (a) (b) FIG.4: (Color online)Finalfracturepatternwithf0 =1(a) FIG.3: (Coloronline)Fractionofbrokenbeamsduetoimme- and f0 = 100 (b) obtained by computer simulations at the diate breaking (squares) and memory accumulation (circles) sameload σ/σc =0.8andτ →∞. Atthehighervalueofthe asafunctionofσ,forf0 =100,τ =500,togetherwiththeto- memory factor f0 more dispersed cracking occurs (b), while talfractionofbrokenbeams(diamonds). Theinsetsshowthe forlowf0 correlatedcrackgrowthcanbeobtainedduetothe final fracture pattern for σ=0.2 (left) and σ=0.8 (right). dominance of immediate breaking. in details the evolution of the crack pattern during the role of the memory factor f is to control the relevant 0 simulations. Figure 2 presents snapshots of the evolving timescaleofthesystem,i.e. thelifetimeofthespecimen system for a disk of diameter 20 cm (approximately 80 t has a dependence t ∼1/f . f f 0 polygons)andσ/σc =0.8. Intheinitialstageoffracture, Inordertoanalyzetheimportanceoftheaccumulation shown in Fig. 2(a), a number of cracks appears in the of damage in the fracture process, we have investigated regions situated near the border of the platens contacts. the final fractions of the beams that have been broken Becauseofthe finite widthofthe platens,the stressfield duringthecourseofthesimulationsduetotheimmediate inducedinthesampleisnotpurelyuniaxialtensilestress breakingrule(firstterminEq.(1b))andthosethathave along the load direction as one might expect [40], so the been broken due to damage accumulation (second term cracks do not initiate from the middle of the specimen. in Eq. (1b)) for different applied loads. Figure 3 shows In fact, the cracks initiate from the border and coalesce the fractions of the broken beams due to the different to form wedges, as shown in Figure 2(b). The wedges breakingmodesasafunctionoftheappliedload. Forlow extend up to approximately0.25D,exhibiting a seriesof loadvaluesdamagedominatesthefailureofbeams,while parallelverticalcracks,characteristicofcompressivetests immediate breaking plays only a minor role. As the ex- [25]. As time evolves the wedges penetrate the specimen ternalloadisincreased,thefractionofbrokenbeamsdue (see Fig. 2(c)) originating a fracture in the center of the to damage accumulation decreases monotonically, while specimen, in a narrow region between the load platens. that due to immediate breaking increases and saturates. This mechanism leads to the sudden and catastrophic For σ ≥ 0.65σ the immediate breaking mode becomes c failure of the specimen (as shown in Figure 2 (d)). This dominant. Itisinterestingtonotethatthetotalfraction failure behavior agrees qualitatively well with previous of broken elements is a decreasing function of the load, experimental observations [14, 26]. implying that correlated crack growth dominates when It is important to emphasize that without the time σ →σc dependent damage accumulation, the system would not Oneofthemostimportantfeaturesofourmodelingap- have macroscopic failure under a load σ < σ . After proach is that it accounts for the inhomogeneous stress c some cracking events the disc would attain equilibrium field in the specimen and provides also information on resulting in an infinite lifetime. However, according to the microstructure and spatial distribution of damage. thefailurecriterionEq.(1b),intactelementsaccumulate The ageingofbeams is less sensitiveto the details ofthe damagewhichresultsinbeambreakingafterafinitetime stress field in the specimen. That is why when damage even under a constant local load. These beam breakings accumulation dominates the breaking of beams microc- due to damage accumulation give rise to load redistri- racks are more dispersed in the stripe of the specimen bution and stress enhancements on intact elements trig- between the loading platens. The insets of Fig. 3 show gering additional breakings. In this way, if healing is thefinalfracturepatterncorrespondingtoσ =0.2σ and c neglected, or analogously τ → ∞ is set in Eq. (1b), the σ =0.8σ . For σ =0.2 (left), where the effect of fatigue c system would fail after a finite time at any small load is more important, we observe more disordered cracks, value. One can conclude from Eq. (1b) that the main while for σ =0.8, where the immediate breakingprocess 5 FIG. 5: Deformation ε of the disk as a function of time t for FIG. 6: Fit of the experimental results of Ref. [14]. Using f0 =100andτ =∞varyingtheexternalloadσ. Macroscopic the parameter values f0 = 100 and τ = 1300 a good quality failure is preceeded by an acceleration of the deformation. agreement is obtained with theexperiments. Theverticaldashedlineindicatesthelifetimet ofthesystem. f Note that t decreases with increasing σ. f that, by increasing σ, the lifetime t rapidly decreases f is moreimportant, longercrackscanbe identified due to andactuallygoestozeroastheexternalloadapproaches crack growth. the tensile strength of the material σ . c In our model the importance of damage accumulation Forthepurposeoffittingthe experimentalresults,the forthebeambreakingiscontrolledbythememoryfactor most important parameters of the model that can be f withrespecttotheimmediatefailureofelements. Fig- 0 tuned are the two threshold values ε and θ for the ure 4 demonstrates crack patterns obtained at the same th th stretchingandbendingdeformationofthebeams,aswell load without healing τ → ∞ varying the value of f . 0 as the memory factor f and the range of memory τ. As It can be observed that for large values of f , even at 0 0 we have pointed out, f only sets the time scale of the a relatively high load σ/σ = 0.85, a large amount of 0 c system, while the other three parameters have a strong distributed cracking occurs in the specimen. Correlated effect on the qualitative form of ε(t) and t (σ). crackgrowthcanonlybeobservedforlowf (seeFig.4a) f 0 where immediate breaking has dominance. These crack In the experiments of Ref. [14] the asphalt specimen patternsareverysimilartothoseoriginatedduetomem- wassubjectedto periodicloadingwitha constantampli- ory effects in thermal fuse networks [41]. The effect of tudemonitoringthemaximumdeformationasafunction increasing the memory factor f0 resembles the effect of ofthe numberofloadingcyclesandthenumberofcycles decreasing the applied load, in the sense that increasing to failure. In our model a constant load or a periodic f0 has the same effect of enlarging the lifetime of the load with a constant amplitude would have the same ef- sample, as can be inferred from Eq. (1b). fect,hence,wedirectlycomparethesimulationresultsto the experimental findings. Figure6 showsa fit ofthe experimentalresults ofRef. B. Macroscopic time evolution [14]tothesimulations,obtainedbyrelatingtimetonum- ber of cycles and an appropriate choice of geometrical Onthemacroscopiclevelthetimeevolutionofthesys- andphysicalparametersforthemesoscopicelementsand tem is characterized by the overall deformation ε of the beams. We seethatthe simulationandthe experimental disk and the lifetime t of the specimen. The total de- f curves show excellent agreement. formation of the specimen, ε, is obtained directly from the simulation by monitoring the position of the loading In order to examine the consequences of the predom- platens. Figure 5 shows the behavior of ε as a func- inance of stretching or bending breaking modes we per- tion of time for different values of σ with respect to the formed a number of simulations with different threshold tensile strength σ of the simulated specimen. Due to values. The dependence of the lifetime on the external c accumulationof damage, ε increasesmonotonically until applied load is shown in Figure 7, for τ = ∞ and dif- catastrophic failure of the material occurs. Some fluctu- ferent values of ε and θ . The simulations show that th th ationscanbeseeninεbeforeitstartsincreasingrapidly, for external loads σ → σ , for all curves the lifetime de- c due to the natural randomness of the beam breaking. creases rapidly and approaches the immediate failure of The lifetime t of the specimen is defined as the time the specimen. For intermediate values of applied load, f at which the deformation ε diverges. It can be observed t exhibits a powerlaw behavior,in accordancewith the f 6 a) θ = 3o th experiment θ = 5o th τ=∞ 104 θth = 7o τ=21150 −1.6 θth = 10o 104 τ=16450 θ = 15o τ=11750 th tf θ = 20o tf th 103 −1 2 10 ε = 0.03 th 10−1 10−1 σ/σ 100 σ c b) εth = 0.01 FIG. 8: Lifetime as a function of the load σ/σc. The square ε = 0.015 filled symbols correspond to experimental results, the open th circles correspond to τ =∞, and the ∗, + and × symbols to ε = 0.02 104 −1.4 εth = 0.03 τ =21150, 16450 and 11750, respectively. th ε = 0.04 th tf stretching breaking mode. −1.8 In order to obtain the measured value γ = 2.0±0.1 the breaking parameters had to be set such that the stretching mode of Eq. (1a) dominates the breaking, i.e. ε = 0.01 and θ = 20o were used which implies that th th θ = 15o thebendingmodehasonlyaminorroleinbreaking. The th upturn of the curves of t for small loads in Fig. 8 indi- f 10−1 cates the emergence of the fatigue limit σ in the system σ l belowwhichσ <σ thespecimenonlysufferspartialfail- l ure and has an infinite lifetime. Figure 8 also illustrates thatthevalueofσ isdeterminedbytherangeofmemory l FIG. 7: Lifetime as a function of the load σ for different τ, since healing has only a dominating effect on the time breaking thresholds. In (a) the stretching threshold is kept evolutionofthe systemwhenτ is comparableto the life- fixed εth = 0.03 and we vary the bending angle threshold timeofthespecimenmeasuredwithouthealing(τ →∞). θth, while in (b) θth = 15o and the lifetime is calculated for Hence, it follows from the Basquin law, Eq.2, that the different values of εth . fatigue limit σl scales with τ as σl ∼ τ−γ1. The good quantitative agreement between simulations and experi- mentobtainedforτ =21150indicatesthatthestretching Basquin law of fatigue [15] breakingmechanismismoreimportantforthisparticular material. This is reasonable since the polymer binding −γ σ in asphalt can not sustain torsion. t ∼ , (2) f (cid:18)σ (cid:19) c where the value of the exponent γ can be controlled by IV. CONCLUSIONS the breaking parameters ε ,θ . We can see in Fig. 7a th th thatforafixedvalueofε =0.03the Basquinexponent th variesfromγ =−1.1±0.1forθ =3o(bendingbreaking We carried out a computational study of the fatigue th mode dominance) to γ = −1.6±0.1 as we increase θ failureofdisorderedmaterialsoccurringunderaconstant th to 20o, hence increasing the influence of the stretching externalload. Atwo-dimensionaldiscreteelementmodel breaking mode. Obviously, the value of σ also increases was extended to capture microscopic failure mechanisms c withincreasingθ ,since oneofthe breakingmodes gets relevantfortheprocessoffatigue. Asaspecificexample, th completelysuppressed. SimilarlyinFig. 7b,withafixed we considered experiments on asphalt where damage re- value of θ =15o the Basquinexponent variesfromγ = covery in the form of healing is known to play a crucial th −1.4±0.1 for ε = 0.04 to γ = −1.8±0.1 for ε = role for the long term performance of the material. th th 0.01, in which case there is a clear predominance of the The breaking of beams is caused by two mechanisms: 7 immediate breaking occurs when the failure thresholds the model is controlledby the relative importance of the of stretching and bending are exceeded. Intact beams stretching and bending modes in beam breaking. The undergoadamageaccumulationprocesswhichislimited fatigue limit σ is obtained as a threshold value of the l by the finite range of memory (healing). The relevant external load below which the specimen does not suffer parameters of the model to fit the experimental results macroscopic failure and has an infinite lifetime. In the arethebreakingthresholdsofstretchingε andbending framework of our model σ is solely determined by the th l θ and the range of memory τ. We carried out a large rangememory τ overwhichlocal loadscontribute to the th amount of computer simulations to study the time evo- total accumulated damage in the system. lution of the system at the macro- and micro-level. We demonstratedthatatthemicro-levelfailureofbeamsdue to damage is responsible for the diffuse crack pattern, while correlated crack growth is obtained when immedi- V. ACKNOWLEDGMENTS ate breakinghasdominance. Themodelprovidesagood quantitative agreementwith previous experimental find- ings on the fatigue failure of asphalt, i.e. varying only WethanktheBrazilianagenciesCNPq,FUNCAPand threeparametersε ,θ ,andτ ofthemodelwecouldre- CAPESforfinancialsupport. H.J.Herrmannisgrateful th th producethe deformation-time diagramandthe Basquin- to Alexander von Humboldt Stiftung and for the Max law of asphalt obtained in experiments. Very interest- PlanckPrize. F.KunwassupportedbyOTKAM041537, inglywefoundthatthe valueofthe Basquinexponentof and T049209. [1] M. J. Alava, P. K. V. V. Nukala, and S. Zapperi, Adv. [21] Y.Q.LiandJ.B.Metcalf, J.Mater. CivelEng. 14,303 Phys.55, 349?476 (2006). (2002). [2] H.J.HerrmannandS.Roux,eds., Statistical Models for [22] F. Kun and H. J. Herrmann, Comput. Methods Appl. theFracture ofDisordered Media (North-Holland,1990). Mech. Engrg. 138, 3 (1996). [3] R.Shcherbakovand D. L.Turcotte, Theor. Appl.Fract. [23] F.KunandH.J.Herrmann,Int.J.Mod.Phys.C7,837 Mec. pp. 245–258 (2003). (1996). [4] D.L. Turcotte, Fractals and Chaos in Geology and Geo- [24] F. Kun and H. J. Herrmann, Phys. Rev. E 59, 2623 physics (Cambridge Univ.Press, 1992). (1999). [5] L. de Arcangelis, A. Hansen, H. J. Herrmann, and [25] F. Kun, G. A. D’Addetta, H. J. Herrmann, and S.Roux,Phys. Rev.B 40, 877 (1989). E. Ramm,Comp. Ass. Mech. Eng. Sci. 6, 385 (1999). [6] A.Carpinteri and N.Pugno, NatureMaterials 4 (2005). [26] G. Lilliu and J. van Mier, in H. W. Reinhardt 60 years [7] H.J. Herrmann,J. Kert´esz,and L. deArcangelis, Euro- (1999). phys.Lett.10, 147 (1989). [27] W.C.ZhuandC.A.Tang,Int.J.Rock.Mech.Min.Sci. [8] J. M. Matthews, C. L. Monismith, and J. Craus, J. 43, 236 (2006). Transp. Eng.-ASCE 119, 634 (1993). [28] A. Brara, F. Camborde, J. Klepaczko, and C. Mariotti, [9] J.S.DanielandY.R.Kim,J.Mater.CivilEng.13,434 Mech. Mater. 33, 33 (2001). (2001). [29] C. Thornton, M. T. Ciomocos, and M. J. Adams, Pow. [10] S. Zhiming, D. N. Little, and R. L. Lytton, J. Mater. Technol. 140, 258 (2004). Civil Eng. 14, 461 (2002). [30] B. van de Steen and A. V. ad J. A. L. Napier, Int. J. [11] Y.-R.Kim,D.N.Little,andR.L.Lytton,J.Mater.Civil Fract. 131, 35 (2005). Eng. 15, 75 (2003). [31] H. J. Herrmann, A. Hansen, and S. Roux, Phys. Rev. B [12] R. Lundstrom, H. D. Benedetto, and U. Isacsson, J. 39, 637 (1989). Mater. Civil Eng. 16 (2004). [32] B.Skjetne,T.Helle,andA.Hansen,Phys.Rev.Lett.87 [13] M.CastroandJ.A.S´anchez,J.Transp.Eng.-ASCE132, (2001). 168 (2006). [33] J.A.Astr¨om,R.P.Linna,J.Timonen,P.F.Moller,and [14] F. Kun, M. H. S. Costa, R. N. C. Filho, J. S. A. Jr., L. Oddershede,Phys. Rev.E 70, 026104 (2004). J. B. Soares, S. Zapperi, and H. J. Herrmann, cond- [34] R. P. Linna, J. A. Astr¨om, and J. Timonen, Phys. Rev. mat/0607606 (2006). E 72, 015601(R) (2005). [15] O. H. Basquin, in Am. Soc. Test. Mater. Proc. (1910), [35] G.-A.D’AddettaandE.Ramm,GranularMatter8,159 vol. 10, pp. 625–630. (2006). [16] D.Sornette,T.Magnin,andY.Br´echet,Europhys.Lett. [36] J. A.Astr¨om, Adv.Phys. 55, 247?278 (2006). 20, 433 (1992). [37] G. D’Addetta, F. Kun, E. Ramm, and H. J. Herrmann, [17] H.J.Lee, Y.R.Kim,and S.H.Kim,Struc.Eng. Mech. inContinuous andDiscontinuous ModellingofCohesive- 7, 225 (1999). Frictional Materials (Berlin: Springer,2001), vol.568 of [18] Y.-R.KimandD.N.Little,J.Mater.CivilEng.16,122 Springer Lecture Notes in Physics, pp.231–258. (2004). [38] G.A.D’Addetta,Ph.D.thesis,Institutfu¨rBaustatikder [19] Y. R. Kim, D. N. Little, and R. C. Burghardt, J. Mat. Universit¨at Stuttgart (2004). Civil Engrg. 3, 140 (1991). [39] B. Behera, F. Kun, S. McNamara, and H. J. Herrmann, [20] S.Prager and M. Tirrel, J. Chem. Phys. 75 (1981). J. Phys.: Condens. Matter (2005). 8 [40] O. Tsoungui, D. Vallet, J.-C. Charmet, and S. Roux, [41] D. Sornetteand C. Vanneste,PRL 68, 612 (1992). Gran. Matter 2, 19 (1999).