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

Authors: Xiaofeng Gao; Georg Koval; Cyrille Chazallon
Journal: Advances in Materials Science and Engineering
 (2019)

---

## Abstract

The 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. The proposed model achieves the implementation of a continuum damage model in a discrete element code, which is a helpful enrichment of this numerical method. The evaluation method of the stress intensity range during the damage evolution provides a way to couple both failure mechanisms. This 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. The 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.


---

## Body

## 1. Introduction

The size and boundary effects of quasi-brittle fractures have been studied systematically in the past decades [1–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–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–17]. The 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. The 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].The fatigue behaviour can be studied experimentally and numerically. The 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. The nonlocal continuum damage models [19–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]. Therefore, the residual strength cannot be determined using these models. Moës et al. proposed an alternative nonlocal damage model called the Thick Level Set (TLS) approach [23–26] and implemented it in an extended finite element code to model the damage growth in solids. The 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. The damage growth is expressed as a level set propagation [23]. This 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 Thick 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. The Thick 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. The cohesive zone models represent quasi-brittle behaviours with good accuracy but require extra equations to determine the crack path. These models were widely used to predict the fatigue crack growth [29–31] and also improved to study the damage initiation and evolution in the interconnections [15, 32–34].Both continuous and discrete numerical methods have their own advantages and shortcomings in the modelling of fracture and damage phenomena. The finite element method (FEM) is cumbersome because the mesh requires to be updated to match the geometry of the discontinuity [35]. The mesoscopic models developed for studying the influence of material composition on the overall behaviour [36–38] have extended the capacity and quality of the finite element modelling but not yet overcome its shortcomings. The 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]. The 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. The 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. The continuum damage model is implemented in a discrete element code and verifies the theoretical prediction. In Section5, the continuum damage model is coupled with a fatigue crack growth model. This 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.


## 2. Damage and Fatigue Crack Growth Models

### 2.1. 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 variableD. 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. The stress-strain relation can be written as follows:(1)σij=Cijklεkl=1−DCijkl0εkl,where σij and εkl are the elastic stress and strain tensor components and Cijkl and Cijkl0 are the damaged and initial (elastic) secant stiffness of the material. Cijkl0 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 quasi-brittle 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. The 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)ε˜=∑i=1i=3σi+E1−D2.In the nonlocal version of Bodin’s model, the “local” equivalent strainε˜ is replaced by its weighted average strain. The expressions for weighted average strain calculation can be found in [20, 21].The rate of damage growth is defined as a function of the local equivalent strain rate:(3)D˙=fDε˜βε˜˙,where ε˜˙=ε˜i+1−ε˜i is the increment of local equivalent strain from t=i to t=i+1 and fD is a function of damage. The exponent β is a material parameter, which relates to the slope (−1−β) of the S-N curve in the log-log diagram.The function of damagefD 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)fD=CDα.


### 2.2. Fatigue Crack Growth Models (FCGMs)

The 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.” The 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.


#### 2.2.1. Paris’ Law

