Meso-Scale Simulation of Concrete Based on Fracture and Interaction Behavior

: Based on the cohesive zone model, a meso-scale model is developed for numerical studies of three-phase concrete under tension and compression. The model is characterized by adopting mixed-mode fracture and interaction behavior to describe fracture, friction and collision in tension and compression processes. The simulation results match satisfactorily with the experimental results in both mechanical characteristics and failure mode. Whole deformation and crack propagation process analyses are conducted to reveal damage evolution of concrete. The analyses also set a foundation for the following parametric studies in which mode II fracture energy, material parameter, frictional angle and aggregates’ mechanical characteristics are considered as variables. It shows that the mixed-mode fracture accounts for a considerable proportion, even in tension failure. Under compression, the frictional stress can constrain crack propagation at the beginning of the damage and reestablish loading path during the softening stage. Aggregates’ mechanical characteristics mainly a ﬀ ect concrete’s performance in the mid-and-late softening stage.


Introduction
It is widely accepted that the dominant failure mode of concrete is quasi-brittle fracture, which is mainly due to the expanding and connecting of initial defects and microcracks.To explain the behavior of concrete's fracture, concepts of fracture mechanics were introduced more than 50 years ago [1].Although it is suitable for applying fracture mechanics theories to the analysis of mass concrete structures, it is inappropriate for application when microcrack and subcritical crack are not ignorable [2].The sizes of microcracks and subcritical cracks of concrete are several orders of magnitude larger than that of metal.Meanwhile, concrete is a kind of multiphase and heterogeneous material, so the COD (crack opening displacement) and J-integral in elastic-plastic fracture mechanics are inappropriate for application [3].
The macro-mechanic response of concrete is the reflection of its meso-mechanical and micro-mechanic characteristics [2].According to this idea, scholars have attempted to use finite element methods to research fractures in concrete, similar to the research of metal fractures [4].Due to the advance of multi-processors' parallel processing capability, finite element simulation of composite material's meso-mechanical characteristics has progressed [5].Through simulation, it is possible to predict the macro performance of concrete with its specific meso-scale structure, material and internal property [6].
Based on the macroscopic experimental phenomena of concrete, several concrete crack models have been proposed to simulate the fracture behavior of concrete.In the 1960s, two basic models were proposed: the discrete crack model [7] and smeared crack model [8].The former simulates cracking by the definition of nodes' split criterion, which can describe the geometric discontinuity at cracks.The latter simulates cracking by the definition of elements' damage, and the model is maintained to be geometrically continuous.Representative discontinuous models include the lattice model [9], rigid bodies-spring model [10], and interface constitutive model [11].Although the previously mentioned discrete crack model can provide intuitionistic and realistic simulations of crack behavior, their calculation results are sensitive to mesh distortion.In recent years, simulating discrete crack by cohesive zone model [12] and extended finite element method [13] are widely recognized as reliable approaches for their low mesh sensitivity and satisfactory convergent property [14].
The cohesive zone model simulates fractures by inserting cohesive elements at the potential crack propagation paths.The method breaks through the limit of using stress fields and stress intensity factor to describe concrete fracture.Unlike traditional investigation of fracture mechanics, the method uses fracture energy to define the constitutive relation between element deformation and damage on a much grander scale.As mentioned above, the special fracture behavior caused by microcracks and subcritical cracking can be described by elements' fracture energy.Ongoing studies have been being conducted on concrete fracture energy [15][16][17][18][19][20], and corresponding testing methods of different fracture types have been proposed [21,22].
The cohesive zone model has been frequently adopted in the simulation of the concrete failure problem and proved to be feasible in concrete tension fracture simulation [23][24][25][26][27][28].However, very few studies report its application in concrete compression simulation, which has a more complex damage evolution in the fracture process [29].
Although it is widely accepted that the main concrete failure process is fracturing [30] and many concrete crack models have been proposed to describe this process, there is still a need of a unified approach to simulate the mechanical behavior and appearance of concrete under different loading conditions.The objective of this research is to extend the application range of the cohesive zone model, making the model suitable for concrete fracture simulation both under tension and compression.Quantifying the effects of the different material and interaction parameters studied in this paper contributes to revealing the relationship between meso-scale fracture and macroscopic behavior of concrete.

Finite Element Modeling Method
The explicit method in the Abaqus/Explicit module is used in this paper to simulate concrete tension and the compression loading process.Since friction and collision of particles play a significant role in concrete's compression softening stage, interaction behavior is defined to describe the contact behavior.The general contact algorithm in Abaqus/Explicit can only be used with three-dimensional surfaces [31], so three-dimensional models are used for simulation.
The models are built to simulate concrete uniaxial tension and compression experiments conducted by Ren et al. [32].The specimens used for experiments are cuboids with dimensions of 150 × 150 × 50 mm, which were cut from large plate specimens with dimensions of 500 × 500 × 50 mm.The mean tension and compression elastic modules are 42.0 GPa and 37.5 GPa, respectively.The mean tension and compression strengths are 3.24 MPa and 37.5 MPa, respectively.The volume proportion of aggregates calculated by the mixture ratio is 45.1%.For computing efficiency, aggregates whose particle diameters are greater than 2.5 mm are modeled.The fine aggregates and cement matrix are treated as mortar whose mechanical behavior is simulated by mortar particle elements and internal interface elements.
As incomputable microcracks occur during the fracture process, it is difficult to apply the fracture criterion to simulate concrete fracture at the micro-scale.However, it is appropriate to use the traction-separation law to describe the meso-scale damage and fracture caused by the microcracks.For this reason, cohesive elements are used to simulate damage evolution and energy dissipation of concrete.In the main study of this paper (Section 6.4), the cohesive elements are only inserted between the elements of aggregates and mortar (interfacial transitional zones, ITZs), and inside the mortar elements (mortar internal interface, MII), as shown in Figure 1.In the study of Section 6.4, the cohesive elements are also inserted inside the aggregate elements (aggregate internal interface, AII), to study the effects of aggregates' mechanical characteristics.In this paper, the study adopts the simplifying assumption that the interface properties of ITZs, MII, AII are the same.
Appl.Sci.2019, 9, x FOR PEER REVIEW 3 of 22 interface, AII), to study the effects of aggregates' mechanical characteristics.In this paper, the study adopts the simplifying assumption that the interface properties of ITZs, MII, AII are the same.The models are meshed by irregular triangle elements in the plane (Figure 1).Compared with other polygon elements, triangle elements provide the maximum number of interfaces with the same nodes in the plane.Cohesive elements are inserted into the interfaces to simulate damage and fracture behaviors.The displacement loading control is adopted to obtain the complete tension and compression process data.The vertical motions of the nodes at the bottom boundary are constrained.The load is applied by controlling the motions of the nodes at the top boundary.All nodes' motions in the direction normal to the plane are constrained.In the tension and compression simulations, the final vertical displacements of the top nodes are 0.09 mm and 0.9 mm, respectively.The total time of loading is 100 s.In the first 20 s and last 80 s, the vertical motions of the top nodes are uniformly accelerated and uniform, respectively.

