Meso-mechanical modelling of damage in concrete using discrete element method with porous ITZs of defined width around aggregates

The article deals with two-dimensional mesoscopic numerical results of fracture in concrete at the aggregate scale in a beam with a notch under bending. The focus was on the effects of both initial micro-porosity in ITZs and aggregate roughness on the load-deflection curve and propagation and location of microand a macro-damage. An improved discrete element method (DEM) was used for concrete to predict its strength, brittleness and fracture. A phase concrete description was used to explicitly take its heterogeneity into account. Concrete included aggregate, mortar, interfacial transition zones (ITZs) and macro-pores. In contrast to existing approaches, ITZs were simulated as porous mortar zones around aggregates with a defined width (without reducing mechanical properties). The real shape and place of aggregates in the beam were established in numerical calculations upon on x-ray micro-CT scans. The findings presented in this paper offer a new perspective as to the understanding of micro-cracking formation in concrete

The paper continues our previous research works to develop a realistic 3D multi-scale mechanical model for concrete in order to replace costly experimental tests and to establish recommendations for the most advantageous concrete mixes as to the amount of the aggregate, its shape and roughness and porosity of the cement matrix (for receiving the optimal concrete in strength and ductility) [12,13]. The model might also investigate the effect of the presence of fibres and additional grains of different syntactic foams (as e.g. glass and polymer micro-spheres, and hollow spheres) on concrete behaviour. In addition, realistic meso-mechanical models for concrete allow also for a better understanding of micro-and macro-cracking phenomena at the aggregate level that have a dominant impact on the global concrete behaviour.
In the current paper, the plain concrete behaviour (stiffness, strength, brittleness and micro-and macro-cracking) was investigated in mesoscopic simulations of a concrete beam under three-point bending using the improved discrete element method (DEM) in the particle-based open-source code YADE, formulated at the University of Grenoble [38,39]. A phase concrete DEM model was again assumed, composed of aggregate, cement matrix with ITZs and macro-voids [12,13]. However, ITZs were simulated now as porous zones around aggregates of a defined width with the same mechanical properties and higher porosity than the mortar, which is an innovative element as compared to existing mesoscopic approaches for concrete. Thus, no special mechanical properties were assumed in ITZs with respect to stiffness and strength (which are usually difficult to be determined in contrast to porosity). Considering ITZs in such a way is more realistic for understanding not only micro-cracking phenomena in concrete at the neighbourhood of aggregates with the different roughness but also humidity flow and heat transport phenomena [10]. The meso-structure of the concrete beam was obtained in numerical calculations with the x-ray tomography system (micro-CT) [12,13,40,41].
The main focus was on investigations of the quantitative effect of initial micro-porosity in ITZs of a defined with around aggregates on the fracturing pattern and load-deflection curve in concrete. The DEM results were compared with the experiment. Moreover, two other aspects were qualitatively studied: (1) the effect of aggregate micro-roughness on the load-deflection curve and micro-damage evolution and (2) the effect of a loading type (bending and uniaxial compression) on the formation of micro-cracks around aggregates. Both, the processes of micro-cracking initiation in the vicinity of aggregates and the ITZ/aggregate interaction are still not fully understood. Research works in [42] on different ITZ micro-cracking scenarios provided the motivation for investigations of those two aspects. To cut a huge computer effort related to the number of particles, four simplifications were assumed in the present computations: (1) 2D simulations instead 3D ones, (2) a few times larger width of ITZs than in experiments, (3) larger particles in ITZs and (4) constant porosity in ITZs. The 2D DEM simulations were carried out with about 40,000 particles. The computation time of one numerical test (being almost proportional to the number of particles) was about 7 days using CPU 2.8 GHz.
With the more powerful computer and parallelization of DEM in the code YADE (in the course of the realization), the full 3D computations will be performed in a reasonable time with the real both the width of ITZs and particle sizes in ITZs.
The paper is arranged as follows. A brief summary of related past work is given in Section 2. In Section 3, we describe our experiments on concrete, before detailing our proposed numerical framework in Section 4 and input data in Section 5. Results on the effect of initial porosity of ITZs and grain roughness on cracking during bending with the key findings are discussed in Sections 6 and 7. We summarize and conclude in Section 8.

