Micromechanical Numerical Modelling on Compressive Failure of Recycled Concrete using Discrete Element Method (DEM)

This paper investigates the failure processes of recycled aggregate concrete by a model test and numerical simulations. A micromechanical numerical modeling approach to simulate the progressive cracking behavior of the modeled recycled aggregate concrete, considering its actual meso-structures, is established based on the discrete element method (DEM). The determination procedure of contact microparameters is analyzed, and a series of microscopic contact parameters for different components of modeled recycled aggregate concrete (MRAC) is calibrated using nanoindentation test results. The complete stress–strain curves, cracking process, and failure pattern of the numerical model are verified by the experimental results, proving their accuracy and validation. The initiation, growth, interaction, coalescence of microcracks, and subsequent macroscopic failure of the MRAC specimen are captured through DEM numerical simulations and compared with digital image correlation (DIC) results. The typical cracking modes controlled by meso-structures of MRAC are concluded according to numerical observations. A parameter study indicates the dominant influence of the macroscopic mechanical behaviors from the shear strength of the interfacial transition zones (ITZs).


Introduction
Concrete materials show complex mechanical behaviors because of their internal structures induced by the distribution of different components and weak interfaces at the mesoscopic level [1][2][3]. Recycled aggregate concretes contain various components including natural aggregates, new or old cement mortars, and numerous interfacial transition zones, which extremely increase the complexity of the meso-structure. The meso-structure of recycled aggregate concrete significantly influences the microcrack initiation, propagation and, consequently, macroscopic failure.
A lot of studies have been conducted to investigate the failure behaviors of recycled aggregate concretes using model tests. Buyukozturk et al. [4] firstly proposed an imitated model for plain concrete specimens composed of coarse aggregates and cement mortar to investigate the deformation and failure behaviors. Choi and Shah [5] used modeled concretes to study how the aggregate spacing influences the fracture process of concretes containing different amounts of embedded aggregates. Li et al. [6] used the digital image correlation method to track the cracking process of modeled recycled

Nanoindentation Test
It is well recognized that the cracking and failure process of concretes is strongly influenced by the interfacial transition zones (ITZ). Therefore, the mechanical and microstructural properties of ITZs in RAC were investigated by using nanoindentation tests and scanning electron microscopy (SEM). Two indentation areas with a size of 150 μm × 100 μm were located within the old ITZ and new ITZ region, as shown in Figure 2. The contour maps of the indentation modulus of old and new ITZs are obtained from nanoindentation tests. The darker areas mean a lower modulus region, while the brighter areas mean higher modulus regions. Figure 2 clearly shows the modulus difference in different portions, which indicates that the microstructure in the ITZ region is quite different from that of the mortar matrix. Figure 2a shows that the indentation modulus of old ITZs is around 70%~80% of the old mortar matrix, and Figure 2b shows that the indentation modulus of new ITZs is around 80%~90% of the new mortar matrix. The thickness of ITZs is around 40~70 μm. There are many more voids that can be observed in the ITZ region from the SEM images (Figure 3), which leads to the lower stiffness and strength of the ITZs. The indentation modulus distributions of different components of the MRAC specimen provide the reference for the quantitative analysis of the microscopic parameters in numerical modeling.

Nanoindentation Test
It is well recognized that the cracking and failure process of concretes is strongly influenced by the interfacial transition zones (ITZ). Therefore, the mechanical and microstructural properties of ITZs in RAC were investigated by using nanoindentation tests and scanning electron microscopy (SEM). Two indentation areas with a size of 150 µm × 100 µm were located within the old ITZ and new ITZ region, as shown in Figure 2. The contour maps of the indentation modulus of old and new ITZs are obtained from nanoindentation tests. The darker areas mean a lower modulus region, while the brighter areas mean higher modulus regions. Figure 2 clearly shows the modulus difference in different portions, which indicates that the microstructure in the ITZ region is quite different from that of the mortar matrix. Figure 2a shows that the indentation modulus of old ITZs is around 70%~80% of the old mortar matrix, and Figure 2b shows that the indentation modulus of new ITZs is around 80%~90% of the new mortar matrix. The thickness of ITZs is around 40~70 µm. There are many more voids that can be observed in the ITZ region from the SEM images (Figure 3), which leads to the lower stiffness and strength of the ITZs. The indentation modulus distributions of different components of the MRAC specimen provide the reference for the quantitative analysis of the microscopic parameters in numerical modeling.

Nanoindentation Test
It is well recognized that the cracking and failure process of concretes is strongly influenced by the interfacial transition zones (ITZ). Therefore, the mechanical and microstructural properties of ITZs in RAC were investigated by using nanoindentation tests and scanning electron microscopy (SEM). Two indentation areas with a size of 150 μm × 100 μm were located within the old ITZ and new ITZ region, as shown in Figure 2. The contour maps of the indentation modulus of old and new ITZs are obtained from nanoindentation tests. The darker areas mean a lower modulus region, while the brighter areas mean higher modulus regions. Figure 2 clearly shows the modulus difference in different portions, which indicates that the microstructure in the ITZ region is quite different from that of the mortar matrix. Figure 2a shows that the indentation modulus of old ITZs is around 70%~80% of the old mortar matrix, and Figure 2b shows that the indentation modulus of new ITZs is around 80%~90% of the new mortar matrix. The thickness of ITZs is around 40~70 μm. There are many more voids that can be observed in the ITZ region from the SEM images ( Figure 3), which leads to the lower stiffness and strength of the ITZs. The indentation modulus distributions of different components of the MRAC specimen provide the reference for the quantitative analysis of the microscopic parameters in numerical modeling.