Mortar Particle and Aggregate Material Model
The linear elastic model is adopted for modeling the behavior of the mortar particles and aggregates.The mortar particles and aggregates are modeled by the six-node linear triangular prism (C3D6) elements.Although the mortar particles and aggregates are assumed to be isotropic in the two-dimensional plane, they need to be defined as anisotropic to eliminate the Poisson effect in the direction normal to the plane.The Poisson ratios normal to the plane are set to zero.
Referring to previous material test results, different elastic moduli are set to the mortar particles and aggregates to simulate the crack propagation caused by local deformation difference.The macroscopic modulus ( m E ) of specimens can be estimated by the following formula: where a E and m E are the moduli of aggregates and mortar, and a P and m P are the volume proportions of aggregates and mortar.The formula is obtained through a short derivation.In the derivation, it is assumed that mortar and aggregates are connected in series.Based on the volume proportions and moduli, Equation (1) was derived.

Aggregate Generation
The aggregates are generated through modeling arbitrary polygon aggregates based on a pregenerated mesh [33].The plane is meshed by triangular grids.The nodes' coordinate information is The models are meshed by irregular triangle elements in the plane (Figure 1).Compared with other polygon elements, triangle elements provide the maximum number of interfaces with the same nodes in the plane.Cohesive elements are inserted into the interfaces to simulate damage and fracture behaviors.
The displacement loading control is adopted to obtain the complete tension and compression process data.The vertical motions of the nodes at the bottom boundary are constrained.The load is applied by controlling the motions of the nodes at the top boundary.All nodes' motions in the direction normal to the plane are constrained.In the tension and compression simulations, the final vertical displacements of the top nodes are 0.09 mm and 0.9 mm, respectively.The total time of loading is 100 s.In the first 20 s and last 80 s, the vertical motions of the top nodes are uniformly accelerated and uniform, respectively.

Mortar Particle and Aggregate Material Model
The linear elastic model is adopted for modeling the behavior of the mortar particles and aggregates.The mortar particles and aggregates are modeled by the six-node linear triangular prism (C3D6) elements.Although the mortar particles and aggregates are assumed to be isotropic in the two-dimensional plane, they need to be defined as anisotropic to eliminate the Poisson effect in the direction normal to the plane.The Poisson ratios normal to the plane are set to zero.
Referring to previous material test results, different elastic moduli are set to the mortar particles and aggregates to simulate the crack propagation caused by local deformation difference.The macroscopic modulus (E m ) of specimens can be estimated by the following formula: where E a and E m are the moduli of aggregates and mortar, and P a and P m are the volume proportions of aggregates and mortar.The formula is obtained through a short derivation.In the derivation, it is assumed that mortar and aggregates are connected in series.Based on the volume proportions and moduli, Equation (1) was derived.

Aggregate Generation
The aggregates are generated through modeling arbitrary polygon aggregates based on a pre-generated mesh [33].The plane is meshed by triangular grids.The nodes' coordinate information is extracted as the basis to model aggregates, and check conflicts and overlaps.The distributions of aggregates accord with Fuller grading cumulative distribution.The function of Fuller grading curve can be decomposed into two dimensions [34] as: where P c (D < D 0 ) denotes the proportion of aggregates whose particle diameter D are smaller than the threshold value D 0 , P k is the ratio of the total volume of aggregates to the concrete volume, and D max is the diameter of the largest aggregate particle.
The geometrical shapes of aggregates are formed by the polygons' vertices.A random function expressed in terms of the polar coordinates r and θ is used to generate the vertices, as shown in Figure 2. To enforce randomicity and ergodicity, random numbers are introduced while generating the radius and angle.The vertices location (r, θ) is defined as follows: where rand(−1,1) is a random number less than 1 and greater than −1, α is the edge number, and f r and f θ are the radius and angle fluctuations respectively.f r and f θ are used to control the irregularity of aggregates.The greater f r and f θ , the more irregular the generated aggregates will be.If f r = 0 and f θ = 0, the generated aggregates are regular polygons.The geometrical shapes of aggregates are formed by the polygons' vertices.A random function expressed in terms of the polar coordinates r and  is used to generate the vertices, as shown in where ( )

Interfacial Transitional Zones (ITZs) and Mortar Internal Interface (MII) Model
The ITZs and MII are modeled by eight-node three-dimensional cohesive elements (CH3D8).Cohesive elements only allow separation in one specific direction (Figure 3).Therefore, cracks can only propagate along the boundaries of particle elements, which may lead to mesh sensitivity [35].In the case of control computing cost, the fine and irregular mesh is generated to provide more path

Interfacial Transitional Zones (ITZs) and Mortar Internal Interface (MII) Model
The ITZs and MII are modeled by eight-node three-dimensional cohesive elements (CH3D8).Cohesive elements only allow separation in one specific direction (Figure 3).Therefore, cracks can only propagate along the boundaries of particle elements, which may lead to mesh sensitivity [35].In the case of control computing cost, the fine and irregular mesh is generated to provide more path choices for the crack propagation and to reduce the mesh sensitivity to some extent [36].The geometrical (constitutive) thickness of cohesive elements is set to zero to maintain the original geometrical partition.The quadratic nominal stress criterion is adopted to define the beginning of an element's degradation, which is presented as follows: In this paper, linear damage evolution is selected to describe the traction-separation law, for its excellent computing efficiency and high accuracy [37].Damage evolution is defined based on the mode I and II fracture energies, which are obtainable from the experiments.Since the model is two-dimensional, the out-of-plane shear failure mode will not happen, and the fracture energy in mode III is not taken into account.Different types of concrete shear experiments have shown that the mode II fracture energy of concrete is about 20 to 25 times larger than mode I fracture energy [22,38,39].The findings are used in the present research.
Since crack propagation paths are meandering curves at the meso-scale, mixed-mode (I + II) fractures are more general than pure mode I or II fractures in concrete failure.The Benzeggagh-Kenane (B-K) fracture criterion [40] is applied to define the mixed-mode fracture energy.It is given by: ( ) where C G is mixed-mode fracture energy, C n G and C s G are mode I and II fracture energies, s G is the sum of shear strain energies in two directions, n G is strain energy in the normal direction, and  is the material parameter.Since tearing mode fracture will not happen in the 2D models, s G is equal to in-plane shear energy.n G and s G of a cohesive element are computed based on the deformation history at integration points in computation, which also relates to the damage evolution of the element [31].The quadratic nominal stress criterion and mixed-mode fracture criterion are shown in Figure 4.In this paper, linear damage evolution is selected to describe the traction-separation law, for its excellent computing efficiency and high accuracy [37].
Damage evolution is defined based on the mode I and II fracture energies, which are obtainable from the experiments.Since the model is two-dimensional, the out-of-plane shear failure mode will not happen, and the fracture energy in mode III is not taken into account.Different types of concrete shear experiments have shown that the mode II fracture energy of concrete is about 20 to 25 times larger than mode I fracture energy [22,38,39].The findings are used in the present research.
Since crack propagation paths are meandering curves at the meso-scale, mixed-mode (I + II) fractures are more general than pure mode I or II fractures in concrete failure.The Benzeggagh-Kenane (B-K) fracture criterion [40] is applied to define the mixed-mode fracture energy.It is given by: where G C is mixed-mode fracture energy, G C n and G C s are mode I and II fracture energies, G s is the sum of shear strain energies in two directions, G n is strain energy in the normal direction, and η is the material parameter.Since tearing mode fracture will not happen in the 2D models, G s is equal to in-plane shear energy.G n and G s of a cohesive element are computed based on the deformation history at integration points in computation, which also relates to the damage evolution of the element [31].The quadratic nominal stress criterion and mixed-mode fracture criterion are shown in Figure 4.

