A Discrete Element Model for Damage and Fatigue Crack Growth of Quasi-Brittle Materials

.e repeatedly applied low-intensity loads would lead to the damage and fatigue crack growth of mechanical structures made of quasi-brittle materials. In numerical modelling, these two mechanisms are normally treated differently and separately; the damage is usually associated with nonlocal approaches, while the fatigue crack growth is related to the local stress intensity range at the crack tip. In this study, a discrete element model for damage and fatigue crack growth of quasi-brittle materials is proposed, which is able to model the damage and fatigue crack growth simultaneously in one single model. .e proposed model achieves the implementation of a continuum damage model in a discrete element code, which is a helpful enrichment of this numerical method. .e evaluation method of the stress intensity range during the damage evolution provides a way to couple both failure mechanisms. .is feature allows crack initiation to be induced by localized damage and a progressive transition to a fracture behaviour with the crack propagation. Independent parameters for the fatigue damage model and fatigue crack growth model are admitted without any previous calibration. .e numerical results are in good agreement with the theoretical predictions of damage and fracture mechanics, and intact and precracked samples are analysed under fatigue loading to show the consistent coexistence of fractured and damaged zones in a single model.


Introduction
e size and boundary effects of quasi-brittle fractures have been studied systematically in the past decades [1][2][3][4][5][6][7][8][9], and some significant progresses have been achieved in recent years by taking the aggregate or grain size into account in the analytical solutions [10][11][12][13][14].Besides the complexities of failures under monotonic loading condition, fatigue is also one of the most common yet complicated failures that can cause damage to mechanical structures [15][16][17].e fatigue loading (repeated moving loads, cycles of temperature, etc.) applied to structures made of quasi-brittle materials (concrete, asphalt concrete, masonry, etc.) generates efforts that are far below the tensile strength or the fracture toughness of the material.However, they are responsible for the degradation of the material properties and fatigue crack growth, which may lead to the final failure state of these structures.
e number of cycles, of a specified character that a specimen sustains before failure of a specified nature occurs, is defined as the fatigue life.Fatigue life is very important for the structure design; however, its prediction is still an empirical science rather than a theoretical one [18].
e fatigue behaviour can be studied experimentally and numerically.e experimental study is commonly expensive and time-consuming and sometimes impossible in the case of huge structures, while the numerical study is time-and cost-efficient and can effectively enable researchers to optimize the experimental effort required [15].It is admitted that the continuum damage models are able to effectively characterize the fatigue life of the specimen without large cracks, including the stages of crack nucleation and short crack propagation.Once a large crack appears, the existence of a stress singularity at the crack tip leads the stress-based or strain-based continuum damage model to a fast and unreal propagation of the damaged zone, a representation of the crack.
e nonlocal continuum damage models [19][20][21] intend to solve this problem; however, these models are not able to indicate the damage or the evolution of the damage during the fatigue loading [22].
erefore, the residual strength cannot be determined using these models.Moës et al. proposed an alternative nonlocal damage model called the ick Level Set (TLS) approach [23][24][25][26] and implemented it in an extended finite element code to model the damage growth in solids.
e undamaged zone and the damaged zone are separated by a level set, and in the damaged zone, the damage variable is an explicit function of the level set.e damage growth is expressed as a level set propagation [23].
is model can be considered as a continuous transition from damage to fracture [26], which is able to handle initiation, growth, branching, and coalescence of crack-like patterns.Latif et al. [27] developed an interface damage model based on the ick Level Set approach, for simulating fatigue-driven delamination in composite laminates using the finite element method.It provides a link between damage mechanics and fracture mechanics through the nonlocal evaluation of the energy release rate, which is able to predict the fatigue crack growth rate and the delamination growth pattern accurately.e ick Level Set approach was compared with cohesive zone models in [28], which are the crack-based models that deal with crack evolution in elastic materials.
e cohesive zone models represent quasi-brittle behaviours with good accuracy but require extra equations to determine the crack path.ese models were widely used to predict the fatigue crack growth [29][30][31] and also improved to study the damage initiation and evolution in the interconnections [15,[32][33][34].
Both continuous and discrete numerical methods have their own advantages and shortcomings in the modelling of fracture and damage phenomena.e finite element method (FEM) is cumbersome because the mesh requires to be updated to match the geometry of the discontinuity [35].e mesoscopic models developed for studying the influence of material composition on the overall behaviour [36][37][38] have extended the capacity and quality of the finite element modelling but not yet overcome its shortcomings.e extended finite element method (XFEM) was developed in 1999 by Belytschko and Black [39] and improved by Moës et al. [40] to help alleviate shortcomings of FEM and has been used to model the propagation of various discontinuities.XFEM has been successfully adopted for simulation of various engineering problems, which have been reviewed in detail in [41].
e discrete element method (DEM) is particularly attractive for modelling quasi-brittle materials due to its ability to construct a mesh that is not completely continuous and homogeneous.Since the mesh is constructed by rigid elements that interact with each other at points of contact, it is much easier to construct the mesoscopic DEM models with voids, imperfections, and heterogeneities than the FEM and XFEM models.e XFEM and DEM were compared in [42] for the simulations of crack initiation, propagation, and coalescence in brittle materials.It shows that DEM yields the better results compared to the XFEM, which is able to predict the crack propagation and coalescence in good agreement with the test results, while XFEM failed to model the shear cracks and has difficulty to generate the crack coalescence [42].
In this study, a discrete element approach is proposed based on a local description of damage and fracture.e continuum damage model is implemented in a discrete element code and verifies the theoretical prediction.In Section 5, the continuum damage model is coupled with a fatigue crack growth model.is feature allows crack initiation to be induced by localized damage and a progressive transition to a fracture behaviour with the crack propagation.Furthermore, independent parameters for the damage and the crack propagation laws are admitted without any previous calibration.Intact and precracked samples are analysed under fatigue loading to show the consistent coexistence of fractured and damaged zones in a single model.Finally, the numerical results are compared to theoretical predictions of fracture mechanics.

Continuum Damage Models (CDMs).
CDMs are capable of predicting the fatigue life of the specimen without large cracks.During the fatigue tests, as the material is deformed, the initiation, growth, and coalescence of microdefects decrease the stiffness (degradation of material properties), which is represented by the growth of the damage variable D. D is a scalar variable ranging from 0 to 1.For the virgin, undamaged material, D is 0. D � 1 corresponds to a completely damaged material with zero stiffness.
e stressstrain relation can be written as follows: where σ ij and ε kl are the elastic stress and strain tensor components and C ijkl and C 0 ijkl are the damaged and initial (elastic) secant stiffness of the material.C 0 ijkl is a function of Young's modulus E and Poisson's ratio ].
Several stress-and strain-based continuum damage models have been established by the researchers, including the models established by Castro and Sanchez [43], Di Benedetto [44], Lee [45], and Bodin et al. [19,21].All the models are capable of predicting the fatigue life of quasibrittle materials.However, the damage evolution law for D and stress and strain adopted in the damage models might be different.In this study, the local version of Bodin's L2R continuum damage model [19] will be utilized, which is widely used for quasi-brittle materials [19,46,47].
Bodin et al. proposed an elastic isotropic continuum damage model for fatigue, which characterizes the decrease in stiffness with cyclic loading.e damage growth criterion is based on a modified Rankine criterion with zero threshold damage growth.Evolution of damage (local version of the model) is controlled by the strain state of the material by a scalar equivalent strain, which can be written as follows: 2 Advances in Materials Science and Engineering In the nonlocal version of Bodin's model, the "local" equivalent strain  ε is replaced by its weighted average strain.
e expressions for weighted average strain calculation can be found in [20,21].
e rate of damage growth is defined as a function of the local equivalent strain rate: where  ε _ �  ε i+1 −  ε i is the increment of local equivalent strain from t � i to t � i + 1 and f(D) is a function of damage.e exponent β is a material parameter, which relates to the slope (−1 − β) of the S-N curve in the log-log diagram.e function of damage f(D) given by Paas (Equation ( 4)) [48] is adopted by Bodin for his L2R [19] damage model.Bodin's L2R damage model can capture the stiffness degradation process before the fatigue crack becomes the dominant factor.It should be noticed that, in order to model all the stages of stiffness degradation, Bodin proposed an L3R damage model with a new function of damage with three adjustable model parameters [19].With Bodin's L3R damage model, the fatigue life can be well predicted.However, the numerical results of failure patterns are quite different in comparison with the test results [22]: (4)