Literature overview
ITZs originates from a wall effect of cement grain packing against the aggregate surfaces which disrupts the packing of cement grains and contributes to their higher porosity [3]. Another origin of ITZs is a micro-bleeding phenomenon -ITZs accumulate free water, especially under large aggregates what also increases their local porosity [43]. ITZs are believed to be affected by several factors, including the particle diameter, water/cement ratio, aggregate roughness, aggregate size distribution and curing age. Insight Fig. 1. Experiments with concrete under bending: (a) geometry and boundary conditions of concrete beams subjected to three-point bending [12,13] and (b) view on some aggregate particles with different diameter used for concrete.
into the micro-structural details of micro-cracking processes around aggregate surfaces is still very limited due to a low imaging resolution. Usually, two microscopic experimental methods are mainly applied for examining the properties of ITZs: backscattered electron (BSE) imaging [3,44,45] and nano-indentation [45][46][47][48][49]. Other applied techniques are mercury intrusion porosimetry and transport tests. According to [44], the ITZ porosity and modulus of elasticity of ITZs decrease in a linear/parabolic way towards the mortar. Concrete specimens cast with a cement of higher flowability have a higher filling density and lower porosity in the ITZ area. Concrete specimens cast with a higher water to binder (w/b) ratio have wider and more porous ITZs. Larger aggregates usually result in wider more porous ITZs and lower bond strength. A decreasing aggregate size causes weaker and more porous ITZs in rock-filled concrete [43]. The ITZ thickness grows with increasing aggregate size and w/c ratio and decreasing curing time [46,50,51]. The curing age has an effect on both ITZ thickness and porosity while the w/b-ratio and aggregate content influence the ITZ porosity only [52]. The elastic modulus of ITZs is much lower than of the mortar [51,53]. The thickness of ITZ is comparable to the median cement grain size used to make a cementitious composite [3]. The thickness of ITZs may be close to zero if the aggregate surface is extremely smooth [24]. A more porous ITZ micro-structure is created around rough aggregates [5]. In general, an uneven distribution of both the thickness and porosity of ITZs happens [5]. Due to gravity, more porous ITZ regions form underneath the aggregates than above them [1,5]. The most effective way to densify ITZs is to add fine particles (e.g. silica fume) which pack closer to aggregates [3,54]. A reduction of the water/cement ratio by applying plasticizers may be also applied. Two different crack initiation scenarios in ITZs were proposed in [42] based on experiments and mechanical analytical considerations in the framework of continuum micro-mechanics, namely ITZ aggregate separation/delamination and ITZ failure. Those considerations were restricted to the identification of the onset of micro-cracking based on tensile elastic limit criteria. The initial micro-cracks occurred directly at the aggregates if the strength of ITZs was larger than the separation strength or were created in ITZs at a certain distance from aggregates if the ITZ strength was smaller than the strength of interfaces between aggregates and ITZs. Under uniaxial tension-dominated loading, both the ITZ failure and ITZ aggregate separation mainly appeared while under uniaxial compression-dominated loading, the ITZ failure happened as the more likely mechanism [42].
Mechanical properties of concretes (related to the interface behaviour in ITZs around aggregates) also depend on micro-structural characteristics of aggregates; they vary with the different roughness of coarse aggregate surfaces [55][56][57][58][59]. Concretes with rougher aggregate surfaces produced higher splitting tensile strength, compressive strength, elastic modulus, fracture energy and Poisson's ratio due to a growth of the interface strength [59]. The tensile and shear bond strength increased in a parabolic way as the roughness of the aggregate surface increased [59]. The mechanical parameters increased slowly after the aggregate surface becomes rougher [59]. During splitting tests, the failure cracks passed through ITZs independently of the aggregate roughness [59]. Moreover, the concrete mechanical behaviour was affected by the aggregate type [58] due to differences in chemical reactions with the cement paste. The relationship between the ITZ porosity and aggregate roughness was not studied in experiments.