Uniaxial Compression Test and DIC Measurement
The MRAC and CM specimens were uniaxially compressed under a minimal loading rate by the serve-control machine. During the complete loading process, the cracking process of the specimen was captured by digital image correlation (DIC) measurement [8,25] (Figure 4). The microcracks can be recognized according to the image analysis on the displacement or strain distribution images taken from the specimen. The axial strain of the specimen was measured by both the DIC method and conventional strain gages, which are compared and shown in the complete stress-strain curves in Figure 4b. The other DIC measurement results including displacement and strain distribution contours and failure patterns are compared with the numerical results in the following section.

Uniaxial Compression Test and DIC Measurement
The MRAC and CM specimens were uniaxially compressed under a minimal loading rate by the serve-control machine. During the complete loading process, the cracking process of the specimen was captured by digital image correlation (DIC) measurement [8,25] (Figure 4). The microcracks can be recognized according to the image analysis on the displacement or strain distribution images taken from the specimen. The axial strain of the specimen was measured by both the DIC method and conventional strain gages, which are compared and shown in the complete stress-strain curves in Figure 4b. The other DIC measurement results including displacement and strain distribution contours and failure patterns are compared with the numerical results in the following section.

DEM and UDEC Program
Cundal and Strack [18] introduced the original idea of Discrete Element Method (DEM) to calculate the motion of a block system by the explicit time-marching scheme. Since then, DEM approaches have been widely adopted to simulate fracture processes in mechanical rock problems [26]. In contrast to continuum-based numerical models, which need hypothetical constitutive relations to simulate material behaviors, DEM models use discrete particles or blocks represent actual physical objects like grains in granular rocks or aggregates in concretes. All the discontinuities inside the materials, such as fractures, cracks, and pores, are inherent parts of the DEM model. Therefore, DEM models are particularly favorable for investigating the cracking processes of brittle heterogeneous materials.
The universal distinct element code (UDEC) program is a 2D numerical code under the framework of the DEM concept. UDEC can simulate either static or dynamic responses of a discontinuous system composed of either rigid or deformable discrete blocks. As shown in Figure 5, a plane concrete specimen can be modeled as the assembly of discrete blocks with random sizes and shapes in UDEC. The Voronoi tessellation procedure is introduced to subdivide the concrete specimen into randomly sized polygons. The Voronoi blocks can be either rigid or subdivided into elastic elements with Young's modulus and Poisson's ratio . The contacts between the blocks are deformable and breakable, which is represented by the Coulomb slip model with tensile strength, as shown in Figure 5. The following contact parameters are involved to calculate the force-displacement relationship of a contact: normal stiffness ( ), shear stiffness ( ), friction angle ( ), cohesion ( ), and tensile strength ( ). A microcrack is generated when the contact breakage

DEM and UDEC Program
Cundal and Strack [18] introduced the original idea of Discrete Element Method (DEM) to calculate the motion of a block system by the explicit time-marching scheme. Since then, DEM approaches have been widely adopted to simulate fracture processes in mechanical rock problems [26]. In contrast to continuum-based numerical models, which need hypothetical constitutive relations to simulate material behaviors, DEM models use discrete particles or blocks represent actual physical objects like grains in granular rocks or aggregates in concretes. All the discontinuities inside the materials, such as fractures, cracks, and pores, are inherent parts of the DEM model. Therefore, DEM models are particularly favorable for investigating the cracking processes of brittle heterogeneous materials.
The universal distinct element code (UDEC) program is a 2D numerical code under the framework of the DEM concept. UDEC can simulate either static or dynamic responses of a discontinuous system composed of either rigid or deformable discrete blocks. As shown in Figure 5, a plane concrete specimen can be modeled as the assembly of discrete blocks with random sizes and shapes in UDEC. The Voronoi tessellation procedure is introduced to subdivide the concrete specimen into randomly sized polygons. The Voronoi blocks can be either rigid or subdivided into elastic elements with Young's modulus E g and Poisson's ratio ν g . The contacts between the blocks are deformable and breakable, which is represented by the Coulomb slip model with tensile strength, as shown in Figure 5. The following contact parameters are involved to calculate the force-displacement relationship of a contact: normal stiffness (k n ), shear stiffness (k s ), friction angle (φ c ), cohesion (c c ), and tensile strength (σ c t ). A microcrack is generated when the contact breakage happens by shear or tension. Therefore, the block boundaries provide numerous potential internal paths for the growth of microcracks. Macroscopic fractures can eventually be formed due to the initiation, propagation, and connection of microcracks. Tan et al. [24] presented a numerical analysis of the failure process of modeled recycled aggregate concrete using the UDEC program, whose results showed very similar observations as the experimental results.
Materials 2020, 13, x FOR PEER REVIEW 5 of 18 happens by shear or tension. Therefore, the block boundaries provide numerous potential internal paths for the growth of microcracks. Macroscopic fractures can eventually be formed due to the initiation, propagation, and connection of microcracks. Tan et al. [24] presented a numerical analysis of the failure process of modeled recycled aggregate concrete using the UDEC program, whose results showed very similar observations as the experimental results.