Fatigue Crack Growth Models (FCGMs).
e existence of cracks can significantly reduce the fatigue life of a component or the whole structure.Some physical evidence is quite intuitive: at higher stresses, a crack tends to propagate "faster" (with respect of the number of cycles), and at similar stress levels, a bigger crack tends also to propagate "faster."e propagation of localized cracks depends on many parameters and is not well characterized by damage models.Fatigue crack growth models (FCGMs) take into account the variation of the stress intensity factor during the cycles to describe the crack propagation [49].By simplicity, only Paris' law is discussed.

Paris' Law.
Fatigue crack growth in a wide variety of brittle and quasi-brittle materials [50][51][52] is described well by the well-known Paris (or Paris-Ergodan) law [49], which relates the stress intensity factor range ΔK to the crack growth rate da/dN c via a power law, with ΔK � K max − K min .
e basic formula is as follows: where da/dN c is the crack growth rate, a is the crack length, and N c is the number of load cycles; c and m are the material parameters dependent on the material composition, environment, stress ratio, etc. [53].ΔK th and K c are the fatigue threshold and fracture toughness of the material, respectively.Paris' law works for sufficiently large cracks, where ΔK is higher than the fatigue threshold ΔK th , and the maximum value of the stress intensity factor K max remains below the material fracture toughness K c .

