Two scale models for fracture behaviours of cementitious materials subjected to static and cyclic loadings

This study employs a lattice fracture model to simulate static and fatigue fracture behaviour of Interfacial Transition Zone (ITZ) at microscale and mortar at mesoscale. The heterogeneous microstructure of ITZ and mesostructure of mortar are explicitly considered in the models. The initial step involves calibrating and validating the microscopic model of the ITZ through micro-cantilever bending tests. Subsequently, this validated ITZ model serves as a constitutive law to simulate the fracture behavior of mortar at the mesoscale using an uncoupled upscaling method. The influence of microstructural features, such as w/c ratio and microscopic roughness, on the fracture behaviour of ITZ is investigated. Moreover, the effect of ITZ properties and stress level on the fracture performance and fatigue damage evolution of mortar is also studied. The simulation results for both the ITZ and mortar demonstrate good agreement with experimental results. The proposed two models provide insights into the fracture mechanisms and fatigue damage evolution in cementitious materials subjected to static and cyclic loadings.


Introduction
Concrete is generally considered as a three-phase composite at the mesoscale (millimeter to centimeter scale), comprising aggregate, cement paste matrix and the interfacial transition zone (ITZ) -a highly porous zone between the aggregate particles and the cement paste matrix.The ITZ plays a pivotal role in determining the mechanical performance and durability of concrete structures.Understanding and simulating the ITZ have become imperative in the field of concrete research and engineering.
Accurate simulation of the ITZ behaviour necessitates a thorough understanding of its distinct microstructure compared to the bulk of cement paste and aggregates.Numerous studies have been conducted to investigate the microstructure development of ITZ [1][2][3][4][5][6].The ITZ is generally characterized by high porosity, deposition of crystals with preferential orientations and fewer cement particles than the bulk cement matrix.The major origin of this porous ITZ is the so-called 'wall effect', which disrupts the packing of the cement grains against the relatively flat aggregate surface [1,7].A reduced amount of cement grains leads to a more porous microstructure during the hydration process.Moreover, a water film which forms around the aggregate particles in fresh concrete also results in a higher local w/c ratio in the vicinity of aggregate particles.Other mechanisms, such as the microbleeding, one-side growth and filtration effect, may also contribute to the porous ITZ [8][9][10].Another important microstructural feature of the ITZ is the excess content of large ettringite (AFt) and calcium hydroxide (CH) crystals with preferential orientations deposited in the vicinity of aggregate, which forms the porous framework of ITZ at the early age of hydration [6].With ongoing hydration, the remaining space is further filled with calcium silicate hydrate (C-S-H) gel and smaller crystals of AFt and CH.It is believed that incorporating as many microstructural characteristics of the ITZ as possible in the model is beneficial for accurately simulating its microscopic mechanical properties.
Several attempts have been made to directly simulate the micromechanical properties of ITZ [11][12][13][14][15][16], which are difficult to measure experimentally.Sun et al. [12] used the Virtual Cement and Concrete Testing Laboratory (VCCTL) cement hydration model to generate the virtual microstructures of ITZ at different hydration ages.The finite element method is then applied to calculate the elastic properties of ITZ based on the simulated ITZ microstructures.The estimated ratios of ITZ elastic modulus to the bulk paste modulus are in good agreement with experimental results.On the other hand, Zhang et al. [17] employed X-ray CT scanning tests to acquire realistic ITZ microstructures and mapped them to a lattice fracture model.The experimentally informed model can effectively simulate the strength and elastic modulus of ITZ.In addition, Xi et al. [18] used of meso-scale concrete fracture models to simulate the fracture behaviour of concrete considering the aggregate randomness and ITZ fracture properties.The simulated outcomes were employed to train an artificial neural network (ANN) model.Afterwards, they applied the trained ANN model to inversely determine the fracture properties of ITZ based on the macroscopic testing results.Despite notable progresses in simulating the mechanical behaviour of the ITZ, there are still significant knowledge gaps that require further research to fully understand this complex region.Considerable efforts remain to address many aspects that have not yet been explored through modelling, including the effects of mixture composition, microscopic roughness of aggregate and the mechanical response under the cyclic loading [19][20][21].
In addition to the microscopic simulation of the ITZ mechanical behaviours, numerous studies have been dedicated to evaluating the contribution of ITZ to the overall mechanical performance of concrete through mesoscopic models [22][23][24][25][26][27].Yang [26] formulated a double-inclusion analytical model to investigate the effect of ITZ thickness on the elastic modulus of mortar.In this model, spherical aggregates were assumed to be surrounded by a concentric layer of ITZ.Sun et al. [12] also used a differential effective medium theory to homogenise the elastic properties of the concrete composite, accounting for the thickness and the elastic modulus of ITZ.Meanwhile, others [22][23][24][25] studied the effects of ITZ properties using either discrete element or finite element models.For instance, Zhao et al. [24] adopted a meso-scale discrete element model to simulate the static fracture behaviour of concrete considering the ITZ strength.Mirshad et al. [25] utilised a finite element model to investigate the effects of parameters including ITZ thickness, aggregate shape, volume fraction, and size on the macroscopic mechanical properties of concrete.In exploring cyclic fracture behaviour, Matsumoto et al. [28] and Gong et al. [29] employed a mesoscopic Rigid Body Spring Model to examine the comprehensive fatigue properties of concrete.These investigations took into account the cyclic constitutive laws of the ITZ.Simon and Kishen [30] proposed a multiscale approach to simulate the fatigue crack growth in concrete.By establishing a correlation between the crack opening displacements associated with aggregate bridging on a macroscale and the microcracking occurring at the microscale, a modified stress intensity factor (SIF) is derived.This modified SIF is further employed in developing an analytical model, which predicts the complete crack growth curve for plain concrete.Huang and Zhang [31] developed a degradation model for micromechanical properties and a cycle-dependent mesoscopic crack bridging model to simulate the degradation behaviour of Engineered Cementitious Composites (ECC) under fatigue tensile loading.The macroscopic degradation behaviour of the bridging stress with presence of multiple cracks is then modelled using the extended finite element method.However, the complex effects of microscopic and mesoscopic characteristics cannot be incorporated into the model because of the simplified assumptions made at these scales.A micromechanics-based model explicitly considering the microcrack growth is developed by Dutta and Kishen [32].The state of damage at the macroscale is established through the application of an energy-based criterion.By using the Mori-Tanaka method of homogenization, the proposed model can effectively investigate the influence of several microstructural properties on the fatigue of plain concrete.It needs to be noted that the key to accurately simulating the mortar or concrete samples are the constitutive behaviours of the ITZ, which should be derived from lower-scale simulations.However, in most studies, ITZ properties are determined through inverse analysis method based on mesoscopic models or simplified assumptions.Directly deriving a microscopic model is challenging, primarily due to the absence of experimental calibration and validation.Moreover, there is a missing link between the microscopic model of ITZ and the mesoscopic model of concrete.Recently, an uncoupled parameter-passing upscaling method has emerged as a valid tool to bridge the two scales [33] and has been successfully used to simulate the static fracture behaviours of cementitious materials.It is worthwhile to explore the feasibility of this upscaling method to the fatigue fracture simulation.
In this paper, the flexural and uniaxial tensile fracture behaviours of the ITZ under both static and cyclic loadings are examined using the lattice fracture model.The model explicitly incorporates microstructures of the ITZ obtained from X-ray computed tomography (XCT).Local mechanical and fatigue properties of distinct phases within the ITZ are calibrated and validated using experimental results at the same length scale.The validated model is subsequently employed to predict the uniaxial tensile fracture behavior of mortar under static and cyclic loadings.Digital mesostructures of mortar are established, and the mechanical and fatigue properties of the ITZ and cement paste from microscopic modelling were used as inputs for the mesoscopic model.The investigation focuses on the influence of various ITZ properties on the fracture performance of mortar, with particular attention to the effect of ITZ properties on fatigue damage evolution in the mortar.