Interaction
Since particles bump and grind during compression, it is necessary to define the interaction characteristics to simulate the collision and friction.It is assumed that when the particle surfaces are in contact, any contact pressure can be transmitted.The 'hard' contact is adopted to describe the normal behavior of the interaction.According to the Coulomb friction mechanism, the friction behavior influences ultimate compressive strength and concrete performance in the softening stage [41].The Coulomb friction model is adopted to describe the tangential behavior of the interaction.As with geomaterial, the friction behavior of concrete can be depicted by the friction angle ϕ [42].The angle can be measured by the direct shear test and uniaxial compression test [43].With the increase of concrete strength, the angle grows slightly.The friction angle ranges from 50 to 60 degrees.In the modeling, the friction angle ϕ is transformed into the friction coefficient µ.The transformation formula is shown as follows:

Interaction
Since particles bump and grind during compression, it is necessary to define the interaction characteristics to simulate the collision and friction.It is assumed that when the particle surfaces are in contact, any contact pressure can be transmitted.The 'hard' contact is adopted to describe the normal behavior of the interaction.According to the Coulomb friction mechanism, the friction behavior influences ultimate compressive strength and concrete performance in the softening stage [41].The Coulomb friction model is adopted to describe the tangential behavior of the interaction.As with geomaterial, the friction behavior of concrete can be depicted by the friction angle  [42].The angle can be measured by the direct shear test and uniaxial compression test [43].With the increase of concrete strength, the angle grows slightly.The friction angle ranges from 50 to 60 degrees.In the modeling, the friction angle  is transformed into the friction coefficient  .The transformation formula is shown as follows: =tan

Calibration of Numerical Simulation
Modeling concrete, which contains multiple materials and interfaces at the meso-scale, requires a large number of parameters.Few experiments can provide the corresponding parameters that can one-to-one match with those used in the numerical simulation presented in this paper.Moreover, the cohesive zone model has a mesh sensitivity problem, so calibration is necessary for obtaining a suitable set of parameters and mesh precision.
In the calibration, it is found that the specimens with the same grading cumulative distribution but different aggregate distributions present different compressive properties, especially during the softening stage.Therefore, it is inadequate to include a single specimen to study the entire concrete performances.In order to conduct universal studies, ten specimens with different aggregate distributions are modeled randomly according to Equation ( 2) and (3).The number and proportion of aggregates included in the specimens are summarized in Table 1.
Illustration of quadratic nominal stress criterion and mixed-mode fracture criterion.

Calibration of Numerical Simulation
Modeling concrete, which contains multiple materials and interfaces at the meso-scale, requires a large number of parameters.Few experiments can provide the corresponding parameters that can one-to-one match with those used in the numerical simulation presented in this paper.Moreover, the cohesive zone model has a mesh sensitivity problem, so calibration is necessary for obtaining a suitable set of parameters and mesh precision.
In the calibration, it is found that the specimens with the same grading cumulative distribution but different aggregate distributions present different compressive properties, especially during the softening stage.Therefore, it is inadequate to include a single specimen to study the entire concrete performances.In order to conduct universal studies, ten specimens with different aggregate distributions are modeled randomly according to Equation ( 2) and (3).The number and proportion of aggregates included in the specimens are summarized in Table 1.To check the mesh dependence and computing efficiency, three models with the same aggregate distribution and different mesh precisions are tested.The approximate element sizes of the models are 1.5 mm, 1.0 mm, 0.5 mm (Figure 5).The number of bulk elements and cohesive elements are listed in Table 2.The time consumption of the models with element size of 1.0 mm and 0.5 mm are 2.5 and 11.6 times that of the model with 1.5 mm element size, respectively.To check the mesh dependence and computing efficiency, three models with the same aggregate distribution and different mesh precisions are tested.The approximate element sizes of the models are 1.5 mm, 1.0 mm, 0.5 mm (Figure 5).The number of bulk elements and cohesive elements are listed in Table 2.The time consumption of the models with element size of 1.0 mm and 0.5 mm are 2.5 and 11.6 times that of the model with 1.5 mm element size, respectively.Figure 6 shows the stress-strain curves of the specimens with different element sizes.Comparing the tension and compression simulation results, it is shown that the compressive performance is more sensitive to mesh precision.It is indicated that the damage and fracture of cohesive elements lead to stress redistribution.With the increase of plastic deformation ability in the local area, the stress distribution becomes more uniform, which can make the particles stay at a higher average stress level.Hence, the increase of cohesive element number can contribute to higher ultimate compressive strength and higher residual strength in the softening stage.The mesh size of 1.0 mm is chosen for further experiments due to its good simulation accuracy and acceptable time consumption.
Based on the outcomes of previous research and the simulation technique discussed in Section 2, the models are calibrated by considering interaction relations, material and interface parameters.Through the calibration, the frictional angle of 54.4° and the material parameters listed in Table 3 can make the models achieve satisfactory results, which are in close agreement with experimental results in both mechanical response and failure mode.In the following, the standard models used for analysis and parametric study are the models with ten different aggregate distributions and the calibrated parameters.Figure 6 shows the stress-strain curves of the specimens with different element sizes.Comparing the tension and compression simulation results, it is shown that the compressive performance is more sensitive to mesh precision.It is indicated that the damage and fracture of cohesive elements lead to stress redistribution.With the increase of plastic deformation ability in the local area, the stress distribution becomes more uniform, which can make the particles stay at a higher average stress level.Hence, the increase of cohesive element number can contribute to higher ultimate compressive strength and higher residual strength in the softening stage.The mesh size of 1.0 mm is chosen for further experiments due to its good simulation accuracy and acceptable time consumption.Based on the outcomes of previous research and the simulation technique discussed in Section 2, the models are calibrated by considering interaction relations, material and interface parameters.Through the calibration, the frictional angle of 54.4 • and the material parameters listed in Table 3 can make the models achieve satisfactory results, which are in close agreement with experimental results in both mechanical response and failure mode.In the following, the standard models used for analysis and parametric study are the models with ten different aggregate distributions and the calibrated parameters.