Discrete Element Method
e discrete element method (DEM) was originally developed by Cundall for modelling granular and particulate systems [54].e method is based on the use of a numerical scheme in which the interaction of the particles is monitored at every contact and the motion of the particles is modelled for every particle.
e constitutive parameters for the contacts between the discrete elements, such as stiffness and strength, influence the behaviour of the model at the macroscopic scale.ese parameters are usually calibrated in order to reproduce experimental results if the random packing arrangement is adopted, which is able to simulate real material behaviour including Poisson's effect and random crack paths.
Challenges related to calibration between the macro-and micromaterial parameters can be avoided if a close-packed assembly (regular hexagonal packing) is adopted, as shown in Figure 1(b).e close-packed assembly is composed of particles with identical size.e mechanical properties of the specimen (Figure 1(a)) depend on the mechanical parameters of the particles.e close-packed assembly has almost the same benefits as the random packing arrangement, and it is much easier to assign and update the contact stiffness during the damage calculation.e heterogeneous model can be constructed by representing different materials (aggregates, mortar, etc.) as clusters of discrete elements with different parameters [55]. is close-packed assembly has been used by  and Liu et al. [60] for quasi-brittle failure of asphalt concrete, rocks, and crystals.

Contact Model.
DEM discretizes a material using elements of simple shapes (circles, spheres, or blocks) that interact with neighbouring elements according to laws of interaction that are applied at points of contact.At each time step, the computation of all contact forces is followed by the application of Newton's second law to the particles.Each contact force has normal and tangential components, N and T, respectively (Figure 2).e contact behaviour follows a standard linear spring and dash-pot model.When the values of the damping parameters, c n and c t , are a sufficiently small fraction of ���� mk n  (where m is the particle mass), the inelastic effect is negligible [59].In this study, only the elastic contribution of the contact force and the relative displacements, δ n and δ t (Figure 2), will be considered.
Young's modulus E and Poisson's ratio ] are the two elastic constants used to characterize the macroscopic linear elastic behaviour of isotropic materials.A direct relationship between these macroscopic parameters and discrete elastic parameters (normal and tangential stiffness, k n and k t , Advances in Materials Science and Engineering respectively) has been established by the researchers [56,61], which for plane stress is as follows: where t is the thickness of the discrete element model.For plane strain condition, this direct relationship can also be found in [56,61].

