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.

1. Introduction

The size and boundary effects of quasi-brittle fractures have been studied systematically in the past decades [19], and some significant progresses have been achieved in recent years by taking the aggregate or grain size into account in the analytical solutions [1014]. 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 [1517]. 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 [1921] 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 [2326] 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 [2931] and also improved to study the damage initiation and evolution in the interconnections [15, 3234].

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 [3638] 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 Section 5, 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 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. The stress-strain relation can be written as follows:where and are the elastic stress and strain tensor components and and are the damaged and initial (elastic) secant stiffness of the material. 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:

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:where is the increment of local equivalent strain from to and is a function of damage. The exponent β is a material parameter, which relates to the slope () of the S-N curve in the log-log diagram.

The function of damage 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]:

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 [5052] is described well by the well-known Paris (or Paris–Ergodan) law [49], which relates the stress intensity factor range to the crack growth rate via a power law, with . The basic formula is as follows:where is the crack growth rate, a is the crack length, and is the number of load cycles; c and m are the material parameters dependent on the material composition, environment, stress ratio, etc. [53]. and are the fatigue threshold and fracture toughness of the material, respectively. Paris’ law works for sufficiently large cracks, where is higher than the fatigue threshold , and the maximum value of the stress intensity factor remains below the material fracture toughness .

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 Figure 1(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. [5559] and Liu et al. [60] for quasi-brittle failure of asphalt concrete, rocks, and crystals.

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, and , are a sufficiently small fraction of (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, and (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, and , 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].

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 ( and ) 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 and (, and , , respectively, as shown in Figure 3(a)) give rise to mean strain values, and :

The normal and tangential components of each contact (, and , , 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 (), mean stresses (Figure 3(c)) may be associated with these forces:

3.2.2. Principal Stresses

The stress tensor (in two dimensions) can be defined by the values of the principal stresses and 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:

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:where .

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 modulus E and Poisson’s ratio ν) and discrete elastic parameters (normal and tangential stiffness, and , 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 (Section 3.2): the contact forces N and T and contact displacements and 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 and (in two dimensions) can be defined by the value of the principal stresses and 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 defined as a function of the local equivalent strain rate (Equation (3)). The damage functions 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 () is updated, and finally, the values of the normal and tangential stiffness are updated (Equation (6)).

4.2. Model Verification

An uncracked plate with dimensions and subjected to sinusoidal fatigue loading with the amplitude (Figure 5), the frequency , and the corresponding angular frequency is considered. The material presents Young’s modulus and Poisson’s ratio . The damage function (Equation (4)) proposed by Paas [48] is adopted in this analysis. The model parameters are , , and . The higher value of parameter C is adopted in order to ensure the fatigue life of applied loading condition is around 100 cycles ( [62]), avoiding the jump-in-cycle procedure in this calculation.

The discrete element model is shown in Figure 6. The particle radius is , 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 and . A time step and a low viscous damping 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 is applied at the upper boundary. In order to apply the strain rate during all cycles, the amplitude of the displacement of the upper boundary is constant and equal to .

The 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 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 is identified.

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-6700K processor with 4 cores and 8 threads running on a Windows 8.1 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 (without jump-in-cycle procedure). More than 72 hours and around 12 hours were needed, respectively, for and . 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) (Figure 9). 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.

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 (Section 2.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 (Section 2.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 (Figure 10).

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 cycle , 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 in Figure 11 is located at the crack tip.

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 , 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

A damage value 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 reaches 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 in contact , 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.

The principal components of the contact forces N and T and contact displacements and 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: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, and 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 and naturally vary from 0 to and , respectively. In this case, the damage of a contact induces a maximum energy release rate [64] . The minimum energy release rate is equal to zero; consequently, the variation of the energy release rate is simply defined aswhere is the number of cycles to reach (total release of the contact energy) and is the surface of the triangle formed by the points , , and , as shown in Figure 14.

The second contact is not identified as a crack tip. The degradation of the elastic properties of the contact induces an increase of the contact force in , as shown in Figure 15. After contact is totally released, in the case of a crack propagation, contact becomes the new crack tip. The value of the damage of is set back to zero , and the value of the number of cycles starts to be incremented. Once for contact , 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 defined as .

5.1.3. Crack Initiation and Propagation

The evolution of the damage variable D characterizes the weakening of a contact before rupture. The rupture associated with the propagation of a crack is defined by the value of (Equation (5)):

A contact which presents is definitely broken (and ceases to exist) only ifwhich represents a crack growth equals to for one contact break. Otherwise, D is set to 0, and the fatigue loading continues, until the number of cycles 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

Figure 16 shows a center cracked plate with dimensions and , which is tested under sinusoidal fatigue loading with an amplitude and a frequency . The material presents Young’s modulus and Poisson’s ratio . A viscous damping and a time step are adopted in the simulations. The damage model parameters are , , and ; the fatigue crack growth model parameters are and . Negligible effects of the frequency indicate very low dynamic effects (quasi-static regime) on the presented results.

5.2.2. Evaluation of the Stress Intensity Range

In Figure 17, the numerical evaluation of the stress intensity range is compared to the theoretical result for different particle diameters d. When the crack length a is below , the discretization levels that we have chosen can all provide an acceptable stress intensity range . Some deviations can be observed when the crack becomes larger (). These deviations tend to decrease for smaller values of d, which indicates a convergence of the data.

The effect of the particle size d tends to be neglectful of the definition of the stress intensity range (Figure 17) and also of the crack growth rate (Equation (15)). However, d defines the propagation step of an existing crack and, according to Equation (16), affects proportionally the number of cycles . 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 cycles .

Figure 18 shows the stiffness degradation process obtained for plates with different initial crack lengths under fatigue loading. The numerical results are compared to the theoretical predictions of the damage model (plate with ) and the fatigue crack growth model (plates with ). 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., , 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 () 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.

5.2.4. Damage Field

Figure 19(a) illustrates the damage distribution of a precracked plate with an initial crack length after 166 fatigue loading cycles. The red values correspond to 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 () 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.

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 Figure 19. 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.

Data Availability

The datasets generated and/or analysed during the current study are available from the corresponding author on reasonable request.

Conflicts of Interest

The authors declare no conflicts of interest.

Acknowledgments

This manuscript is based on Xiaofeng Gao’s doctoral thesis. This research was supported by the program of the China Scholarship Council and China Postdoctoral Science Foundation (Grant no. 2018M631478).