Own experiments on concrete
The experimental data of a three-point bending concrete beam with a notch was described in details in our papers [12,13]. The geometry of the used beam was 320 × 80 × 40 mm 3 (length, height and width) [12]. The distance between supports was 240 mm (3 × H) (Fig. 1a). The ordinary Portland cement, the round-shaped gravel (Fig. 1b) and sand particles were used for concrete mixture, with maximum aggregate d max = 16 mm and mean aggregate d 50 = 2 mm [12]. The aggregate greater or equal to 2 mm in diameter filled the specimen in 37% (aggregate and sand filled in together 75% the entire volume) [12]. The aggregate surfaces were smooth, i.e. the maximum normalized roughness R n(max) , which was expressed as the ratio between the maximum mean distance between the peaks and valleys along the aggregate surface, and the aggregate diameter was low (R n(max) = 0.045). The micro-porosity of the mortar was established at the 5% level in the x-ray micro-CT measurements [41]. The compressive strength of concrete was f c = 51.81 MPa, Young's modulus was E = 36.1 GPa, the Poisson's ratio was v = 0.22 and the flexural tensile strength was f t,flex = 4.0 MPa [12]. The maximum vertical force F max obtained during three-point bending tests was varied between 2.15 and 2.25 kN (f t,flex = 3.73-3.91 MPa) [12].
The size, shape and positions of the aggregate and macro-voids were established with micro-CT images of the specimen's mid-part (80 × 50 × 40 mm 3 ), cut out from each beam, just after the test [12]. The new generation x-ray micro-tomography system SkyScan 1173 was used (the pixel size varied from 3 µm up to 90 µm and large objects with a diameter up to 200 mm might be used) [41,60]. Due to a random distribution of aggregates, the macro-crack above the notch was pronounced curved in the 3D space [12,13]. It occurred by bridging the interfacial cracks [12] (micro-damage first occurred in ITZs). The width of porous ITZs on the concrete surface was between 20 and 50 μm (Fig. 2a) that was measured with the scanning electron-microscope (SEM) Hitachi TM303 [12]. ITZs were visible around about 90% of the aggregate circumference (in all aggregates (d a ≥ 2 mm)) due to the creation of water lenses beneath aggregates during concrete mixing [3]. An image binarization technique was next used to find the porosity distribution in ITZs around aggregates (Fig. 2). From the digital image of ITZ (Fig. 2a), a white-black bitmap was first created, wherein the black colour corresponded to porous and white colour to aggregate and cement matrix. Next, the width of ITZ was divided into strips of the aggregate shape 5 μm wide (Fig. 2b), and the black and white pixels were calculated. The porosity of each strip was determined as the percentage of the area of pores to the total strip area. The experimental porosity of ITZs changed between 25% (directly at the aggregates) down to 5% (in the cement matrix) (Fig. 2c). Fig. 3 shows an exemplary experimental micro-crack above the main macro-crack which was created inside of ITZ (it was observed on the concrete surface again with SEM). Finally, some nanoindentation tests [47] were carried out on small concrete specimens to investigate the variation of the modulus of elasticity in three different concrete phases (aggregate, mortar and ITZ) [45,47]. The results of one exemplary test for the region of one aggregate are shown in Fig. 4. They demonstrate that the modulus of elasticity of ITZs is smaller by 35% than the modulus of elasticity of the mortar and by 80-85% than the modulus of elasticity of the aggregate (Fig. 4). Figs. 2-4 clearly indicate on some specific properties of ITZs as compared to other phases in concrete. It has been shown in [49] that it is impossible to probe representative volumes of highly heterogeneous materials like mortar/concrete based on nanoindentation tests. The results of nanoindentation tests (Fig. 4) were not used in the DEM calculations.