Mean Strains and Stresses.
e mean values of the components of the tensors of stress and strain are based on the behaviour of one pair of contacts (ki and kj) associated with three particles (i, j, and k).A local coordinate system (n, t) is de ned, where t virtually connects both contacts for which n is an orthogonal axis (Figure 3(a)).e normal and tangential (relative) displacements associated with contacts ki and kj (δ nik , δ tik and δ njk , δ tjk , respectively, as shown in Figure 3(a)) give rise to mean strain values, ε nn and ε tt : e normal and tangential components of each contact (N ik , T ik and N jk , T jk , respectively, as indicated in Figure 3(b)) can be projected over the directions (n, t) giving rise to the resultant forces.Considering a particle with diameter d and thickness t (t d), mean stresses (Figure 3(c)) may be associated with these forces:

Principal Stresses.
e stress tensor (in two dimensions) can be de ned by the values of the principal stresses σ I and σ II and their orientations.Hence, ψ is de ned as the angle between (n, t) and the coordinate system associated with the principal stresses (Figure 3(c)).Consequently, the principal stresses may be written as follows: e value of ψ is determined based on the information from mean strains and mean stresses at a pair of contacts, which can be written as follows: Unit cell 4 Advances in Materials Science and Engineering where

Discrete Element Fatigue Damage Model
Most of the fatigue laboratory tests of a high number of cycles consist of applying a sinusoidal displacement (or force) with constant amplitude at the boundary of the sample.During testing, the variation of global sti ness is monitored, which is de ned as the ratio between the amplitudes of the force and the displacement.In order to numerically model these time-consuming laboratory tests, Bodin et al. [19,21] proposed a nonlocal damage model to predict the material behaviour (Section 2).e model was implemented in a nite element code, along with a selfadaptive jump-in-cycle procedure for high cycle fatigue computations.
In this section, a local version of Bodin's L2R damage model is implemented in a discrete element code.A closepacked assembly (regular hexagonal packing) is adopted due to the direct relationship between the macroscopic parameters (Young's modulus E and Poisson's ratio ]) and discrete elastic parameters (normal and tangential sti ness, k n and k t , respectively).During damage modelling, Young's modulus of the material decreases as the evolution of the damage value D, which can be characterized by the decrease of local contact sti ness.