Tension Simulation Result
The complete tensile stress-strain curves of the specimens are shown in Figure 7.A mean curve is plotted to show the integral trend of the curves.Assume that σ i (ε) is the stress of i-th specimen at the global strain ε.The points (σ, ε) on the mean curve can be calculated as follows: As shown in Figure 8, there are two typical tension failure modes: the specimen fracture with one or two dominant cracks perpendicular to the tension direction.In the failure mode with only one dominant crack, the crack extends through the whole cross-section as the experimental result.The majority of the tension failures are this mode.In the other failure mode, the two dominant cracks at different horizontal positions extend from both sides.With the crack propagation, the horizontal projections of the two cracks intersect each other.Afterward, one of the crack's propagation rates declines, and the other crack extends through the whole cross-section.Figure 9 shows the strain-stress curves of the specimens with one or two dominant cracks.The strength and energy dissipation of the specimen with two dominant cracks is higher than that of the specimen with one dominant crack Through the calibration, the mean curve agrees well with the curve of the experimental result [32].The curves of standard deviation and standard deviation coefficient that are calculated with simulation data are also plotted in Figure 7.The curves represent the absolute and relative errors of the simulation data, respectively.It is worth noting that the standard deviation almost remains constant in the softening stage.The stable deviation indicates that the specimens have a close damage evolution rate at the same global deformation level in spite of their different aggregate distributions.
As shown in Figure 8, there are two typical tension failure modes: the specimen fracture with one or two dominant cracks perpendicular to the tension direction.In the failure mode with only one dominant crack, the crack extends through the whole cross-section as the experimental result.The majority of the tension failures are this mode.In the other failure mode, the two dominant cracks at different horizontal positions extend from both sides.With the crack propagation, the horizontal projections of the two cracks intersect each other.Afterward, one of the crack's propagation rates declines, and the other crack extends through the whole cross-section.Figure 9 shows the strain-stress curves of the specimens with one or two dominant cracks.The strength and energy dissipation of the specimen with two dominant cracks is higher than that of the specimen with one dominant crack As shown in Figure 8, there are two typical tension failure modes: the specimen fracture with one or two dominant cracks perpendicular to the tension direction.In the failure mode with only one dominant crack, the crack extends through the whole cross-section as the experimental result.The majority of the tension failures are this mode.In the other failure mode, the two dominant cracks at different horizontal positions extend from both sides.With the crack propagation, the horizontal projections of the two cracks intersect each other.Afterward, one of the crack's propagation rates declines, and the other crack extends through the whole cross-section.Figure 9 shows the strain-stress curves of the specimens with one or two dominant cracks.The strength and energy dissipation of the specimen with two dominant cracks is higher than that of the specimen with one dominant crack The specimen with one dominant crack in Figure 8 is selected for the deformation process analysis.In order to observe the whole deformation process, the deformation of the specimen is enlarged geometrically, as shown in Figure 10.The specimen's original deformation in the direction of tension is multiplied by a deformation scale factor v D to make the aspect ratio become 2:1 from 1:1.Eight characteristic states are selected to describe the deformation process.The corresponding points of the states are plotted on the stress-strain curve in Figure 9.
After the specimen is loaded, before state A, the deformation evenly distributes at the aggregates and mortar, as shown in Figure 10a.The obvious plastic deformation appears since state B. Although the global stress hasn't reached the damage initial stress of the ITZ elements, a few ITZ elements The specimen with one dominant crack in Figure 8 is selected for the deformation process analysis.In order to observe the whole deformation process, the deformation of the specimen is enlarged geometrically, as shown in Figure 10.The specimen's original deformation in the direction of tension is multiplied by a deformation scale factor D v to make the aspect ratio become 2:1 from 1:1.Eight characteristic states are selected to describe the deformation process.The corresponding points of the states are plotted on the stress-strain curve in Figure 9.
After the specimen is loaded, before state A, the deformation evenly distributes at the aggregates and mortar, as shown in Figure 10a.The obvious plastic deformation appears since state B. Although the global stress hasn't reached the damage initial stress of the ITZ elements, a few ITZ elements begin sustaining damage due to the stress concentration caused by the deformation difference between mortar and aggregates.
It is obvious from Figure 10c that the deformation concentrates at ITZ elements and nearby MII elements at state C.Many deformation concentrated areas are formed as spindle-shaped, which are comprised by ITZ elements at mid-piece and MII elements at both ends.When the specimen attained the ultimate tensile strength state, two potential fracture zones (FZ1 and FZ2) are formed by the deformation concentrated areas, as shown in Figure 10d.
the ultimate tensile strength state, two potential fracture zones (FZ1 and FZ2) are formed by the deformation concentrated areas, as shown in Figure 10d.
At state E, one of the potential fracture zones (FZ1) degrades, and another one (FZ2) extends unstably.At state F, the first fracture initiates at the center-left.The deformation of MII elements is very close to that of the nearby ITZ elements.The original spindle-shaped characteristic of the deformation areas has almost disappeared.Although the cracks develop and extend quickly afterward, the drop rate of bearing capacity slows down, and the curve gradually flattens.At state G, a large number of ITZ elements are eliminated.As can be noticed, a new deformation concentrated area (FZ4) has been generated at the bottom right corner of the potential fracture zone during the stage between state F and G, while the potential fracture zone (FZ3), which was generated since state D, begins contracting.At state H, fracture happens in this new deformation concentrated area (FZ4).This illustrates that although the tension crack generally propagates through the main deformation concentrated areas, the crack would bypass parts of areas and propagate through a new fracture path due to the force redistribution.

Compression Simulation Result
The aggregate distribution has a substantial influence on the concrete compressive property, especially in the softening stage.Compared with the tension simulation result, the result of compression simulation shows greater discreteness.The tensile behavior is determined by the weakest section, while compressive behavior is determined by the entire specimen [41].
The crack initiation is caused by the localized shear failure, as shown in Figure 11.There are two to four main shear bands crossing the whole section in the middle region accompanied by several small bands on both sides.The angles between the shear bands and the direction of compression range from 5 to 15 degrees.The specimens are divided into several wedge-shaped blocks by the shear At state E, one of the potential fracture zones (FZ1) degrades, and another one (FZ2) extends unstably.At state F, the first fracture initiates at the center-left.The deformation of MII elements is very close to that of the nearby ITZ elements.The original spindle-shaped characteristic of the deformation areas has almost disappeared.
Although the cracks develop and extend quickly afterward, the drop rate of bearing capacity slows down, and the curve gradually flattens.At state G, a large number of ITZ elements are eliminated.As can be noticed, a new deformation concentrated area (FZ4) has been generated at the bottom right corner of the potential fracture zone during the stage between state F and G, while the potential fracture zone (FZ3), which was generated since state D, begins contracting.At state H, fracture happens in this new deformation concentrated area (FZ4).This illustrates that although the tension crack generally propagates through the main deformation concentrated areas, the crack would bypass parts of areas and propagate through a new fracture path due to the force redistribution.