Formulation of discrete element method (DEM) for concrete
Meso-scale models directly simulate the material-structure and therefore may be used to comprehensively study the mechanism of the initiation, growth and formation of micro-and macro-cracks [13]. They may also be applied for studying different local phenomena at the aggregate level (e.g. force chains, vortex structures, bottlenecks) in order to early identify a final crack location [11,62]. The disadvantages of meso-scale models are the huge computational cost.
The DEM calculations were performed with the three-dimensional spherical explicit discrete element open code YADE [38,39], allowing for a small overlap to be formed between two contacted bodies (soft-particle model). In DEM, particles interact with each other during translational and rotational motions through a contact law and Newton's 2nd law of motion using an explicit timestepping scheme [63]. The ease of access to information at the particle scale makes DEM a very useful tool to study the dynamics of systems involving particles. Our DEM model was able to faithfully reproduce shear localization in granular materials [64][65][66][67][68] and fracture in concrete [12,13,24,29,30,61]. A cohesive bond was assumed at the grain contact exhibiting brittle failure under the critical normal tensile load [39,69]. The shear cohesion failure initiated contact slip and sliding obeying the Coulomb friction law under normal compression [39]. The tensile failure initiated contact separation. The DEM model for concrete was described in detail in [12,13]. It may be summarized as follows (Eqs. (1)-(7), Fig. 5): where F n -the normal contact force, U -the overlap between discrete elements, N -the unit normal vector at each contact point, F sthe tangential contact force, F s prev , -the tangential contact force from the previous iteration, X s -the relative tangential displacement increment, K n -the normal contact stiffness, K s -the tangential contact stiffness, E c , -Young's modulus of the particle contact, ν c -the Poisson's ratio of particle contact, R -the particle radius, R A and R B -the contacting particle radii, μ -the Coulomb inter-particle friction angle, F max s -the critical cohesive contact force, F min n -the minimum tensile force, C -the cohesive contact stress (maximum shear stress at the pressure equal to zero) [69] and T -the tensile normal contact stress [69], F damp k -the damped contact force, F k and v p k -the k th components of the residual force and translational particle velocity v p and α d -the positive damping coefficient smaller than 1 (sgn(•) returns the sign of the k th component of velocity) [70]. If any contact between spheres after failure re-appeared, the cohesion will not turn up (Eq. (5)). The damage was created if a cohesive joint between spheres (Eq. (6)) disappeared after reaching a critical threshold. A simple local non-viscous damping was involved [70] to accelerate convergence towards quasi-static equilibrium (Eq. (7)). Note that material softening was not considered in the DEM model. Albeit bonds can break by shear, the essential micro-scale mechanism for damage in the pre-failure regime was the bond breakage in tension [61]. Five main local material parameters are required for discrete simulations: E c , υ c , μ, C and T. In addition, the particle radii R, particle mass density ρ and damping parameters α d have to be known. The detailed calibration procedure was described in [12,13,29], based on simple standard  [38,39]. laboratory tests of concrete specimens due to a data lack regarding mechanical properties of ITZs (in contrast, the properties of the mortar and aggregates can be easily determined since they can be individually tested). The DEM model is also able to realistically reproduce the concrete behaviour under high confinement [71]. The crack was not allowed to propagate through aggregate, i.e. the particle breakage was not been taken into account.
Our combined DEM/x-ray μ-CT images mesoscopic approach proved to be an extremely appealing computational tool for investigating fracture in concrete under bending (2D and 3D analyses) [12,13], uniaxial compression (2D and 3D simulations) [24,29] and splitting tension (2D analyses) [30]. In those calculations, ITZs had no physical width and were always simulated by weaker contacts between aggregates and mortar paste [12]. Such an approach was effective enough to capture mechanical fracture. However, the approach is not appropriate to study coupled mechanical-hydro-thermal problems in concrete that requires to take porous ITZs into account since the porosity favours the heat and fluid transport [10]. In the current paper, the mechanical problems were studied only (the coupled mechanical-hydro-thermal problems were not discussed).