Model Implementation.
e damage of all the contacts in the assembly is calculated in the same way and updated at the same time.e evaluation of the damage per contact can be summarized by the following operations: (1) DEM elastic analysis and identi cation of the stress and strain elds (Section 3.2): the contact forces N and T and contact displacements δ n and δ t are the direct values that can be obtained from a discrete element analysis; based on this information, the stress and strain of each contact pair in the close-packed assembly (Figure 3) can be calculated by Equations ( 7) and ( 8). ( 2) Evaluation of the principal stresses: the stress tensor σ n and σ t (in two dimensions) can be de ned by the value of the principal stresses σ I and σ II and their orientation by Equations ( 9) and ( 10).(3) Calculation of the equivalent strain (Equation ( 2)) for each contact pair: the evolution of damage is controlled by the strain state of the material by a scalar equivalent strain, which can be written for a 2D structure as follows: (4) Local equivalent strain of the contact: the equivalent strains of the contact pairs are averaged to contact points.Since damage is associated with the single contact, the mean equivalent strain of the contact pairs around a certain contact (calculated during the previous operation) is adopted as the local equivalent strain for the contact.Only the existing contact pairs of the scheme shown in Figure 4 are considered in the averaging.(5) Evaluation of the damage growth: the damage growth rate is de ned as a function of the local equivalent strain rate ε _ (Equation ( 3)). e damage functions f(D) given by Paas [48] is adopted in Bodin's L2R damage model (Equation ( 4)) and used in this study.(6) Evaluation of the damage and update: the value of Young's modulus (E E 0 (1 − D)) is updated, and nally, the values of the normal and tangential sti ness are updated (Equation ( 6)).Advances in Materials Science and Engineering (Figure 5), the frequency f 25 Hz, and the corresponding angular frequency ω 2πf is considered.
e material presents Young's modulus E 0 30 GPa and Poisson's ratio ] 1/3.e damage function (Equation ( 4)) proposed by Paas [48] is adopted in this analysis.e model parameters are α −2.25, β 4.0, and C 5.0 × 10 16 .e higher value of parameter C is adopted in order to ensure the fatigue life of applied loading condition is around 100 cycles (E/E 0 0.5 [62]), avoiding the jump-in-cycle procedure in this calculation.
e discrete element model is shown in Figure 6. e particle radius is r 1 mm, each row has 39 or 40 particles according to its position, and there are 69 rows in total.ere are 2726 particles and 7961 contacts.For this discretization level, boundary e ects are highly reduced; therefore, accurate stress and stain can be obtained in the whole geometry [56,57].According to Equation ( 6), the contact sti ness is calculated as k n 5.2 × 10 7 N/m and k t 5.2 × 10 4 N/m.A time step Δt 5 × 10 −7 s and a low viscous damping c n c t 0.0295 Ns/m are adopted in the simulations, to ensure stable calculations.e lower boundary is xed at each particle center in the vertical direction, and the sinusoidal velocity v ε 0 hωcos(ωt) is applied at the upper boundary.In order to apply the strain rate ε 0 1.25 × 10 −4 during all cycles, the amplitude of the displacement of the upper boundary is constant and equal to 14.7 mm.
e numerical results are compared with Bodin's L2R damage model predictions, as shown in Figures 7 and 8.For an imposed strain condition, the global damage increase and sti ness decrease are calculated from the decrease of the reaction force of the xed support, which are in good agreement with the theoretical L2R damage model (Equation (4)) predictions.No e ects of the particle diameter are observed under homogeneous conditions, which means the stress and strain of the contacts or contact pairs are calculated accurately for these discretization levels, and the equations of the L2R damage model have been successfully implemented.For the case of the bending beam with a stress gradient, the stress values of the contacts can be accurately calculated, when an appropriate ratio of beam size to particle size is adopted.
us, the numerical results of damage evolution and sti ness degradation for the beam case will be consistent with the theoretical predictions of the implemented continuum damage model.According to the European fatigue standards [62], when the sti ness ratio reaches its critical value 0.5, the corresponding fatigue life N f 105 cycles is identi ed. e calculation speed of the discrete element modelling depends on the discretization level, the discrete element   6 Advances in Materials Science and Engineering code, computing power, etc. is calculation was performed on an Intel Core i7-6700K 4.0 GHz processor with 4 cores and 8 threads running on a Windows 8.1 64 − bit operation system.e calculation required around 48 hours, using parallel calculation with all 8 CPU threads, to run 350 fatigue loading cycles when the particle diameter is d 2 mm (without jump-in-cycle procedure).More than 72 hours and around 12 hours were needed, respectively, for d 1 mm and d 4 mm.Hence, it is necessary to enhance the eciency of the discrete element fatigue modelling, in order to perform a large number of cycles with an appropriate discretization level.