Compression Simulation Result
The aggregate distribution has a substantial influence on the concrete compressive property, especially in the softening stage.Compared with the tension simulation result, the result of compression simulation shows greater discreteness.The tensile behavior is determined by the weakest section, while compressive behavior is determined by the entire specimen [41].
The crack initiation is caused by the localized shear failure, as shown in Figure 11.There are two to four main shear bands crossing the whole section in the middle region accompanied by several small bands on both sides.The angles between the shear bands and the direction of compression range from 5 to 15 degrees.The specimens are divided into several wedge-shaped blocks by the shear bands as the experimental result [32].
bottom right corner of the potential fracture zone during the stage between state F and G, while the potential fracture zone (FZ3), which was generated since state D, begins contracting.At state H, fracture happens in this new deformation concentrated area (FZ4).This illustrates that although the tension crack generally propagates through the main deformation concentrated areas, the crack would bypass parts of areas and propagate through a new fracture path due to the force redistribution.

Compression Simulation Result
The aggregate distribution has a substantial influence on the concrete compressive property, especially in the softening stage.Compared with the tension simulation result, the result of compression simulation shows greater discreteness.The tensile behavior is determined by the weakest section, while compressive behavior is determined by the entire specimen [41].
The crack initiation is caused by the localized shear failure, as shown in Figure 11.There are two to four main shear bands crossing the whole section in the middle region accompanied by several small bands on both sides.The angles between the shear bands and the direction of compression range from 5 to 15 degrees.The specimens are divided into several wedge-shaped blocks by the shear bands as the experimental result [32].The complete compressive stress-strain curves are shown in Figure 12.As with the statistical characteristics of tension stress-curves, the datum maintains a low dispersion before the mean compressive stress reaches its peak.Remarkably, both the curves of standard deviation and standard deviation coefficient have a 'valley' at the specimens' softening stage.In the descending branch, the stress-strain curves of the specimens are pinched together at a certain state.The 'valley' symbolizes that there is a relatively steady state in the softening stage, at which the specimens have a close residual bearing capacity.
Figure 13 shows the complete compressive stress-strain curve of the model with the same aggregate distribution as the one analyzed in Section 4. Seven states labeled from 'A' to 'G' on the curve are selected to describe the crack propagation process as shown in Figure 14.The cracks newly appearing in each state are differentiated by colors.
The sums of the lengths of ITZ and MII elements with different damage degrees are accumulated, as shown in Figure 15.SDEG (overall value of the scalar damage variable) values of the elements are extracted from the result file to determine the damage degree.When an element's SDEG The complete compressive stress-strain curves are shown in Figure 12.As with the statistical characteristics of tension stress-curves, the datum maintains a low dispersion before the mean compressive stress reaches its peak.Remarkably, both the curves of standard deviation and standard deviation coefficient have a 'valley' at the specimens' softening stage.In the descending branch, the stress-strain curves of the specimens are pinched together at a certain state.The 'valley' symbolizes that there is a relatively steady state in the softening stage, at which the specimens have a close residual bearing capacity. Figure 13 shows the complete compressive stress-strain curve of the model with the same aggregate distribution as the one analyzed in Section 4. Seven states labeled from 'A' to 'G' on the curve are selected to describe the crack propagation process as shown in Figure 14.The cracks newly appearing in each state are differentiated by colors.The sums of the lengths of ITZ and MII elements with different damage degrees are accumulated, as shown in Figure 15.SDEG (overall value of the scalar damage variable) values of the elements are extracted from the result file to determine the damage degree.When an element's SDEG value equals 1, the element would be eliminated to simulate visible fracture.The lengths of the elements are calculated by the coordinates of the elements' nodes.Since the nonlinear behavior of the specimen is simulated by the damage of ITZ and MII elements, the total length of the elements represents a specimen's damage condition.It is interesting that the valleys of the standard deviation curve and standard deviation coefficient curve occur when the specimen is at state E, which indicates that the relatively steady state is attributed to the formation of main shear bands.At this state, the residual capacity depends on the friction action of the blocks divided by the shear bands.In the following loading process, the crack propagation rate slows down.The specimen's deformation gradually evolves into the slipping of the wedge-shaped blocks.

The Effect of Mode II Fracture Energy
Since mode I fracture was considered as the dominant form in tension damage, single fracture energy was adopted to describe the interface fracture behavior in the previous concrete tension simulation.In this method, the fracture energies of the interface elements are fixed values regardless of fracture type.Moreover, since the compression failure of concrete is mostly due to shear fracture [41], modeling with one specific fracture energy is unable to simulate the fracture behavior of concrete both under tension and compression.Due to the multiphase and heterogeneous property of concrete, the mixed-mode fracture would occur inevitably in concrete failure, even under uniaxial tension conditions.Mixed-mode fracture energy is determined by mode I and mode II fracture energy, and the two pure-mode fracture energies both have influences on concrete tensile and compressive behavior.In this section, the effect of mode II fracture energy on the performance of concrete under tension and is studied.
As mentioned in Section 2, the experiment results show that there is a relationship between mode I and II fracture energy.In this section, mode I fracture energy is fixed.The ratio of mode II and mode I fracture energy (  ) is taken as the variable to investigate the effect of mode II fracture energy.The ratio  is calculated by the following equation: The tension simulation result of varying  is shown in Figure 16a.The influence on the ultimate tensile strength is negligible.It is noteworthy that the mode II fracture energy influences the specimen's tensile behavior in the softening stage.The mixed-mode fracture accounts for a large proportion.In the case that the mode II fracture energy equals the mode I fracture energy ( =1  ), the total fracture energy is 48.2% lower than that of the standard model ( =25


).It can be inferred that the mixed-mode fracture energy caused by the difference of the two pure-mode fracture energy accounts for about half of the total fracture energy.The specimen is linear-elastic until the first damage occurs.Damage firstly occurs at ITZ elements, and then MII elements.The damage evolution rate of the specimen gradually accelerates and achieves a constant speed with the deformation increased.Although the ultimate bearing capacity hasn't been reached, the crack initiates at state A. The crack nucleates at the ITZ elements parallel to the loading direction and propagates slowly before the state B.
Afterwards, several small cracks appear from the bottom middle-right to the top-right corner during the stage between state B and C. The bearing capacity gradually declines at this stage.A rudiment of a shear band seems to be formed by these small cracks.
The residual bearing capacity descends rapidly during the stage between states C and E, accompanied by the appearance of main shear bands.The MII elements' fracture leads to the branching and coalescence of the crack.In this stage, the damage evolution rate is fast but steady.
It is interesting that the valleys of the standard deviation curve and standard deviation coefficient curve occur when the specimen is at state E, which indicates that the relatively steady state is attributed to the formation of main shear bands.At this state, the residual capacity depends on the friction action of the blocks divided by the shear bands.In the following loading process, the crack propagation rate slows down.The specimen's deformation gradually evolves into the slipping of the wedge-shaped blocks.