Fatigue crack growth in a wide variety of brittle and quasi-brittle materials [50–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/dNc via a power law, with ΔK=Kmax−Kmin. The basic formula is as follows:(5)dadNc=cΔKm,ΔK>ΔKth,Kmax<Kc,where da/dNc is the crack growth rate, a is the crack length, and Nc is the number of load cycles; c and m are the material parameters dependent on the material composition, environment, stress ratio, etc. [53]. ΔKth and Kc 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 ΔKth, and the maximum value of the stress intensity factor Kmax remains below the material fracture toughness Kc.


## 2.1. 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 variableD. 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. The stress-strain relation can be written as follows:(1)σij=Cijklεkl=1−DCijkl0εkl,where σij and εkl are the elastic stress and strain tensor components and Cijkl and Cijkl0 are the damaged and initial (elastic) secant stiffness of the material. Cijkl0 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 quasi-brittle 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. The 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)ε˜=∑i=1i=3σi+E1−D2.In the nonlocal version of Bodin’s model, the “local” equivalent strainε˜ is replaced by its weighted average strain. The expressions for weighted average strain calculation can be found in [20, 21].The rate of damage growth is defined as a function of the local equivalent strain rate:(3)D˙=fDε˜βε˜˙,where ε˜˙=ε˜i+1−ε˜i is the increment of local equivalent strain from t=i to t=i+1 and fD is a function of damage. The exponent β is a material parameter, which relates to the slope (−1−β) of the S-N curve in the log-log diagram.The function of damagefD 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)fD=CDα.


## 2.2. Fatigue Crack Growth Models (FCGMs)

The 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.” The 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.


### 2.2.1. Paris’ Law

Fatigue crack growth in a wide variety of brittle and quasi-brittle materials [50–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/dNc via a power law, with ΔK=Kmax−Kmin. The basic formula is as follows:(5)dadNc=cΔKm,ΔK>ΔKth,Kmax<Kc,where da/dNc is the crack growth rate, a is the crack length, and Nc is the number of load cycles; c and m are the material parameters dependent on the material composition, environment, stress ratio, etc. [53]. ΔKth and Kc 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 ΔKth, and the maximum value of the stress intensity factor Kmax remains below the material fracture toughness Kc.


## 2.2.1. Paris’ Law

Fatigue crack growth in a wide variety of brittle and quasi-brittle materials [50–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/dNc via a power law, with ΔK=Kmax−Kmin. The basic formula is as follows:(5)dadNc=cΔKm,ΔK>ΔKth,Kmax<Kc,where da/dNc is the crack growth rate, a is the crack length, and Nc is the number of load cycles; c and m are the material parameters dependent on the material composition, environment, stress ratio, etc. [53]. ΔKth and Kc 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 ΔKth, and the maximum value of the stress intensity factor Kmax remains below the material fracture toughness Kc.


## 3. Discrete Element Method

The discrete element method (DEM) was originally developed by Cundall for modelling granular and particulate systems [54]. The 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. The constitutive parameters for the contacts between the discrete elements, such as stiffness and strength, influence the behaviour of the model at the macroscopic scale. These 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 Figure1(b). The close-packed assembly is composed of particles with identical size. The mechanical properties of the specimen (Figure 1(a)) depend on the mechanical parameters of the particles. The 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. The heterogeneous model can be constructed by representing different materials (aggregates, mortar, etc.) as clusters of discrete elements with different parameters [55]. This close-packed assembly has been used by Le et al. [55–59] and Liu et al. [60] for quasi-brittle failure of asphalt concrete, rocks, and crystals.Figure 1
Equivalent (a) continuous and (b) discrete media (close-packed assembly).
(a)
(b)


### 3.1. 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). The contact behaviour follows a standard linear spring and dash-pot model. When the values of the damping parameters, cn and ct, are a sufficiently small fraction of mkn (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.Figure 2
(a) Unit cell and (b) contact law (modified from [56]).
(a)
(b)Young’s modulusE 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, kn and kt, respectively) has been established by the researchers [56, 61], which for plane stress is as follows:(6)kn=Et31−ν,kt=1−3ν1+νkn=1−3ν31−ν2Et,where t is the thickness of the discrete element model. For plane strain condition, this direct relationship can also be found in [56, 61].


### 3.2. Definition of Strain and Stress in a Discrete Medium

#### 3.2.1. Mean Strains and Stresses

The 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 defined, where t virtually connects both contacts for which n is an orthogonal axis (Figure 3(a)). The 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:(7)ε¯nn=14δnik3+δnjk3+δtik−δtjk23d,ε¯tt=12δnik+δnjk−δtik3+δtjk31d.Figure 3
(a) Contact displacements and (b) contact forces of respective adjacent particles. (c) Mean stresses and the orientation of their principal values [56].
(a)
(b)
(c)The normal and tangential components of each contact (Nik, Tik and Njk, Tjk, 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:(8)σ¯n=Nik3+Tik+Njk3−Tjk2dt,σ¯t=−Nik+Tik3+Njk+Tjk32dt.


#### 3.2.2. Principal Stresses

The stress tensor (in two dimensions) can be defined by the values of the principal stressesσ¯I and σ¯II and their orientations. Hence, ψ is defined 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:(9)σ¯I=σ¯n+σ¯ttanψ,σ¯II=σ¯n−σ¯ttanψ.The 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:(10)ψ=−12arctan2σ¯tA,where A=E/1−νε¯nn+ε¯tt−2σ¯n.


## 3.1. 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). The contact behaviour follows a standard linear spring and dash-pot model. When the values of the damping parameters, cn and ct, are a sufficiently small fraction of mkn (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.Figure 2
(a) Unit cell and (b) contact law (modified from [56]).
(a)
(b)Young’s modulusE 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, kn and kt, respectively) has been established by the researchers [56, 61], which for plane stress is as follows:(6)kn=Et31−ν,kt=1−3ν1+νkn=1−3ν31−ν2Et,where t is the thickness of the discrete element model. For plane strain condition, this direct relationship can also be found in [56, 61].


## 3.2. Definition of Strain and Stress in a Discrete Medium

### 3.2.1. Mean Strains and Stresses

The 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 defined, where t virtually connects both contacts for which n is an orthogonal axis (Figure 3(a)). The 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:(7)ε¯nn=14δnik3+δnjk3+δtik−δtjk23d,ε¯tt=12δnik+δnjk−δtik3+δtjk31d.Figure 3
(a) Contact displacements and (b) contact forces of respective adjacent particles. (c) Mean stresses and the orientation of their principal values [56].
(a)
(b)
(c)The normal and tangential components of each contact (Nik, Tik and Njk, Tjk, 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:(8)σ¯n=Nik3+Tik+Njk3−Tjk2dt,σ¯t=−Nik+Tik3+Njk+Tjk32dt.


### 3.2.2. Principal Stresses

The stress tensor (in two dimensions) can be defined by the values of the principal stressesσ¯I and σ¯II and their orientations. Hence, ψ is defined 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:(9)σ¯I=σ¯n+σ¯ttanψ,σ¯II=σ¯n−σ¯ttanψ.The 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:(10)ψ=−12arctan2σ¯tA,where A=E/1−νε¯nn+ε¯tt−2σ¯n.


## 3.2.1. Mean Strains and Stresses

The 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 defined, where t virtually connects both contacts for which n is an orthogonal axis (Figure 3(a)). The 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:(7)ε¯nn=14δnik3+δnjk3+δtik−δtjk23d,ε¯tt=12δnik+δnjk−δtik3+δtjk31d.Figure 3
(a) Contact displacements and (b) contact forces of respective adjacent particles. (c) Mean stresses and the orientation of their principal values [56].
(a)
(b)
(c)The normal and tangential components of each contact (Nik, Tik and Njk, Tjk, 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:(8)σ¯n=Nik3+Tik+Njk3−Tjk2dt,σ¯t=−Nik+Tik3+Njk+Tjk32dt.


## 3.2.2. Principal Stresses

The stress tensor (in two dimensions) can be defined by the values of the principal stressesσ¯I and σ¯II and their orientations. Hence, ψ is defined 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:(9)σ¯I=σ¯n+σ¯ttanψ,σ¯II=σ¯n−σ¯ttanψ.The 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:(10)ψ=−12arctan2σ¯tA,where A=E/1−νε¯nn+ε¯tt−2σ¯n.


## 4. 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 stiffness is monitored, which is defined 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). The model was implemented in a finite element code, along with a self-adaptive 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 close-packed assembly (regular hexagonal packing) is adopted due to the direct relationship between the macroscopic parameters (Young’s modulusE and Poisson’s ratio ν) and discrete elastic parameters (normal and tangential stiffness, kn and kt, 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 stiffness.


### 4.1. Model Implementation

The damage of all the contacts in the assembly is calculated in the same way and updated at the same time. The evaluation of the damage per contact can be summarized by the following operations:(1)
DEM elastic analysis and identification of the stress and strain fields (Section3.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 defined 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:(11)ε˜=σI+E01−D2+σII+E01−D2.(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 Figure4 are considered in the averaging.(5)
Evaluation of the damage growth: the damage growth rate is defined as a function of the local equivalent strain rateε˜˙ (Equation (3)). The damage functions fD 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=E01−D) is updated, and finally, the values of the normal and tangential stiffness are updated (Equation (6)).Figure 4
Contact pairs associated with the average of the equivalent strain of the contactC.


### 4.2. Model Verification

An uncracked plate with dimensionsw=80mm and h=117.8mm subjected to sinusoidal fatigue loading with the amplitude ε0=1.25×10−4 (Figure 5), the frequency f=25Hz, and the corresponding angular frequency ω=2πf is considered. The material presents Young’s modulus E0=30GPa and Poisson’s ratio ν=1/3. The damage function (Equation (4)) proposed by Paas [48] is adopted in this analysis. The model parameters are α=−2.25, β=4.0, and C=5.0×1016. The higher value of parameter C is adopted in order to ensure the fatigue life of applied loading condition is around 100 cycles (E/E0=0.5 [62]), avoiding the jump-in-cycle procedure in this calculation.Figure 5
An uncracked plate subjected to imposed sinusoidal strain.The discrete element model is shown in Figure6. The particle radius is r=1mm, each row has 39 or 40 particles according to its position, and there are 69 rows in total. There are 2726 particles and 7961 contacts. For this discretization level, boundary effects are highly reduced; therefore, accurate stress and stain can be obtained in the whole geometry [56, 57]. According to Equation (6), the contact stiffness is calculated as kn=5.2×107N/m and kt=5.2×104N/m. A time step Δt=5×10−7s and a low viscous damping cn=ct=0.0295Ns/m are adopted in the simulations, to ensure stable calculations. The lower boundary is fixed at each particle center in the vertical direction, and the sinusoidal velocity v=ε0hω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.7mm.Figure 6
The discrete element model of an uncracked plate subjected to imposed sinusoidal strain.The numerical results are compared with Bodin’s L2R damage model predictions, as shown in Figures7 and 8. For an imposed strain condition, the global damage increase and stiffness decrease are calculated from the decrease of the reaction force of the fixed support, which are in good agreement with the theoretical L2R damage model (Equation (4)) predictions. No effects 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. Thus, the numerical results of damage evolution and stiffness 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 stiffness ratio reaches its critical value 0.5, the corresponding fatigue life Nf=105cycles is identified.Figure 7
The theoretical and numerical predictions of the increase of damage for an uncracked plate subjected to fatigue loading.Figure 8
The theoretical and numerical predictions of the decrease of global stiffness for an uncracked plate subjected to fatigue loading.The calculation speed of the discrete element modelling depends on the discretization level, the discrete element code, computing power, etc. This calculation was performed on an Intel Core i7-6700K4.0GHz processor with 4 cores and 8 threads running on a Windows 8.1 64−bit operation system. The 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=2mm (without jump-in-cycle procedure). More than 72 hours and around 12 hours were needed, respectively, for d=1mm and d=4mm. Hence, it is necessary to enhance the efficiency of the discrete element fatigue modelling, in order to perform a large number of cycles with an appropriate discretization level.


## 4.1. Model Implementation

The damage of all the contacts in the assembly is calculated in the same way and updated at the same time. The evaluation of the damage per contact can be summarized by the following operations:(1)
DEM elastic analysis and identification of the stress and strain fields (Section3.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 defined 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:(11)ε˜=σI+E01−D2+σII+E01−D2.(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 Figure4 are considered in the averaging.(5)
Evaluation of the damage growth: the damage growth rate is defined as a function of the local equivalent strain rateε˜˙ (Equation (3)). The damage functions fD 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=E01−D) is updated, and finally, the values of the normal and tangential stiffness are updated (Equation (6)).Figure 4
Contact pairs associated with the average of the equivalent strain of the contactC.


## 4.2. Model Verification

An uncracked plate with dimensionsw=80mm and h=117.8mm subjected to sinusoidal fatigue loading with the amplitude ε0=1.25×10−4 (Figure 5), the frequency f=25Hz, and the corresponding angular frequency ω=2πf is considered. The material presents Young’s modulus E0=30GPa and Poisson’s ratio ν=1/3. The damage function (Equation (4)) proposed by Paas [48] is adopted in this analysis. The model parameters are α=−2.25, β=4.0, and C=5.0×1016. The higher value of parameter C is adopted in order to ensure the fatigue life of applied loading condition is around 100 cycles (E/E0=0.5 [62]), avoiding the jump-in-cycle procedure in this calculation.Figure 5
An uncracked plate subjected to imposed sinusoidal strain.The discrete element model is shown in Figure6. The particle radius is r=1mm, each row has 39 or 40 particles according to its position, and there are 69 rows in total. There are 2726 particles and 7961 contacts. For this discretization level, boundary effects are highly reduced; therefore, accurate stress and stain can be obtained in the whole geometry [56, 57]. According to Equation (6), the contact stiffness is calculated as kn=5.2×107N/m and kt=5.2×104N/m. A time step Δt=5×10−7s and a low viscous damping cn=ct=0.0295Ns/m are adopted in the simulations, to ensure stable calculations. The lower boundary is fixed at each particle center in the vertical direction, and the sinusoidal velocity v=ε0hω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.7mm.Figure 6
The discrete element model of an uncracked plate subjected to imposed sinusoidal strain.The numerical results are compared with Bodin’s L2R damage model predictions, as shown in Figures7 and 8. For an imposed strain condition, the global damage increase and stiffness decrease are calculated from the decrease of the reaction force of the fixed support, which are in good agreement with the theoretical L2R damage model (Equation (4)) predictions. No effects 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. Thus, the numerical results of damage evolution and stiffness 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 stiffness ratio reaches its critical value 0.5, the corresponding fatigue life Nf=105cycles is identified.Figure 7
The theoretical and numerical predictions of the increase of damage for an uncracked plate subjected to fatigue loading.Figure 8
The theoretical and numerical predictions of the decrease of global stiffness for an uncracked plate subjected to fatigue loading.The calculation speed of the discrete element modelling depends on the discretization level, the discrete element code, computing power, etc. This calculation was performed on an Intel Core i7-6700K4.0GHz processor with 4 cores and 8 threads running on a Windows 8.1 64−bit operation system. The 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=2mm (without jump-in-cycle procedure). More than 72 hours and around 12 hours were needed, respectively, for d=1mm and d=4mm. Hence, it is necessary to enhance the efficiency of the discrete element fatigue modelling, in order to perform a large number of cycles with an appropriate discretization level.


## 5. 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) (Figure9). In the beginning of the lifetime, the material presents only intrinsic defects (microcracks, 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 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.Figure 9
The four stages of the fatigue failure process [16].The 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 stiffness 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.


### 5.1. Model Implementation

Damage models (Section2.1) describe the effect of distributed (micro-) defects over a certain region. The rupture is caused by the coalescence of these defects, giving rise to cracks which subsequently propagate. This phenomenon is not well described by standard fatigue models, which suffer from discretization effects due to the stress gradients.In the present work, a damage approach (Section2.1) is adopted to describe the behaviour before contact rupture. The rupture of a contact resulting in crack propagation is limited by a crack growth criterion (Section 2.2).


#### 5.1.1. Local Identification of the Position of Crack Tips

The 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 different existing criteria of crack propagation, roughly a crack may be created or propagated where the stress (and/or strain) is maximized (Figure10).Figure 10
Localization of the local maxima of the stress/strain for (a) a cracked plate and (b) a simply supported beam.
(a)
(b)In an elastic system, a simple verification 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 cycledD/dNc, which depends on the stress value but also on damage itself, and can be calculated numerically [19]. Figure 11 shows an example of the damage increment per cycle obtained from DEM simulation for a center cracked plate subjected to sinusoidal fatigue loading. The maximum dD/dNc in Figure 11 is located at the crack tip.Figure 11
Damage increment per cycle of a center cracked plate.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 Figure11. For these two contacts at the crack tip, when the damage reaches D=1, indicating the total degradation of the stiffness, the energy release rate is calculated, as shown in Section 5.1.2.


#### 5.1.2. Evaluation of the Range of the Energy Release RateΔG

A damage valueD=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. This 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 reachesD=1 in a contact not identified 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 defined as the effect 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 C1, a crack tip. The energy release rate is evaluated, and a potential crack extension is identified. 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]. The real crack propagation direction can be modelled directly after the different material properties (e.g., aggregates, mortar, etc.) are given to the groups of contacts in the discrete element model. This improvement of the proposed model will be realized in the future study.Figure 12
(a) Scale effect on the damage value and (b) contact points in the potential crack path.
(a)
(b)The principal components of the contact forcesN and T and contact displacements δn and δt can be written as Equations (12) and (13). For a certain contact, when it is the first contact (defined 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:(12)FI=Nsinπ/3+ψ−Tcosπ/3+ψ,δI=δnsinπ/3+ψ−δtcosπ/3+ψ,(13)FI=Nsinπ/3−ψ+Tcosπ/3−ψ,δI=δnsinπ/3−ψ+δtcosπ/3−ψ,where ψ is defined as the angle between (n, t) and the coordinate system associated with the principal stresses (Figure 3(c)). During a fatigue test, FI 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 FI and δI naturally vary from 0 to maxFI and maxδI, respectively. In this case, the damage of a contact induces a maximum energy release rate [64] Gmax. The minimum energy release rate Gmin is equal to zero; consequently, the variation of the energy release rate is simply defined as(14)ΔG=∑1NcDgid2cosψ/2,where NcD is the number of cycles to reach D=1 (total release of the contact energy) and gi is the surface of the triangle formed by the points 0,0, maxδIi−1,maxFIi−1, and maxδIi,maxFIi, as shown in Figure 14.Figure 13
Principal components of force and displacement for a certain contact being (a) the first contact (defined clockwise) in a pair and (b) the second contact in a pair.
(a)
(b)Figure 14
Evaluation of the energy release at the crack tip (contactC1) during a fatigue test.The second contactC2 is not identified as a crack tip. The degradation of the elastic properties of the contact C1 induces an increase of the contact force in C2, as shown in Figure 15. After contact C1 is totally released, in the case of a crack propagation, contact C2 becomes the new crack tip. The value of the damage of C2 is set back to zero D=0, and the value of the number of cycles NcD starts to be incremented. Once D=1 for contact C2, 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).Figure 15
Evaluation of the energy release at the second contact close to the crack tip during a fatigue test.Based on the relation between the energy release rate and the stress intensity factor in plane stress [64], the stress intensity range is simply defined as ΔK=EΔG.


#### 5.1.3. Crack Initiation and Propagation

The evolution of the damage variableD characterizes the weakening of a contact before rupture. The rupture associated with the propagation of a crack is defined by the value of da/dNc (Equation (5)):(15)dadNc=cΔKm.A contact which presentsD=1 is definitely broken (and ceases to exist) only if(16)dadNcNc≥dcosψ2,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 Nc increases sufficiently to fulfill the rupture condition. This 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.


### 5.2. Numerical Modelling Results for a Center Cracked Plate

#### 5.2.1. Geometry, Loading, and Material Properties

Figure16 shows a center cracked plate with dimensions w=80mm and L=120mm, which is tested under sinusoidal fatigue loading with an amplitude σ0=1.25MPa and a frequency f=25Hz. The material presents Young’s modulus E0=30GPa and Poisson’s ratio ν=1/3. A viscous damping cn=ct=0.0295Ns/m and a time step Δt=5×10−7s are adopted in the simulations. The damage model parameters are α=−2.25, β=4.0, and γ=5.0×1016; the fatigue crack growth model parameters are c=3.0×10−10m/cyclePamm and m=1.03. Negligible effects of the frequency indicate very low dynamic effects (quasi-static regime) on the presented results.Figure 16
Center cracked plate subjected to imposed sinusoidal stress.


#### 5.2.2. Evaluation of the Stress Intensity RangeΔK

In Figure17, the numerical evaluation of the stress intensity range ΔK is compared to the theoretical result for different particle diameters d. When the crack length a is below 25mm, 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>25mm). These deviations tend to decrease for smaller values of d, which indicates a convergence of the data.Figure 17
Comparison of the range of stress intensity factorΔK calculated using numerical and theoretical methods for different crack lengths a (a0=5mm) and particle diameters d.The effect of the particle sized tends to be neglectful of the definition of the stress intensity range ΔK (Figure 17) and also of the crack growth rate da/dNc (Equation (15)). However, d defines the propagation step of an existing crack and, according to Equation (16), affects proportionally the number of cycles Nc. It means that the choice of smaller d induces only more precision on the numerical calculation.


#### 5.2.3. Stiffness 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. This behaviour is due to the decrease of the global stiffness of the sample caused by damage of the material and the propagation of the crack. The process of the stiffness degradation can be quantified 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 cyclesNc.Figure18 shows the stiffness degradation process obtained for plates with different initial crack lengths a0 under fatigue loading. The numerical results are compared to the theoretical predictions of the damage model (plate with a0=0mm) and the fatigue crack growth model (plates with a0>0mm). 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., a0=10mm, 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 stiffness decrease. However, a smaller crack (a0=2mm) deviates from fatigue crack growth because of the additional contribution of the material damage to the stiffness degradation. The 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 affect the sensitivity of the material to cracks.Figure 18
Stiffness degradation as a function of the number of cyclesNc for different initial crack sizes.


#### 5.2.4. Damage Field

Figure19(a) illustrates the damage distribution of a precracked plate with an initial crack length a0=5mm after 166 fatigue loading cycles. The red values correspond to D=1 and indicate the propagation of the crack. The 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 effects. The 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. The 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.Figure 19
Damage distribution and crack propagation of a precracked plate with the initial crack lengtha0=5mm after (a) 166, (b) 179, (c) 187, and (d) 215 fatigue loading cycles.
(a)
(b)
(c)
(d)The 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 Figure19. The size of the fracture process zone seems to depend on the crack size. The reduction of the global stiffness of the structure (Figure 18) is indeed controlled by the propagation of the crack associated with the damage of the uncracked zone, as shown in Figure 19.


## 5.1. Model Implementation

Damage models (Section2.1) describe the effect of distributed (micro-) defects over a certain region. The rupture is caused by the coalescence of these defects, giving rise to cracks which subsequently propagate. This phenomenon is not well described by standard fatigue models, which suffer from discretization effects due to the stress gradients.In the present work, a damage approach (Section2.1) is adopted to describe the behaviour before contact rupture. The rupture of a contact resulting in crack propagation is limited by a crack growth criterion (Section 2.2).


### 5.1.1. Local Identification of the Position of Crack Tips

The 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 different existing criteria of crack propagation, roughly a crack may be created or propagated where the stress (and/or strain) is maximized (Figure10).Figure 10
Localization of the local maxima of the stress/strain for (a) a cracked plate and (b) a simply supported beam.
(a)
(b)In an elastic system, a simple verification 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 cycledD/dNc, which depends on the stress value but also on damage itself, and can be calculated numerically [19]. Figure 11 shows an example of the damage increment per cycle obtained from DEM simulation for a center cracked plate subjected to sinusoidal fatigue loading. The maximum dD/dNc in Figure 11 is located at the crack tip.Figure 11
Damage increment per cycle of a center cracked plate.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 Figure11. For these two contacts at the crack tip, when the damage reaches D=1, indicating the total degradation of the stiffness, the energy release rate is calculated, as shown in Section 5.1.2.


### 5.1.2. Evaluation of the Range of the Energy Release RateΔG

A damage valueD=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. This 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 reachesD=1 in a contact not identified 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 defined as the effect 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 C1, a crack tip. The energy release rate is evaluated, and a potential crack extension is identified. 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]. The real crack propagation direction can be modelled directly after the different material properties (e.g., aggregates, mortar, etc.) are given to the groups of contacts in the discrete element model. This improvement of the proposed model will be realized in the future study.Figure 12
(a) Scale effect on the damage value and (b) contact points in the potential crack path.
(a)
(b)The principal components of the contact forcesN and T and contact displacements δn and δt can be written as Equations (12) and (13). For a certain contact, when it is the first contact (defined 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:(12)FI=Nsinπ/3+ψ−Tcosπ/3+ψ,δI=δnsinπ/3+ψ−δtcosπ/3+ψ,(13)FI=Nsinπ/3−ψ+Tcosπ/3−ψ,δI=δnsinπ/3−ψ+δtcosπ/3−ψ,where ψ is defined as the angle between (n, t) and the coordinate system associated with the principal stresses (Figure 3(c)). During a fatigue test, FI 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 FI and δI naturally vary from 0 to maxFI and maxδI, respectively. In this case, the damage of a contact induces a maximum energy release rate [64] Gmax. The minimum energy release rate Gmin is equal to zero; consequently, the variation of the energy release rate is simply defined as(14)ΔG=∑1NcDgid2cosψ/2,where NcD is the number of cycles to reach D=1 (total release of the contact energy) and gi is the surface of the triangle formed by the points 0,0, maxδIi−1,maxFIi−1, and maxδIi,maxFIi, as shown in Figure 14.Figure 13
Principal components of force and displacement for a certain contact being (a) the first contact (defined clockwise) in a pair and (b) the second contact in a pair.
(a)
(b)Figure 14
Evaluation of the energy release at the crack tip (contactC1) during a fatigue test.The second contactC2 is not identified as a crack tip. The degradation of the elastic properties of the contact C1 induces an increase of the contact force in C2, as shown in Figure 15. After contact C1 is totally released, in the case of a crack propagation, contact C2 becomes the new crack tip. The value of the damage of C2 is set back to zero D=0, and the value of the number of cycles NcD starts to be incremented. Once D=1 for contact C2, 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).Figure 15
Evaluation of the energy release at the second contact close to the crack tip during a fatigue test.Based on the relation between the energy release rate and the stress intensity factor in plane stress [64], the stress intensity range is simply defined as ΔK=EΔG.


### 5.1.3. Crack Initiation and Propagation

The evolution of the damage variableD characterizes the weakening of a contact before rupture. The rupture associated with the propagation of a crack is defined by the value of da/dNc (Equation (5)):(15)dadNc=cΔKm.A contact which presentsD=1 is definitely broken (and ceases to exist) only if(16)dadNcNc≥dcosψ2,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 Nc increases sufficiently to fulfill the rupture condition. This 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.


## 5.1.1. Local Identification of the Position of Crack Tips

The 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 different existing criteria of crack propagation, roughly a crack may be created or propagated where the stress (and/or strain) is maximized (Figure10).Figure 10
Localization of the local maxima of the stress/strain for (a) a cracked plate and (b) a simply supported beam.
(a)
(b)In an elastic system, a simple verification 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 cycledD/dNc, which depends on the stress value but also on damage itself, and can be calculated numerically [19]. Figure 11 shows an example of the damage increment per cycle obtained from DEM simulation for a center cracked plate subjected to sinusoidal fatigue loading. The maximum dD/dNc in Figure 11 is located at the crack tip.Figure 11
Damage increment per cycle of a center cracked plate.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 Figure11. For these two contacts at the crack tip, when the damage reaches D=1, indicating the total degradation of the stiffness, the energy release rate is calculated, as shown in Section 5.1.2.


## 5.1.2. Evaluation of the Range of the Energy Release RateΔG

A damage valueD=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. This 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 reachesD=1 in a contact not identified 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 defined as the effect 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 C1, a crack tip. The energy release rate is evaluated, and a potential crack extension is identified. 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]. The real crack propagation direction can be modelled directly after the different material properties (e.g., aggregates, mortar, etc.) are given to the groups of contacts in the discrete element model. This improvement of the proposed model will be realized in the future study.Figure 12
(a) Scale effect on the damage value and (b) contact points in the potential crack path.
(a)
(b)The principal components of the contact forcesN and T and contact displacements δn and δt can be written as Equations (12) and (13). For a certain contact, when it is the first contact (defined 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:(12)FI=Nsinπ/3+ψ−Tcosπ/3+ψ,δI=δnsinπ/3+ψ−δtcosπ/3+ψ,(13)FI=Nsinπ/3−ψ+Tcosπ/3−ψ,δI=δnsinπ/3−ψ+δtcosπ/3−ψ,where ψ is defined as the angle between (n, t) and the coordinate system associated with the principal stresses (Figure 3(c)). During a fatigue test, FI 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 FI and δI naturally vary from 0 to maxFI and maxδI, respectively. In this case, the damage of a contact induces a maximum energy release rate [64] Gmax. The minimum energy release rate Gmin is equal to zero; consequently, the variation of the energy release rate is simply defined as(14)ΔG=∑1NcDgid2cosψ/2,where NcD is the number of cycles to reach D=1 (total release of the contact energy) and gi is the surface of the triangle formed by the points 0,0, maxδIi−1,maxFIi−1, and maxδIi,maxFIi, as shown in Figure 14.Figure 13
Principal components of force and displacement for a certain contact being (a) the first contact (defined clockwise) in a pair and (b) the second contact in a pair.
(a)
(b)Figure 14
Evaluation of the energy release at the crack tip (contactC1) during a fatigue test.The second contactC2 is not identified as a crack tip. The degradation of the elastic properties of the contact C1 induces an increase of the contact force in C2, as shown in Figure 15. After contact C1 is totally released, in the case of a crack propagation, contact C2 becomes the new crack tip. The value of the damage of C2 is set back to zero D=0, and the value of the number of cycles NcD starts to be incremented. Once D=1 for contact C2, 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).Figure 15
Evaluation of the energy release at the second contact close to the crack tip during a fatigue test.Based on the relation between the energy release rate and the stress intensity factor in plane stress [64], the stress intensity range is simply defined as ΔK=EΔG.


## 5.1.3. Crack Initiation and Propagation

The evolution of the damage variableD characterizes the weakening of a contact before rupture. The rupture associated with the propagation of a crack is defined by the value of da/dNc (Equation (5)):(15)dadNc=cΔKm.A contact which presentsD=1 is definitely broken (and ceases to exist) only if(16)dadNcNc≥dcosψ2,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 Nc increases sufficiently to fulfill the rupture condition. This 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.


## 5.2. Numerical Modelling Results for a Center Cracked Plate

### 5.2.1. Geometry, Loading, and Material Properties

Figure16 shows a center cracked plate with dimensions w=80mm and L=120mm, which is tested under sinusoidal fatigue loading with an amplitude σ0=1.25MPa and a frequency f=25Hz. The material presents Young’s modulus E0=30GPa and Poisson’s ratio ν=1/3. A viscous damping cn=ct=0.0295Ns/m and a time step Δt=5×10−7s are adopted in the simulations. The damage model parameters are α=−2.25, β=4.0, and γ=5.0×1016; the fatigue crack growth model parameters are c=3.0×10−10m/cyclePamm and m=1.03. Negligible effects of the frequency indicate very low dynamic effects (quasi-static regime) on the presented results.Figure 16
Center cracked plate subjected to imposed sinusoidal stress.


### 5.2.2. Evaluation of the Stress Intensity RangeΔK

In Figure17, the numerical evaluation of the stress intensity range ΔK is compared to the theoretical result for different particle diameters d. When the crack length a is below 25mm, 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>25mm). These deviations tend to decrease for smaller values of d, which indicates a convergence of the data.Figure 17
Comparison of the range of stress intensity factorΔK calculated using numerical and theoretical methods for different crack lengths a (a0=5mm) and particle diameters d.The effect of the particle sized tends to be neglectful of the definition of the stress intensity range ΔK (Figure 17) and also of the crack growth rate da/dNc (Equation (15)). However, d defines the propagation step of an existing crack and, according to Equation (16), affects proportionally the number of cycles Nc. It means that the choice of smaller d induces only more precision on the numerical calculation.


### 5.2.3. Stiffness 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. This behaviour is due to the decrease of the global stiffness of the sample caused by damage of the material and the propagation of the crack. The process of the stiffness degradation can be quantified 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 cyclesNc.Figure18 shows the stiffness degradation process obtained for plates with different initial crack lengths a0 under fatigue loading. The numerical results are compared to the theoretical predictions of the damage model (plate with a0=0mm) and the fatigue crack growth model (plates with a0>0mm). 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., a0=10mm, 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 stiffness decrease. However, a smaller crack (a0=2mm) deviates from fatigue crack growth because of the additional contribution of the material damage to the stiffness degradation. The 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 affect the sensitivity of the material to cracks.Figure 18
Stiffness degradation as a function of the number of cyclesNc for different initial crack sizes.


### 5.2.4. Damage Field

Figure19(a) illustrates the damage distribution of a precracked plate with an initial crack length a0=5mm after 166 fatigue loading cycles. The red values correspond to D=1 and indicate the propagation of the crack. The 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 effects. The 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. The 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.Figure 19
Damage distribution and crack propagation of a precracked plate with the initial crack lengtha0=5mm after (a) 166, (b) 179, (c) 187, and (d) 215 fatigue loading cycles.
(a)
(b)
(c)
(d)The 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 Figure19. The size of the fracture process zone seems to depend on the crack size. The reduction of the global stiffness of the structure (Figure 18) is indeed controlled by the propagation of the crack associated with the damage of the uncracked zone, as shown in Figure 19.


## 5.2.1. Geometry, Loading, and Material Properties

Figure16 shows a center cracked plate with dimensions w=80mm and L=120mm, which is tested under sinusoidal fatigue loading with an amplitude σ0=1.25MPa and a frequency f=25Hz. The material presents Young’s modulus E0=30GPa and Poisson’s ratio ν=1/3. A viscous damping cn=ct=0.0295Ns/m and a time step Δt=5×10−7s are adopted in the simulations. The damage model parameters are α=−2.25, β=4.0, and γ=5.0×1016; the fatigue crack growth model parameters are c=3.0×10−10m/cyclePamm and m=1.03. Negligible effects of the frequency indicate very low dynamic effects (quasi-static regime) on the presented results.Figure 16
Center cracked plate subjected to imposed sinusoidal stress.


## 5.2.2. Evaluation of the Stress Intensity RangeΔK

In Figure17, the numerical evaluation of the stress intensity range ΔK is compared to the theoretical result for different particle diameters d. When the crack length a is below 25mm, 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>25mm). These deviations tend to decrease for smaller values of d, which indicates a convergence of the data.Figure 17
Comparison of the range of stress intensity factorΔK calculated using numerical and theoretical methods for different crack lengths a (a0=5mm) and particle diameters d.The effect of the particle sized tends to be neglectful of the definition of the stress intensity range ΔK (Figure 17) and also of the crack growth rate da/dNc (Equation (15)). However, d defines the propagation step of an existing crack and, according to Equation (16), affects proportionally the number of cycles Nc. It means that the choice of smaller d induces only more precision on the numerical calculation.