Microscopic simulation of ITZ
To properly simulate the mesoscopic fatigue behaviour of mortar, a constitutive model of ITZ at the microscale is needed.In this section, a lattice fracture model is employed to simulate the mechanical and fatigue behaviour of ITZ at the microscale.The realistic ITZ microstructures obtained from XCT results are mapped to the lattice network.The effect of microscale roughness of aggregate on the overall mechanical and fatigue performance of ITZ is investigated.Fig. 1 2

.1.1. Two-dimensional virtual ITZ specimens
To obtain the realistic microstructure, miniature ITZ specimens are fabricated and scanned by a Micro CT-Scanner (Phoenix Nanotom, Boston, MA, USA).Standard grade CEM I 42.5 N Portland cement (ENCI, the Netherlands) and deionized water are first mixed to prepare the fresh paste (w/c ratio 0.3 and 0.4).Two types of quartzite aggregate surfaces, namely artificial flat surface and natural rough surfaces, are prepared to investigate the effect of micro-roughness on the overall mechanical performance of the ITZ.For the flat aggregate, the aggregate is cut using a Minitom low-speed cutting machine, and then ground to obtain a flat surface using a Struers Labopol-5 thin sectioning machine.The other type of aggregate, featuring a naturally rough surface, undergoes cleaning without additional grinding.Following these processes, all aggregate samples are then dried in an oven (60 • C) for 48 hours.After cooling down to the room temperature, the fresh paste is then cast onto the aggregate surface.After sealed curing for 28 days, a precision microdicing machine (MicroAce Series 3 Dicing Saw) is employed to cut the cement paste-aggregate composite, generating micro-cantilever ITZ beams with a square cross-section measuring 150 × 150 μm 2 .The cantilevered length of the beam is around 750 μm.For more experimental details regarding the preparation process of miniatured ITZ specimens, the reader is referred to previous work of the authors [34].Besides, two large ITZ prisms with the size of around 1 × 1 × 1.5 mm 3 (w/c 0.4) containing flat and rough aggregate surfaces are prepared using the similar cutting procedure.Both the ITZ cantilever beams and ITZ prisms are scanned by XCT to obtain greyscale-based 2D images, see Fig. 2. The X-ray source tube used for ITZ cantilever beams is set at 90 kV/170 μA during scanning, which results in a voxel resolution of 0.5 × 0.5 × 0.5 μm 3 .For ITZ prisms, the X-ray source tube is set at 90 kV/110 μA, resulting in a voxel resolution of 1.5 × 1.5 × 1.5 μm 3 .In general, a smaller mesh size is preferred as it allows for more detailed consideration of microstructural information.However, the choice of mesh size is also influenced by factors such as XCT image resolutions and computational efficiency.Given that the recognized width of the ITZ falls within the range of 30-100 μm [35], an element size of 2.5 μm is chosen to capture the microstructural features in the interfacial zone.The spatial resolution of XCT images is subsequently interpolated to 2.5 μm/pixel using the bilinear interpolation algorithm.Consequently, the element size for ITZ simulation is set at 2.5 μm.