The Effect of Mode II Fracture Energy
Since mode I fracture was considered as the dominant form in tension damage, single fracture energy was adopted to describe the interface fracture behavior in the previous concrete tension simulation.In this method, the fracture energies of the interface elements are fixed values regardless of fracture type.Moreover, since the compression failure of concrete is mostly due to shear fracture [41], modeling with one specific fracture energy is unable to simulate the fracture behavior of concrete both under tension and compression.Due to the multiphase and heterogeneous property of concrete, the mixed-mode fracture would occur inevitably in concrete failure, even under uniaxial tension conditions.Mixed-mode fracture energy is determined by mode I and mode II fracture energy, and the two pure-mode fracture energies both have influences on concrete tensile and compressive behavior.
In this section, the effect of mode II fracture energy on the performance of concrete under tension and compression is studied.
As mentioned in Section 2, the experiment results show that there is a relationship between mode I and II fracture energy.In this section, mode I fracture energy is fixed.The ratio of mode II and mode I fracture energy (ξ) is taken as the variable to investigate the effect of mode II fracture energy.The ratio ξ is calculated by the following equation: The tension simulation result of varying ξ is shown in Figure 16a.The influence on the ultimate tensile strength is negligible.It is noteworthy that the mode II fracture energy influences the specimen's tensile behavior in the softening stage.The mixed-mode fracture accounts for a large proportion.In the case that the mode II fracture energy equals the mode I fracture energy (ξ = 1), the total fracture energy is 48.2% lower than that of the standard model (ξ = 25).It can be inferred that the mixed-mode fracture energy caused by the difference of the two pure-mode fracture energy accounts for about half of the total fracture energy.The compression simulation result of varying the multiple is shown in Figure 16b.The ultimate compressive strength is significantly influenced by the mode II fracture energy, which is different from tensile strength.When the mode II fracture energy equals the mode I fracture energy ( =1  ), the compressive strength and total fracture energy is 44.7% and 82.2% lower than that of the standard model ( =25  ), respectively.As described in Sections 3 and 4, the crack initiates before and after the specimen reaches the ultimate bearing capacity state under compression and tension, respectively.Low mode II fracture energy will limit interface elements' shear deformation under compression, while it has a relatively low influence on the interface elements' tensile deformation under tension.It can be noticed that with the increase of mode II fracture energy, the growth of compressive strength decreases.The limit of compressive strength gradually translates from deformation ability to the shear strength of interface elements.

The Effect of Mixed-Mode Fracture Criterion
The concrete's fracture is generally mixed-mode, especially under compression.The Benzeggagh-Kenane (B-K) fracture criterion with the material parameter  is adopted to define the damage evolution of mixed-mode fracture.It is accepted that the B-K criterion can adequately specify the mixed-mode fracture energy for a wide range of composite material [44].
In this section, the material parameter  is taken as the variable to investigate the effect of mixed-mode fracture criterion.The other parameters are held the same, as listed in Table 3.Since there are limited studies of the material parameter  of concrete, the parameter value of brittle epoxy resin is taken as a reference, which ranges from 1.0 to 2.0.Through calibration, it seems that 1.2 is an appropriate value for the parameter  under both compression and tension conditions.The compression simulation result of varying the multiple is shown in Figure 16b.The ultimate compressive strength is significantly influenced by the mode II fracture energy, which is different from tensile strength.When the mode II fracture energy equals the mode I fracture energy (ξ = 1), the compressive strength and total fracture energy is 44.7% and 82.2% lower than that of the standard model (ξ = 25), respectively.As described in Sections 3 and 4, the crack initiates before and after the specimen reaches the ultimate bearing capacity state under compression and tension, respectively.Low mode II fracture energy will limit interface elements' shear deformation under compression, while it has a relatively low influence on the interface elements' tensile deformation under tension.It can be noticed that with the increase of mode II fracture energy, the growth of compressive strength decreases.The limit of compressive strength gradually translates from deformation ability to the shear strength of interface elements.

The Effect of Mixed-Mode Fracture Criterion
The concrete's fracture is generally mixed-mode, especially under compression.The Benzeggagh-Kenane (B-K) fracture criterion with the material parameter η is adopted to define the damage evolution of mixed-mode fracture.It is accepted that the B-K criterion can adequately specify the mixed-mode fracture energy for a wide range of composite material [44].
In this section, the material parameter η is taken as the variable to investigate the effect of mixed-mode fracture criterion.The other parameters are held the same, as listed in Table 3.Since there are limited studies of the material parameter η of concrete, the parameter value of brittle epoxy resin is taken as a reference, which ranges from 1.0 to 2.0.Through calibration, it seems that 1.2 is an appropriate value for the parameter η under both compression and tension conditions.
Figure 17 shows the complete tensile and compressive stress-strain mean curves for varying the parameter η.The results indicate that the material parameter mainly affects the concrete performance in the softening stage under both tension and compression.The ultimate tensile and compressive strengths are slightly affected.It seems that the curves separate at the peak and gradually tend to decrease in parallel with the increase of global strain.For the specimens with the same aggregate distribution, the crack propagation paths are similar under tension, as are the failure modes and the sum of crack length under compression, as shown in Figure 8.According to the B-K criterion expression (Equation ( 5)), if the ratio of strain energy at a tangent to that in the normal direction is constant, a decrease of material parameter η leads to the increase of total fracture energy G C .Due to the similarity in failure modes, the bearing capacities of the specimens are mainly determined by the fracture energy.

The Effect of Frictional Behavior
Coulombic friction is considered a key factor that influences the ultimate strength and residual bearing capacity of concrete.The results of the direct shear test and uniaxial compression test showed that the concrete with higher compressive strength has a greater internal friction angle.Although a desirable strength of concrete specimens can be achieved in experiments by adjusting the material types and aggregate gradation, it is hard to control the maximum allowable frictional stress.Therefore, the frictional angle α is taken as the variable to investigate the effect of frictional behavior.
As shown in Figure 18, the internal frictional behavior of concrete affects the final failure mode.With the increase of frictional angle, the inclination angles of main shear bands decrease.The local cracks tend to propagate dispersedly into concrete blocks rather than in the direction of the shear bands.The stress-strain curves of varying frictional angle α are shown in Figure 19.As expected, the increase of frictional angle value has a positive influence on the concrete bearing capacity in the whole loading process.The influence strengthens once the plastic deformation appears and weakens when the specimens' bearing capacity declines by about a third.The sum of the crack length of the specimens is calculated as shown in Figure 20.
When the tangential interaction is frictionless (α = 0), the compressive strength and energy dissipations are 35.9% and 77.1% lower than that of the standard model (α = 54.4°),respectively.It is noticed that the crack propagates more sufficiently and the inclination angle of the shear bands decreases with the increase of frictional angle.Although the cracks' local propagation paths are different, the main shear bands of the specimens are similar.The friction stress could constrain the crack propagation at the beginning of the damage, which benefits the ultimate bearing capacity.In the softening stage, the friction stress could help reestablish the loading path and make more blocks participate in bearing load, which may cause more cracks.