Beam geometry, concrete meso-structure and construction of porous ITZs
The concrete beam was assumed as a phase material. It included aggregate, cement matrix with porous ITZs and macro-pores in the mid-region (50 × 80 mm 2 ) of the beam (Fig. 6) [12]. The place and shape of aggregates in this region were the same as on the micro-CT images [12,13]. In 2D calculations, the beam thickness included one row of aggregates and mortar particles [12]. All aggregate grains (2 mm ≤ d a ≤ 16 mm) were simulated as grain clusters composed of infinitely rigid bodies [12,13]. The surface of aggregates was smooth as in the experiments (Section 2). All aggregate particles with d a ≥ 2 mm included porous ITZs [12]. The spherical elements without ITZs were used to describe the cement matrix (0.25 mm ≤ d cm < 0.75 mm with d min = 0.25 mm). The presence of the hydrate-foam matrix in the mortar was no taken into account [27]. ITZs were solely composed of mortar spheres of d ITZ = 0.25-0.35 mm [12]. Their micro-porosity p (ITZ)  -the total volume). The width of ITZs was chosen to be 0.75 mm, i.e. ITZs were composed of 3 rows of spheres (Fig. 7). Thus, it was about 10-times larger than in reality. The 10-times decrease of the width of ITZs (connected with a strong decrease of the minimum particle diameters) would be linked with too huge computation time. The DEM calculations were carried out with the different initial porosity of ITZs: p (ITZ) = 5% (as in the cement matrix), 7%, 10%, 15%, 20% and 25% (maximum experimental value in Fig. 2c).
To construct porous ITZs, the spheres of the cement matrix were inserted into the box of the concrete specimen size (with μ = 0°) until the desired porosity was obtained (e.g. 20% for the specimens with ITZs) [12]. The specimen was kept until its total kinetic energy dropped under the required value. The particles at the distance of 0.25-0.75 mm (3 layers) from aggregates were assumed as ITZs with the prescribed constant porosity. The aggregates and ITZs around them were kept blocked while the other aggregate particles were removed. Then, the cement matrix particles were added to the specimen with the porosity of 5% that corresponded to the real concrete compactness. Next, all spheres were again allowed to settle down below a certain kinetic energy level [12]. Finally, all contact forces due to the particle penetration U were deleted and μ was set on the target value [12]. To reproduce the starting configuration, the initial overlap was subtracted in each calculation step when determining the normal contact forces , where U o -the initial overlap and U n -the overlap in the calculation step n) [61]. Macro-pores larger than 1.0 mm in diameter were simulated as empty regions with a real shape [12]. Initial stresses in concrete due to shrinkage were not taken into account since the experimental concrete beams were carefully prepared to avoid the surface evaporation and autogenous shrinkage [12]. The beam with 2D particle clusters (Fig. 6) included in total about 40 000 spheres (38 000 in the meso-region: 31 000 in the mortar matrix and ITZs and 7000 in the aggregate) with the minimum diameter of d min = 0.25 mm. Fig. 6. Geometry of concrete beam model in 2D DEM with meso-scale mid-region (4-phase model for concrete): front side of the entire beam (midregion: dark grey colour corresponds to aggregates, grey colour to cement matrix, light grey colour to ITZs and white colour to macro-voids) [12,13].