Lattice modelling approach
Lattice models are extensively employed for simulating the fracture behavior of cement-based materials, primarily because the simulated cracks closely resemble those observed in both laboratory tests and practical applications [36][37][38].In these models, the material is discretized as a network of lattice beam elements.A group of cells is generated and each cell coincides with a single pixel.Inside each cell, a sub-cell is created, as shown in Fig. 3.By using a pseudo-random number generator, a node is randomly created inside each sub-cell.The length ratio of the sub-cell to the cell is defined as the randomness of the lattice mesh.The value of randomness affects the simulated fracture behaviour of material as well as the Poisson's ratio of the lattice system, as suggested in [39,40].In order to avoid large variations in the length of elements and introduce geometry disorder of texture, a randomness of 0.5 is adopted herein.This value has been proved to produce more realistic and consistent results that closely align with experimental data [41].Delaunay tessellation of the set of nodes was then performed to connect nodes with adjacent Voronoi cells by Timoshenko beam elements [38].
By assigning these beam elements with different properties based on the pixel images obtained from XCT tests, the microstructure of the material can be mapped on the lattice network, see Fig. 3.Besides the aggregate, there are four main phases in the paste matrix: unhydrated cement (UHC), high-density C-S-H (HD), low-density C-S-H (LD) and pores.For the ITZ simulation, the global thresholding approach [37,42] is used to segment these four phases in 2D images randomly selected from XCT results.Three threshold values are defined based on the grey-scale histogram.These threshold values are determined using either the inflection point of the cumulative grey-scale value curve or the change of the tangent slope.The hydration degree and the w/c ratio are also taken into account to determine the LD/HD CSH ratio.Note that other hydration products, e.g.calcium hydroxide (CH), ettringite (AFt) and monosulfate (AFm), are not considered in current study.More  Y. Gan et al. details on the segmentation procedure can be found in [37,42].After specifying the aggregate boundary, hydration that connecting the aggregate surface with the paste matrix are defined as interfacial elements.Therefore, two interfacial elements can be detected, i.e.Aggregate-LD C-S-H and Aggregate-HD C-S-H.Note that the interfacial bonding between the aggregate and UHC is assumed to be extremely weak and is neglected.Elements linked to pores are likewise eliminated from the mesh and treated as initial flaws.
It should be emphasized that conducting a 2D analysis represents a significant simplification of the behaviors observed in 3D specimens.In this study, it is assumed that properties along the third direction remain constant, and displacements, strains, and stresses within the plane are uniformly distributed across the thickness of the specimen.Several distinctions between 2D and 3D fracture simulations can be expected.For example, realistic 3D connectivity and tortuosity of pore structure, which significantly influence crack patterns, cannot be accounted for in 2D simulations.Consequently, there is a decrease in fracture energy in 2D compared to 3D simulations.Moreover, Huang et al. [43] found that the 3D model exhibits a higher modulus (approximately 5.6%) than 2D models, which fail to simulate out-of-plane interlock effects of stiff inclusions.In contrast, the simulations conducted by Luković et al. [44] suggest that 2D and 3D simulations of cementitious materials may produce similar average pre-peak behaviors.These minor differences in elastic behaviors between 2D and 3D simulations may support the use of a 2D model to interpret 3D elastic deformation measurements.
Nevertheless, it is always critical and necessary to use the 3D models to reproduce the experimental results, which inherently exist in 3D.The preference for the 2D model in this study primarily stems from its numerical convenience and computational efficiency.Fig. 4 shows the generated lattice models for two types of aggregate surface.
To quantify the surface texture, the arithmetic average roughness (R a ) is calculated for each ITZ model and summarized in Table 4.To evaluate the sole effect of aggregate surface roughness, a homogenized paste using only LD C-S-H phase is generated for the ITZ model, see Fig. 4. Therefore, the variation of overall strength of these ITZ models can be directly attributed to the differences in surface roughness.
For the static fracture simulation, a set of linear elastic analyses is performed by calculating the comparative stress within each element using the following equation: where A is the cross-section area of the beam element; W is the crosssectional moment of resistance; N and M are the local normal force and bending moment, respectively.α N and α M denote the normal force influence factor and the bending influence factor, which are herein chosen as 1.0 and 0.05, respectively [45,46].The choice of these values is to ensure that the failure mode is controlled by the tensile stress [47,48].The stress level in each element is calculated by dividing the comparative stress with the pre-defined tensile strength f t of the element.The beam element with the highest stress-to-strength ratio is identified and removed from the lattice network in a step-wise manner.Details about the numerical algorithms can be found in [45,46].For the fatigue simulation, a cyclic constitutive law that incorporates the fatigue damage evolution is adopted, see Fig. 5.In each analysis step, the maximum fatigue load is imposed on the system.The fatigue damage in each element is calculated based on the local stress and number of  cycles.A fatigue damage index D, which describes the degree of strength reduction with respect to the initial strength, is introduced to consider the fatigue damage evolution.The remaining strength of the damaged element can be calculated as (1-D)f t .When D = 1 the element is assumed to have failed and will be removed from the mesh.It is clear that different phases in the paste-aggregate composite should exhibit different fatigue properties.Therefore, each phase is assigned with a fatigue property following the well-known S-N relationship [49,50], which reads as: where S i is the stress level of the element at the i-th analysis step; N i is the corresponding fatigue life; a and b are two major parameters deciding the local fatigue properties.Since the stiffness matrix is consistently updated due to the accumulated fatigue damage, the internal stress will be redistributed for each loading cycle.Therefore, the fatigue damage D i at the i-th step is calculated according to the Miner's rule [51][52][53][54][55]: where n c is defined as the cycle number for each analysis step; N i is the fatigue life calculated from Eq.( 2) based on the i-th step of stress level and D i-1 is the fatigue damage at previous step.To reduce the computational cost, a method referred to as the "block cycle jump" technique is employed [56].This technique involves bundling a thousand cycles into a single loading cycle, thereby eliminating the need to simulate each individual loading cycle.Consequently, the fatigue damage for 1000 cycles is taken into account in a single analysis step.
If the calculated comparative stress of an element is higher than its remaining strength (i.e.σ sf /(1-D i-1 )f t >1), a post-peak cyclic softening behaviour is assumed for this element, see Fig. 5.In this case, the remaining strength of the element is gradually reduced following a linear softening curve [52,57,58].In total, four softening steps are introduced to represent the relationship between crack opening and bridging stress during the post-peak stage.At each softening step, a certain percentage of reduction d soft in strength is assumed for this element.Therefore, the fatigue damage index of this element at the i-th step should be ).Following the approach described in [29], a constant strain (ε c =-400 µε) at the compression branch is introduced for all phases to determine the reduction percentage of elastic moduli and the residual strains (e.g.ε p1 and ε p2 ).ε 0 is the strain at the original strength f t and ε a is the strain at the reduced strength f a .
When the maximum strain ε max is known, the residual strain ε p1 and the corresponding damaged elastic modulus E d1 of the element can be calculated using the following equations: Note that the UHC and aggregate phase are assumed to have infinite fatigue life without any fatigue damage accumulation.Only the two C-S-  H phases in the cement paste and the corresponding interfacial elements are assumed to experience fatigue damage.After one loading cycle, the system will be updated due to the new stiffness matrix and local strength.If the global stiffness of the system falls below 20% of its original stiffness, the simulation is stopped and the total number of cycles is assumed to be the final fatigue life.Otherwise, another analysis step is executed.