Coupled DEM Model for Damage and Fatigue Crack Growth
In view of the crack initiation and propagation, the failure modes of quasi-brittle materials 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) (Figure 9).In the beginning of the lifetime, the material presents only intrinsic defects (microcracks, voids, etc.).Due to the e ect of the cyclic loading, these small defects tend to grow in size and quantity which damage the material, reducing its sti ness.
When the intrinsic 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 characterize the fracture behaviour.e crack nucleation stage and most of the short crack growth stage were shown to be well described by the continuum damage model.However, as the crack length increases, the decrease in the global sti ness becomes dominated by crack propagation.At this point, the continuum damage model failed, resulting in fast propagation due to the stress singularity at the crack tip.Hence, it is necessary to adopt a fatigue crack growth model (e.g., Paris' law) to estimate better the fracture behaviour during the end of stage II and stage III.

Model Implementation. Damage models (Section 2.1)
describe the e ect of distributed (micro-) defects over a certain region.e rupture is caused by the coalescence of these defects, giving rise to cracks which subsequently propagate.
is phenomenon is not well described by standard fatigue models, which su er from discretization e ects due to the stress gradients.
In the present work, a damage approach (Section 2.1) is adopted to describe the behaviour before contact rupture.
e rupture of a contact resulting in crack propagation is limited by a crack growth criterion (Section 2.2).

Local Identi cation of the Position of Crack Tips.
e propagation of a crack can be analysed as the creation or extension of the boundaries of a given geometry.In fracture mechanics, this transformation is usually controlled by the energy release during the process.Despite the di erent existing criteria of crack propagation, roughly a crack may be created or propagated where the stress (and/or strain) is maximized (Figure 10).In an elastic system, a simple veri cation of the local maximum value of the principal stress may be enough to identify potential localization of crack tips.However, during fatigue, the value of the stress tends to decrease due to the degradation of elastic properties of the material where the stress is concentrated.A better indicator, in this case, is shown to be the damage increment per cycle dD/dN c , which Advances in Materials Science and Engineering depends on the stress value but also on damage itself, and can be calculated numerically [19].Figure 11  In the present model, the possibility of crack propagation will only be considered on contacts, which locally maximizes the damage increment per cycle, as shown in Figure 11.For these two contacts at the crack tip, when the damage reaches D 1, indicating the total degradation of the sti ness, the energy release rate is calculated, as shown in Section 5.1.2.

Evaluation of the Range of the Energy Release Rate ΔG.
A damage value D 1 indicates the possibility of crack propagation, if it happens at a crack tip (as described in Section 5.1.1),or indicates an overestimated damage value.
is is usually the case at the neighbourhood of crack tips, for example, where the damage grows until unrealistic values due to the stress concentration.
In crack growth models, the energy release is considered to be localized exclusively at the crack tip.An overdamaged zone near a crack tip leads automatically to inconsistent evaluations of the energy release rates at the crack tip.In order to avoid this disturbance due to unphysical damage values, if damage reaches D 1 in a contact not identi ed as a crack tip, the value of D is automatically set to zero until this point eventually becomes a crack tip.Mathematically, this point is treated as an intact point, considering that it recovers its initial elastic properties.Physically, it indicates a scale decrease on the rupture process due to the proximity with a crack.As suggested in Figure 12(a), damage is de ned as the e ect of a certain number of defects inside a certain zone; however, in a smaller scale, there is only an intact material.In Figure 12(b), the damage value reaches D 1 in contact C 1 , a crack tip.e energy release rate is evaluated, and a potential crack extension is identi ed.It should be noted that the fracture process zone with a certain width around the crack tip can appear naturally in this model, thanks to the contribution of the continuum damage model  implemented in the discrete element codes.However, due to the heterogeneity of engineering materials, the crack will not follow strictly the crack path presented in Figure 12(b) [63].e real crack propagation direction can be modelled directly after the di erent material properties (e.g., aggregates, mortar, etc.) are given to the groups of contacts in the discrete element model.is improvement of the proposed model will be realized in the future study.
e principal components of the contact forces N and T and contact displacements δ n and δ t can be written as Equations ( 12) and (13).For a certain contact, when it is the rst contact (de ned clockwise) in a pair of contacts (Figure 13(a)), the principal components can be obtained by Equation (12), and when it is the second contact in a pair (Figure 13(b)), Equation ( 13) should be used to calculate the principal components: where ψ is de ned as the angle between (n, t) and the coordinate system associated with the principal stresses (Figure 3(c)).During a fatigue test, F I and δ I oscillate between a minimum level 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 I and δ I naturally vary from 0 to max F I and max δ I , respectively.In this case, the damage of a contact induces a maximum energy release rate [64] G max .e minimum energy release rate G min is equal to zero; consequently, the variation of the energy release rate is simply de ned as where N c D is the number of cycles to reach D 1 (total release of the contact energy) and g i is the surface of the triangle formed by the points (0, 0), (max δ I(i−1) , max F I(i−1) ), and (max δ I(i) , max F I(i) ), as shown in Figure 14. e second contact C 2 is not identi ed as a crack tip.e degradation of the elastic properties of the contact C 1 induces an increase of the contact force in C 2 , as shown in Figure 15.After contact C 1 is totally released, in the case of a crack propagation, contact C 2 becomes the new crack tip.
e value of the damage of C 2 is set back to zero D 0, and the value of the number of cycles N cD starts to be incremented.Once D 1 for contact C 2 , the surface of the full triangle (Figure 15) can then be computed, which is substituted into Equation ( 14) to estimate the range of the energy release rate.Following the same principle, the range of the energy release rate can be obtained systematically at the crack tip during crack propagation (discussed in Section 5.1.3).
Based on the relation between the energy release rate and the stress intensity factor in plane stress [64], the stress intensity range is simply de ned as ΔK E ΔG √ .

Crack Initiation and Propagation.
e evolution of the damage variable D characterizes the weakening of a contact before rupture.e rupture associated with the propagation of a crack is de ned by the value of da/dN c (Equation ( 5)): A contact which presents D 1 is de nitely broken (and ceases to exist) only if which represents a crack growth da equals to dcosψ/2 for one contact break.Otherwise, D is set to 0, and the fatigue loading continues, until the number of cycles N c increases su ciently to ful ll the rupture condition.is rupture criterion prevents deviations from the crack growth criterion induced by the damage model (Section 2.1) in conditions of a stress singularity at the crack tip.

Evaluation of the Stress Intensity Range ΔK.
In Figure 17, the numerical evaluation of the stress intensity range ΔK is compared to the theoretical result for di erent particle diameters d.When the crack length a is below 25 mm, the discretization levels that we have chosen can all provide an acceptable stress intensity range ΔK.Some deviations can be observed when the crack becomes larger (a > 25 mm).ese deviations tend to decrease for smaller values of d, which indicates a convergence of the data.e e ect of the particle size d tends to be neglectful of the de nition of the stress intensity range ΔK (Figure 17) and also of the crack growth rate da/dN c (Equation ( 15)).However, d de nes the propagation step of an existing crack and, according to Equation ( 16), a ects proportionally the number of cycles N c .It means that the choice of smaller d induces only more precision on the numerical calculation.

Sti ness Degradation.
For the condition of imposed stress with constant amplitude, the amplitude of the displacement at the edge of the sample tends to increase during the fatigue test. is behaviour is due to the decrease of the global sti ness of the sample caused by damage of the material and the propagation of the crack.e process of the sti ness degradation can be quanti ed by the ratio of the initial displacement amplitude and its instantaneous value at a time during the test, which is a function of the number of cycles N c .  Figure 18 shows the sti ness degradation process obtained for plates with di erent initial crack lengths a 0 under fatigue loading.e numerical results are compared to the theoretical predictions of the damage model (plate with a 0 0 mm) and the fatigue crack growth model (plates with a 0 > 0 mm).For the uncracked plate, the numerical result accurately follows the theoretical prediction of the continuum damage model.For a plate with a long crack, e.g., a 0 10 mm, the numerical result behaves according to the prediction of the fatigue crack growth model (Paris' law), as the crack growth is the dominant factor for sti ness decrease.However, a smaller crack (a 0 2 mm) deviates from fatigue crack growth because of the additional contribution of the material damage to the sti ness degradation.
e competition between the two mechanisms depends on the choice of the material parameters for the damage and fatigue crack growth models.Beyond the fatigue life, different parameters may a ect the sensitivity of the material to cracks.

Damage Field.
Figure 19(a) illustrates the damage distribution of a precracked plate with an initial crack length a 0 5 mm after 166 fatigue loading cycles.e red values correspond to D 1 and indicate the propagation of the crack.e deep blue values correspond to very low values of damage, which happens near the initial crack surface, where the stresses (and strains) are low during the whole test.Far from the crack, the damage value tends to be homogeneous, except two points on the top of the sample, due to boundary e ects.e damage tends naturally to be higher near the crack path, induced by the increase of the stresses (and strains).However, very close to the crack tip, an undamaged zone (D 0) is observed.e value itself should not be considered as a damage quantity, but it indicates where the damage model fails to describe the material behaviour, leading to an unrealistic fast rupture.It can be associated with a fracture process zone (FPZ) [65,66], acting as a bridging zone between the cracked region and uncracked region.
e increase of the number of cycles causes the evolution of the high damage zone due to the extension of the fatigue crack, as shown in Figure 19.e size of the fracture process zone seems to depend on the crack size.e reduction of the global sti ness of the structure (Figure 18) is Advances in Materials Science and Engineering indeed controlled by the propagation of the crack associated with the damage of the uncracked zone, as shown in Figure 19.

Conclusions
e damage and fatigue crack growth are normally modelled separately by di erent approaches.e damage is usually modelled by nonlocal approaches implemented in a nite element code, which may produce reasonable global sample behaviour, based on unrealistic material behaviour.e fatigue crack growth model is related to the local stress intensity range at the crack tip, which can produce reasonable rupture patterns but fails to predict the correct fatigue life and global behaviour when there is no crack or only small cracks in the sample.In order to reduce these limitations of two di erent approaches, a simple numerical scheme coupling damage and fracture mechanics in a discrete element environment was proposed.
e association of these di erent mechanical formulations allows the reproduction of experimental evidence: before material rupture by damage models and during crack propagation by crack growth models.In parallel, important drawbacks of each approach are avoided.
In this study, the local version of the continuum damage model was successfully implemented in the discrete element code and compared to the theoretical prediction showing good agreement under homogeneous stress conditions. is discrete element fatigue model would be a helpful enrichment for the discrete element method and shares as well the advantages of this numerical method in terms of the construction of the models with voids, imperfections, or heterogeneities and the simulation of crack initiation and propagation.e evaluation method of the stress intensity range during the damage evolution provides a way to couple both failure mechanisms. is feature allows crack initiation to be induced by localized damage and a progressive transition to a fracture behaviour with the crack propagation.Independent parameters for the fatigue damage model and fatigue crack growth model are admitted without any

4. 2 .Figure 3 :
Figure 3: (a) Contact displacements and (b) contact forces of respective adjacent particles.(c) Mean stresses and the orientation of their principal values [56].

Figure 6 :Figure 4 :
Figure 6: e discrete element model of an uncracked plate subjected to imposed sinusoidal strain.

Figure 8 :Figure 7 :
Figure 8: e theoretical and numerical predictions of the decrease of global sti ness for an uncracked plate subjected to fatigue loading.
shows an example of the damage increment per cycle obtained from DEM simulation for a center cracked plate subjected to sinusoidal fatigue loading.e maximum dD/dN c in Figure 11 is located at the crack tip.

Figure 10 :
Figure 10: Localization of the local maxima of the stress/strain for (a) a cracked plate and (b) a simply supported beam.

Figure 16 CrackFigure 12 :
Figure 12: (a) Scale e ect on the damage value and (b) contact points in the potential crack path.

Figure 13 :
Figure 13: Principal components of force and displacement for a certain contact being (a) the rst contact (de ned clockwise) in a pair and (b) the second contact in a pair.

Figure 14 :Figure 15 :
Figure 14: Evaluation of the energy release at the crack tip (contact C 1 ) during a fatigue test.

Figure 19 :
Figure 19: Damage distribution and crack propagation of a precracked plate with the initial crack length a 0 5 mm after (a) 166, (b) 179, (c) 187, and (d) 215 fatigue loading cycles.
Results for a Center Cracked Plate 5.2.1.Geometry, Loading, and Material Properties.