The Effect of Frictional Behavior
Coulombic friction is considered a key factor that influences the ultimate strength and residual bearing capacity of concrete.The results of the direct shear test and uniaxial compression test showed that the concrete with higher compressive strength has a greater internal friction angle.Although a desirable strength of concrete specimens can be achieved in experiments by adjusting the material types and aggregate gradation, it is hard to control the maximum allowable frictional stress.Therefore, the frictional angle α is taken as the variable to investigate the effect of frictional behavior.
As shown in Figure 18, the internal frictional behavior of concrete affects the final failure mode.With the increase of frictional angle, the inclination angles of main shear bands decrease.The local cracks tend to propagate dispersedly into concrete blocks rather than in the direction of the shear bands.The stress-strain curves of varying frictional angle α are shown in Figure 19.As expected, the increase of frictional angle value has a positive influence on the concrete bearing capacity in the whole loading process.The influence strengthens once the plastic deformation appears and weakens when the specimens' bearing capacity declines by about a third.The sum of the crack length of the specimens is calculated as shown in Figure 20.
noticed that the crack propagates more sufficiently and the inclination angle of the shear bands decreases with the increase of frictional angle.Although the cracks' local propagation paths are different, the main shear bands of the specimens are similar.The friction stress could constrain the crack propagation at the beginning of the damage, which benefits the ultimate bearing capacity.In the softening stage, the friction stress could help reestablish the loading path and make more blocks participate in bearing load, which may cause more cracks.

The Effect of Aggregates' Mechanical Characteristics
Although the tensile behavior of concrete is dominated by the properties of ITZs and mortar, damage and fracture may happen in aggregates during compression [45].Different kinds of

The Effect of Aggregates' Mechanical Characteristics
Although the tensile behavior of concrete is dominated by the properties of ITZs and mortar, damage and fracture may happen in aggregates during compression [45].Different kinds of

The Effect of Aggregates' Mechanical Characteristics
Although the tensile behavior of concrete is dominated by the properties of ITZs and mortar, damage and fracture may happen in aggregates during compression [45].Different kinds of aggregates have differences in mechanical characteristics, geometrical shape, and corresponding ITZ properties [46].To study the effect of aggregates' mechanical characteristics on concrete compressive When the tangential interaction is frictionless (α = 0), the compressive strength and energy dissipations are 35.9% and 77.1% lower than that of the standard model (α = 54.4 • ), respectively.It is noticed that the crack propagates more sufficiently and the inclination angle of the shear bands decreases with the increase of frictional angle.Although the cracks' local propagation paths are different, the main shear bands of the specimens are similar.The friction stress could constrain the crack propagation at the beginning of the damage, which benefits the ultimate bearing capacity.In the softening stage, the friction stress could help reestablish the loading path and make more blocks participate in bearing load, which may cause more cracks.

The Effect of Aggregates' Mechanical Characteristics
Although the tensile behavior of concrete is dominated by the properties of ITZs and mortar, damage and fracture may happen in aggregates during compression [45].Different kinds of aggregates have differences in mechanical characteristics, geometrical shape, and corresponding ITZ properties [46].To study the effect of aggregates' mechanical characteristics on concrete compressive performance, based on standard models, cohesive elements are inserted into aggregates to simulate aggregate internal interfaces (AII).The material parameters of AII elements are set by referring to the material tests of the following common aggregate rocks: limestone [47], marble [48], granite [49], and basalt [50].The material parameters of AII elements are listed in Table 4.The stress-strain curves of varying aggregates' mechanical characteristics are shown in Figure 21.Compared with the ultimate bearing capacity of the models without considering aggregate damage and fracture, the ultimate bearing capacities of the models with limestone, marble, granite, basalt aggregates fall by averages of 1.1%, 1.0%, 0.23%, 0.18%, respectively.The effect of aggregate type on the stress-strain relationship is mainly reflected in the softening stage.As shown in Figure 22, damage evolution of aggregates accelerates in the early softening stage.Since marble aggregates have the lowest tensile and shear strength, they are more likely to be damaged during the compression process.The energy dissipation of aggregates can lead to stress redistribution at the local area and enhance residual strength.
Figure 23 shows the failure modes of the models with the same aggregate distribution and different aggregates' mechanical characteristics.It can be seen that the fracture of aggregates only constitutes a minority, but the crack propagation paths at the models' left half are obviously influenced by local crack propagation differences caused by aggregates' fracture.Since the shear bands at the models' right half form earlier than those at left half, the models have very similar failure modes (FZ5-FZ9) at the right half.In view of these observations, it can be concluded that the influence of aggregates' mechanical characteristics on concrete's mechanical characteristics and failure modes is mainly concentrated in the mid-and-late softening stage.
different aggregates' mechanical characteristics.It can be seen that the fracture of aggregates only constitutes a minority, but the crack propagation paths at the models' left half are obviously influenced by local crack propagation differences caused by aggregates' fracture.Since the shear bands at the models' right half form earlier than those at left half, the models have very similar failure modes (FZ5-FZ9) at the right half.In view of these observations, it can be concluded that the influence of aggregates' mechanical characteristics on concrete's mechanical characteristics and failure modes is mainly concentrated in the mid-and-late softening stage.

Conclusions
The application range of the cohesive zone model for simulation of three-phase concrete (i.e., aggregate, mortar, and ITZs) has been extended.Through the definition of mixed-mode fracture and interaction behavior, the model is suitable for simulation both under tension and compression.The modeling method and relevant material parameters are discussed based on the simulation technique and previous material experiment results.A balance between mesh sensitivity and computing efficiency has been achieved through the calibration of mesh precision.The simulation results have

Conclusions
The application range of the cohesive zone model for simulation of three-phase concrete (i.e., aggregate, mortar, and ITZs) has been extended.Through the definition of mixed-mode fracture and