Model calibration
Calibration of the ITZ model is performed by adjusting the parameters based on experimental results [34].Basically, the static mechanical properties of different phases in paste matrix and ITZ are first calibrated based on the authors' previous experimental works [34,[59][60][61].Subsequently, after determining these properties, the fatigue parameters are calibrated based on the fatigue experimental results in [34,61].For local mechanical properties, static flexural bending tests using a KLA Nano indenter G200 are conducted [34].A brief description of the micro-scale experiment is given as follows: firstly, the cantilever ITZ beam is fabricated by using a precision micro-dicing machine (MicroAce Series 3 Dicing Saw) equipped with a resin blade.Two perpendicular cutting directions with the same cutting spacing were applied on the cement paste-aggregate composite sample using the micro-dicing machine, as shown in Fig. 6(a).In this way, multiple rows of cement paste cantilever beams with a square cross section of the beam of 150 × 150 μm 2 can be fabricated.For the testing configuration, one side of the cantilever ITZ beam is fixed on the metal surface and the free end is loaded by the Nano-indenter equipped with a cylindrical wedge indenter tip, as shown in Fig. 6(b).The load is progressively increased in a monotonic or cyclic manner until the failure of beam.The load-displacement responses of

Table 1
Local properties for different interfacial phases [65].Fig. 9. Simulated S-N results using 10 different virtual samples comparison with experimental results in [61].