Model calibration (based on uniaxial compression test)
In the current paper, the following parameters were used for the mortar in the mid-region during all 2D analyses (Fig. 6): E c,cm = 16.8 GPa, ν c = 0.2, C cm = 180 MPa and T cm = 28 MPa (C cm /T cm = 5.5), based on a 2D uniaxial compression tests for concrete (see Fig. 8). The values of E c,cm , C cm and T cm were slightly higher than in the 2D analyses during bending wherein ITZs were merely simulated as weak contacts between the cement paste and aggregates with the reduced stiffness and strength of the cement matrix (E c,cm = 11.2 GPa, ν c = 0.2, C cm = 140 MPa and T cm = 25 MPa [12,13]). The remaining DEM parameters were equal to μ = 18°(inter-particle friction angle) and ρ = 2.6 kG/m 3 (mass density) [12]. The inter-particle friction angle μ was assumed based on triaxial compression tests with granulates in the form of clumps [64]. The damping parameter (equal to α d = 0.08) did not affect the results [65]. Note that due to a tensile failure mode in the beam, the effect of μ on fracture was negligible. The effect of d min on DEM results was comprehensively studied during uniaxial tension [29,30] and uniaxial compression [24]. If d min was small enough, i.e. d min ≤ 0.25-0.50 mm, the effect of d min on the stress-strain curve was insignificant. Since the post-peak characteristics in accompanying laboratory uniaxial compression experiments (following standard regulations) was not measured, the concrete DEM behaviour in the softening regime might not be compared with the laboratory test outcomes. The calculation results obviously depend on the width of ITZs [19]; the concrete strength diminishes with growing the width of weak ITZ zones (Fig. 10B). With the decreasing width of ITZs (connected to the smaller d min ), the material constants for the mortar E c,cm , C cm and T cm become lower. Thus, they remain in the range: E c,cm = 11.2-16.8 GPa, C cm = 140-180 MPa and T cm = 25-28 MPa for the ITZ width between 0 and 0.75 mm. The remaining beam parts consisted of larger spheres (Fig. 6) and were modelled with the same material parameters as mortar, i.e. E c,cm = 16.8 GPa, ν c = 0.2, C cm = 180 MPa and T cm = 28 MPa (C cm /T cm = 5.5) [13].
To calibrate the DEM model, some numerical compression tests were carried out on the concrete specimen 50 × 50 mm 2 , cut out from the beam mid-region (50 × 80 mm 2 ).  (ITZ) in ITZs, both the concrete strength and concrete stiffness decreased but concrete ductility increased at the same time during uniaxial compression. The computed failure mode in specimens indicated the occurrence of a few vertical and skew macro-cracks (Fig. 8B).
Based on the evolution of broken normal contacts (equivalent to micro-cracking) in individual concrete phases of the specimen with p (ITZ) = 20% (Fig. 9), it can be observed that most of the contacts were damaged in the mortar (curve 'd') and sequentially in ITZs (curve 'a'), directly at the aggregate (curve 'c') and in the mortar/ITZ layers (curve 'b') in the entire post-peak region. The number of broken contacts in ITZs for ε = 0.3% was about 3 times higher than directly at the aggregates and in the mortar/ITZ layers (3000 against 900 for ε = 0.3%). Thus, under compression-dominated deformation, the ITZ failure dominated as in the experiments [42].
The 3D simulations provide obviously more realistic outcomes than the 2D ones [13]. The stress differences between 2D and 3D DEM simulations for concrete were found to be minimum at the peak load [13]. Those differences increase with growing post-peak deformation when assuming the same material parameters [13]. However, the 2D analyses may be also used for obtaining realistic results of fracture (location, shape and width) and a post-peak stress-strain response when: (a) the 2D DEM model is properly calibrated for concrete with laboratory tests [13,30] and (b) the 2D calculations are performed for the given concrete meso-structure in a cross-section. In addition, the 2D analyses strongly cut the computation time of DEM simulations. Note that the model extension to 3D conditions is numerically straightforward [13]. The smaller spherical elements in ITZs can be also easily incorporated into the model.

Force-displacement curves
The numerical curves of the resultant vertical force F versus CMOD with the different initial micro-porosity p (ITZ) of ITZs (p (ITZ) = 5-25%) are shown in Fig. 10A. They were compared to one experimental outcome for the vertical cross-section of the beam mid-region at the depth of 10 mm from the front side [12,13]. In addition, the effect of the number of spherical layers in ITZs (1)(2)(3) was also demonstrated for p (ITZ) = 20% (Fig. 10B).
The effect of p (ITZ) proved to be strong on the beam strength and stiffness, and medium on the beam brittleness. The calculated maximum vertical force from 2D DEM obviously decreased with increasing porosity of ITZs p (ITZ) (Fig. 10A). It was equal to F max = 2.85 kN (CMOD = 0.020 mm) for p (ITZ) = 5% and F max = 1.9 kN (CMOD = 0.018 mm) for p (ITZ) = 25% (decrease by 30%). The concrete brittleness decreased with increasing p (ITZ). The best agreement with the experiment was achieved for p (ITZ) = 20% with respect to strength and brittleness. The reduction of the width of porous ITZs obviously increased the beam strength (Fig. 10B). In order to get the smoother curves F = f(CMOD) and better accordance with experiments, the 3D DEM analyses are required that use a significantly larger number of particles [13]. With the much lower experimental width of ITZs, the effect of p (ITZ) would be smaller [19].

Beam fracture
The final crack path in the region above the notch (in red) is described in Fig. 11. The curvature of the final crack became larger with increasing p (ITZ) . For p (ITZ) = 15-25%, the shape of the macro-crack was similar as in the experiment. The height h cr and width w cr of the macro-crack were also similar and equal to h cr = 0.65D = 4.5 cm and to w cr = 0.1 mm above the notch. The length of the crack varied from 4.7 up to 5.2 cm; it grew with increasing p (ITZ) .