Image Based Meso-Structures Model and Microscopic Contact Parameter Determination
The UDEC program was adopted to simulate the MRAC specimen and analyze the cracking evolvement and progressive failure process of the specimen, considering meso-structures. As described in the last section, the MRAC specimen was simulated by discrete Voronoi blocks. The actual MRAC specimen consists of several different components at the meso-level scale that may influence its overall mechanical behaviors, including natural aggregates, new and old cement mortar, and new and old ITZs. In particular, the ITZs with the lowest physical strength may be the critical factor in controlling the microcracking process of the MRAC specimen. To model the meso-structures of the MRAC specimen, Digital Image Processing (DIP) is introduced, as illustrated in Figure 6. Binarization processing was used to transform the digital image of the MRAC specimen into a black and white image. The contrast between the old mortar matrix and the natural aggregate can be easily observed from the black and white image. The boundaries of each natural aggregate can then be extracted as polylines from the image by the outline identification technique. The concrete specimen was subdivided into abundant convex-shaped blocks of random size by the Voronoi tessellation procedure. The Voronoi blocks in the numerical specimen were then divided into three groups according to the boundaries recognized from the digital image, which are the new mortar matrix, the old mortar matrix, and the natural aggregate. The boundaries between different groups are therefore the ITZs.

Image Based Meso-Structures Model and Microscopic Contact Parameter Determination
The UDEC program was adopted to simulate the MRAC specimen and analyze the cracking evolvement and progressive failure process of the specimen, considering meso-structures. As described in the last section, the MRAC specimen was simulated by discrete Voronoi blocks. The actual MRAC specimen consists of several different components at the meso-level scale that may influence its overall mechanical behaviors, including natural aggregates, new and old cement mortar, and new and old ITZs. In particular, the ITZs with the lowest physical strength may be the critical factor in controlling the microcracking process of the MRAC specimen. To model the meso-structures of the MRAC specimen, Digital Image Processing (DIP) is introduced, as illustrated in Figure 6. Binarization processing was used to transform the digital image of the MRAC specimen into a black and white image. The contrast between the old mortar matrix and the natural aggregate can be easily observed from the black and white image. The boundaries of each natural aggregate can then be extracted as polylines from the image by the outline identification technique. The concrete specimen was subdivided into abundant convex-shaped blocks of random size by the Voronoi tessellation procedure. The Voronoi blocks in the numerical specimen were then divided into three groups according to the boundaries recognized from the digital image, which are the new mortar matrix, the old mortar matrix, and the natural aggregate. The boundaries between different groups are therefore the ITZs.
The ITZs were modeled in the DEM model as zero-thickness contacts and endowed with corresponding contact parameters according to the coordinates of the extracted polylines, while, in the continuum numerical models (Anuruddha et al. 2018 [13]), the thickness of ITZs was considered from 0.5~2.0 mm, which is obviously much higher than the thickness of 40~60 µm measured from nanoindentation tests [27,28]. The old mortar and the aggregate portions can also be distinguished depending on which one is internal and which is external of the ITZ polylines. Rigid blocks are assumed in the proposed model. Thus, the mechanical behaviors of the MRAC specimen are completely controlled by its microscopic contact parameters. Three different groups of contact parameters are needed to model the MRAC specimen, which are natural aggregates, cement mortar, and ITZs, respectively. To simplify the parameter calibration procedure, the mortar and ITZs were not distinguished into old or new groups.  The ITZs were modeled in the DEM model as zero-thickness contacts and endowed with corresponding contact parameters according to the coordinates of the extracted polylines, while, in the continuum numerical models (Anuruddha et al. 2018 [13]), the thickness of ITZs was considered from 0.5~2.0 mm, which is obviously much higher than the thickness of 40~60 μm measured from nanoindentation tests [27,28]. The old mortar and the aggregate portions can also be distinguished depending on which one is internal and which is external of the ITZ polylines. Rigid blocks are assumed in the proposed model. Thus, the mechanical behaviors of the MRAC specimen are completely controlled by its microscopic contact parameters. Three different groups of contact parameters are needed to model the MRAC specimen, which are natural aggregates, cement mortar, and ITZs, respectively. To simplify the parameter calibration procedure, the mortar and ITZs were not distinguished into old or new groups.
The contact parameters for different components in the MRAC specimen must be determined firstly on the condition that the macroscopic mechanical behaviors of the simulated material such as deformability, uniaxial compression strength, and failure pattern can be generally reproduced. Nevertheless, the microscopic contact parameters are diverse and different from the regular parameters obtained from experiments [21,22]. A determination procedure has to be performed in advance to determine proper contact parameters according to the macroscopic mechanical behaviors of the DEM model.
There is only one group of contact parameters needed for the CM specimen model; therefore, the contact parameters of mortar can be determined firstly by comparison with the uniaxial compression results of the CM specimen. Tan et al. [29] have provided empirical relations to calculate the contact stiffnesses and based on the elastic modulus and Poisson's ratio of the Voronoi-based DEM model. The contact parameters ( , and ) of the CM specimen can be predefined according to their corresponding macroscopic strength parameters. The peak strength and corresponding axial strain in the compressive stress-strain curve of the CM specimen ( Figure 8) are the target parameters to calibrate the contact parameters Then, the contact parameters can be further modified according to several trials and errors. The determination flow diagram is illustrated in Figure 7. The complete stress-strain curves and failure patterns of the numerical CM specimen with calibrated contact parameters are shown in  The contact parameters for different components in the MRAC specimen must be determined firstly on the condition that the macroscopic mechanical behaviors of the simulated material such as deformability, uniaxial compression strength, and failure pattern can be generally reproduced. Nevertheless, the microscopic contact parameters are diverse and different from the regular parameters obtained from experiments [21,22]. A determination procedure has to be performed in advance to determine proper contact parameters according to the macroscopic mechanical behaviors of the DEM model.
There is only one group of contact parameters needed for the CM specimen model; therefore, the contact parameters of mortar can be determined firstly by comparison with the uniaxial compression results of the CM specimen. Tan et al. [29] have provided empirical relations to calculate the contact stiffnesses k n and k s based on the elastic modulus and Poisson's ratio of the Voronoi-based DEM model. The contact parameters (φ c , c c and σ c t ) of the CM specimen can be predefined according to their corresponding macroscopic strength parameters. The peak strength and corresponding axial strain in the compressive stress-strain curve of the CM specimen ( Figure 8) are the target parameters to calibrate the contact parameters.
Then, the contact parameters can be further modified according to several trials and errors. The determination flow diagram is illustrated in Figure 7. The complete stress-strain curves and failure patterns of the numerical CM specimen with calibrated contact parameters are shown in Figures 8 and 9 in comparison with the experimental results.
The numerical result in Figure 9a shows more cracks and main fractures compared to the experimental results in Figure 9b,c. The main reason for these discrepancies is that the numerical CM specimens are almost homogeneous, which leads to a multi-fracture pattern at the failure state. In reality, the location where the main fracture appears depends on the material flaws inside the CM specimen. Moreover, the numerical model ( Figure 9a) is able to show much more details compared to the experimental observations, including micro-cracks with an extremely small crack width. However, the real photo in Figure 9c can only show the main fractures with relatively wider apertures. The DIC picture in (Figure 9b) can show the micro-cracks to a certain degree, but not precisely. Nevertheless, the typical cracking behaviors are actually captured by the numerical model.       The numerical result in Figure 9a shows more cracks and main fractures compared to the experimental results in Figure 9b,c. The main reason for these discrepancies is that the numerical CM specimens are almost homogeneous, which leads to a multi-fracture pattern at the failure state. In reality, the location where the main fracture appears depends on the material flaws inside the CM specimen. Moreover, the numerical model ( Figure 9a) is able to show much more details compared to the experimental observations, including micro-cracks with an extremely small crack width. However, the real photo in Figure 9c can only show the main fractures with relatively wider apertures. The DIC picture in (Figure 9b) can show the micro-cracks to a certain degree, but not precisely. Nevertheless, the typical cracking behaviors are actually captured by the numerical model.
The contact parameters of other components can then be predetermined by taking the contact parameters of the CM specimen as a standard. The microhardness obtained from the nanoindentation test provides a quantitative scale between the contact stiffness and strength for different components in concrete [30,31]. Still, trials and errors are always necessary before optimum contact parameters are obtained. An entire set of contact parameters for the MRAC specimen has been determined and is listed in Table 1. The loading plate was modeled using the linear elastic constitutive law. A sufficient number of discrete blocks should be contained in a DEM model to avoid significant uncertainty in the macroscopic mechanical behaviors. The average block edge length is set to 1.0 mm, as the DEM specimen contains 12,859 blocks with an average block area 1.8 × 10 −6 m 2 and 62,500 contacts. The axial stress acting on the specimen was controlled by moving the two loading plates in opposite directions. The loading velocity was kept at a small value so that the cracking process could be well recorded during the complete loading process.