Table 3
Local fatigue properties of different phases in paste-aggregate composite [65].the beam are recorded and used to calculate the strength and elastic modulus.The experimentally obtained S-N curves are also used to calibrate the fatigue parameters in the model.In the simulation, the same boundary condition as in the experiments is applied on virtual samples, as is illustrated in Fig. 7.The local mechanical properties of three components (i.e.LD C-S-H, HD C-S-H and UHC) derived in previous studies [37] are directly used herein for the hydration phases within the paste.To calibrate the two interfacial phases (aggregate-LD C-S-H and aggregate-HD C-S-H), the ratio of local mechanical properties among different phases is initially established based on the ratio of the measured indentation moduli from the conventional nanoindentation testing [62].Subsequently, the values for two interfacial phases are fine-tuned through fitting with experimental results obtained in previous studies [59].To examine the effect of local properties of interfacial elements on the fracture performance, three sets of parameters as listed in Table 1 are evaluated.Investigation on this aspect is interesting because real concrete presents varying interfacial properties, as seen in lightweight concrete [63] and ultra-high performance concrete [64].Fig. 7 shows the simulated load-displacement curves and the corresponding fracture patterns at certain analysis steps.Both the interfacial strengths assumed in Case 1 and 2 are much lower than those of the hydration products (i.e.52.2 MPa).Therefore, a similar fracture pattern is observed in both cases, where the critical crack is generated along the aggregate surface.The strength assumed in Case 3 is relatively close to the hydration products.A strong interfacial bonding allows the crack to penetrate into the interfacial transition zone rather than along the aggregate surface.As a result, the maximum load calculated in Case 3 simulation is slightly higher than that in Case 2.
By simulating ten virtual samples for each w/c ratio, it is found that the simulated load-displacement curves using the local properties in   Case 2 agree well with the experimental results, see Fig. 8. Therefore, the local elastic moduli E n and tensile strengths f t of all element types can be finalized and are summarized in Table 2.The local mechanical properties will serve as the initial properties in the following fatigue simulations.
Fatigue testing results of ITZ beams reported in [61] are used to calibrate the parameters related to the fatigue damage evolution.The local fatigue properties of main hydration products in harden cement paste have been previously determined by the authors [37].In total, 13 sets of local fatigue properties for two interfacial element types are evaluated.With the changing of local fatigue properties, the slopes of simulated S-N curves vary significantly.The parameter combination, i.e. a(LD)=0.09,a(HD)=0.07,b=1.05, which results in a S-N curve close to experimental results, is adopted in subsequent simulations.Fig. 9 presents the simulated S-N curves of ten virtual samples and the comparison with experimental results.The local fatigue properties for different phases are summarized in Table 3.Assuming elastic-brittle behaviour of hydration products at the microscale, the maximum strain ε max for each phase is assumed to be slightly larger than the value of strength/modulus ratio.It is important to note that a part of large scatter observed in fatigue tests may arise from the stochastic nature of the cementitious materials.It has been shown that the stochastic mechanical properties of cement paste can be predicted by lattice fracture model using microstructures randomly extracted from XCT images [46].By explicitly considering the heterogeneous microstructures, deterministic values of micromechanical parameters are also used here as inputs for simplicity.Fig. 10 illustrates the simulated S-N curves for ten virtual ITZ beams with varying w/c ratios.Owing to the substantial scatter, distinguishing fatigue strengths between the two w/c ratios is challenging.In Fig. 10 (b), when plotting fatigue life against the applied stress magnitude, it becomes evident that a higher stress level is required to induce failure in samples with lower w/c ratios for a given number of cycles.In general, the simulated fatigue life for both w/c ratios fall within a reasonable range compared to experimental results.The simulation attributes fatigue life scatter solely to the microstructure variation, which is a major contributor to the variability observed in fatigue test results.

Uniaxial tensile simulation of ITZ
ITZ microstructures generated in Section 2.1.2are used to perform computational uniaxial tension tests.This step aims to acquire homogenized mechanical and fatigue properties of the ITZ at the microscale.The outcome of these simulations will serve as input for the meso-scale simulation of mortar (Section 3.1).In the uniaxial tensile simulation, one end of the aggregate is fixed and the nodal load is applied on the top of the cement paste specimen.The calibrated local properties of different phases in the ITZ microstructure from Table 2 and Table 3 are used.
The influence of surface roughness on the overall properties is investigated.The simulated stress-strain curves are presented in Fig. 11.For both surface types, the highest uniaxial strength can be found in the ITZ models with the homogenized paste matrix (zero porosity).This is not supervising as the strength typically decreases with the increasing porosity.By comparing the results using two types of aggregate surfaces (see Table 4), it is found that the rough surface seems to exhibit a higher tensile strength, especially when the heterogeneity of ITZ microstructure is not considered.However, when comparing the simulation results from ITZ-R-6 to ITZ-R-10, the relationship between surface roughness and uniaxial tensile strength does not exhibit a notably linear correlation.It should be noted that, on the one hand, irregularities on the rough surface may potentially act as initiation points for cracks, thereby reducing strength.On the other hand, they may also provide more nucleation sites to enhance the hydration and adhesion [66].Moreover, the mechanical interlocking of rough surface also leads to a stronger interfacial bonding strength.
Moreover, Fig. 11 illustrates that a greater fracture energy (quantified as the area beneath the post-peak stress-strain curve) is always achieved in the model featuring the rough surface.This outcome stems from the formation of more tortuous cracking paths along the aggregate surfaces, as evidenced in Fig. 12.In general, the crack initiate at the interface element and propagate along the aggregate surface.Despite the effect of aggregate roughness, the local distribution of pores generated during hydration also plays an important role in fracture process, as depicted in ITZ_R_5 in Fig. 12.The combination of weak bonding strength and increased local porosity in the ITZ contributes to its low tensile strength.
For the fatigue simulations under uniaxial tensile loading, two upper stress levels, namely 70% and 75% of the tensile strength, are applied to the ITZ models, utilizing the calibrated fatigue parameters from Table 3.The lower stress levels are fixed to 0.1% of the strength.The simulated relationships between the uniaxial tensile fatigue life and the stress level along with the simulation results of cement paste are presented in Fig. 13.Apparently, the fatigue strength of cement paste is much higher than that of the ITZ indicated by the number of cycles to failure at a given stress level.In Fig. 13(a), no significant difference between the two types of ITZ models is observed.However, when the results are plotted against the magnitude of stress, see Fig. 13(b), the model with the rough surface has longer fatigue life at a given stress magnitude.This is in accordance with the findings of Lee et al. [67] that the tensile strength and fatigue resistance are positively correlated.

Mesoscopic simulation of Mortar
In this section, a parameter-passing upscaling scheme developed in [68] is adopted to establish a connection between the micro-and meso-scales.This method has demonstrated its efficacy as a reliable tool for bridging these two scales and has been successfully applied in simulating the static fracture behaviors of cementitious materials [33].At the mesoscale, the mortar specimen is assumed to have a dimension of 50 × 50 mm 2 with mesh size of 500 µm.Given that both the side lengths of the cement paste and ITZ virtual specimens are 500 µm in the microscopic model, the length of 500 µm serves as the crucial connection point that links the two scales.The previously simulated mechanical and fatigue properties of ITZ in Section 2.1 and paste matrix in previous work [37] at this length scale are used as local inputs in the mesoscopic simulation.Integrating this information with the mesostructure of the mortar enables simulating the global fracture behaviour of mortar.Additionally, a parametric study is conducted and the simulation results are compared with experimental results in literature.

