A discrete element model for damage and fracture of geomaterials under fatigue loading

Failure processes in geomaterials (concrete, asphalt concrete, masonry, etc.) under fatigue loading (repeated moving loads, cycles of temperature, etc.) are responsible for most of the dysfunctions in pavements, brick structures, etc. In the beginning of the lifetime of a structure, the material presents only inner defects (micro cracks, voids, etc.). Due to the effect of the cyclic loading, these small defects tend to grow in size and quantity which damage the material, reducing its stiffness. With a relatively high number of cycles, these growing micro cracks become large cracks, which characterizes the fracture behavior. From a theoretical point of view, both mechanisms are treated differently. Fracture is usually described locally, with the propagation of cracks defined by the energy release rate at the crack tip; damage is usually associated to non-local approaches. In the present work, damage and fracture mechanics are combined in a local discrete element approach.


Introduction
Fatigue is the weakening of a material caused by repeatedly applied loads of low intensity.The failure processes of geomaterials subjected to fatigue loading can be described by four stages, including crack nucleation (Stage I), short crack growth (Stage II), large crack growth (Stage III), and ultimate failure (Stage IV) [1].In the beginning of the lifetime, the material presents only inner defects (micro cracks, voids, etc.).Due to the effect of the cyclic loading, these small defects tend to grow in size and quantity which damage the material, reducing its stiffness.When these defects become short cracks, the failure process turns into its second stage of short crack growth.With a relatively high number of cycles, these growing short cracks become large cracks, which characterizes the fracture behavior.
It is still difficult today to capture all the four stages with one single fatigue model.In the present work, a discrete element approach is proposed, based on a local description of damage and fracture which allows a smooth interaction between both formulations.The basic concepts of damage and crack growth models are presented in Sec. 2. The discrete element modeling in elasticity is discussed in Sec. 3, while the fatigue extension is analyzed in Sec. 4. Finally, the numerical results are compared to theoretical predictions in Sec. 5, followed by the conclusions.e-mail: xiaofeng.gao@insa-strasbourg.fr e-mail: georg.koval@insa-strasbourg.fr e-mail: cyrille.chazallon@insa-strasbourg.fr 2 Continuum damage and fatigue crack growth models

Continuum damage model
The elastic isotropic continuum damage model characterizes the stiffness decrease with cyclic loading.As the material is deformed, the initiation, growth and coalescence of micro defects decrease the stiffness (degradation of material properties), which is represented by the growth of the damage variable D (0 ≤ D ≤ 1).For undamaged material, D is 0. While D = 1 corresponds to a completely damaged with zero stiffness.The stress strain relation can be written as [1]: where σ i j and ε kl are the stress and strain tensor components, C i jkl and C 0 i jkl are damaged and initial (elastic) se-cant stiffness of the material respectively.C 0 i jkl is a function of Young's modulus E and Poisson ratio ν.
The evolution of the damage proposed by [2] is controlled by the strain state of the material, defined by the scalar equivalent strain ε = where σ i + = σ i if the principal stresses σ i ≥ 0 or σ i + = 0 otherwise.The Young's modulus of the damaged material E is a function of its initial value and the damage variable: The rate of damage growth is defined as a function of local equivalent strain rate (local model): where β is a constant parameter and f (D) is function of the damage.The expression given by Paas [1] is adopted in the following developments: where γ and α are constant parameters.In practice, the values of α, β and γ are calibrated according to the experimental results.

Fatigue crack growth model
The presence of a crack significantly reduce the fatigue life of a component or structure.However, the stress singularity at the crack tip induce unreal estimation of the propagation of cracks in damage models.
Based on fracture mechanics, Paris' law [3] is the most popular fatigue crack growth model used in materials science, which relates the stress intensity factor range ΔK to the crack growth rate under a fatigue stress regime: where N c is the number of loading cycles, da/dN c is the crack growth rate, a is the crack length, B and M are material constants.Eq. 6 is well adapted for relatively long cracks, nevertheless it produces unrealistic small crack growth rates da/dN c , if a → 0.
3 Discrete element elastic model

Numerical scheme and contact model
The simulations are based on the discrete element method (DEM), which uses an explicit numerical scheme to evaluate the contact interactions of the particles and the their motion [4].Each contact force presents normal and tangential components, respectively where δ n and δ t are the normal and tangential relative displacements at the contact, k n and k t are the normal and tangential stiffness, c n and c t are the normal and tangential viscous damping, respectively (see Fig. 2).
The state of the model is updated at each time step Δt after the following operations [5]: (1) position and velocity of each particle is calculated according to Newton's laws of motion; (2) contact detection; (3) contact forces, calculated by Eq. 7.