Stress-Strain Curve and Contact State Evolvement
The complete stress-strain curve of the MRAC specimen predicted by the numerical simulation is shown in Figure 10, which is in good agreement with the laboratory data. In particular, in the pre-peak loading stages, a very similar curve shape and almost the same macroscopic strength can The contact parameters of other components can then be predetermined by taking the contact parameters of the CM specimen as a standard. The microhardness obtained from the nanoindentation test provides a quantitative scale between the contact stiffness and strength for different components in concrete [30,31]. Still, trials and errors are always necessary before optimum contact parameters are obtained. An entire set of contact parameters for the MRAC specimen has been determined and is listed in Table 1. The loading plate was modeled using the linear elastic constitutive law. A sufficient number of discrete blocks should be contained in a DEM model to avoid significant uncertainty in the macroscopic mechanical behaviors. The average block edge length is set to 1.0 mm, as the DEM specimen contains 12,859 blocks with an average block area 1.8 × 10 −6 m 2 and 62,500 contacts. The axial stress acting on the specimen was controlled by moving the two loading plates in opposite directions. The loading velocity was kept at a small value so that the cracking process could be well recorded during the complete loading process.

Stress-Strain Curve and Contact State Evolvement
The complete stress-strain curve of the MRAC specimen predicted by the numerical simulation is shown in Figure 10, which is in good agreement with the laboratory data. In particular, in the pre-peak loading stages, a very similar curve shape and almost the same macroscopic strength can be observed. The numbers of contacts in tensile opening and shear sliding during the complete loading process of the simulation are shown in Figure 11, which represents the microcrack evolution inside the MRAC specimen. Four typical stages in the cracking process are labeled with Roman numerals in Figure 10  be observed. The numbers of contacts in tensile opening and shear sliding during the complete loading process of the simulation are shown in Figure 11, which represents the microcrack evolution inside the MRAC specimen. Four typical stages in the cracking process are labeled with Roman numerals in Figure 10