Two-dimensional virtual mortar specimens
To generate the mesostructure of mortar considering the realistic aggregate shape, the Anm material model proposed by Qian in [68] is used.The Anm material model is programmed on the basis of the 3D mathematical analysis of particle shape as described in [69].In this model, irregular shape aggregate particles described by the spherical harmonic expansion are placed one after another into a pre-defined domain following a parking algorithm.The aggregates are selected from a dataset according to the given size distribution of particles in terms of sieve range.More details about this model are available in [68].Herein, the volumetric fraction of sand is set to be 30.0%and the size of individual sand is randomly selected from 5 mm to 10 mm.Furthermore, spherical air voids are also considered in a similar way.The air content is set as 3.5% and the diameter range of pores is from 1 mm to 12 mm.After packing the sands and the air voids, the vector-based composite geometrical structure is obtained and digitalized into the digital specimens with a resolution of 500 μm/voxel, resulting in a 100 × 100 ×100 voxels mesostructure.The generated 3D digital mortar sample is shown in Fig. 14.Afterwards, 2D slices are extracted and converted to the lattice model containing four element types: paste matrix, sand, ITZ and air voids.It is crucial to recognize the inherent limitations of 2D simulations, particularly in their inability to account for realistic 3D microstructural features.Adopting 3D simulations has the potential to provide more accurate representations of crack patterns, fracture energy, and fatigue damage evolution.

Constitutive law for mortar simulation
The constitutive stress-strain curves for both cement paste and ITZ obtained in microscopic simulations should be properly incorporated in the mesoscopic model.Therefore, a so-called step-wise approach as described by Qian in [68] is used to formulate the constitutive relations by approximating the original non-linear stress-strain curve with a multi-linear curve.Fig. 15 shows the generated multi-linear curves for both cement paste and ITZ.For each component, five points are selected from the stress-strain diagram to form a multi-linear curve.These points should be chosen in a way that fully represents the major behaviour of each component in terms of elastic modulus, tensile strength and fracture energy.With respect to the cement paste element, its constitutive relation is taken from the authors' previous work [37] in which uniaxial tension simulations were conducted on the 500 × 500 μm 2 cement paste specimens.The previously simulated stress-strain curve on the ITZ virtual specimen with the same size is also approximated using the multi-linear curve.To investigate the effect of ITZ properties on the global performance of mortar, two different surface roughness, i.e.ITZ_F_4 and ITZ_R_4, are used as inputs for mesoscale modelling.The choice of the properties of ITZ_F_4 sample is because that it has the lowest strength in flat ITZ samples and can be considered as the lower bound.Accordingly, the ITZ_R_4 is the higher bound, which is the strongest in rough ITZ samples.There is an apparent distinction in their stress-strain curves, as is shown in Fig. 15, and they can be labelled as 'weak' and 'strong' ITZ.The sand in the model is assumed to experience no degradation subjected to the static and fatigue loadings.Therefore, a linear elastic brittle constitutive law is used for the sand.Its elastic modulus is taken from the literature [17], i.e. 70 GPa.Its tensile strength is assumed to be 1/1000 of the modulus.For the fatigue simulation, the post-peak cyclic behaviour is also described by the chosen points, which are introduced in cyclic constitutive laws using Eqs.( 4)- (5).In addition, the simulated S-N relationships for cement paste and ITZ are used as parameters that control the local fatigue damage evolution in mesoscopic modelling.It is generally acknowledged that there is no fatigue limit for cementitious materials as very subtle and minor fatigue damage (e.g.nano-scale crack) is always accumulated during each loading cycle.However, the fatigue damage is negligible when the local stress is small.The calculation of fatigue damage in every cement paste element seems to be unnecessary.To accelerate this computationally expensive operation, the fatigue damage in one paste element is accumulated only when its stress level exceeds 10% of the remaining strength.Similar to the static fracture simulation, two interfacial element types with different mechanical and fatigue properties are implemented in fatigue simulations.The mechanical and fatigue properties of components in the mortar simulation are summarised in Table 5.The fatigue parameter a is a decisive factor in determining the damage evolution and a higher value of this parameter generally indicates a lower fatigue resistance.In this case, it is evident that the cement paste displays superior fatigue resistance compared to the ITZ in terms of parameter a.Note that the overall fatigue performance depends on both the stress amplitude and the remaining strength of the material.Given that stress concentrations tend to be in proximity to the ITZ in mortar or concrete composites, high stress levels are anticipated in these areas.Consequently, the relatively weak ITZ in concrete is prone to fail under fatigue loading.