Broken contacts
The micro-damage evolution in the beam meso-region (based on the evolution of broken contacts) in individual concrete phases versus CMOD with the different initial porosity of ITZs p (ITZ) (5%, 10% and 20%) is shown in Fig. 12.
The initial contact number was about 75,000-85,000. The particle contact damage started far before strength peak. It started obviously earlier with the higher p (ITZ) : for CMOD = 0.013 mm with p (ITZ) = 5%, for CMOD = 0.009 mm with p (ITZ) = 10% and for CMOD = 0.0045 mm with p (ITZ) = 20% (Fig. 12). The number of all broken contacts (CMOD = 0.1 mm) decreased with increasing p (ITZ): 600 (p (ITZ) = 5%), 515 (p (ITZ) = 10%) and 450 (p (ITZ) = 20%). The majority of all contacts was broken in the cement matrix and in ITZs for p (ITZ) = 10-20%. The number of broken contacts in ITZs increased with increasing p (ITZ) and was 3.5 times higher than both in the mortar/ITZ layer and directly at the aggregates for CMOD = 0.1 mm. Thus, the ITZ failure again dominated, similarly as in our tests (Fig. 3) and splitting tests [59]. However, this outcome was opposite to direct tensile tests in [42]. As compared with uniaxial compression for p (ITZ) = 20% (Fig. 9), the relatively more contacts were broken in ITZs and in the mortar/ITZ layers, the relatively fewer contacts were damaged in the mortar and directly at the aggregates (Fig. 12).

Global micro-porosity change
The evolution of the global micro-porosity increment p-p 0 in the entire 2D meso-scale beam region due to deformation is shown in Fig. 13. The clear dilatancy occurred from the beginning of deformation (due to micro-cracking) that decreased with increasing p (ITZ) . After the peak force, (CMOD = 0.02 mm), the dilatancy rate was similar independently of p (ITZ) . The maximum final value of the increment p-p 0 was 1-1.25%.

Effect of grain roughness during bending
Some comparative numerical computations were carried out with the same location, shape, size of aggregates and global microporosity of the cement matrix as in Section 5, but with the increased micro-roughness of all aggregates. The higher micro-roughness was artificially obtained by pulling out some spheres from the aggregate surfaces (Fig. 14). The highest distance between the peaks and valleys along a rough aggregate surface was on average 0.4 mm and the mean distance between peaks was on average 0.75 mm. Thus, the numerical results were qualitatively compared with the experiments only [59]. Two different DEM analyses were performed. In the first case, both the tensile and shear bond strength around aggregates remained the same as for smooth aggregates (E c,cm = 16.8 GPa, ν c = 0.2, C cm = 180 MPa, T cm = 28 MPa, μ = 18°and ρ = 2.6 kG/m 3 ). In the second case, the tensile bond strength around aggregates was increased by 25% (T cm = 35 MPa) (based on experiments [59]) whereas the remaining parameters E c,cm , ν c , C cm , μ and ρ were the same. Fig. 15 shows the numerical curves of the resultant vertical force F versus CMOD and final crack paths with the different aggregate roughness for p (ITZ) = 5% and p (ITZ) = 20% (CMOD = 0.15 mm), although both the initial porosity of ITZs and the width of ITZs may be higher for rough aggregates after concrete mixing than for smooth aggregates [5].

Force-displacement curves
The growth of the aggregate roughness contributed to a clear increase of the stiffness, strength and fracture energy of the concrete beam (Fig. 15A) due to stronger particle interlocking. This outcome is in agreement with splitting tests in [59]. The strength increased by about 25% for p (ITZ) = 20%. The residual stress was also higher with the rough aggregates by about 25% as in [60]. The trajectory of the final crack (Fig. 15B) in concrete with rough aggregates was similar. However, with the rough aggregates, the branching phenomenon of the macro-crack also occurred at the beam mid-height.

Crack displacements
The evolution of the average normal and shear displacement along the main macro-crack just above the notch for smooth and rough aggregates is demonstrated in Fig. 16 (initial porosity of ITZs was p (ITZ) = 20%). The normal crack displacement was correlated   with CMOD and linearly increased during bending independently of the aggregate roughness. The tangential crack displacement was several times smaller than the normal one due to bending. It was as twice as high for rough aggregates.