The Complete Failure Processes
The deformation field and crack distribution of the MRAC specimen from both DIC measurements and numerical simulations at the selected loading stages (marked in Figure 10) are compared to each other. Figures 12 and 13 show the macroscopic deformation field of the MRAC specimen at different loading stages ( ). Figure 12 shows the lateral displacement distribution of the MRAC specimen obtained by DIC measurements. Several macroscopic fractures were observed for MRAC specimens after failure according to the discontinuity in the lateral displacement contour. The fracture directions are almost vertical. Figure 13 shows the displacement vectors of the block corners of the MRAC specimen obtained by numerical simulations. Similar macroscopic fractures of the MRAC specimen at stage are marked in Figure 13 according to the discontinuity in the displacement. loading process of the simulation are shown in Figure 11, which represents the microcrack evolution inside the MRAC specimen. Four typical stages in the cracking process are labeled with Roman numerals in Figure 10

The Complete Failure Processes
The deformation field and crack distribution of the MRAC specimen from both DIC measurements and numerical simulations at the selected loading stages (marked in Figure 10) are compared to each other. Figures 12 and 13 show the macroscopic deformation field of the MRAC specimen at different loading stages ( ). Figure 12 shows the lateral displacement distribution of the MRAC specimen obtained by DIC measurements. Several macroscopic fractures were observed for MRAC specimens after failure according to the discontinuity in the lateral displacement contour. The fracture directions are almost vertical. Figure 13 shows the displacement vectors of the block corners of the MRAC specimen obtained by numerical simulations. Similar macroscopic fractures of the MRAC specimen at stage are marked in Figure 13 according to the discontinuity in the displacement.