Conclusions
The application range of the cohesive zone model for simulation of three-phase concrete (i.e., aggregate, mortar, and ITZs) has been extended.Through the definition of mixed-mode fracture and interaction behavior, the model is suitable for simulation both under tension and compression.The modeling method and relevant material parameters are discussed based on the simulation technique and previous material experiment results.A balance between mesh sensitivity and computing efficiency has been achieved through the calibration of mesh precision.The simulation results have good agreement with the experiments in both mechanical characteristics and failure modes.The whole process analyses and parametric studies are conducted to investigate the failure mechanism of concrete under tension and compression.The following conclusions can be drawn from this study.
The ITZ and MII elements modeled by cohesive elements are appropriate to simulate concrete meso-scale fracture behavior and damage accumulation at a smaller scale.The mesh precision level in which the approximate element size is 1/150 of the size of the model can provide a relatively high simulation accuracy and calculation efficiency.
The discreteness of compressive behavior caused by aggregate distribution on the stress-strain curve is larger than that of tension behavior.In the tension simulation, the specimens with different aggregate distributions have similar damage evolution rates at the same global deformation level.The failure modes of concrete can be separated into two types: one or two dominant cracks perpendicular to the tension direction.
In the concrete compression simulation, the inclination angles of shear bands range from 5 to 15 degrees.A valley in the standard deviation coefficient curve is formed, which symbolizes a relatively steady state in the softening stage when main shear bands form.
In the parametric study, it is found that the mixed-mode fracture accounts for a considerable proportion both in tension and compression failure.The fracture energy caused by the difference of the two pure-mode fracture energy occupies 48.2% and 82.2% of the total fracture energy under tension and compression, respectively.Mode II fracture energy has a significant impact on both ultimate and residual compressive strength of concrete.In the case that mode II fracture energy is equal to mode I fracture energy, the ultimate compressive strength reduces by 44.7%.
The mixed-mode fracture criterion mainly affects concrete performance in the softening stage under both tension and compression, while the effects on failure mode, tensile and compressive strength are negligible.
Frictional behavior mainly affects concrete compressive performance.In the case that there is no internal frictional behavior, the ultimate compressive strength and energy dissipation reduce by 35.9% and 77.1%, respectively.The frictional stress can constrain crack propagation at the beginning of damage and reestablish the loading path in the softening stage.
If the influence of aggregates' geometrical shapes and ITZ properties are subtracted, aggregates' mechanical characteristics affect concrete's mechanical characteristics and failure modes, mainly concentrated in the mid-and-late softening stage.

Figure 1 .
Figure 1.Element mesh of model: (a) bulk elements of aggregate (grey colored) and mortar particle (white colored); (b) cohesive elements of ITZs (interfacial transitional zones, red colored) and mortar internal interface (blue coloured).

Figure 1 .
Figure 1.Element mesh of model: (a) bulk elements of aggregate (grey colored) and mortar particle (white colored); (b) cohesive elements of ITZs (interfacial transitional zones, red colored) and mortar internal interface (blue coloured).
 denotes the proportion of aggregates whose particle diameter D are smaller than the threshold value 0 D , k P is the ratio of the total volume of aggregates to the concrete volume, and max D is the diameter of the largest aggregate particle.

Figure 2 .
Figure 2. To enforce randomicity and ergodicity, random numbers are introduced while generating the radius and angle.The vertices location ( , ) r  is defined as follows:

r f and f 
are used to control the irregularity of aggregates.The greater r f and f  , the more irregular the generated aggregates will be.If 0 r f = and 0 f  = , the generated aggregates are regular polygons.

Figure 4 .
Figure 4. Illustration of quadratic nominal stress criterion and mixed-mode fracture criterion.

Figure 5 .
Figure 5.The meshes with different approximate element sizes: (a) boundary constraint condition marked by the blue triangles and the zoom area marked by the red square; (b) approximate size of 1.5 mm; (c) approximate size of 1.0 mm; (d) approximate size of 0.5 mm.

Figure 5 .
Figure 5.The meshes with different approximate element sizes: (a) boundary constraint condition marked by the blue triangles and the zoom area marked by the red square; (b) approximate size of 1.5 mm; (c) approximate size of 1.0 mm; (d) approximate size of 0.5 mm.

Figure 6 .
Figure 6.Effect of element size on global stress-strain relation under (a) tension and (b) compression.

Figure 6 .
Figure 6.Effect of element size on global stress-strain relation under (a) tension and (b) compression.

)Figure 7 .
Figure 7. Statistical characteristic of global stress-strain curve under tension.

Figure 7 .
Figure 7. Statistical characteristic of global stress-strain curve under tension.

Figure 8 .
Figure 8. Failure modes of the specimens under tension: (a) simulation result (with one dominant crack); (b) simulation result (with two dominant cracks).

Figure 10 .
Figure 10.Deformation magnification in loading direction of a specimen in eight states.Cohesive elements of ITZs (blue coloured) and mortar internal interface (red coloured).

FZ4Figure 10 .
Figure 10.Deformation magnification in loading direction of a specimen in eight states.Cohesive elements of ITZs (blue coloured) and mortar internal interface (red coloured).

Figure 11 .
Figure 11.Failure modes of the specimens with different aggregate distributions under compression

Figure 11 .
Figure 11.Failure modes of the specimens with different aggregate distributions under compression

Figure 14 .
Figure 14.Crack propagation process of a specimen: (a) State A, B, C; (b) State D, E; (c) State F, G.

Figure 15 .
Figure 15.The relationship between global strain and the sum of length of (a) ITZ elements and (b) MII elements with different damage degree.

Figure 15 .
Figure 15.The relationship between global strain and the sum of length of (a) ITZ elements and (b) MII elements with different damage degree.

Figure 16 .
Figure 16.Effect of mode II fracture energy on global stress-strain relationship under tension (a) and compression (b).

Figure 16 .
Figure 16.Effect of mode II fracture energy on global stress-strain relationship under tension (a) and compression (b).

Figure 17 .
Figure 17.Effect of material parameter on global stress-strain relationship under (a) tension and (b) compression.

Figure 17 .
Figure 17.Effect of material parameter on global stress-strain relationship under (a) tension and (b) compression.

Figure 18 .Figure 19 .Figure 20 .
Figure 18.Failure modes of a specimen under the different frictional angles.

Figure 20 .
Figure 20.Effect of frictional angle on crack propagation under compression.

Figure 21 .
Figure 21.Effect of aggregates' mechanical characteristics on global stress-strain relationship under compression.

Figure 23 .
Figure 23.Failure modes of a specimen with different aggregate types; in the zoom areas of Figure 23b, the cracks through aggregates are marked by red, the others are marked by black.

Figure 22 .Figure 22 .
Figure 22.The relationship between strain and the sum of length of AII elements whose SDEG (overall value of the scalar damage variable) values are higher than 0.95.

Figure 23 .
Figure 23.Failure modes of a specimen with different aggregate types; in the zoom areas of Figure 23b, the cracks through aggregates are marked by red, the others are marked by black.

Figure 23 .
Figure 23.Failure modes of a specimen with different aggregate types; in the zoom areas of Figure 23b, the cracks through aggregates are marked by red, the others are marked by black.
where t 0 n , t 0 s and t 0 t represent the maximum allowable nominal stresses in normal and two shear directions, respectively.Likewise, t n , t s and t t represent the values of the nominal tensile stress and nominal shear stresses in two directions, respectively.

Table 1 .
Aggregate gradation and particle number.

Table 1 .
Aggregate gradation and particle number.

Table 2 .
The bulk and cohesive elements number of the model with different mesh precision.

Table 2 .
The bulk and cohesive elements number of the model with different mesh precision.