Uniaxial tensile simulation of mortar
A total of ten virtual mortar samples are modelled to determine their mechanical properties.Static uniaxial tension tests are performed by applying unit nodal load at one end and fixing the other.The simulated stress-strain curves are given in Fig. 16.The average strength and elastic modulus for mortar samples with weak ITZ are 2.30 ± 0.12 MPa and 19.76 ± 0.92 GPa, while the average results for strong ITZ are 6.02 ± 0.57 MPa and 21.84 ± 1.16 GPa.It is clear that the overall strength is increased due to the stronger ITZ.The simulated mechanical properties are found to be within a reasonable range when compared to the experimental results reported in authors' previous work [41], where the mortars with similar volume fraction of aggregate but lower w/c ratio (0.3) were tested.The experimentally obtained strength and elastic modulus are 5.74 MPa and 26.36 GPa, respectively.It is important to note that the 2D simulations do not consider the realistic 3D connectivity and tortuosity of the pore structure, potentially resulting in a lower fracture energy and possibly higher strength.The preference for the 2D model in this study is due to its computational efficiency.
In fatigue simulations, the applied upper stress levels are 60%, 70% and 80% of corresponding static strength, while the lower stress levels are fixed to 0.1% of the strength.Two types of ITZ elements are used in fatigue simulations of mortar samples.The simulated S-N results are shown in Fig. 17, indicating that the stronger ITZ also results in a higher fatigue resistance.It is clear that the simulated fatigue properties of mortar are much lower than those of cement paste at the same stress level due to the inclusion of ITZ.To validate the numerical results, experimental data on the tensile fatigue of mortar at a similar length scale are normally needed.However, limited experimental investigations on the tensile fatigue of mortar are available in the literature.Alternatively, a comparison is made with existing experimental results focusing on the flexural fatigue of mortar or concrete to evaluate the numerical results.The collected S-N data from [70][71][72][73][74] are presented in Fig. 17    ranging from 0.4 to 0.5, and are cured for a duration of 28 days.The simulated fatigue life of numerical mortar based on inputs from micromechanical models lies within a reasonable range compared to the experimental results.This provides positive affirmation for the multiscale modelling procedures utilized herein, as well as the validity of the fatigue parameters applied in this model.
As an important indicator of fatigue damage, the simulated variation of the stiffness during the fatigue loading is plotted in Fig. 18.The experimental results from fatigue tests conducted on mortar and concrete, as reported in the literature [75][76][77][78], are also provided for comparison.Two distinct stages can be discerned: the initial stage, encompassing nearly 80% of the fatigue life, exhibits a gradual reduction in stiffness.The second stage is marked by a sudden and pronounced decrease in stiffness.This phenomenon occurs as accumulated damage reaches a critical level, leading to the initiation of unstable cracks before eventual failure.Similar experimental observations supporting these stages have been reported in [34,79].
In general, various damage indicators, such as the secant stiffness and residual strain of the sample, are employed in most experimental investigations to indirectly assess fatigue damage evolution.One advantage of the current modeling scheme is its capability to visualize the entire process of damage evolution by mapping the variation of the damage index D in each element.This allows for a direct analysis of fatigue damage development within the mortar sample.It is important to note that fatigue damage primarily relies on the local stress field and the element type.Meanwhile, the stress field will be updated for each cycle due to the degradation of stiffness.It is plausible that a severely damaged element may experience less damage in the subsequent cycle due to stress redistribution.This phenomenon can influence the fatigue cracking path and the ultimate fatigue life.Therefore, the fatigue fracture process is more complicated compared to the static fracture, as different degrees of degradation accumulate simultaneously in the sample.
Fig. 19 shows the fatigue damage contour map at different number of cycles for a mortar sample loaded at the stress level of 60%.The results using two ITZ element types are presented.For the mortar sample with the weak ITZ in Fig. 19(a), damage is mainly accumulated in the ITZ at 10% of fatigue life N f , and some minor damages are also observed in the paste matrix near the air void.With increasing loading cycles, fatigue   damage is generated much faster in ITZ elements than in the paste matrix.Some of ITZ elements are already broken at around half the fatigue life as their damage indices reaches 1.Before the end of the simulation, the damage development in the ITZ continues and multiple fatigue cracks surrounding the aggregates are gradually formed.However, only mild damages is detected in most cement paste elements even when the fatigue failure criterion is fulfilled.For the case of mortar with the strong ITZ, when the simulation is at the initial stage (25% of N f ), most cement paste elements have already been damaged, which is different from the sample with the weak ITZ.This means that the paste matrix may suffer more damage due to the stronger ITZ.However, when the major fatigue crack localizes, the damage in cement paste barely increases with the loading cycles.Another significant difference between the results of the mortar with two ITZ types is the distribution and the number of fatigue cracks.The current mortar sample with the weak ITZ tends to generate more fatigue cracks on the sites of ITZ, while much less interfacial cracks are detected in the strong ITZ sample.Moreover, the main reason for the fatigue failure of the weak ITZ sample is the massive breakdowns of ITZ elements, but the critical crack formed in strong ITZ sample is due to the coalescence of cracks in both the ITZ and the paste matrix.
To gain more insights into the fatigue damage evolution, the contour maps of fatigue damage at the final fatigue life under different stress levels are analysed.The results of two virtual mortar samples, i.e., Mortar_R_4 and Mortar_F_7, are presented in Fig. 20.When the strong ITZ is applied in the mortar simulation, the distributions of final cracks are slightly different for these two stress levels, as is shown in Fig. 20(a).Nevertheless, fatigue damages generated in the paste matrix under two stress levels, in terms of number of damaged elements and damage degree (D), are very similar.For instance, around 52.0% and 56.9% of cement paste elements have been damaged to a certain degree (i.e., 0 < D ≤ 0.2) for the two stress levels, respectively.In the case of Mortar_F_7 with the weak ITZ, the major difference is the fatigue damage evolution in paste matrix.When the applied stress level is 80%, there is almost no damage developed in paste matrix.The mortar fails immediately after the formation of the critical crack in the ITZ.This phenomenon is similar to the crack localization in static fracture.However, when the stress level is low, some damages do accumulate in the paste matrix and the fatigue life is considerably decreased.
In the current example, when the applied stress is the same and the difference in fatigue properties is small, the 244% increase of ITZ strength leads to an increase of the fatigue life by approximately one to two orders of magnitude.It can be concluded that the mechanical properties of ITZ can significantly affect the fatigue properties and eventually the final fatigue life of mortar by means of altering the damage evolution process.Improvement of the ITZ is essential for increasing the fatigue resistance of cementitious composites.