The Complete Failure Processes
The deformation field and crack distribution of the MRAC specimen from both DIC measurements and numerical simulations at the selected loading stages (marked in Figure 10) are compared to each other. Figures 12 and 13 show the macroscopic deformation field of the MRAC specimen at different loading stages (b ∼ e). Figure 12 shows the lateral displacement distribution of the MRAC specimen obtained by DIC measurements. Several macroscopic fractures were observed for MRAC specimens after failure according to the discontinuity in the lateral displacement contour. The fracture directions are almost vertical. Figure 13 shows the displacement vectors of the block corners of the MRAC specimen obtained by numerical simulations. Similar macroscopic fractures of the MRAC specimen at stage e are marked in Figure 13 according to the discontinuity in the displacement. Figure 14 shows the lateral strain distribution of the MRAC specimen obtained by DIC measurements and a photo of the final fractured pattern. The strain distribution contour shows much more evident localization of deformations than the displacement contour, which reflects the microcracking evolvement of the MRAC specimen during the complete loading process. It is observed that microcracks appear mostly in the ITZs and the mortar region at stage b. Vertical cracks mainly originated from the ITZs, as the loading increases at stage c the cracks propagated in the mortar region. A large number of observable microcracks were generated and formed a connected crack-band at peak loading stage d. The crack-band is slightly different from the final fractured pattern at the post-peak loading stage e, which means some macro-fractures are generated after the peak strength of the specimen. The deformation localization can also be observed in the numerical results according to the velocity vectors of block centroid in Figure 15. Very similar fracture patterns are found in both experimental and numerical results at stage e, which indicates the mesoscopic heterogeneities of the MRAC specimen is well represented by the proposed numerical model.  Figure 14 shows the lateral strain distribution of the MRAC specimen obtained by DIC measurements and a photo of the final fractured pattern. The strain distribution contour shows much more evident localization of deformations than the displacement contour, which reflects the microcracking evolvement of the MRAC specimen during the complete loading process. It is observed that microcracks appear mostly in the ITZs and the mortar region at stage . Vertical cracks mainly originated from the ITZs, as the loading increases at stage the cracks propagated in the mortar region. A large number of observable microcracks were generated and formed a connected crack-band at peak loading stage . The crack-band is slightly different from the final fractured pattern at the post-peak loading stage , which means some macro-fractures are generated after the peak strength of the specimen. The deformation localization can also be observed in the numerical results according to the velocity vectors of block centroid in Figure 15. Very similar fracture patterns are found in both experimental and numerical results at stage , which indicates the mesoscopic heterogeneities of the MRAC specimen is well represented by the proposed numerical model.  Figure 14 shows the lateral strain distribution of the MRAC specimen obtained by DIC measurements and a photo of the final fractured pattern. The strain distribution contour shows much more evident localization of deformations than the displacement contour, which reflects the microcracking evolvement of the MRAC specimen during the complete loading process. It is observed that microcracks appear mostly in the ITZs and the mortar region at stage . Vertical cracks mainly originated from the ITZs, as the loading increases at stage the cracks propagated in the mortar region. A large number of observable microcracks were generated and formed a connected crack-band at peak loading stage . The crack-band is slightly different from the final fractured pattern at the post-peak loading stage , which means some macro-fractures are generated after the peak strength of the specimen. The deformation localization can also be observed in the numerical results according to the velocity vectors of block centroid in Figure 15. Very similar fracture patterns are found in both experimental and numerical results at stage , which indicates the mesoscopic heterogeneities of the MRAC specimen is well represented by the proposed numerical model.   The DIC measurement captured the lateral displacement and lateral strain on the MRAC specimen surface, which reflects the cracking propagation of the MRAC specimen during the loading process. The DEM model can provide further information about the cracking mode of each microcrack since either the tensile and shear contact failure mechanism is in consideration as shown in Figure 5. The contact breakages caused by the tensile failure and the shear failure are plotted in Figures 16 and 17. The width of contacts indicates the separation and sliding displacement of contacts, respectively. There are some horizontal shear cracks in Figure 17, which represents the shear displacement between the specimen and the loading plates due to the Poisson effect.  The DIC measurement captured the lateral displacement and lateral strain on the MRAC specimen surface, which reflects the cracking propagation of the MRAC specimen during the loading process. The DEM model can provide further information about the cracking mode of each microcrack since either the tensile and shear contact failure mechanism is in consideration as shown in Figure 5. The contact breakages caused by the tensile failure and the shear failure are plotted in Figures 16 and 17. The width of contacts indicates the separation and sliding displacement of contacts, respectively. There are some horizontal shear cracks in Figure 17, which represents the shear displacement between the specimen and the loading plates due to the Poisson effect. The DIC measurement captured the lateral displacement and lateral strain on the MRAC specimen surface, which reflects the cracking propagation of the MRAC specimen during the loading process. The DEM model can provide further information about the cracking mode of each microcrack since either the tensile and shear contact failure mechanism is in consideration as shown in Figure 5. The contact breakages caused by the tensile failure and the shear failure are plotted in Figures 16 and 17. The width of contacts indicates the separation and sliding displacement of contacts, respectively. There are some horizontal shear cracks in Figure 17, which represents the shear displacement between the specimen and the loading plates due to the Poisson effect.  A great number of tensile microcracks are found at stage b, which is far away from the peak loading stage. A lot of tensile cracks randomly initiate in the mortar region in the loading direction (stage in Figure 18). Moreover, some tensile cracks are generated along the ITZs, which have small angles (<10°) relative to the loading direction. The shear cracks seem to be more dependent on the meso-structures compared to the tensile ones. Almost every shear crack initiate along the ITZs A great number of tensile microcracks are found at stage b, which is far away from the peak loading stage. A lot of tensile cracks randomly initiate in the mortar region in the loading direction (stage b in Figure 18). Moreover, some tensile cracks are generated along the ITZs, which have small angles (<10 • ) relative to the loading direction. The shear cracks seem to be more dependent on the meso-structures compared to the tensile ones. Almost every shear crack initiate along the ITZs (stage b in Figure 17), which have a certain angle (30 • < 60 • ) relative to the vertical direction.
The fracture keeps growing due to the connecting of microcracks during the loading process. Finally, a predominant macroscopic fracture penetrates through the whole MRAC specimen at stage , and peak loading is then achieved. The loading may suddenly distribute on different fragments of the ruptured specimen, which creates more macroscopic fractures after the peak load, as shown in stage e in Figures 16 and 17. Moreover, the overlays of Figures 16 and 17 of stage e are compared with the failure specimen, as shown in Figure 18. Due to the limitation of the 2D model, the real 3D meso-structure of the concrete specimen cannot be fully modeled in the proposed model. The spatial distribution of aggregates inside the concrete specimen was not captured. Therefore, the numerical model is not able to capture all the cracks or fractures observed from experimental results in every detail. Nevertheless, the numerical results show quite good correlation with experimental results, based on both the stress-strain curves and the crack distributions. As the loading increases, more and more tensile cracks grow and propagate in the loading direction. A few tensile cracks concentrate along the vertical ITZs, and in the mortar regions where between an upper and a lower aggregate (stage c in Figure 16). Some tensile cracks and shear cracks begin to be connected as they are continuously forming and propagating at stage c. Some macroscopic fractures, therefore, appear and lead to the localization of deformation observed in Figure 14. These fractures do not penetrate the whole MRAC specimen; therefore, the stress-strain curve in Figure 10 does not lose the load-carrying capability but shows a nonlinear shape at stage c.
The fracture keeps growing due to the connecting of microcracks during the loading process. Finally, a predominant macroscopic fracture penetrates through the whole MRAC specimen at stage d, and peak loading is then achieved. The loading may suddenly distribute on different fragments of the ruptured specimen, which creates more macroscopic fractures after the peak load, as shown in stage e in Figures 16 and 17.
Moreover, the overlays of Figures 16 and 17 of stage e are compared with the failure specimen, as shown in Figure 18. Due to the limitation of the 2D model, the real 3D meso-structure of the concrete specimen cannot be fully modeled in the proposed model. The spatial distribution of aggregates inside the concrete specimen was not captured. Therefore, the numerical model is not able to capture all the cracks or fractures observed from experimental results in every detail. Nevertheless, the numerical results show quite good correlation with experimental results, based on both the stress-strain curves and the crack distributions.

