Accepted Manuscript A discrete element modelling approach for fatigue damage growth in cemented materials Nhu H.T. Nguyen, Ha H. Bui, J. Kodikara, S. Arooran, F. Darve PII: S0749-6419(18)30246-8 DOI: 10.1016/j.ijplas.2018.08.007 Reference: INTPLA 2397 To appear in: International Journal of Plasticity Received Date: 20 April 2018 Revised Date: 14 August 2018 Accepted Date: 14 August 2018 Please cite this article as: Nguyen, N.H.T., Bui, H.H., Kodikara, J., Arooran, S., Darve, F., A discrete element modelling approach for fatigue damage growth in cemented materials, International Journal of Plasticity (2018), doi: 10.1016/j.ijplas.2018.08.007. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. ACCEPTED MANUSCRIPT A discrete element modelling approach for fatigue damage growth in cemented materials RI PT AUTHORS: AUTHOR’S AFFILIATIONS M AN U SC Nhu H.T. Nguyen1, Ha H. Bui1*, J. Kodikara 1, S. Arooran1 and F. Darve2 Department of Civil Engineering, Monash University, Melbourne, Australia 2 Grenoble CNRS, UMR 5521, 3SR, Grenoble Alpes University, Grenoble, France EP TE D 1 AC C *CORRESPONDING AUTHOR Dr. Ha H. Bui, Senior Lecturer Department of Civil Engineering, Monash University Room 126, Building 60, Monash University, Victoria 3800, Australia Email: ha.bui@monash.edu Tel : +61 3 9905 2599 1 ACCEPTED MANUSCRIPT Abstract This study proposes a modelling approach capable of capturing both the fatigue response and localisation of failure in cemented materials. This modelling approach takes advantage of Discrete Element Method (DEM) to reproduce the heterogeneous microstructure and crack development in cemented materials. In conjunction with this, a new constitutive model is RI PT developed to characterise the fatigue behaviour of the materials at the grain scale. The model formulation is based on coupling damage mechanics and plasticity theory and combining with a fatigue damage evolution law to describe the degrading response of cemented materials subjected to cyclic loading. The proposed model is employed to govern the explicit SC behaviours of DEM bonding contacts representing cement bridges between aggregates in the physical materials. The macro-behaviour is then obtained in DEM simulations as the M AN U collective response of all contacts and particles in the material domain. Through numerical experiments, the proposed modelling approach is shown to capture well the fatigue behaviour and cracking process in cemented materials subjected to cyclic loading. The microstructural effects on the fatigue response of the materials are naturally reproduced in simulations thanks to the discrete nature of DEM. These results demonstrate the capability of the proposed modelling approach as well as its potential to be a faithfully numerical technique for TE D modelling and investigating fatigue damage in cemented materials. Keywords: Discrete element method, fatigue damage, cemented materials, couple damage AC C EP plasticity, fatigue model 2 ACCEPTED MANUSCRIPT Introduction Fatigue cracking due to repetitive loads has been recognised as one of the main distress modes in many structures made by cemented materials such as highway pavements, airport pavements and concrete foundations of wind power plans (Dattoma et al., 2006; Fathima & Kishen, 2013). When these structures are subjected to cyclic loading, plastic deformation is RI PT observed to develop even though the applied nominal stress is below the material yielding limit. If the applied nominal stress is sufficiently small, plastic deformation gradually stabilises and eventually unchanged after a number of load cycles. This behaviour is called shakedown (Weichert & Maier, 2014; Zhang et al., 2017), and the magnitude of cyclic loads SC under which the plastic deformation can stabilise is the shakedown limit or endurance limit of the materials. Nevertheless, if the applied stress is greater than the shakedown limit, plastic damage of cemented structures. M AN U deformation keeps accelerating with the increase of cyclic loading, leading to the fatigue At the microstructural scale, fatigue damage in cemented materials is considered as a progressive process of microstructural changes with increasing the number of load cycles. This process is associated with the propagation and coalescence of micro-cracks in the materials (Liang & Zhou, 1997; Fathima & Kishen, 2013). The development of micro-cracks TE D in cemented materials subjected to cyclic loading is mainly attributed to their heterogeneous microstructure, which consists of aggregates and cement mortar with pre-existing microcracks at the mortar matrix or at the interfaces between mortar and aggregates (Zhang, 2001; Ray & Kishen, 2010). This heterogeneous feature of cemented materials causes stress to EP concentrate in areas near micro-crack tips, making local stresses in these areas to be higher than the applied nominal stress. As a consequence, the micro-cracks can grow although the applied nominal stress has not reached the material yielding limit. In this context, the fatigue AC C 1. behaviour of cemented materials can be considered to depend on the rate of micro-cracking per load cycle and material microstructural features. Modelling of fatigue damage in cemented materials thus requires appropriate descriptions of the fatigue crack growth rate and microstructural effects. Existing modelling approach for fatigue damage in cemented materials can be classified into two main methods, based on fracture mechanics or damage mechanics. The fracture mechanics-based method describes the fatigue behaviour of cemented materials as the relationship between the crack growth rate per load cycle and the fracture mechanics terms such as stress intensity factor and fracture energy (Ray & Kishen, 2010; Simon & Kishen, 3 ACCEPTED MANUSCRIPT 2017), or based on the cohesive fracture theory (Nguyen et al., 2001; Yang et al., 2001). On the other hand, the damage mechanics-based method focuses on developing the evolution law of damage variable with respect to load cycles (Papa & Taliercio, 1996; Alliche, 2004). The evolution of damage variable, in this sense, represents the loss of material load-carrying capacity when undergoing cyclic loading. As the fatigue damage of cemented materials is RI PT considered to be an irreversible process (Horii et al., 1992; Fathima & Kishen, 2013), the theory of plasticity has also been adopted to describe the increment of plastic deformation in the material fatigue response (Voyiadjis & Abu-Lebdeh, 1994; Ristinmaa, 1995; Moreo et al., 2007; Al-Rub & Darabi, 2012; Darabi et al., 2012; Carrara & De Lorenzis, 2015; Cheng et al., SC 2016). Although a great effort has been dedicated to developing constitutive models for the fatigue damage of cemented materials, current modelling approaches still face challenges in capturing both the fatigue response and localised failure in the materials subjected to cyclic M AN U loading. This can be attributed to the fact that modelling fatigue damage is more complicated than other types of damage (Xiao et al., 1998), as it requires the appropriate characterisations of different stages of crack growth (i.e. stable crack growth and unstable crack growth (Gao & Hsu, 1998)) as well as the material heterogeneous microstructure. The microstructural effects on the macro-behaviour of cemented materials can be addressed by micromechanics- TE D based models (Krajcinovic & Fanella, 1986; Mohamed & Hansen, 1999; Shen et al., 2012; Shen & Shao, 2016). Besides, most existing fatigue models have been developed for continuum modelling and implemented into continuum computational methods such as finite element method (FEM) for numerical simulations. The continuum methods, nonetheless, EP suffer from over-complicated model formulation and ad-hoc treatments to deal with the heterogeneous response and localised failure occurred in the fatigue damage process of AC C cemented materials. Therefore, it is necessary to have an alternative approach that can effectively handle the heterogeneous response and localised failure in cemented materials, hence easing the complexity of constitutive model formulation. In recent decades, Discrete Element Method (DEM) (Cundall & Strack, 1979) has emerged as an effective numerical tool to model the damage behaviour of cemented materials and reproduce the material microstructural responses. Different from continuum methods, which consider materials as a single continuous body, DEM discretises material domains as an assemblage of particles, carrying their own physical properties and motions. These particles interact with each other through contacts and these interactions are governed by contact models. This discrete nature enables DEM to naturally replicate the heterogeneous 4 ACCEPTED MANUSCRIPT microstructure of cemented materials. Furthermore, as DEM allows large displacement and full detachment between particles, this method can handle the localised failure and discontinuity in the materials with no need for ad-hoc treatments on constitutive contact models. Such advanced features make DEM a potential method to model the fatigue damage of cemented materials. In fact, several bonding contact models have been developed to RI PT characterise the damage behaviour of cemented materials under monotonic loadings (Scholtès & Donzé, 2013; Gui et al., 2015; Gui et al. 2016; Oñate et al., 2015; Nguyen et al., 2017a; Nguyen et al., 2017c). As the fatigue damage is due to the micro-cracking induced from the heterogeneous microstructure of the materials, these bonding contact models might be used to SC characterise fatigue damage in cemented materials if the DEM particles and contacts are modelled at the molecular scale and the true microstructure of cemented materials is accurately replicated. Nonetheless, this modelling approach would remarkably increase the M AN U overall computational cost as millions of particles will be generated in simulations. For more computational efficiency, particles in DEM simulations are generally modelled at larger scales (e.g. grain scale or aggregate scale). At these scales, the above bonding contact models are unable to describe the fatigue damage growth in cemented materials owing to the lack of intrinsic fatigue damage mechanisms in the cement bridges. As a result, it is necessary to TE D have bonding contact models for DEM incorporating the grain-scale fatigue damage mechanisms to capture the fatigue responses of cemented materials. However, to the authors’ knowledge, there is no existing contact model for DEM capable of capturing the fatigue behaviour of cemented materials. This indeed prevents the application of DEM for modelling EP fatigue damage in the materials. The aim of this paper is to propose a modelling approach based on DEM to characterise AC C fatigue damage in cemented materials. In this approach, DEM plays the role in modelling the heterogeneous response and localised failure in the materials, while a new fatigue bonding model is developed to incorporate physical fatigue damage mechanisms at cement bridges in the materials. The unique feature of this approach is its capability in representing the microstructural processes of cemented materials and easily capturing crack initiation and propagation. The new fatigue model remains simple as it is formulated in two principal directions (i.e. normal and shear) of the DEM contact and only characterises the grain-scale fatigue behaviour of cement bridges, while different stages of fatigue crack growth exhibited at the macroscopic scale of cemented materials is captured as the collective response of all particles and contacts in the material domain. The capability of this proposed modelling 5 ACCEPTED MANUSCRIPT approach in capturing and predicting the fatigue behaviour of cemented materials is first compared with the experimental data. Subsequently, its ability in reproducing the microstructural effect on the fatigue damage response of cemented materials is demonstrated in further simulations. This article is organised as follows. The basis of DEM computation is briefly explained in RI PT Section 2. Details of the model formulation and implementation algorithm are then presented in Section 3. This is followed by the demonstrations of model behaviour at the constitutive level and parametric studies in Section 4 in order to illustrate the influences of the model parameters on contact behaviour at the constitutive level. Next, Section 5 presents the SC application of the proposed modelling approach to capture the fatigue damage responses and microstructural behaviour of cemented materials. Finally, Section 6 summarises conclusions Numerical background Particle motion equations The DEM is adopted in this study as the numerical framework to simulate the behaviour of cemented materials. In DEM, the computational domain is discretised into a number of rigid second law as follows: − = where + EP − = TE D particles in different sizes. These particles can move and rotate in accordance with Newton’s is the translational acceleration; AC C 2. M AN U drawn from this study. is the moment of inertia of the particle. (1) (2) is the rotational acceleration; and is the mass and are the total force and moment acting on the particles as illustrated in Fig. 1 and are computed as the vectorial summation of all applied and contact forces as follows: = = where + (3) × + is the applied force; (4) is the contact force, which is obtained based on the force- displacement relationship determined in contact models; particle centre to the contact point and is the vector connecting the is the applied moment. On the other hand, 6 and ACCEPTED MANUSCRIPT are the damping force and moment, incorporated into Eqs. (1) and (2) to damp out the particle accelerations. This is essential to setup a quasi-static condition in simulations as the DEM formulations are fully dynamics. In this study, the damping force and moment are computed using the global damping formulations proposed in Potyondy & Cundall, (2004). The global damping force and moment are proportional to the total force and moment acting on the particle, but imposed in opposite directions as follows: + = RI PT =α (5) (6) where α is the global damping coefficient. To achieve a quasi-static condition in DEM SC simulations, Potyondy & Cundall, (2004) suggested using α = 0.7 so that particle’s M AN U acceleration is heavily damped. In addition, the loading rate in simulations should be chosen small enough to ensure numerical results at the structural level are not affected by changes in loading rate. In this sense, the effect of kinetic energy on the macro-behaviour in DEM simulations can be negligible. Subsequently, the motion equations of particles are integrated in time. In DEM computation, the standard Leap-Frog algorithm is normally used to integrate the particle motion equations. TE D In particular, once particle accelerations at a time step “ ” are known, the particle velocities and displacements at the next time step are updated as follows: !∆ /$ !∆ = = %∆ /$ + + %∆ /$ + ∆ ∆ !∆ /$ =& + ∆ !∆ /$ (7) (8) (9) ∆ AC C & !∆ = EP !∆ /$ (10) As the DEM computation uses the explicit time integration scheme, it is necessary to impose a constraint for the computational time step as below to ensure numerical stability: ∆ ≤∆ () in which ∆ ∆ () = C, (11) () is the critical time step determined as: -. ⁄/ (12) 7 where C, ACCEPTED MANUSCRIPT - ≤ 1 is the stability coefficient and / is the stiffness parameter. In the modelling of samples with numerous particles, the critical time step is normally attained as the minimum of all critical time steps computed on all particles in the assembly. DEM modelling of intact materials DEM was conventionally developed to model the non-cohesive behaviour of granular media. RI PT This method has later been applied to simulate the behaviour of intact materials with the introduction of a bonding contact model (Potyondy & Cundall, 2004). In this model, particles are allowed to be bonded when they are not touching each other, as long as the distance between their centres is smaller than a distance tolerance. The bonded particles can interact even when they are non-overlapped and the bonding contact can resist all tension, shear and links the contact force to contact stress as follows: M AN U = 3 2 = 34 2 + 3, 2 SC compression stresses. The bonding contact is assigned a specific shape and area 2 , which (13) where 34 and 3, are the normal and shear stress components acting in the bonding contact (see Fig. 1), respectively. The normal and shear stresses in the bonding contact are induced from the relative displacements of the bonded particles based on the bonding contact models. TE D In DEM modelling of intact materials, the bonding contact model plays an important role in characterising the local behaviour of materials while the macro-behaviour of materials is AC C EP captured as the combination of all contact and particle responses in numerical samples. Fig. 1: The total force and moment acting on particles and the normal and shear stress vectors acting on the bonding contact 3. A fatigue damage contact model for DEM 8 ACCEPTED MANUSCRIPT As mentioned earlier, DEM is able to reproduce the heterogeneity of cemented materials. However, the uses of circular particles to represent aggregates and contact bonds to represent mortar matrix in the general DEM simulation have significantly simplified the true microstructural details of the materials (e.g. morphologies of aggregates and cohesive mortar matrix). Furthermore, existing bonding contact models in DEM are unable to capture the RI PT fatigue damage growth in cemented materials owing to the lack of fatigue crack growth mechanisms at the grain scale. These issues have indeed limited the applicability of DEM for modelling fatigue cracking in cemented materials. To overcome these limitations, a new fatigue bonding contact model for cemented materials is developed in this study. This SC bonding contact model phenomenologically takes into account those microstructural details of cemented materials missing in DEM simulations and more importantly captures the fatigue miro-cracking mechanism of the materials. The proposed model is formulated on the basis of M AN U the cohesive damage-plastic model developed by the authors to model the grain-scale damage behaviour of cemented materials under monotonic loads (Nguyen et al., 2017b). In the proposed model, new fatigue variables and fatigue evolution law are incorporated into the original cohesive damage-plastic law to account for the micro-crack growth and irreversible displacement occurred in the cement bridge when subjected to cyclic loading. Details of the 3.1. Model formulation TE D fatigue bonding contact model formulation are described as below. The expression of the relative displacement between two particles connected by a bonding contact in the original model is first redefined to incorporate the fatigue displacement EP variable. This variable represents the irreversible displacement at cement bridges caused by fatigue cracking, which distinguishes itself from the plastic displacement developed when AC C stress reaches the yielding limit. The separation of fatigue and plastic displacements is necessary to reflect fundamental differences between the underlying mechanisms of these two irreversible processes. Specifically, fatigue occurs due to the degradation of elastic modulus produced by micro-cracking. The micro-cracking is an irreversible process associated with the local micro-mechanisms of failure. Micro-cracks can grow even when the macro-behaviour of the materials is quasi-static and cyclic stress state is still below the plastic yield surface of the materials. On the other hand, once micro-cracks grow into macro-cracks, the plastic displacement occurs as the result of frictional interaction between two sides of the macro-cracks. At this stage, the stress state of the material lies on its plastic yield surface. As 9 ACCEPTED MANUSCRIPT a consequence, the relative displacement in a bonding contact is now equal to the sum of elastic, plastic and fatigue displacement components as follows: = 5 + 6 + 7 (14) Accordingly, the contact stress-displacement relationships in normal and shear directions are 84 = /49 1 − : ;<4 − <46 − <4 = − :/49 〈−;<4 − <46 − <4 =〉 7 7 8, = /,9 1 − : ;<, − <, − <, = 7 6 RI PT rewritten as: (15) (16) SC where /49 and /,9 are the initial normal and shear stiffness of the bonding contact; : is the damage variable. The Macauley brackets embraced in Eq. (15) is to impose the full recovery of contact normal stiffness in compression modes. M AN U The expression of damage variable : is also revised to take into account the contribution of fatigue damage. In this model, it is considered that within a load cycle (i.e. loading and then unloading), if the stress state in bonding contact is below its yielding limit, only fatigue damage can develop in the contact. On the other hand, if stress in the bonding contact reaches the yielding limit, fatigue damage is lumped into the plastic damage due to yielding. In this sense, the damage variable is expressed as the sum of damage increments in all load cycles, := A )BC @: = A )BC D 1 − E) @:6 + E) @:7 F (17) EP (@:7 ): TE D each of which consists of an incremental damage due to either yielding (@:6 ) or fatigue AC C where E) is a shift factor in each load cycle, which equals 0 if the stress state in the contact reaches its yield surface and 1 if the stress state in the contact is below its yield surface within a load cycle. It can be seen from Eqs. (15), (16) and (17), the normal and shear stresses in the bonding contact are influenced by two sets of internal variables, including plastic variables (<46 , <,6 and :6 ,) and fatigue variables (<4 , <, and :7 ). The plastic variables represent the 7 7 plastic damage state of the bonding contact, while the fatigue variables manifest its fatigue damage state. The yielding and fatigue processes in the bonding contact are thus characterised by the increments of these internal variables. 6 6 Increments of plastic variables: @<4 , @<, and @:6 10 ACCEPTED MANUSCRIPT The increments of plastic variables are calculated based on the original damage-plasticity model. Details of the original model formulations can be found in Nguyen et al., (2017b). Herein, only the key formulations and modifications made to this model for the calculations of incremental plastic displacements and damage variable are presented. Based on the original damage-plasticity model, the plastic displacements and damage RI PT variable can develop when stress in the bonding contact reaches a yield criterion which takes the following form: FH = I $ 8,$ − JK 9 1 − : − 84 L MN$ + JK 9 1 − : − 8 9 1 − : L MN$ $ SC .2K 9 8 9 L M − 8 9 L M where I = K9 (18) In the above equations, 8 9 is the contact tensile strength, K 9 is the contact cohesion and M is M AN U the friction angle. When the bonding contact is yielding, the increments of normal and shear plastic displacements are updated based on the standard plastic flow rule as follows: @<46 = 2@λ L ψJK 9 1 − : − 84 L ψN @<, = 2@λ8, 6 (19) (20) TE D where ψ is the dilatancy angle and @λ is the plastic multiplier. The evolution of plastic damage variable is then derived from the formulation of damage variable evolution proposed in the original model, yet taking into account the reformulation of total displacement and damage variables in Eqs. (14) and (17). This results in the following function of incremental V V ST ST %R WU ! WX Y TU TX AC C @:6 = χ P1 − Q EP plastic damage variable: Z (21) In Eq. (21), χ is the term associated with the deformation state of the bonding contact and is expressed as follows: χ=Q V ] V ] ∑ STU !∑ STU ∑ STX !∑ STX ^ %[ ! W W TU TX (22) where <4 and <, are the softening parameters as in the original model. Eqs. (21) and (22) can facilitate coupling between damage mechanics and plasticity theory to fully describe both cracking and irreversible displacement occurred in the damage process of cemented materials. In addition, the power law obtained for the incremental plastic damage variable can provide a 11 ACCEPTED MANUSCRIPT smooth transition between the cohesive and non-cohesive behaviour of a contact bond when damage variable is close to one (Nguyen et al., 2017b). Increments of fatigue variables: @<4 , @<, and @:7 7 7 RI PT Different from the plastic variables, the fatigue variables are assumed to develop when the maximum stress in the bonding contact within a load cycle satisfies a fatigue criterion, i.e. _7 ` 0. In this study, the fatigue criterion is derived from the yield criterion by incorporating a fatigue threshold parameter, a, as follows: SC Fb = I $ 8,$ − JaK 9 1 − : − 84 L MN$ + JaK 9 1 − : − a8 9 1 − : L MN$ (23) From Eq. (23), it can be seen that a determines the proportion of the fatigue limit to the M AN U plastic yielding limit of the bonding contact. The value of a has to be chosen between 0 and 1 to ensure the fatigue criterion is always below the yield criterion, hence satisfying the predetermined condition that fatigue damage is able to evolve under the yield criterion. Different from the yield criterion, the fatigue criterion only plays the role in determining the stress states in which fatigue damage can evolve in bonding contacts, while the increment of fatigue TE D damage variable is defined independently later on. The shape of the fatigue criterion is displayed in Fig. 2 along with the yield criterion. The shaded area represents the range of AC C EP maximum cyclic stress state in the bonding contact wherein fatigue damage can develop. Fig. 2: The yield and fatigue criteria of the proposed model. The shaded area represents the range of maximum cyclic stress states wherein fatigue damage is allowed to develop 12 ACCEPTED MANUSCRIPT Once the above fatigue conditions are met, the fatigue damage variable will increase. Different from the plastic damage variables whose increments are calculated cycle-by-cycle, the increment of the fatigue damage variable is only updated once at the onset of unloading path in a load cycle. This updating method is considered to be more computationally efficient than the cycle-by-cycle update, especially in the modelling of high-cycle fatigue damage RI PT (Turon et al., 2007). Moreover, it sets a precedence for embedding the “cycle jump” technique to reduce computational time in the simulations of cyclic tests. As a consequence of this fatigue damage updating method, the incremental fatigue damage variable is formulated as the function of the number of load cycles in accordance with Paris’ law (Paris SC & Erdogan, 1963). Originally, Paris’s law defined the fatigue crack growth rate per load cycle as a power function of the stress intensity factor as follows: e M AN U @L7 = K ∆d @c where L7 is the fatigue crack length; c is the number of load cycles; K and (24) are the material fatigue parameters; and ∆d is the stress intensity factor range in a load cycle. The fatigue crack length L7 is then linked to the fatigue damage variable based on the geometrical := L f TE D expression of the damage variable for rectangular DEM contacts as follows (see Fig. 3): (25) (25) leads to (26) AC C @L7 = f@:7 EP where L is the total crack length. In case the crack growth is due to fatigue, i.e. @L = @L7 , Eq. Fig. 3: The shape of a bonding contact and presumed crack length By substituting Eq. (26) into Eq. (24), the incremental fatigue damage variable is obtained as: 13 @:7 K ∆d = @c f ACCEPTED MANUSCRIPT e (27) In Eq. (27), the stress intensity factor range, ∆d, represents the effect of cyclic stress state on the fatigue behaviour of cement bonds. According to this equation, ∆d = ∆dg when the bonding contact is under mode I loading, while ∆d = ∆dgg when the bonding contact is under RI PT mode II loading. The contributions of both modes I and II are also taken into account in a mixed-mode loading based on an energy release rate analysis (Paris & Sih, 1965). As a result, the stress intensity factor is considered to take the following form: ∆d = . ∆dg $ + ∆dgg $ (28) SC In Eq. (28), mode I and mode II stress intensity factors are calculated based on the fracture mechanics theory (Irwin, 1957): k ∆dgg = hgg √jL∆8, = hgg √jL;8,e − 84e)4 = k M AN U ∆dg = hg √jL∆84 = hg √jL;84e − 8,e)4 = (29) (30) where hg and hgg are the geometrical factors calculated by Eq. (31). This equation is derived by Tada et al., (1973) based on the analysis results of mode I and mode II tests on rectangular plates with a horizontal centre notch, which is similar to the shape of the DEM contacts. TE D L $ L p 1 − 0.025 m n + 0.06 m n f f hg = hgg = h = jL qrs m n 2f (31) EP Subsequently, substituting Eq. (25) into (29), (30) and (31) yields: ∆dg = h√jf:;84e h= k − 84e)4 = (32) − 8,e)4 = AC C ∆dgg = h√jf:;8,e k (33) 1 − 0.025:$ + 0.06:p qrs mj: n 2 (34) In Eq. (32) and (33), the value of damage variable has to be greater than 0 so that Paris’ law (i.e. Eq. (27)) can be valid. To satisfy this condition, the damage variable is assumed to take a small initial value :)4) in numerical simulations. This assumption is reasonable in the sense that :)4) represents the pre-existing micro-cracks in cement bonds after the sample generation process. 14 ACCEPTED MANUSCRIPT Along with the evolution of the fatigue damage variable, the fatigue displacements in normal and shear directions increase. The incremental fatigue displacements are related to the incremental fatigue damage variable by Eq. (35). This equation is derived in the same way as the formulation of incremental plastic damage variable in Eq. (21): ] Z (35) RI PT @:7 = χ P1 − Q ] ST ST %[ WU ! WX ^ TU TX This leads to @:7 @<4 @<, + = −t R1 − Y <4 <, χ 7 (36) SC 7 As Eq. (36) is insufficient to solve for the normal and shear incremental fatigue M AN U displacements in mixed-mode conditions, an additional constraint between these two variables is required. In this work, the additional constraint is defined as the ratio between @<4 and @<, . This constraint represents the proportion of incremental fatigue displacements 7 7 in normal and shear directions in a load cycle and can be derived following the concept of plastic flow rule used to calculate the normal and shear plastic displacements in Eqs. (19) and (20) as follows: 7 @<, 7 = L ψJaK 9 1 − : − 84 L ψN 8, TE D @<4 (37) Eqs. (36) and (37) make a set of equations with two unknowns, @<4 and @<, . By solving this 7 7 EP set of equations, the values of @<4 and @<, can be obtained. Accordingly, the plastic and 7 7 fatigue damage variables are calculated and the total damage variable is updated using Eq. AC C (17). During the failure process, the total damage variable increases from :)4) to 1. In parallel to this, both the yield and fatigue criterions are allowed to shrink in order to represent the loss of load-carrying capacity of the bonding contact (see Fig. 4). The yield criterion shrinks faster than the fatigue criterion, and these two criteria meet when the damage variable reaches one, which corresponds to the final damage stage of the bonding contact. 15 RI PT ACCEPTED MANUSCRIPT 3.2. The model implementation M AN U damage variable SC Fig. 4: The evolutions of the yield and fatigue criteria with respect to the increment of For DEM numerical experiments, the proposed constitutive model is implemented into PFC 5.0 (Itasca, 2014). The implementation algorithm of the proposed fatigue model is illustrated in Fig. 5. The computational process involves updating the internal variables and updating stresses. Once all internal variables are updated, the stresses in the bonding contact are then TE D rectified from the stress predictors, 846(5 and 8,6(5 . Specifically, if the bonding contact is yielding, the stresses are updated using the stress return algorithm given in the original model (Nguyen et al., 2017b). On the other hand, if the contact is subjected to fatigue damage, the EP increments of the fatigue displacements and fatigue damage variable are updated at the onset of unloading paths. The stresses at this position are then rectified using Eq. (38) and (39) to account for the change of contact state caused by the increments of fatigue displacements and AC C fatigue damage variable. 84A = 846(5 8,A = 8,6(5 1 − :A 7 − @<4 1 − :A /49 1 − :A%C (38) 1 − :A 7 − @<, 1 − :A /,9 A%C 1−: (39) where c and c − 1 denotes the current and previous number of load cycles, respectively. 16 EP TE D M AN U SC RI PT ACCEPTED MANUSCRIPT AC C Fig. 5: Flow chart of the implementation algorithm of the proposed fatigue model 3.3. Cycle jump technique The high-cycle fatigue damage simulations of cemented materials involve thousands of load cycles, which requires extensive computational time to model cycle by cycle the entire fatigue damage process. In the simulations of regular cyclic loading scenarios with the same loading range, a “cycle-jump” technique can be employed to reduce the computational cost. This technique was first proposed in Chaboche, (1988) and Chaboche & Lesne, (1988) based on the assumption that damage evolution in the bonding contacts is steady in a small range of load cycles. In this sense, the simulation can “jump” from the first load cycle to the final load cycle of the range and the total damage increment in the bonding contact can be obtained by 17 ACCEPTED MANUSCRIPT extrapolating the damage increment in the first load cycle. This approach can considerably reduce the actual number of simulated load cycles, hence saving the overall computational time. In this study, the number of omitted cycles in a “cycle jump”, cuTe6 , is determined using the following equation: cuTe6 = @:uTe6 @:e k (40) RI PT where @:uTe6 is an estimated damage increment after a cycle jump and @:e k is the maximum damage increment in all bonding contacts at the first load cycle in the cycle jump range. As the stress range in the bonding contact is considered to be unchanged in a cycle SC jump, the total increment of fatigue damage variable after the cycle jump can be derived as @:7 = AvwxV uBC R hu √: u hC √:C e Y @:7C M AN U follows: (41) where : u = :C + ∑zBC @:7y , which is the total damage variable of the bonding contact at the u%C beginning of load cycle { | in the cycle jump; hu is the geometrical factor corresponding to TE D : u and @:7C is the incremental fatigue damage at the first load cycle in the cycle jump. Details of the derivation of Eq. (41) is presented in Appendix A. Model behaviour and the influences of fatigue parameters In this section, the model behaviour at the constitutive level is examined by performing EP simulations of two particles connected by a bonding contact. The particle diameter and contact width are both equal to 0.1 m. The model parameters assigned to the bonding contacts are given in Table 1. The initial value of damage variable :)4) is assumed to be 1%. The AC C 4. constitutive responses of the model are illustrated in two fundamental stress modes (pure tension and pure shear) as shown in Fig. 6. For each mode, both monotonic and cyclic tests are conducted to show the correlations between the damage and fatigue behaviour at the constitutive level. In the monotonic test, the fatigue damage is deactivated in the model by setting the fatigue threshold to be greater than one. As a result, the fatigue variables remain zero and the constitutive model formulation is transformed to the original model formulation in Nguyen et al., (2017b). In the cyclic tests, a cyclic loading regime with the maximum stress level of 60 kPa and minimum stress level of 20 kPa is applied to the contact until the damage variable reaches 0.4 (see Fig. 6). Subsequently, the contact is loaded monotonically 18 ACCEPTED MANUSCRIPT until failure. As the maximum cyclic stress can reach to 60% of the contact cohesion, the fatigue threshold has to be chosen less than 0.6 so that fatigue damage can develop in the bonding contact. Table 1: The values of the model parameters used in the pure tension and pure shear tests SC 20.0 10.0 120 100 2.0×10-5 2.0×10-5 10 5 RI PT Selected value M AN U Model parameters Elasto-plastic parameters /49 (GPa/m) /,9 (GPa/m) 8}9 (kPa) K 9 (kPa) <4 (m) <, (m) ϕ (0) ψ (0) Fatigue parameters K EP TE D 1.0 1.0×10-8 Fig. 6: Illustrations of the constitutive tests in pure tension and pure shear modes and the AC C cyclic loading regime 4.1. Pure tension and pure shear test results Fig. 7 illustrates the results of the pure tension tests. The normal stress-displacement curve obtained in the monotonic test follows the same pattern as in the original model (Nguyen et al., 2017b), including a linear elastic stage followed by a softening stage when the bonding contact is yielding (see Fig. 7a). In the cyclic tests, the fatigue damage variable is shown to evolve with increasing the number of load cycles (see Fig. 7b). This is also observed from the gradual reduction of the slope of loading/unloading paths as shown in Fig. 7a. Along with the development of fatigue damage, the fatigue displacement in the normal direction is shown to increase. The increments of normal fatigue displacement and fatigue damage variable are 19 ACCEPTED MANUSCRIPT induced from the hysteresis loops of loading/unloading paths, representing the micro-crack growth and deforming mechanisms occurring in the fatigue damage process of cemented materials. When damage variable in the contact is equal to 0.4, normal stress in the bonding contact is still under the yielding limit. Thus, under a continuous strain increment, the normal stress increases linearly until reaching the yield criterion, then softening in the same path as RI PT in the monotonic test. Fig. 8 presents the results obtained in the pure shear test. It shows that when the damage variable of the bonding contact is equal to 0.4, shear stress in the bonding contact has just reached the yielding limit. As a result, when the monotonic load is applied, the plastic damage immediately develops in the bonding contact and the shear stress is TE D M AN U SC softened without any elastic stress increment as in the pure tension test. Fig. 7: (a) The stress-displacement curves obtained in the pure tension tests; (b) the evolution AC C EP of damage variable in the pure tension cyclic test Fig. 8: (a) The stress-displacement curves obtained in the pure shear tests; (b) the evolution of damage variable in the pure shear cyclic test 20 ACCEPTED MANUSCRIPT To examine the effectiveness of the proposed cycle jump technique, the pure tension and shear cyclic tests are repeated in which the cycle jump technique is employed. The value of @:uTe6 is chosen equal to 0.005 and 0.01. Fig. 9 shows the stress-displacement curves for normal and shear tests obtained in the simulations using the cycle jump technique. The actual number of simulated load cycles is shown to significantly reduce when using the cycle jump technique. Specifically, the numbers of simulated load cycles until the damage variable RI PT reaches 0.4 in both the pure tension and shear tests reduce by 388.8% when @:uTe6 = 0.005 and by 797.4% when @:uTe6 = 0.01, as compared to the simulation without the cycle jump technique. Fig. 10 shows a comparison of the fatigue damage evolution results between the SC original simulation without the cycle jump technique and those using the cycle jump technique. The curves of fatigue damage evolution versus the number of load cycles obtained M AN U in the simulations using the cycle jump technique excellently match with those obtained from the original simulations. These results demonstrate the effectiveness of the cycle jump technique for reducing the overall computational time in fatigue simulations while still AC C EP TE D ensuring the accuracy of numerical results. 21 ACCEPTED MANUSCRIPT Fig. 9: The stress-displacement results in the pure tension and pure shear cyclic tests using M AN U SC RI PT the cycle jump technique with (a) @:uTe6 = 0.005 and (b) @:uTe6 = 0.01 Fig. 10: The influences of cycle jump technique on the results of fatigue damage evolution in (a) the pure tension cyclic test and (b) the pure shear cyclic test 4.2. The effects of fatigue parameters on the constitutive behaviour The influences of fatigue parameters on the constitutive behaviour of the bonding contact are TE D investigated in this section to provide a basis for calibrations of model parameters in numerical experiments. The fatigue parameters of the proposed model consist of and K, which represent the fatigue properties of the bonding contact. Their effects are investigated by repeating the pure tension and pure shear cyclic tests. However, the values of parameters of parameter EP and K are varied from 0.8 to 1.2 and from 0.5×10-8 to 1.5×10-8, respectively. The influence on the results of damage evolution against the number of cyclic loads in the AC C pure tension cyclic test is depicted in Fig. 11a. The increment of fatigue damage in the bonding contact is shown to be sensitive to the value of . More specifically, it requires 2008 load cycles for the fatigue damage variable to increase from 0.01 to 0.4 when the number of load cycles reduces to only 50 when = 0.8, but = 1.2 for the same damage increment. This can be explained based on Eq. (27), which shows the increment of fatigue damage variable is proportional to the stress intensity factor power of . Therefore, a small change of would significantly affect the evolution of fatigue damage variable. Fig. 11b displays the results of fatigue damage evolution with respect to the number of load cycles when varying K. The increment of fatigue damage variable is shown to be proportional to the value of K, which is consistent with Eq. (27). Note that as the same loading magnitude is used in both the 22 ACCEPTED MANUSCRIPT pure tension and pure shear cyclic tests, the results of fatigue damage evolution against the SC RI PT number of load cycles in these tests exactly match with each other. M AN U Fig. 11: Influences of fatigue parameters on the damage evolution in the pure tension cyclic tests. The results of the pure shear cyclic test exactly match with those in the pure tension cyclic tests 5. Applications of the proposed modelling approach The proposed modelling approach is now applied to simulate the experiments of cemented TE D materials. The four-point bending test of cemented-treated samples subjected to both monotonic and cyclic loads carried out by Sounthararajah et al., (2018) are chosen for the simulation. This four-point bending test has been used in the pavement industry to develop practical design models and measure the material fatigue properties for inputs in pavement EP designs. The numerical setup of this test is first described, followed by the calibration of model parameters. Through the calibration process, the capability of the proposed modelling AC C method in capturing the fatigue behaviour of cemented materials is evaluated by comparing numerical results with experimental data. To further examine the predictive capability of the proposed approach, the calibrated model parameters are applied to simulate the fatigue test at different stress levels and numerical results of the S-N curve are then compared with the experimental results. 5.1. Simulation setup The experimental and numerical setups of the four-point bending tests are illustrated in Fig. 12. The beam samples are 200 mm in length, 50 mm in height and 50 mm in width, placed on two supporting rollers at the span of 150 mm and loaded by two loading rollers placed on the 23 ACCEPTED MANUSCRIPT top surface of the sample at the distance of 50 mm. A numerical sample in the same dimensions is generated in DEM by an assemblage of poly-dispersed particles with radii ranging from 1 mm to 19 mm. The particle size distribution of the assemblage is taken similar to the aggregate size distribution of the physical material (see Fig. 13). Grains with diameters less than 1 mm are considered to not participate well in the force transmission of RI PT the aggregate skeleton. Thus, they are simulated by particles of 1 mm diameter for more computational efficiency. All of these particles are then randomly distributed in the sample domain and are allowed to rearrange until reaching an equilibrium condition. The supporting and loading rollers are modelled in DEM by wall elements. Only the intersecting areas SC between the rollers and samples, instead of the entire rollers, are modelled in order to reduce the allocated searching domain in PFC 5.0, thereby saving the computational time. The quasistatic condition as described in Section 2 is set for these simulations. The proposed fatigue M AN U model is subsequently incorporated into all the inter-particle contacts in the sample. On the other hand, the classical cohesionless model developed by Cundall & Strack, (1979) is incorporated into the wall-to-particle contacts and frictionless interactions are set for all of these contacts. Next, the numerical sample is loaded by applying velocities to the loading walls. The TE D monotonic and cyclic loading patterns in the experiments are reproduced in the simulations using a wall servo-control technique. The basis of this technique is to control the velocity of the loading walls in a way that the interactional forces between the walls and samples match the objective loading patterns. During the loading processes, the values of applied load, •, and mid-span deflection, @, are continuously recorded. The flexural stress, flexural strain and equations: εb = Eb = P‚ W„ $ AC C FS = EP flexural modulus of the material are then calculated from • and @ using the following (42) 108„@ 23‚$ (43) 23• ‚ˆ 108@ W„ ˆ (44) where FS is the flexural stress, ‚ is the support span, ‰ is the average beam width, „ is the average beam height, εb is the flexural strain, Eb is the flexural modulus, • is the maximum magnitude of applied cyclic load and @ is the sample deflection corresponding to • . 24 ACCEPTED MANUSCRIPT Similar to the experiments, the initial flexural modulus EbŠ‹Š is determined as the average value of the flexural modulus results obtained in first 50 load cycles in the fatigue test. The fatigue cracking process of the sample can be observed at the macroscopic scale through the curve of Eb /EbŠ‹Š versus the number of load cycles, or so-called the flexural modulus SC RI PT degradation curve. M AN U Fig. 12: The experimental and numerical setups of the four-point bending test. Particles with AC C EP TE D diameters more than 1mm are coloured darker than particles with diameters equal to 1mm Fig. 13: The particle size distributions of experimental and numerical samples 5.2. Model calibration The proposed constitutive model consists of two sets of parameters: elasto-plastic parameters and fatigue parameters. These parameters are calibrated based on a trial and error process by directly comparing the simulation results to those of the experiments. The elasto-plastic parameters are first calibrated with the monotonic test data. Details of the calibration procedure can be found in Nguyen et al., (2017b). The calibrated values of elasto-plastic parameters are then used in the simulation of fatigue test to further calibrate the fatigue 25 ACCEPTED MANUSCRIPT parameters. The calibration procedure of the fatigue parameters is as follows. The fatigue threshold is first calibrated. This parameter controls the number of contacts subjected to fatigue in the sample, hence affecting the sample deformation corresponding to its fatigue damage evolution. The fatigue threshold is calibrated to fit the deformation of experimental samples at the final stage, which is represented by the value of Eb /EbŠ‹Š . Parameters is then calibrated to fit the flexural modulus degradation curve at the initial fatigue stage. Finally, RI PT parameter K is calibrated to fit the fatigue life of a numerical sample with its experimental counterpart. The deformable method proposed in Potyondy & Cundall, (2004) is also adopted in the calibration process to make the constitutive law independent of the contact size. Accordingly, SC the contact normal and shear stiffness are replaced by the effective modulus, Œ 9 , and normal- M AN U to-shear stiffness ratio, /49 ⁄/,9 , in the calibration. The same values of the effective modulus and normal-to-shear stiffness ratio are assigned to all bonding contacts in numerical samples. The normal and shear stiffnesses of each bonding contact are then back-calculated using the relationship given in Nguyen et al., (2017b). The calibrated values of the model parameters for the chosen quasi-static condition (i.e. global damping coefficient of 0.7) are presented in simulations. TE D Table 2. The initial value of damage variable :)4) is assumed to be equal to 1% in all Table 2: The values of calibrated model parameters for the flexural tests AC C EP Model parameters Elasto-plastic parameters Œ 9 (GPa) /49 ⁄/,9 8}9 (kPa) K 9 (kPa) <4 (m) <, (m) ϕ (0) ψ (0) Fatigue parameters K a Calibrated values 1.55 2.4 74 74 2.0×10-5 2.0×10-5 20 10 1.2 1.8×10-11 0.84 Fig. 14 compares the flexural stress-strain curves obtained in the experiments and simulation for the monotonic test. It shows that the numerical curve falls within the experimental range. 26 ACCEPTED MANUSCRIPT During the loading process, the cemented sample undergoes three stages including a linear elasticity stage, a non-linear hardening stage prior to peak stress and a residual stage after the peak. All of these observed responses are captured well in the simulation. In addition, the cracking pattern in the numerical sample agrees with the experimental result as illustrated in Fig. 15. In both the experiment and simulations, the pattern of major cracks can slightly vary RI PT depending on the arrangement of particles, but the major cracks always initiate at the bottom of the beam in the area between the two loading rollers and then propagates upwards along particle interfaces. The localised fracture and crack development are naturally captured in the simulation without any ad-hoc treatment on the contact model. This demonstrates the SC effectiveness of the proposed method in modelling the damage response of cemented EP TE D M AN U materials under monotonic loads. AC C Fig. 14: The flexural stress-strain curves obtained in the experiment and simulation of the monotonic test. Fig. 15: Sample crack patterns in the experiment and simulation of the monotonic test The numerical fatigue test is then performed with the maximum cyclic stress equal to 1029 kPa, which is 77% of the maximum flexural stress obtained in the monotonic test (FS•Ž• = 1336.4 kPa). This corresponds to the maximum cyclic load of 857.5 N according to Eq. (42). 27 ACCEPTED MANUSCRIPT The cyclic load in haversine shape as illustrated in Fig. 16 is continuously applied to the numerical sample until failure. The cycle jump technique is also employed in this simulation with @:uTe6 = 0.003. This chosen value of @:uTe6 was checked and shown to marginally affect the fatigue response of the sample at macro-scale. Fig. 17a shows a comparison of the flexural modulus degradation curves obtained in the RI PT experiment and simulations using the proposed fatigue model and the original model (Nguyen et al., 2017b). The reduction of flexural modulus in the experiment is observed to go through three stages including the initial, middle and final stages. The simulation using the proposed fatigue model can closely capture the gradual reduction of flexural modulus in the SC middle stage and the unstable reduction of flexural modulus at the final stage. However, there stands a deviation in the initial stage in which the flexural modulus in the simulation is M AN U observed to drop faster than the flexural modulus in the experiment. This can be attributed to the difference in microstructural configuration between the experimental and numerical samples. The reduction of flexural modulus during the cyclic loading process is the consequence of damage evolution in the numerical sample, evidenced by the reduction of the slope of loading/unloading curves with increasing the number of load cycles as can be seen in Fig. 17b. Along with the damage evolution, permanent deformation gradually increases in the TE D sample, exhibited by that flexural strain values measured at the end of unloading paths gradually shifts away from the origin (see Fig. 17b). On the contrary to the simulation result of the proposed fatigue model, the simulation using the original cohesive model is unable to capture the reduction of the flexural modulus as well as the sample failure in the cyclic test. EP The flexural modulus in the simulation with the original cohesive model remains constant with increasing the number of load cycles (see Fig. 17a). These results clearly show the AC C capability as well as the necessity of the proposed fatigue model for the DEM simulations of fatigue damage in cemented materials. 28 ACCEPTED MANUSCRIPT Fig. 16: The shape of the applied cyclic load used in the fatigue test Fig. 18 illustrates the damage evolution and crack development in the numerical sample during the cyclic test. In the initial stage, damage starts to evolve in several places in the area between two loading rollers but no clear fracture is observed in the numerical sample. During RI PT the middle stage, damage is shown to localise at the bottom mid-span of the sample, then gradually developing with increasing number of load cycles. At the final stage, a major crack initiates and rapidly propagates from the bottom to the top surface of the sample. The final crack pattern in the numerical sample matches well with the crack patterns observed in the monotonic test and in the experiment. This result indicates the capability of this modelling TE D M AN U SC approach in reproducing fatigue cracking in cemented samples. Fig. 17: (a) The flexural modulus degradation curves obtained in the experiment and EP simulations with the cyclic stress level equal to 77% of the maximum flexural stress; (b) The AC C flexural stress-strain result obtained in the simulation of the monotonic and fatigue tests 29 SC RI PT ACCEPTED MANUSCRIPT Fig. 18: Damage evolution in the numerical sample during the fatigue test M AN U 5.3. Numerical prediction In Section 5.2, the proposed modelling approach has been shown to be able to capture the fatigue behaviour of cemented materials. In this section, additional simulations of the fatigue test at different maximum cyclic stress levels (i.e. 95%, 92%, 84%, 81% and 70% of FS•Ž• ) are carried out to examine its predictive capability. The numerical sample and boundary condition remain the same as in Section 5.1. The same model parameters as given in Table 2 TE D are used in all of these simulations. The results of S-N curves obtained in the simulation are then compared with the S-N curves of the experiments. Similar to the experimental method, the fatigue life of the numerical sample is defined as the number of load cycles applied to that AC C EbŠ‹Š = 50%. EP sample until the flexural modulus is equal to 50% of the initial flexural modulus, i.e. Eb / 30 ACCEPTED MANUSCRIPT Fig. 19: (a) The S-N curve obtained in the experiment and simulation of the flexural fatigue test; (b) The flexural modulus degradation curves obtained in the numerical fatigue tests with different maximum cyclic stress levels Fig. 19a displays the S-N curves obtained in the experiments and simulations. It shows that, RI PT given a single set of model parameters and a microstructural configuration of the numerical sample, the proposed modelling approach can predict the decrease of sample fatigue life with increasing the maximum cyclic stress level. Similar to the experiment, the results of fatigue life obtained in the simulations with the maximum stress level ranging from 77% to 95% of SC FS•Ž• relatively fall to a logarithmic trend line closing to the experimental trend line. The small deviation is considered due to the differences in the microstructural configurations of the numerical and experimental samples. On the other hand, the numerical sample subjected M AN U to the cyclic stress level equal to 70% of FS•Ž• is shown to be self-stabilised after hundreds of load cycles. This can be attributed to the shakedown effect, in which particles and contacts in the sample gradually rearrange to a more stable condition when undergoing cyclic loading. The difference between the experimental and numerical results, in this case, can be attributed to the randomness of sample microstructures in the experiment and simulation. Fig. 19b depicts the flexural modulus degradation curves obtained in the numerical fatigue tests with TE D different stress levels. It can be seen that when increasing the applied cyclic stress level, the sample fails earlier and the value of Eb /EbŠ‹Š at the failure points, where the curves suddenly drop, is also larger. This is due to the fact that the reduction of Eb /EbŠ‹Š is as a result of fatigue EP damage evolution in the sample, while the sample failure occurs when the applied cyclic stress reaches the yield criterion after a number of load cycles. A higher cyclic stress level AC C applied to numerical samples results in greater stress magnitudes in the bonding contacts. Hence, the contacts can reach their yield criterion faster while the evolution of fatigue damage is still small. The sample also fails earlier and the value of Eb /EbŠ‹Š at the failure points is greater than those of the fatigue tests with lower cyclic stress levels. 5.4. The influence of particle arrangement on the numerical results In Section 5.2, the model parameters are calibrated with a specific microstructural configuration of the numerical sample. However, in the experiments of cemented materials, the microstructure configuration of samples would vary owing to the randomness of aggregates and cement bonds arrangement in different samples. This is one of the main causes of highly scattered results observed in the fatigue tests of cemented materials. As the 31 ACCEPTED MANUSCRIPT proposed modelling approach can directly replicate the heterogeneous microstructure of cemented materials through the sample generation process, the influence of particle arrangement on the material behaviour can be naturally reproduced in simulations. To examine the influence of particle arrangement on numerical results, four additional beam samples are generated, namely Samples 1 to 4. The sample dimensions and particle size RI PT distribution remain the same as the original sample given in Section 5.1, yet the particle arrangement is varied differently in each new sample. The same values of model parameters as given in Table 2 are assigned to all bonding contacts in these new samples. Both the monotonic and fatigue tests at the maximum stress level of 1029 kPa (i.e. 77% of FS•Ž• of SC the original sample) are then performed on the new samples and numerical results are compared with those of the original sample. M AN U Fig. 20 compares the flexural stress-strain curves and the flexural modulus degradation curves obtained in numerical tests of the new samples with those of the original sample. It can be seen while the results of the monotonic test are relatively consistent and falling well within the experimental range (Fig.20a), the results of the fatigue test vary substantially among the numerical sample (Fig.20b). The scatter results of the fatigue tests can be attributed to the difference in the microstructure of each numerical sample due to the TE D rearrangement of particles. In fact, the variation of the material microstructure can affect the distribution of local stress in bonding contacts, causing changes of the damage evolution rate per load cycle in each bonding contact. As the fatigue behaviour of the sample at the macroscopic scale is the collective result of all the damage evolution per load cycle, the EP change of damage evolution in every load cycle would result in a considerable deviation in the final fatigue life of the samples. These responses are consistent with the experimental AC C observation of cemented materials, which normally shows high scattered results of fatigue life in different samples (Do et al., 1993; Dyduch et al., 1994). 32 ACCEPTED MANUSCRIPT Fig. 20: The variation of flexural strength results (a) and fatigue results (b) obtained in the simulations when varying the particle size distribution of samples To provide an insight into the scatter results of the fatigue test, the coordination number is employed to explain the relationship between the particle arrangement in numerical samples RI PT and their mechanical responses. In DEM modelling, the coordination number of a sample is a determining factor of the sample packing quality, which is obtained by the following equation: C4 = 2c c6 (45) where c and c6 are the total number of contacts and particles in the sample, respectively. SC Fig. 21 presents the relationships between the coordination numbers C4 and the results of M AN U elastic flexural modulus and fatigue life of the samples, respectively. The elastic flexural modulus and fatigue life may increase with the increments of the coordination number in the numerical samples. These observations are reasonable in the sense that the increment of the coordination number represents the increase of the number of contacts to encounter loads in numerical samples, which improves the sample stiffness and load-carrying capacity. However, it is acknowledged that the fatigue response of numerical samples also depends on how stress TE D transmits and concentrates in the samples during the cyclic loading process. The magnitude of stress concentration can vary from one sample to another due to the differences in their microstructure configuration. As the evolution of fatigue damage relies on the magnitude of stress in accordance with Paris’ law, the variation of local stress causes the rate of fatigue EP damage evolution to be varied in different samples. This microscopic mechanism contributes to the scattered results observed in the numerical fatigue tests of cemented materials. As a AC C result of such microstructural effect, the fatigue life of samples 2 and 3 is shown to be quite similar although their coordination numbers are relatively different (see Fig. 21). This behaviour can be further explained by having a closer look at the stress distribution in these numerical samples. Fig. 22 illustrated the distribution of tensile stress in the numerical samples under the application of flexural stress equal to 400 kPa. It shows that even though the same stress is applied to these samples, the magnitude of tensile stress in bonding contacts at the location where cracks are initiated (zoomed area) of sample 3 is greater than those in other samples. This suggests that the rate of fatigue damage evolution in sample 3 is higher, and thus this sample can reach its fatigue life faster than the others. 33 RI PT ACCEPTED MANUSCRIPT SC Fig. 21: The relationships between the coordination number and the elastic flexural modulus and fatigue life of the samples M AN U Fig. 23 depicts the crack patterns of the experimental and numerical samples obtained in the fatigue tests. In all experiments and simulations, cracks initiate in the area between the two loading rollers and develop from the bottom to the top surface of the samples. However, the location and shape of cracks are varied in each sample, and this can also be attributed to the difference in the microstructural configuration of each sample due to the rearrangement of particles. The difference of sample microstructure induces the change of stress transmission TE D and concentration in the sample, thus changing the crack location. Besides, as the crack propagates along particle interfaces, the direction of crack propagation highly depends on how particles arrange near the crack tip. The variation of crack patterns can be well captured AC C behaviour. EP by the proposed approach thanks to it physical reproduction of the material microstructural 34 ACCEPTED MANUSCRIPT Fig. 22: The distribution of tensile stress in numerical samples under the application of M AN U SC RI PT flexural stress of 400 kPa Fig. 23: Crack patterns of experimental and numerical samples obtained in the fatigue tests Conclusion A discrete-based approach is proposed in this study to model the fatigue damage of cemented TE D materials. This modelling approach takes advantages of DEM to reproduce the heterogeneous microstructure of cemented materials and capturing real crack development. Alongside with the application of DEM, a new fatigue bonding contact model is developed to characterise the fatigue response of cemented materials at the grain scale. This model is derived on the basis EP of the cohesive damage-plasticity model proposed in Nguyen et al., (2017b), along with fatigue damage variables and a fatigue evolution law, to describe both cracking and irreversible displacements in the fatigue process of cemented materials. Fatigue damage is AC C 6. allowed to evolve when stress is inside the yield surface through the introduction of a fatigue criterion under the yield criterion. The evolution of fatigue damage is then described by a single fatigue law of a scalar fatigue damage variable (i.e. Paris’ law), while different stages of fatigue crack growth exhibited in the materials at the macroscopic scale are captured as the collective response of all bonding contacts and particles in numerical samples. Through numerical flexural tests, the capability of the proposed modelling approach in capturing and predicting the fatigue behaviour of cemented materials under flexural boundary condition is demonstrated. This is evidenced by the good agreements in the flexural modulus degradation curves, the S-N curves and sample crack patterns between the experiments and 35 ACCEPTED MANUSCRIPT simulations. In addition, thanks to the capability of DEM in reproducing physical material microstructure and the explicit characterisation of contact behaviour in the new fatigue model, the proposed modelling approach is able to reproduce the microstructural effects on the fatigue behaviour of cemented materials. The variations of stress-strain results, fatigue life and crack patterns in the samples with changing particle arrangement can be naturally RI PT captured in the DEM simulations. This demonstrates the advantages of the proposed approach over existing continuum methods for modelling the fatigue response and localised fracture in cemented materials, while still retaining the simplicity of contact model formulation. Moreover, the incorporation of fatigue damage mechanisms of cement bridges SC into the DEM bonding contact model is shown to be necessary given the existing damageplastic contact model is unable to capture these mechanisms in the DEM simulations at the M AN U grain scale. Acknowledgement Funding supported by Australian Research Council (grant number LP130100884 and DP 160100775) is gratefully acknowledged. Appendix A: Derivation of the total increment of fatigue damage variable in a cycle jump TE D ∆d is first derived by substituting Eq. (32) and (33) into (28) as: ∆d = h√jf:.∆84$ + ∆8,$ EP The incremental fatigue damage in a load cycle is then obtained by substituting the above equation into Eq. (27) as: K ∆d @:7 = f K;√jf.∆84$ + ∆8,$ = = ;h√:= f e AC C e e Hence, the incremental fatigue damage in each load cycle of a cycle jump is calculated as: K;√jf.∆84$ + ∆8,$ = @:7C = mhC .:C n f e @:7$ … = mh e K;√jf.∆84$ + ∆8,$ = :$ n f $. e e e h $ √:$ =R Y @:7C C C h √: 36 ACCEPTED MANUSCRIPT u @:7 K;√jf.∆84$ + ∆8,$ = u. u = mh : n f e e =R hu √: u hC √:C e Y @:7C The total increment of fatigue damage variable in a cycle jump is then obtained as follows: uBC u @:7 = AvwxV uBC Reference hu √: u e R Y @:7C C C h √: RI PT @:7 = AvwxV SC Al-Rub, R. K. A. & Darabi, M. K. (2012). A thermodynamic framework for constitutive modeling of time-and rate-dependent materials. Part I: Theory. International Journal of Plasticity, 34, 61-92. Alliche, A. (2004). Damage model for fatigue loading of concrete. International Journal of Fatigue, 26, 915-921. M AN U Carrara, P. & De Lorenzis, L. (2015). A coupled damage-plasticity model for the cyclic behavior of shear-loaded interfaces. Journal of the Mechanics and Physics of Solids, 85, 33-53. Chaboche, J.-L. (1988). Continuum damage mechanics: Part II—Damage growth, crack initiation, and crack growth. Journal of applied mechanics, 55, 65-72. Chaboche, J. L. & Lesne, P. M. (1988). A non‐linear continuous fatigue damage model. Fatigue & fracture of engineering materials & structures, 11, 1-17. TE D Cheng, J., Qian, X. & Zhao, T. (2016). Rheological viscoplastic models of asphalt concrete and rate-dependent numerical implement. International Journal of Plasticity, 81, 209230. Cundall, P. A. & Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. Géotechnique, 29, 18. AC C EP Darabi, M. K., Al-Rub, R. K. A., Masad, E. A., Huang, C.-W. & Little, D. N. (2012). A modified viscoplastic model to predict the permanent deformation of asphaltic materials under cyclic-compression loading at high temperatures. International Journal of Plasticity, 35, 100-134. Dattoma, V., Giancane, S., Nobile, R. & Panella, F. W. (2006). Fatigue life prediction under variable loading based on a new non-linear continuum damage mechanics model. International Journal of Fatigue, 28, 89-95. Do, M.-T., Chaallal, O. & Aïtcin, P.-C. (1993). Fatigue behavior of high-performance concrete. Journal of Materials in civil Engineering, 5, 96-111. Dyduch, K., Szerszeń, M. & Destrebecq, J.-F. (1994). Experimental investigation of the fatigue strength of plain concrete under high compressive loading. Materials and Structures, 27, 505-509. Fathima, K. M. P. & Kishen, J. M. C. (2013). A thermodynamic framework for fatigue crack growth in concrete. International Journal of Fatigue, 54, 17-24. Gao, L. & Hsu, C.-T. T. (1998). Fatigue of concrete under uniaxial compression cyclic loading. Materials Journal, 95, 575-581. 37 ACCEPTED MANUSCRIPT Gui, Y., Bui, H. H. & Kodikara, J. (2015). An application of a cohesive fracture model combining compression, tension and shear in soft rocks. Computers and Geotechnics, 66, 142-157. Gui, Y. L., Zhao, Z. Y., Kodikara, J., Bui, H. H., & Yang, S. Q. (2016). Numerical modelling of laboratory soil desiccation cracking using UDEC with a mix-mode cohesive fracture model. Engineering Geology, 202, 14-23. RI PT Horii, H., Shin, H. C. & Pallewatta, T. M. (1992). Mechanism of fatigue crack growth in concrete. Cement and concrete composites, 14, 83-89. Irwin, G. R. (1957). Analysis of stresses and strains near the end of a crack traversing a plate. Journal of applied mechanics, 24, 361-364. Itasca (2014). PFC 5.0 User Manual. Itasca Consulting group, Inc, Washington DC, United State, p. 210. SC Krajcinovic, D. & Fanella, D. (1986). A micromechanical damage model for concrete. Engineering Fracture Mechanics, 25, 585-596. M AN U Liang, R. Y. & Zhou, J. (1997). Prediction of fatigue life of asphalt concrete beams. International Journal of Fatigue, 19, 117-124. Mohamed, A. R. & Hansen, W. (1999). Micromechanical modeling of concrete response under static loading—part 1: model development and validation. Materials Journal, 96, 196-203. Moreo, P., Garcia-Aznar, J. M. & Doblare, M. (2007). A coupled viscoplastic rate-dependent damage model for the simulation of fatigue failure of cement–bone interfaces. International journal of plasticity, 23, 2058-2084. TE D Nguyen, N., Bui, H. H., Nguyen, G. D., Kodikara, J., Arooran, S. & Jitsangiam, P. (2017a). A thermodynamics-based cohesive model for discrete element modelling of fracture in cemented materials. International Journal of Solids and Structures, 117, 159-176. EP Nguyen, N. H. T., Bui, H. H., Nguyen, G. D. & Kodikara, J. (2017b). A cohesive damageplasticity model for DEM and its application for numerical investigation of soft rock fracture properties. International Journal of Plasticity, 98, 175-196. AC C Nguyen, T. T., Bui, H. H., Ngo, T. D., & Nguyen, G. D. (2017c). Experimental and numerical investigation of influence of air-voids on the compressive behaviour of foamed concrete. Materials & Design, 130, 103-119. Nguyen, O., Repetto, E. A., Ortiz, M. & Radovitzky, R. A. (2001). A cohesive model of fatigue crack growth. International Journal of Fracture, 110, 351-369. Oñate, E., Zárate, F., Miquel, J., Santasusana, M., Celigueta, M., Arrufat, F., Gandikota, R., Valiullin, K. & Ring, L. (2015). A local constitutive model for the discrete element method. Application to geomaterials and concrete. Computational Particle Mechanics, 2, 139-160. Papa, E. & Taliercio, A. (1996). Anisotropic damage model for the multiaxial static and fatigue behaviour of plain concrete. Engineering Fracture Mechanics, 55, 163-179. Paris, P. C. & Erdogan, F. (1963). A critical analysis of crack propagation laws. Journal of basic engineering, 85, 528-533. Paris, P. C. & Sih, G. C. (1965). Stress analysis of cracks. Fracture toughness testing and its applications. ASTM International. 38 ACCEPTED MANUSCRIPT Potyondy, D. O. & Cundall, P. A. (2004). A bonded-particle model for rock. International Journal of Rock Mechanics and Mining Sciences, 41, 1329-1364. Ray, S. & Kishen, J. M. C. (2010). Fatigue crack propagation model for plain concrete–An analogy with population growth. Engineering Fracture Mechanics, 77, 3418-3433. Ristinmaa, M. (1995). Cyclic plasticity model using one yield surface only. International journal of plasticity, 11, 163-181. RI PT Scholtès, L. & Donzé, F.-V. (2013). A DEM model for soft and hard rocks: Role of grain interlocking on strength. Journal of the Mechanics and Physics of Solids, 61, 352-369. Shen, W. Q., Shao, J.-F., Kondo, D. & Gatmiri, B. (2012). A micro–macro model for clayey rocks with a plastic compressible porous matrix. International journal of plasticity, 36, 64-85. SC Shen, W. Q. & Shao, J. F. (2016). An incremental micro-macro model for porous geomaterials with double porosity and inclusion. International Journal of Plasticity, 83, 37-54. M AN U Simon, K. M. & Kishen, J. M. C. (2017). A multiscale approach for modeling fatigue crack growth in concrete. International Journal of Fatigue, 98, 1-13. Sounthararajah, A., Bui, H. H., Nguyen, N. H. T., Jitsangiam, P. & Kodikara, J. (2018). Early-Age Fatigue Damage Assessment of Cement-Treated Bases Under Repetitive Heavy Traffic Loading. Journal of materials in civil engineering, (in press). Tada, H., Paris, P. C. & Irwin, G. R. (1973). The stress analysis of cracks. Handbook, Del Research Corporation. TE D Turon, A., Costa, J., Camanho, P. P. & Dávila, C. G. (2007). Simulation of delamination in composites under high-cycle fatigue. Composites Part A: applied science and manufacturing, 38, 2270-2282. Voyiadjis, G. Z. & Abu-Lebdeh, T. M. (1994). Plasticity model for concrete using the bounding surface concept. International Journal of Plasticity, 10, 1-21. EP Weichert, D. & Maier, G. (2014). Inelastic behaviour of structures under variable repeated loads: direct analysis methods, Springer. Xiao, Y. C., Li, S. & Gao, Z. (1998). A continuum damage mechanics model for high cycle fatigue. International journal of fatigue, 20, 503-508. AC C Yang, B., Mall, S. & Ravi-Chandar, K. (2001). A cohesive zone model for fatigue crack growth in quasibrittle materials. International Journal of Solids and Structures, 38, 3927-3944. Zhang, J. (2001). Size effect on fatigue in bending of concrete. Journal of materials in civil engineering, Vol. 13. Zhang, J., Shen, W. Q., Oueslati, A. & De Saxce, G. (2017). Shakedown of porous materials. International Journal of Plasticity, 95, 123-141. 39 ACCEPTED MANUSCRIPT Highlights AC C EP TE D M AN U SC RI PT • The fatigue behaviour of cemented material can be effectively captured by DEM combining with a new fatigue contact model • The fatigue contact model coupling damage mechanics and plasticity theory can describe well the degrading response of cemented materials subjected to cyclic loading • This proposed modelling method facilitates the characterisation of fatigue damage in cemented materials for better understanding and design of the material behaviour

1/--страниц