Conclusions
In this paper, a 2D lattice model is used to investigate the fracture behaviour microscopic ITZ and mesoscopic mortar under static and cyclic loadings.The influence of w/c ratio, aggregate surface roughness on the fatigue behaviour of ITZ was studied.Additionally, the effect of ITZ properties and stress levels on the fatigue behaviour of mortar was investigated.Based on the presented results, the following conclusions can be drawn: (1) The microscopic model of the ITZ, derived directly from the microscale, is established and verified using experimental data obtained at the microscale, which is subsequently integrated into the mesoscopic mortar model.It was found that the uncoupled upscaling approach is a suitable tool to bridge the two scales in a multiscale modelling scheme.The simulated mechanical properties and S-N curves of mortar lie within a reasonable range compared to experimental results found in literature.(2) An increase in aggregate surface roughness generally leads to a higher tensile strength, fracture energy and fatigue resistance of ITZ at the microscale.However, the relationship between surface roughness and mechanical properties does not exhibit a notably linear correlation.The overall mechanical properties depend on the local microstructural properties, such as porosity and stacking of hydration products, as well as the mechanical interlocking around the aggregate surface.
Two types of models are generated: the first type of model is the ITZ beam model with the size of 150 × 750 μm 2 , employed for the calibration and validation of local mechanical and fatigue properties.The second is the ITZ model with the size of 500 × 500 μm 2 containing the ITZ and aggregate, which is used to predict the uniaxial tensile behaviour of ITZ.Given the assumption that the linear elastic aggregate element remains intact in the model, the portion of the ITZ model containing only the aggregate is omitted.Subsequently, the model is resized to 500 × 250 μm 2 , aiming to enhance computational efficiency.

Fig. 4 .
Fig. 4. Generation lattice model based on 2D images of ITZ microstuctures with different surface textures (black pixel represents UHC, orange represents aggregate, yellow represents LD C-S-H and red represents HD C-S-H).

Fig. 7 .
Fig. 7. Comparison of simulated load-displacement curves and fracture patterns of ITZ beams with the w/c ratio of 0.3 using three different sets of local mechanical properties.

Fig.
Fig. The simulated results of (a) stress level-fatigue life curves and (b) applied stress magnitude-fatigue life curves.

Fig. 11 .
Fig. 11.The simulated stress-strain diagrams of ITZ models with (a) flat aggregate surfaces and (b) natural rough surfaces under uniaxial tension.

Fig. 13 .
Fig. 13.(a) The comparison of S-N curves between cement paste and ITZ; (b) The relationship between the applied stress and number of cycles to failure.
. The collected data mostly have similar mix compositions comprising cement mortar and normal concrete, with a w/c ratio

Fig. 14 .
Fig. 14.The 3D digital mortar sample and the extracted 2D slice (red represents sand, grey represents paste matrix and blue is pore).

Fig. 15 .
Fig.15.Approximations of the non-linear stress-strain responses of paste matrix and ITZ using multi-linear curves.

Fig. 17 .
Fig. 17.The comparison of and numerical S-N curves.

Fig. 18 .
Fig. 18.The between simulated stiffness degradation of mesoscopic mortar model and experimental results.

Fig. 19 .
Fig. 19.The contour maps of fatigue damage at different stages of cyclic loading at the stress level of 60% for virtual mortar samples with (a) the weak ITZ (Mortar_F_2, N f =30486) and (b) the strong ITZ (Mortar_R_2, N f = 660866).

( 3 )
The mechanical properties of the ITZ significantly affect the fracture behaviours of mortar subjected to both static and fatigue loading.Typically, a weak ITZ leads to localized damage during fatigue loading, whereas a strong ITZ promotes a more homogenized mortar, resulting in a more evenly distributed damage pattern.Variations in stress levels also influence the fatigue resistance of mortar by modifying the fatigue damage evolution within the mesostructure.The simulation model presented in this paper is a first attempt to develop a lattice fatigue model at the micro-and mesoscale, and is considered an essential starting point for predicting the fatigue performance of concrete structures at the macroscale.Although the current 2D model offers the opportunity to simulate the fatigue behaviour of cementitious material in a multiscale framework, 3D simulation is still required to consider more realistic material structures.Nevertheless, the current model has successfully demonstrated the effect of heterogeneities at different scales on the fatigue performance of cementitious material in an efficient way and provided some insights into the fatigue damage evolution and fatigue fracture phenomena.CRediT authorship contribution statementBranko Šavija: Writingreview & editing, Supervision, Funding acquisition.Klaas van Breugel: Writingreview & editing, Supervision, Funding acquisition.Erik Schlangen: Writingreview & editing, Supervision, Funding acquisition, Conceptualization.Minfei Liang: Writingreview & editing, Software, Investigation, Data curation, Conceptualization.Yidong Gan: Writingreview & editing, Writingoriginal draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization.

Table 4
The summary of simulation results for different ITZ models.

Table 5
Local fatigue properties of different phases in mortar.