Cracking Modes
According to the analysis in the last section, the initiation and evolution of microcracks are dependent on the cracking modes, which will control the final failure behaviors of the MRAC specimen. The mode I (tensile microcracks) and mode II (shear microcracks) cracks at the peak loading stage are thus plotted in a single figure (Figure 19). The influence on the cracking modes from the mesoscopic structure of the MRAC specimen can be observed. Either mode I or mode II cracks can be found in the new mortar region. However, most of them are mode I cracks. A small amount of mode II cracks was generated after the macroscopic fractures formed. The cracks of mode I widely distribute in the new mortar region and angle less than 10 • from the loading direction. The cracks of mode I, with large openings, are found to appear mainly between an upper and a lower natural aggregate, where the compressional stress concertation always happens. The cracks of mode I can also be observed in the recycled aggregate region, where they are mostly inside the old mortar region and in the old ITZs, with directions nearly parallel to the loading direction. Because of the high strength of the natural aggregates, none of the mode II cracks are observed in the natural aggregate region. Very few cracks of mode I are found, apart from some particular ones with long, narrow shapes and a long axis perpendicular to the loading direction. dependent on the cracking modes, which will control the final failure behaviors of the MRAC specimen. The mode I (tensile microcracks) and mode II (shear microcracks) cracks at the peak loading stage are thus plotted in a single figure (Figure 19). The influence on the cracking modes from the mesoscopic structure of the MRAC specimen can be observed. Either mode I or mode II cracks can be found in the new mortar region. However, most of them are mode I cracks. A small amount of mode II cracks was generated after the macroscopic fractures formed. The cracks of mode I widely distribute in the new mortar region and angle less than 10° from the loading direction. The cracks of mode I, with large openings, are found to appear mainly between an upper and a lower natural aggregate, where the compressional stress concertation always happens. The cracks of mode I can also be observed in the recycled aggregate region, where they are mostly inside the old mortar region and in the old ITZs, with directions nearly parallel to the loading direction. Because of the high strength of the natural aggregates, none of the mode II cracks are observed in the natural aggregate region. Very few cracks of mode I are found, apart from some particular ones with long, narrow shapes and a long axis perpendicular to the loading direction. The mode II cracks appear mostly along the old ITZs inside the recycled aggregate region, which has angles larger than 10° in the loading direction. While outside the recycled aggregate region, both mode I and mode II cracks can be observed along some new ITZs. These cracks, therefore, can be defined as mixed mode cracks.
The typical cracking modes controlled by their corresponding mesoscopic structures revealed in the numerical simulation are schematically illustrated in Figure 20. Although the cracks of mode I are dominant in quantity, as shown in Figure 11, the final macroscopic fractures are generated due to the propagation and the connection of all different modes of cracks. The failure process of the simulated MRAC specimen can be considered as a self-organizing process of the microcracks. The process is dependent on the spatial distribution of different components of the MRAC specimen at the mesoscopic scale. The mode II cracks appear mostly along the old ITZs inside the recycled aggregate region, which has angles larger than 10 • in the loading direction. While outside the recycled aggregate region, both mode I and mode II cracks can be observed along some new ITZs. These cracks, therefore, can be defined as mixed mode cracks.
The typical cracking modes controlled by their corresponding mesoscopic structures revealed in the numerical simulation are schematically illustrated in Figure 20. Although the cracks of mode I are dominant in quantity, as shown in Figure 11, the final macroscopic fractures are generated due to the propagation and the connection of all different modes of cracks. The failure process of the simulated MRAC specimen can be considered as a self-organizing process of the microcracks. The process is dependent on the spatial distribution of different components of the MRAC specimen at the mesoscopic scale.

Parameter Study and Discuss
The macroscopic failure process of the MRAC specimen is closely related to microcrack development during the complete loading process. Therefore, variations in the mechanical properties of microscopic contacts ultimately affect the macroscopic mechanical behaviors. Some key microscopic contact parameters are selected for further parameter analysis.
The microcracks can hardly be observed inside the natural aggregates during the loading process. It is believed that the natural aggregates have enough tensile and shear strength, which will not strongly influence the macroscopic behaviors of the MRAC specimen. Thus, only the normal stiffness of contact of natural aggregates has been tested for parameter analysis. As shown in