Continuum media description
An isotropic elastic behavior can be reproduced by the periodic compact arrangement of particles of diameter d shown on Fig. 2b [6,7].The relation between continuum elastic parameters (Young's modulus E and Poisson ratio ν) and discrete elastic parameters (normal and tangential stiffness, k n and k t , respectively) reads in plane stress [6]: The proposed approach can be generalized to nonperiodic particle packings, however a calibration procedure would be necessary in order to identify the material parameters.

Principal stresses
The mean values of the components of the tensors of stress and strain are based on the behavior of one pair of contacts (ki and k j) associated to three particles (i, j and k) [8].A local coordinate system (n; t) is defined; t virtually connects both contacts and n is an orthogonal axis (Fig. 3a).
The stress tensor (in two dimensions) can be defined by the value of the principal stresses σ I and σ II , and their orientation.β is the angle between (n; t) and the coordinate system associated to the principal stresses.Consequently, the principal stresses may be written as: where ) is the tangential stress, which are calculated based on the normal and tangential components of each contact (respectively N ik and T ik , N jk and T jk , as indicated in Fig. 3b).
The normal and the orthonormal strains, respectively 3 /d are obtained by means of the normal and tangential (relative) displacements associated to the contacts ki and k j (respectively δ nik and δ tik , δ n jk and δ t jk , as shown in Fig. 3a).Hence, the angle β is determined as

Discrete element fatigue model
In the present work, a damage approach (see Sec. 2.1) is adopted to describe the behavior before contact rupture.The rupture of a contact, giving rise to a crack propagation, is limited by a crack growth criterion (see Sec. 2.2), which depends on the stress intensity range ΔK.

Damage evaluation
The evaluation of the damage per contact can be summarized by the following operations: (1) DEM elastic analysis (as in Sec.3.1); (2) evaluation of the principal stresses (Eq.9); (3) calculation of the equivalent strain (Eq.2) and strain rate; (4) evaluation of the damage growth (Eqs.4 and 5); ( 5) evaluation of the damage and update of the value of the Young's modulus (Eq. 3) and the normal and tangential stiffness (Eq.8).

Evaluation of the stress intensity range ΔK
In order to obtain consistent results with fatigue crack growth models for crack propagation, when a value of D = 1 is indicated by the fatigue model at a crack tip, the value of the energy release rate is evaluated and a potential crack extension is identified.
The normal components of the contact forces N and T , and contact displacements δ n and δ t , respectively write During a fatigue test, F and δ oscillates between a minimum and a maximum level, which depends on the shape of the cyclic loading.For a sinusoidal loading centered at zero stress, the positive values of F and δ naturally vary from 0 to max F and max δ, respectively.In this case, the damage of a contact induces a maximum energy release rate [9] G max .The minimal energy release rate G min is equal to zero; consequently, the variation of the energy release rate is simply defined as where N cD is the number of cycles to reach D = 1 (total release of the contact energy), g i is the surface of the triangle formed by the points (0, 0), (max δ i−1 , max F i−1 ) and (max δ i , max F i ) as shown in Fig. 4.  Based on the relation between the energy release rate and the stress intensity factor in plane stress [9], the stress intensity range is simply defined as ΔK = E √ ΔG.

Crack initiation and propagation
The evolution of the damage variable D characterizes the weakening of a contact before rupture.The rupture, associated with the propagation of a crack is defined by the value of da/dN c (Eq.6).A contact which presents D = 1, is definitely broken (and cease to exist) only if which represents a crack growth da equivalent to a particle radius d/2.Otherwise, D is set to 0, and the fatigue loading continues, until the number of cycles N cD is sufficiently incremented to fulfill the rupture condition.This rupture criterion prevent deviations from the crack growth criterion induced by the damage model (Sec.2.1) in conditions of singular stress at a crack tip.

Numerical results
A center cracked plate with dimensions w = 80 mm and L = 120 mm (as indicated in Fig. 5) is tested under sinusoidal fatigue loading with amplitude σ 0 = 1.25 MPa and frequency f = 25 Hz.The material presents Young's  In Fig. 6 the numerical evaluation of the stress intensity range ΔK is compared to the theoretical result for different diameters d.Deviations are observed only for very large cracks (2a > 50 mm) and tend to decrease for smaller values of d, which indicates a convergent result.The stiffness degradation can be quantified by the ratio between the initial amplitude of the displacements and its value during the test (which is a function of the number of cycle N c ).In Fig. 7 the fatigue stiffness degradation obtained for plates with different initial crack sizes a 0 are compared to theoretical predictions.For the uncracked plate, the numerical solution follows naturally the damage theoretical prediction.If the plate present a long crack, such as a 0 ≥ 10 mm the numerical results tend to prediction of Paris' law.However, a smaller crack (a 0 = 2 mm) deviates from Paris' law because of the supplementary contribution of the material damage on the stiffness degradation.

Conclusions
A numerical scheme coupling damage and fracture mechanics in a discrete element environment is presented.The association of these different mechanical formulations allows the reproduction of experimental evidences: before material rupture by damage models; during crack propagation by crack growth models.In parallel, important drawbacks of each approach are avoided, such as discretization effects, as shown by the convergent behavior of the results; and the unrealistic results of crack growth models for very short or simply no cracks.

Figure 1 .
Figure 1.The fours stages of the fatigue failure process.

Figure 2 .
Figure 2. (a) Contact forces and associated displacements.(b) Periodic compact arrangement of particles and the contact law.

Figure 4 .
Figure 4. Evaluation of the energy release at the crack tip during a fatigue test.

Figure 5 .
Figure 5. Center cracked plate subjected to fatigue loading.

Figure 6 .
Figure 6.Comparison between the estimation of the stress intensity range ΔK and the theory for different crack lengths a (a 0 = 5 mm).

Figure 7 .
Figure 7. Stiffness degradation as a function of number of cycles N c for different initial crack sizes.The lines indicate the theoretical solutions, while the symbols, the numerical predictions.