## 5.2.3. Stiffness 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. This behaviour is due to the decrease of the global stiffness of the sample caused by damage of the material and the propagation of the crack. The process of the stiffness degradation can be quantified 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 cyclesNc.Figure18 shows the stiffness degradation process obtained for plates with different initial crack lengths a0 under fatigue loading. The numerical results are compared to the theoretical predictions of the damage model (plate with a0=0mm) and the fatigue crack growth model (plates with a0>0mm). 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., a0=10mm, 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 stiffness decrease. However, a smaller crack (a0=2mm) deviates from fatigue crack growth because of the additional contribution of the material damage to the stiffness degradation. The 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 affect the sensitivity of the material to cracks.Figure 18
Stiffness degradation as a function of the number of cyclesNc for different initial crack sizes.


## 5.2.4. Damage Field

Figure19(a) illustrates the damage distribution of a precracked plate with an initial crack length a0=5mm after 166 fatigue loading cycles. The red values correspond to D=1 and indicate the propagation of the crack. The 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 effects. The 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. The 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.Figure 19
Damage distribution and crack propagation of a precracked plate with the initial crack lengtha0=5mm after (a) 166, (b) 179, (c) 187, and (d) 215 fatigue loading cycles.
(a)
(b)
(c)
(d)The 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 Figure19. The size of the fracture process zone seems to depend on the crack size. The reduction of the global stiffness of the structure (Figure 18) is indeed controlled by the propagation of the crack associated with the damage of the uncracked zone, as shown in Figure 19.


## 6. Conclusions

The damage and fatigue crack growth are normally modelled separately by different approaches. The damage is usually modelled by nonlocal approaches implemented in a finite element code, which may produce reasonable global sample behaviour, based on unrealistic material behaviour. The 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 different approaches, a simple numerical scheme coupling damage and fracture mechanics in a discrete element environment was proposed. The association of these different 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. This 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. The evaluation method of the stress intensity range during the damage evolution provides a way to couple both failure mechanisms. This 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. The 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. The proposed model can be further extended to study the crack nucleation and short crack growth problems for varied microstructures. Such a numerical model can be constructed by representing different materials (aggregates, mortar, etc.) as clusters of discrete elements with different parameters, including the mechanical properties of the contacts and model parameters for the continuum damage model and fatigue crack growth model.The practical applications and precisely comparison with the experimental results will be performed in the future study, and the implementations of the aggregate or grain size in the analytical equations and discrete element code would make this model more realistic and robust.



---

*Source: 6962394-2019-02-26.xml*