Broken contacts
The evolution of broken contacts N in the beam meso-region for individual concrete phases versus CMOD with the initial porosity of ITZs p (ITZ) = 20% for rough aggregates is demonstrated in Fig. 17. Fig. 18 shows the distribution of broken contacts at arbitrary aggregates for p (ITZ) = 20% with smooth and rough aggregates.
The growth of the aggregate roughness strongly affected the micro-damage (Figs. 17 and 18). For the rough aggregates in comparison to the smooth aggregate, the more contacts were broken in ITZs (200 against 175) and the fewer in the mortar (140 against 170) and similarly in the mortar/ITZ layers (75 against 75) and directly at the aggregates (50 against 40). Thus, ITZ failure again dominated.

Contact force network
In Fig. 19, the maps of the normal and tangential contact forces in the post-peak regime are plotted. The red colour corresponds to the large compression normal contact forces and the blue one to large tensile normal contact forces (higher than the mean ones) [13]. The high normal and tangential contact forces occurred along the crack path due to interlocking. The maximum compressive normal contact forces were almost four times higher in the specimen with the rough aggregates (F n max = 4.16 N versus F n max = 1.10 N) (Fig. 19a). The maximum tensile normal contact forces were higher (by factor 2) close to the crack path with the rough aggregates (F n max = 1.70 N against F n max = 0.98 N). The maximum tangential contact forces were 3 times higher for the rough aggregates ( Fig. 19b) (F s max = 6.7 N against F s max = 2.11 N).

Tensile bond strength around aggregates T cm = (35 MPa)
The growth of the tensile bond strength from 28 MPa to 35 MPa for rough aggregates straightened a final macro-crack ( Fig. 20  versus Fig. 15Bd). The macro-crack branching was also weaker. The force-displacement curve was almost the same. The number of broken contacts directly at the aggregates was reduced (Fig. 21) whereas the number of broken contacts in ITZs was similar.
To further extend the knowledge on micro-cracking initiation and propagation in concrete, some progress in experimental techniques for examining ITZ regions is also needed with respect to the imaging resolution. In addition, the bond strength depending upon the aggregate roughness has to be investigated in detail since it is a very influential fracture factor. Our DEM model with thin porous ITZs around aggregates will be also connected soon with the computational fluid dynamics (CFD) to study problems of humidity transport in concrete [10,72].

Summary and conclusions
The detailed 2D analyses were conducted to examine a particle-scale response during the deformation process to elucidate the micro-cracking evolution in ITZs around aggregates responsible for the macroscopic concrete properties. In contrast to existing numerical approaches, ITZs were simulated as porous mortar zones around aggregates with a defined width and without reducing mechanical properties. Our 2D DEM calculations based on micro-structure obtained from x-ray micro-CT images showed good agreement with quasi-static laboratory bending tests in concrete damage and strength behaviour at the macro level. The new approach showed a better and more realistic correspondence between DEM and experimental as compared to other approaches. The approach may be used for humidity and heat transport problems. From our numerical analyses and experimental data, the following -The effect of initial micro-porosity in ITZs on the beam response was pronounced. The beam strength, stiffness and brittleness increased with decreasing micro-porosity of ITZs. The macro-crack curvature increased with increasing micro-porosity of ITZs. The growth of porosity in ITZs induced the growth of broken contacts both in ITZs and cement matrix/ITZ layers and reduced the number of broken contacts both in the cement matrix and directly at the aggregates. -The growth of the aggregate roughness contributed to a pronounced increase of the global strength (peak and residual) and stiffness due to stronger particle interlocking. The normal and tangential contact forces at the macro-crack were much higher for rough aggregates. The tangential crack displacement was also as twice as high. The geometry of the macro-crack was similar although its branching was additionally obtained for rough aggregates. -With the growth of the aggregate roughness, more contacts broke in ITZs and less in the mortar. With the growth of the tensile bond strength around rough aggregates, fewer contacts were damaged directly at the aggregates. -As compared with uniaxial compression, the more %-contacts were broken in ITZs and in the mortar/ITZ layers and the fewer %-contacts directly at the aggregates and in the mortar under bending-dominated deformation.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.