Parameter Study and Discuss
The macroscopic failure process of the MRAC specimen is closely related to microcrack development during the complete loading process. Therefore, variations in the mechanical properties of microscopic contacts ultimately affect the macroscopic mechanical behaviors. Some key microscopic contact parameters are selected for further parameter analysis.
The microcracks can hardly be observed inside the natural aggregates during the loading process. It is believed that the natural aggregates have enough tensile and shear strength, which will not strongly influence the macroscopic behaviors of the MRAC specimen. Thus, only the normal stiffness of contact k n of natural aggregates has been tested for parameter analysis. As shown in Figure 21, the macroscopic stiffness of the MRAC specimen increases with the increase in natural aggregate stiffness, but not evidently. The loading level to start the unstable cracking stage became higher for stiffer natural aggregate cases, which led to higher peak strength.

Parameter Study and Discuss
The macroscopic failure process of the MRAC specimen is closely related to microcrack development during the complete loading process. Therefore, variations in the mechanical properties of microscopic contacts ultimately affect the macroscopic mechanical behaviors. Some key microscopic contact parameters are selected for further parameter analysis.
The microcracks can hardly be observed inside the natural aggregates during the loading process. It is believed that the natural aggregates have enough tensile and shear strength, which will not strongly influence the macroscopic behaviors of the MRAC specimen. Thus, only the normal stiffness of contact of natural aggregates has been tested for parameter analysis. As shown in Figure 21, the macroscopic stiffness of the MRAC specimen increases with the increase in natural aggregate stiffness, but not evidently. The loading level to start the unstable cracking stage became higher for stiffer natural aggregate cases, which led to higher peak strength.    Figure 22 shows that the contact stiffness of ITZs has a greater influence on the macroscopic stiffness of the MRAC specimen than the stiffness of the natural aggregate, while it has little effect on the peak strength.  Axial strain e / % Figure 22. Complete compressive stress-strain curves of MRAC specimen with various ITZ stiffness values. Figure 23 shows that the tensile strength of ITZs seems to have little influence on the macroscopic mechanical behaviors of the MRAC specimen. As we discussed in the last section, the cracks of mode I only appear in the ITZs that are nearly parallel to the loading direction. Tensile failure does not control the cracking process of the ITZs.    Figure 23 shows that the tensile strength of ITZs seems to have little influence on the macroscopic mechanical behaviors of the MRAC specimen. As we discussed in the last section, the cracks of mode I only appear in the ITZs that are nearly parallel to the loading direction. Tensile failure does not control the cracking process of the ITZs. Unlike the tensile strength, the shear strength of the ITZs shows a significant influence on the macroscopic mechanical behaviors of the MRAC specimen, as shown in Figure 24. Because the failure modes of most ITZs are mode II or mixed mode, as shown in Figure 19, the increase in the shear strength of the ITZs can significantly increase the peak strength and the loading levels to start both the stable and the unstable cracking stages. The shear strength of the ITZs needs to be very carefully considered in DEM numerical models. Unlike the tensile strength, the shear strength of the ITZs shows a significant influence on the macroscopic mechanical behaviors of the MRAC specimen, as shown in Figure 24. Because the failure modes of most ITZs are mode II or mixed mode, as shown in Figure 19, the increase in the shear strength of the ITZs can significantly increase the peak strength and the loading levels to start both the stable and the unstable cracking stages. The shear strength of the ITZs needs to be very carefully considered in DEM numerical models.

Conclusions
The failure process of the MRAC specimen under a uniaxial compression test was investigated using both experimental and numerical methods. The actual heterogeneity of the modeled concrete specimen contains different components and ITZs can be well reproduced by combining the DEM model and DIP method. The main conclusions drawn from the proposed study are summarized as follows: (1) The simulated cracking and failure process of the MRAC specimen shows very similar observations to the experimental results. The proposed numerical model provides the possibility to investigate the failure mechanism of the recycled aggregate concrete at a microscopic level, which is beyond the ability of continuum approaches.

Conclusions
The failure process of the MRAC specimen under a uniaxial compression test was investigated using both experimental and numerical methods. The actual heterogeneity of the modeled concrete specimen contains different components and ITZs can be well reproduced by combining the DEM model and DIP method. The main conclusions drawn from the proposed study are summarized as follows: (1) The simulated cracking and failure process of the MRAC specimen shows very similar observations to the experimental results. The proposed numerical model provides the possibility to investigate the failure mechanism of the recycled aggregate concrete at a microscopic level, which is beyond the ability of continuum approaches. (2) The failure process of the MRAC specimen can be considered as a self-organizing process of microcracks and is dependent on the mesoscopic structure of the specimen. The typical cracking modes controlled by the corresponding mesoscopic structures are concluded according to the numerical observations. (3) The parameter study shows that the shear strength of the ITZs has a dominant influence on the macroscopic strength of the recycled aggregate concrete, which needs careful consideration in the DEM models.