Computational technology for analysis of 3D meso-structure effects on damage and failure of concrete

Methodology for analysis of meso-structure effects on longer-scale mechanical response of concrete is developed. Eﬃcient algorithms for particle generation and packing are proposed to represent 3D meso-structures as collections of discrete features distributed randomly in a continuous phase. Specialised to concrete, the continuous phase represents mortar, while the features are aggregates and voids. Intra-and inter-phase co-hesive zones are used for failure initiation and crack propagation. A Monte Carlo approach is proposed to analyse the effects of meso-structure geometrical (volume density, size distribution and shape of features) and physical (strength and energy of cohesive zones) properties, whereas a number of model realisations with identical properties are used for statistical analysis. The results present the relative signiﬁcance of each meso-structure parameter for the emergent load capacity (tensile strength), damage evolution via micro-crack coalescence and macro-crack patterns, and failure energy density (toughness) of concrete. The proposed methodology is an effective tool for meso-structure optimisation in the design of concrete structures with prescribed requirements for strength and toughness.


Introduction
In structural design, concrete is considered as a homogeneous continuum.This first approximation is acceptable as long as components remain in the elastic regime of deformation, i.e. no energy dissipation mechanisms such as plasticity and surface separations (micro-cracking) are activated.When energy dissipation is present, it is localised, guided and governed by sub-continuum structures, the knowledge of which is critical for understanding and explaining inelastic and failure behaviour.
Concrete is a complex composite with sub-continuum structures at multiple length scales.The largest length-scale with observable heterogeneity is traditionally called the meso-scale.Meso-structure of concrete contains voids, aggregates, mortar, and interfaces between them.As closest to the component scale, the meso-structure plays a critical role in the observable macroscopic behaviour of concrete (Wittmann, 1983).In the past few decades, a number of meso-scale models considering both meso-structure and local failure mechanisms have been proposed and developed to understand and quantify their effects on longer-scale response of concrete.
With respect to the meso-structure representation of concrete, the two main approaches are the image-based modelling and the parameterization modelling.Although the meso-structure obtained from images, e.g. using computed tomography, is closest to nature, it is presently very expensive and time-consuming to generate sufficient scanned images from samples to make meaningful statistical analyses (e.g.Jivkov et al., 2013;Roubin et al., 2014;Sharma et al., 2013), especially for three-dimensional (3D) problems.For parameterization modelling, both direct (e.g.Du et al., 2014;López et al., 2007;Wang and Jivkov, 2015) and indirect approaches (e.g.Jivkov and Xiong, 2014;Leite et al., 2004;Reichert, 2009;Xiong et al., 2014) have been widely used to characterize heterogeneous materials.An advantage of the direct approach is that it allows for varying independent key meso-structure parameters, such as volume fractions, shape, size and spatial distribution of pores/voids and aggregates.Because of this feature, the direct approach is particularly suitable for statistical analysis of damage and failure of concrete performed in this study.
For meso-structure models of concrete, identification and generation of unit cell geometry is a vital step, since both shape and size of aggregates have significant influences on stress distribution, crack initiation and damage accumulation up to macroscopic failure within concrete (Wriggers and Moftah, 2006).Regarding two-dimensional cases, particles with regular shapes such as circle, ellipse and polygon are generally used.Beddow and Meloy (1980) proposed a morphological law to acquire rounded aggregates.Angular aggregates were generated respectively as polygons with prescribed elongation ratios (Wang et al., 1999) and convex polygons using Voronoi tessellation method (Galindo-Torres and Pedroso, 2010).More recently, Wang et al. (2015b) developed a "generate-and-place" procedure to simulate circular, elliptical and polygonal aggregates for meso-scale modelling of fracture and damage of concrete in 2D.
Regarding three-dimensional cases, aggregates are usually assumed to be spherical for simplicity, as presented in Bazant et al. (1990), Man and Van Mier (2008) and Wriggers and Moftah (2006).In recent years, considerable attention has been paid to generate various aggregate shapes, e.g., ellipsoidal aggregates using functions with varying parameters (Liu et al., 2014;Qian, 2012), and polyhedral aggregates by random packing systems (He et al., 2012) and Voronoi tessellation methods (Caballero et al., 2006;Galindo-Torres et al., 2012).It is worth noting that most of the existing meso-structure models of concrete only consider random aggregates (Bailakanavar et al., 2012;López et al., 2007) but neglect voids.However, X-ray computed tomography (XCT) images of concrete (Daudeville, 2013;Hunter, 2004;Man, 2010) clearly show that voids exist in concrete at mesoscale.It is therefore imperative to develop a procedure for automatic generation of morphological details of materials with both randomly distributed aggregates and voids of different shapes.
With respect to the numerical models for local material failure (new surface generation, micro-cracking), a number of models have been used to study damage and failure of concrete.The dominant approach at present employs continuum finite element modelling of the constituents with allowance for failure via cohesive interfaces between the continuum solid elements (Unger et al., 2007;Vořechovský and Sadílek, 2009).The key element is the behaviour of the cohesive interfaces, based on cohesive zone model developed by Barenblatt (1959Barenblatt ( ,1962) ) and Dugdale (1960) and later extended by Hillerborg et al. (1976).For cohesive zone model, the cohesive interface elements (CIEs) are pre-inserted or dynamically inserted into the initial finite element mesh so that realistic crack patterns can be simulated (Espinosa and Zavattieri, 2003).Because of its simple formation, easy implementation in the form of CIEs, cohesive zone model becomes more and more popular in modelling crack propagation in concrete and other quasi-brittle materials.It should be noted that other numerical methods, such as continuum strong discontinuity approach (CSDA) (Oliver et al., 2008;Simo et al., 1993), extended finite element method (XFEM) (Patzák and Jirásek, 2003), embedded finite element method (EFEM) (Roubin et al., 2015), etc. have been used for modelling concrete structures and show growing popularity in dealing with strong discontinuities.Other alternative approaches, such as discrete element models (Shiu et al., 2008) and lattice models (Jivkov, 2014;Reichert, 2009), have also been developed.A key issue associated with the class of discrete models is the difficulty in determining the model parameters which provide observed macroscopic behaviour, particularly for non-regular arrangements (Tu and Lu, 2011).This is avoided in the continuum-based modelling on the expense of the need to calibrate cohesive laws.Therefore, cohesive zone model is used in this work to represent failure initiation and propagation in concrete.
The aim of this work is to develop a holistic procedure for analysis of meso-structure and local properties effects on the macroscopic behaviour of concrete.The procedure contains three parts: generation of realistic meso-structures, which extends the work of Wang et al. (2015b) from 2D to 3D; incorporation of cohesive interfaces and description of their properties; Monte Carlo simulations (MCS) with spatially randomised meso-structures and statistical analysis of emergent behaviour of concrete.
Firstly, an in-house program HMG (heterogeneous material generator) consisting of particle generation and packing is developed and applied to generate 3D meso-structures of concrete with randomly distributed spherical, ellipsoidal and polyhedral aggregates and voids, respectively.Secondly, interfaces between all aggregate-Table 1 Particle size distribution of coarse aggretates in concrete (Hirsch, 1962).The proposed methodology allows for informed design of concrete structures.

Meso-structures of concrete
In this section, details of the developed HMG and the proposed intersection and overlap checking algorithms for random distribution of spherical, ellipsoidal and polyhedral aggregates and voids are presented.The proposed algorithms are suitable for generating a variety of heterogeneous materials, such as cement composites, nuclear graphite, graphene nano-composites, etc. Plain concrete with random aggregates and voids is set as an example in this section.

Particle size distribution -aggregate and void
The aggregate size distribution of concrete is often characterised by the Fuller curve that represents the gradation of aggregate particles resulting in optimum density and strength of concrete mixture.Fuller curve can be described by (Wriggers and Moftah, 2006) where P(d) is the cumulative percentage passing a sieve with aperture diameter d, d max is the maximum size of aggregate particles, and n is the exponent of the equation (n = 0.45-0.70).
For the ease of numerical implementation, the gradation curve expressed in Eq. ( 1) can be replaced with a number of segments, where the amount of aggregate (V agg ) within each grading segment [d i , d i + 1 ] can be expressed as (Wang et al., 1999) V where d max and d min denote the maximum and minimum sizes of aggregates, respectively, P agg is the volume fraction of aggregate and V is the sample volume of concrete.
A typical particle size distribution of coarse aggregates in concrete (Hirsch, 1962) as given in Table 1 is used in this study.According to their size, aggregates in concrete are generally classified into two categories: fine aggregate (e.g., sand) and coarse aggregate (e.g., gravel and crushed stone).Only coarse aggregates larger than 2.36 mm are considered in this study, while the large number of fine aggregates together with cement matrix is regarded as mortar (You and Buttlar, 2004), which keeps the model fidelity and avoids the computational difficulty caused by enormous number of elements.For normal strength concrete cast with a mould in the laboratory, the coarse aggregates typically occupy 30-50% of the concrete volume.
Mechanical behaviour of concrete is highly related to aggregate, in particular, the shape of the aggregate particles that depends on aggregate type.It was found that gravel aggregates have spherical or ellipsoidal shape, while crushed aggregates have polyhedral shape (Daudeville, 2013).Pores in concrete can be roughly divided into air voids, capillary pores and gel pores.The capillary and gel pores (nanometre or micrometre scales) are too small and outside the scope of mesoscale modelling.Larger pores of dimensions of up to a few millimetres are the result of the air entrapped during mixing and not removed by vibration of fresh concrete (Bertolini et al., 2013), which could be clearly observed from CT images at meso-scale (Hunter, 2004;Man, 2010).Therefore, it is necessary to include voids in the meso-structure of concrete.In this study, voids with size ranging from 2 to 4 mm as presented in Hordijk (1991) are taken into account.The shape of voids is assumed to be sphere or ellipsoid.The algorithms and procedures to generate individual particle with various shape of sphere, ellipsoid and polyhedron are given in Appendix.Wang et al. (1999) presented a comprehensive procedure using a commonly adopted "taking" and "placing" method to generate a 2D random geometric arrangement of aggregates.A similar procedure was adopted in our previous published work for generating 2D meso-structures of concrete with different shape of aggregates and voids (Wang et al., 2015b).The present study extends it to randomly distribute 3D particles of aggregate and void with various shape into a cube according to the prescribed particle size distribution and density.The flowchart of particle packing composed of "input", "taking" and "placing" processes is shown in detail in Fig. 1.

Particle packing -aggregate and void
The central idea is to create aggregates and voids in the cubic concrete in a repeated manner, until the target volume of aggregates and voids is achieved.The "input" process reads the controlling parameters for generating a meso-structure with random aggregates and voids.The "taking" process generates an individual aggregate or void in accordance with the random size and shape descriptions.The "placing" process subsequently places the aggregates and voids into the predefined area in a random manner, subjected to the prescribed physical constraints.The generation of synthetic 3D meso-structures is implemented in MATLAB.

Intersection and overlap checking algorithms
The last step of particle packing is to check whether the distributed aggregates and voids intersect or overlap with each other.The straightforward, efficient but easily implemented algorithms for intersection and overlap checking of 3D spheres, ellipsoids, polyhedrons and their different combinations are proposed in this section.In order to place an aggregate particle at a free position and orientation within the concrete volume, three conditions need to be satisfied: (1) the whole particle should be completely within the boundaries of the concrete volume; (2) there must not be any overlap with previously placed particles; (3) the minimum distance between the edge of a particle and the boundary of the concrete specimen, and the minimum gap width between two adjacent particles should be limited.
The first condition can be easily met by detecting if the minimum and maximum coordinates of the particles are inside the concrete section.The third condition can be checked concurrently with the checking of first two conditions by enlarging the size of the particle to (1+γ ) its size.The same criterion applies to the distribution of voids in concrete.Given the set of independent and dependent parameters, the intersection and overlap checking algorithms in different situations are given below.

Random spherical particles
For two spherical particles, the intersection and overlap condition can be checked easily by comparing the distance between the particle centres and the sum of two radii.
where x 0 , y 0 , z 0 are the coordinates of centre of the newly generated sphere; r is the radius of the new sphere.

Random ellipsoidal particles
It is neither economical nor necessary to check each newly generated particle against all the existing particles for intersection and overlap checking.In order to reduce the computational cost, intersection and overlap checking is performed by using a less costly hierarchy method.The computational cost is decreased by only considering particles in the near field of the particles being added and excluding all the other particles in the far field.
In Hierarchy I, for two ellipsoids with centres A and B, if the distance between A and B is larger than 2a (a is the length of major semi principal axis), then these two ellipsoids do not intersect or overlap.Thus, all the particles that are at a radial distance greater than 2a from the centre of the particles can be excluded for next intersection and overlap checking.Consequently, no particles or very few particles need to be further checked which results in a significant increase in the computational efficiency.
Hierarchy II could further reduce the computational cost by checking whether the centres of any previously generated ellipsoids filtered in Hierarchy I are located inside or on the surface of new ellipsoid before it goes to Hierarchy III.This could make sure that the new ellipsoid does not totally contain one of existing particles which will be ignored in Hierarchy III.Finally, Hierarchy III checks whether any nodes on new ellipsoids are located inside or on the surface of ellipsoids filtered in Hierarchy II.
The nodes on the surface of an ellipsoid may be parameterized in several ways.One possible choice is spherical coordinate system: where θ is the polar angle measured from a fixed zenith direction, ϕ is the azimuth angle of its orthogonal projection on a reference plane that passes through the origin and is orthogonal to the zenith, measured from a fixed reference direction on that plane.
Considering an arbitrary node N( x, ȳ, z) on an ellipsoid surface, its location relationship to another ellipsoid can be determined by calculating F( x, ȳ, z): where is the ellipsoid, is the boundary of the ellipsoid.If F(x, y, z) > 0, then node N is outside of the ellipsoid; if F( x, ȳ, z) = 0, then node N is on the surface of the ellipsoid; if F( x, ȳ, z) < 0, then node N is inside of the ellipsoid.So if any node on the surface of the newly generated ellipsoid is inside of any previously generated ellipsoids, the particle will not be placed into the model.In such a case, another set of random numbers (a, b, c, α, β, γ , x 0 , y 0 , z 0 ) is generated for an ellipsoid and an attempt of placing the particle at a new location and orientation is made.Two ellipsoids are regarded as separated only when all the nodes on one ellipsoid are outside of the other one.It could be achieved by calculating F 1 (x 0 , y 0 , z 0 ) in Hierarchy II and F 2 ( x, ȳ, z) in Hierarchy III.If F 1 (x 0 , y 0 , z 0 ) and F 2 ( x, ȳ, z) are both positive, then the two ellipsoids are considered to be separated.The coordinates of all the nodes on the ellipsoid surface could be obtained by search algorithms in MATLAB.

Random polyhedral particles
This section discusses the method of intersection and overlap checking for convex polyhedrons.The generated random convex polyhedrons are made of the random nodes extracted from spheres.Likewise, the hierarchy algorithm as described above is used here to reduce the computational cost.So all the particles that are at a radial distance of the original spheres greater than (r 1 + r 2 ) from the centre of the particle being added can be excluded for next intersection and overlap checking.
The method of intersection and overlap checking for polyhedrons is based on determining whether two convex polyhedrons intersect with each other.This is achieved by detecting whether all the vertexes of the existing polyhedrons lie on one side of an arbitrary plane of the new polyhedron, while any point inside the new polyhedron lie on the another side.A 3D example with polyhedron L on the left and polyhedron R on the right is shown in Fig. 2 to illustrate the idea of the algorithm.If all the vertexes 1-6 of polyhedron L lie on the left side of the plane ABC consists of vertexes A, B and C, while an The equation of plane ABC can be obtained with nodal coordinates denoted by A(x 1 , y 1 , z 1 ), B(x 2 , y 2 , z 2 ) and C(x 3 , y 3 , z 3 ) as follows So the general equation of plane ABC can be expressed as: Considering all the vertexes N(x v , y v , z v ) on polyhedron L and the point O(x o , y o , z o ) inside polyhedron R, the location of point N relative to the polyhedron R can be determined by calculating for all the vertexes, then these two polyhedrons are considered to be separate.

Random spherical/ellipsoidal particles and random spherical/ellipsoidal voids
The first hierarchy algorithm that a radial distance of the sphere and ellipsoid greater than (r + a) from the centre of the particles being added can be excluded for the next intersection and overlap checking is used.The general equation of sphere can be expressed by the standard equation (see Eq. (A1) in Appendix): The surface of the sphere is parameterized using spherical coordinates (r, θ , ϕ): where θ and ϕ have the same meaning as those given in Eq. ( 4), r 0 and r 1 are the minimum and maximum radii of particle , η denote the random number uniformly distributed between 0 and 1.
Likewise those for ellipsoids, the intersection and overlap between ellipsoids and spheres can be evaluated by checking whether all the nodes on an ellipsoid/sphere are outside of the other sphere/ellipsoid.Considering all nodes N(x n , y n , z n ) on an sphere/ellipsoid, their location relative to another ellipsoid/sphere can be determined with F(x n , y n , z n ) or H(x n , y n , z n ).If F(x n , y n , z n ) > 0 or H(x n , y n , z n ) > 0 for all the vertexes, then the sphere and ellipsoid are regarded as separate ones.Meanwhile, Hierarchy II is performed to make sure that the new sphere/ellipsoid does not contain any previous ones, the centres of which may be inside or on the surface of the new particle.

Random polyhedral particles and random spherical/ellipsoidal voids
The first hierarchy algorithm similar to the previous one but with different region range, (r + r') for polyhedron and sphere, and (r + a) for polyhedron and ellipsoid, is used.The further checking is carried out by using the similar method as that used for polyhedrons as described above, but all nodes on the spheres/ellipsoids surfaces are regarded as vertexes.Here the spheres or ellipsoids are assumed to be generated before the polyhedrons.Considering all nodes N(x n , y n , z n ) on a sphere/ellipsoid and the point O(x o , y o , z o ) inside of the polyhedron, the sphere/ellipsoid location relative to the polyhedron can be determined with for all nodes, then the polyhedron and sphere/ellipsoid are regarded as separate ones.

Concrete with gravel aggregates
Based on the developed particle generation and packing procedure, a series of cubic concrete samples with dimensions of 50 mm × 50 mm × 50 mm containing spherical and ellipsoidal gravel aggregates and voids are generated, as shown in Fig. 3, where aggregates are in green and voids are in red.The other phase in the cube represents mortar.The used aggregate size distribution is given in Table captions Table 1.The aspect ratios for ellipsoidal aggregates and voids are The minimum thickness of the mortar film around the aggregates is a challenge, because no simple theory is available yet for its evaluation.The thickness was assumed to be 0.1-1 mm in most of the previous numerical studies (Kim and Al-Rub, 2011;Pedersen et al., 2013).Although it has been found to have an insignificant effect on the concrete behaviour (Kim and Al-Rub, 2011), it plays an important role in the spatial distribution of the aggregate particles: a larger value would result in a more uniform spatial distribution but could lead to the limit for reaching a high volume fraction.In present study, the minimum space between the edge of an aggregate and the boundary of the concrete specimen (γ 1 ), and minimum gap width between any two aggregates (γ 2 ) are set to be 0.2 mm and 0.2 mm, respectively.
Concrete specimens composed of spherical gravel aggregates with volume fraction (P agg ) of 30% and spherical/ellipsoidal voids with volume fraction (P void ) of 0% and 2% are shown in Fig. 3(a), (b) and (c), respectively.Concrete specimens made up of ellipsoidal gravel aggregates with volume fraction of 30% and spherical/ellipsoidal voids with volume fraction of 0% and 2% are shown in Fig. 3(d), (e) and (f), respectively.

Concrete with crushed aggregates
Polyhedral aggregates and spherical or ellipsoidal voids are considered for concrete with crushed aggregates.A series of 3D cubic concrete samples (50 mm × 50 mm × 50 mm) consisting of mortar, polyhedral crushed aggregates and spherical/ellipsoidal voids are generated and shown in Fig. 7.The minimum and maximum edges of  the polyhedrons are set to N min = 8, N max = 16.All other parameters are the same as those used for concrete with gravel aggregates above.

Efficiency of the proposed algorithm
In order to investigate the quality of the proposed hierarchy algorithm, the computational times involved in generating typical unit cells (50 mm × 50 mm × 50 mm) with ellipsoidal particles (P agg = 20%) and spherical voids (P void = 1%) by employing the present proposed algorithm and generally used traditional algorithm are shown in Fig. 8.These simulations were performed using an Intel Core i7 @ 3.4 GHz desktop computer with 8 GB RAM.It is evident that the present algorithm significantly reduces computational cost as compared to the traditional algorithm, with mean computational (a)  time 16.2 s for the present algorithm and 515.8 s for the traditional algorithm.

Mesh generation
The generated meso-structures of concrete as illustrated in Section 2 are subsequently meshed for modelling with cohesive ele-ments.In this work, all finite element meshing are performed with the pre-processing functionality in commercial FE package ANSYS.The free meshing technique with mesh adaptation functions was used to avoid the restriction on model geometry.For meshing, different components in concrete, e.g., mortar and aggregates should maintain continuity on their surfaces.Thus, the finite element boundaries are coincident with material surfaces and there are no material discontinuities within the elements.The detailed meshing procedure is described as follows: Firstly, the entire concrete volume V all and the volume representing aggregates V agg are created.The "overlapping" Boolean operation in the ANSYS pre-processor is applied to separate aggregates from entire domain.The remaining volume belongs to mortar and voids, denoted by V mix .Secondly, the volume representing voids V void is created, and the "overlapping" Boolean operation is utilized again to separate V void from V mix .The remaining volume is mortar, denoted as V mortar .Finally, the volume representing voids is deleted, since voids do not contribute to macroscopic mechanical properties of concrete.
The mortar and aggregates are meshed with tetrahedral elements so that more realistic crack paths can be obtained.The ANSYS APDL (parametric design language) programs in combination with ANSYS batch processing are developed to provide a powerful tool of automatic mesh generation for a large number of samples required for Monte Carlo simulation.

Cohesive element insertion
After finite element meshing, the generated mesh will automatically have shared nodes at the interfaces between two elements.A duplicate set of nodes is required at the interfaces in order to simulate the potential micro-crack initiation in cohesive zone model.Here 6-node zero in-plane thickness cohesive interface elements (CIEs) are pre-inserted into the existing element surfaces.This is implemented by using an in-house computer program CIEIN, which is programmed in FORTRAN.
An example of cohesive interface elements inserting into the surface between aggregate and mortar is shown in Fig. 9, where three sets of CIEs with different traction-separation softening laws, namely, CIE_AGG, CIE_MOR and CIE_INT are assigned to aggregates, mortar and interfaces, respectively.As a result, Fig. 9(a) and (b) show the initial FE mesh (4 elements and 7 nodes) is changed into FE mesh with inserted cohesive elements (7 elements, 15 nodes).
The entire procedure for interface elements inserting into finite element mesh is described in detail below: (1) Reading the original mesh file generated from ANSYS including information on element (nodal connectivity), node (nodal coordinate) and material numbering (1 for aggregate and 2 for mortar); (2) Every normal node between each pair of solid elements (i.e., aggregate and mortar) is duplicated so that a separated node is created in the same position for every two adjacent elements.
The nodal connectivity is changed by using these duplicated nodes, which results in each element being not connected with its neighbouring elements; (3) Inserting three different zero-thickness cohesive interface elements (CIEs) with distinct mechanical properties, namely CIE_AGG, CIE_MOR, and CIE_INT into aggregate-aggregate, mortar-mortar, mortar-aggregate interfaces.As a result, all solid elements are connected with each other again; (4) Generating a new input file in INP format for ABAQUS by updating the mesh information of solid elements C3D4_AGG for aggregate, C3D4_MOR for mortar, and three cohesive interface elements CIE_AGG, CIE_MOR, and CIE_INT for interfaces.
A loop algorithm has been developed to generate every INP file including information about material properties, boundary and loading conditions for ABAQUS in combination with the mesh files automatically obtained from ANSYS.
A typical finite element model including cohesive interface elements for cubic concrete specimen with 30% ellipsoidal aggregates and 2% spherical voids is taken as an example and shown in Fig. 10.Different phases made of meso-structure of concrete, i.e., aggregate, void and mortar, are identified and meshed firstly, as shown in Fig. 10(a

Cohesive interface constitutive behaviour
The behaviour of zero-thickness cohesive elements is given with a relation between surface traction (stress) and separation (displacement between originally coinciding surfaces).In 3D cohesive zone model, it is generally assumed that there exist a normal traction t n and two tangential traction (shear cohesion) t s and t t across the crack surfaces, through mechanisms such as material bonding, aggregate interlocking and surface friction in fracture process zone.The traction-separation relation has two distinct phases: (1) surfaces maintain full integrity; and (2) surfaces gradually separate to complete loss of interaction.For quasi-brittle materials, such as the concrete constituents, mortar and aggregate, phase (1) can be assumed linear elastic.The critical point between the two phases, called the damage initiation point, can be defined in terms of either critical traction or critical separation.In this study, critical points are defined in terms of critical tractions for the three modes, t n0 , t s0 and t t0 , respectively, and the element damage initiates when the following condition is met For the separation phase (2), we have selected a linear-softening response, i.e. the traction decreases linearly with the separation beyond the critical.This leads to bilinear cohesive crack model, illustrated in Fig. 11 for normal and tangential separations.In the figure, δ n0 and δ s0 (δ t0 ) denote the critical separations, corresponding to the critical tractions t n0 , t s0 and t t0 , respectively, and δ nf and δ sf (δ tf ) denote the separations at complete failure for normal traction and tangential traction, respectively.The unloading paths are also indicated.The initial stiffness k n0 , k s0 and k t0 must be set high enough to minimize an overly-flexible response in cohesive zone system, but not too high to produce instabilities in the stiffness matrix.This is because the compliance due to the presence of interfaces should be negligible compared with the compliance of continuum elements; they represent the potential crack paths in a physical sense.Meanwhile, it should be noted that excessively high values may lead to illconditioning of the system equations and low values cannot ensure displacement continuity across the interfaces in the elastic range.The effects of initial stiffness on computational results have been investigated previously (Qiang et al., 2000;Yang et al., 2015).The following relation is suggested in Qiang et al. (2000) as a guideline for initial stiffness selection: where E and v are Young's modulus and Poisson's ratio, b is the characteristic size of elements, and c is taken as 10-100 from the experience in Qiang et al. (2000).An appropriate elastic stiffness of the interfaces equal to 10 6 MPa/mm is chosen as used in 2D models (López et al., 2007;Wang et al., 2015b) and 3D studies (Caballero et al., 2006).
The areas under the curves respectively stand for the normal fracture energy and tangential fracture energy G sf (G tf ).These are related to the tensile/shear strength and failure separation according to Hence knowledge of two of the parameters, G, t 0 and δ f , fully determines the traction separation law.In this work we have used knowledge of fracture energies and critical tractions from previous works.
The evolution of damage in phase (2) under the combined normal and tangential separations is described via a scalar index D. To this end an effective relative displacements δ m , introduced as where < > is the Macaulay bracket is used to define the scalar damage where δ m,max is the maximum effective relative displacement obtained during the loading history, δ m0 and δ mf denote the effective relative displacements at damage initiation and final failure, respectively.Notably, the damage variable D evolves monotonically from 0 to 1 during loading.Damage evolution affects the stiffness coefficients k n , k s and k t for unloading and reloading, which change according to This affects the tractions, which change according to where tn and ts are the traction components predicted by the elastic traction-displacement behaviour for the current separation without (a) t n n curve in normal direction (b) t s (t t ) s ( t ) curve in tangential direction Loading Unloading damage.The theoretical framework described above is used for the intra-phase and inter-phase cohesive elements in the meso-structure with different parameters for mortar-mortar, mortar-aggregate and aggregate-aggregate boundaries.

Loading Unloading
The constitutive laws selected for CIEs in this study are primarily for modelling interface fracture, and do not represent fully the coupling between shear and normal responses, particularly the interaction between shear and compression.It is expected that this simplification has minor effect on the analysis of concrete under uniaxial tension in the pre-peak regime dominated by tensile micro-cracking.The effect on the post-peak behaviour, governed by micro-crack coalescence and associated development of shear failures, might be captured less accurately.In order to accommodate general loading conditions, specific cohesive models are required to incorporate the mechanism that the shear strength is highly dependent on the normal stress at the interfaces.Such models, based on finite-thickness cohesive elements are under development.

Models description for Monte Carlo analysis
All models analysed in this work have dimensions 25 mm × 25 mm × 25 mm.For a given set of meso-structure characteristics (aggregate and void content, shape and size distribution), a number of model realisations differing in the spatial distributions of these features, have been generated for the purpose of subsequent statistical analysis of the macroscopic responses.Simulations with the realisations corresponding to a fixed set of meso-structure parameters constitute a Monte Carlo analysis of the set.
The mechanical properties of the two phases in concrete as well as the properties of the interface elements used in this work, obtained from Caballero et al. (2006), are given in Table 2. Regarding the aggregate-aggregate interfaces (CIE_AGG), it is noted that aggregates have much higher strength than mortar-mortar (CIE_MOR) and mortar-aggregate interfaces (CIE_INT) in normal concrete.Therefore, cracks are normally not allowed to initiate inside aggregates unless otherwise specified.This is implemented by defining elastic behaviour without damage for CIE_AGG.Due to lack of experimental data, the shear fracture properties were assumed to be the same as the normal ones (Wang et al., 2015a(Wang et al., ,2015b) ) for both CIE_MOR and CIE_INT.The linear tension/shear softening laws described in Section 3.3 are used to model these interfaces.
The loading on all model specimens has been applied via prescribed displacement, where horizontal displacements were prescribed to all nodes on the left and right surfaces, with a value equal to zero for the left surface, and a uniformly distributed displacement on the right surface.The analyses were terminated at a displacement d = 0.1 mm, corresponding to macroscopic strain of ε = 0.002, the same as those used in our previous 2D study (Wang et al., 2015b).Solutions were performed with the commercial software ABAQUS/Explicit (Dessault Systems, 2012), it determines a solution to the dynamic equilibrium equation without iterating by explicitly advancing the kinematic state from the previous increment, which allows for solving highly nonlinear equation systems with softeninginduced instabilities efficiently.
The Abaqus/Implicit solver with automatic incrementing often fails to perform simulations for highly nonlinear problems with material softening.A very small value of the maximum increment is required to model the whole loading process, while it will lead to a high computational cost.The Dynamic/Explicit solver, with sufficiently long loading time to minimize the dynamic effect, is able to simulate the full response accurately with a reasonable time and thus is used in this study.

Analysis of simulation parameters
Cohesive zone models are known to produce results affected by element size and, in the case of dynamic analysis, such as in Abaqus/Explicit, by loading time.This sub-section presents an advance study of these effects to inform on most appropriate size and time.Further the effect of Monte Carlo sample number, i.e. number of realisations per fixed set of meso-structure parameters, is analysed.The outcomes are summarised as those used for the parametric studies in the following sub-sections.

Effect of mesh density/element size
When the cohesive crack paths are unknown a priori, fine meshes are needed to minimise the dependence of simulation results, e.g., stress-displacement curves and crack patterns, on initial FE mesh.However, this will lead to large-scale nonlinear equation systems, and the computational costs often become prohibitively expensive,  especially for Monte Carlo simulations.Moreover, the effect of the minimum thickness around aggregates may potentially be related to the global mesh size.Therefore, an appropriate mesh density is required to achieve a balance between the accuracy and computational efficiency.Herein, three meshes with different average element length (i.e., L e = 1 mm, L e = 0.8 mm and L e = 0.6 mm) are used, as shown in Fig. 12. Mesh I has 156,241 nodes, 28,432 solid elements and 307,475 cohesive elements.The numbers for the Meshes II and III are 236,260,42,714 and 465,435,and 349,777,62,466 and 690,464, respectively.
The simulated stress-displacement and dissipated fracture energy-displacement curves for each element length are shown in Fig. 13, where the mean stress is calculated as the sum of reaction force of nodes on the left boundary divided by the cross-section area of specimen.It can be observed that the mesh dependence of the stress-displacement and dissipated fracture energy-displacement curves is negligible for the selected element sizes.
The obtained final crack paths for different meshes are shown in Figs.14-16, where plots (a) represent the failure surfaces and plots (b) represent the damage via the separations of the cohesive elements; cohesive elements with damage index D ≥ 0.95, where D = 1 means complete failure, are shown in red.Despite a few differences, the three figures show similarity in terms of crack pattern and failed surface morphology.Considering the excellent agreement between the stress-displacement responses in Fig. 14, and the balance between accuracy and efficiency, the mesh size of 1 mm (Mesh I) has been selected for all meshes in the following simulations.

Effect of loading time
The loading time has a significant effect on the computational efficiency and accuracy of results when the ABAQUS/Explicit solver is used to model quasi-static loading conditions.In principle, the loading time must be long enough to minimise any dynamic effects.However, a long loading time results in greater computational cost.So a balance also has to be made between the computational efficiency and accuracy.Fig. 17 shows a comparison of σ -d curves obtained from 50 random samples with three loading times of 0.001 s, 0.005 s and 0.01 s.It can be seen that a loading rate influence occurs when loading time is 0.001 s (a), and there is an obvious oscillation.For loading rates

Typical tensile behaviour
3D modelling allows for more realistic representation of aggregate and void distribution in concrete compared to previous 2D models.The need for statistical representativeness requires a significant number of model realizations per set of meso-structure parameters to en-sure convergence.This makes the Monte Carlo approach substantially more computationally demanding in 3D and requires some balance between number of realizations and available computer power.
Stress-displacement curves from a Monte Carlo simulation with 50 realisations are shown in Fig. 18 together with the numerical mean curve and experimental data acquired from tension tests of six samples by Hordijk Hordijk (1991).The simulated stress-displacement curves in this work fit well with the experimental scatter.They are similar in terms of clear peak and sharp initial post-peak drop followed by a long shallow tail.The peak stress has a mean of 3.49 MPa   with standard deviation of 0.05 MPa, demonstrating that it is relatively insensitive to the spatial randomness of aggregates and voids.However, notably large differences are observed between the softening responses of the different realisations.This suggests that sufficient number of realisations is necessary to capture a required accuracy.In order to select suitable points for comparison between different numbers of realisations, the standard deviation of stress during loading is plotted in Fig. 19.The standard deviation of stress is relatively low (around 0.05 MPa) before it reaches the peak stress (displacement = 0.01 mm), and the deviation first increases to around 0.5 MPa and then decreases.

Effect of sample number
In order to evaluate the effect of Monte Carlo sample number on simulation results, two special loading points, corresponding to the peak stress and the displacement of 0.03 mm when the highest standard deviation is observed during loading process (Fig. 19), are investigated, the results are shown in Figs.20 and 21, respectively.Both the mean value and standard deviation of stress tend to be stable when the sample number is about 50, which indicates that 50 random samples are enough to achieve statistical convergence.It should also be noted that Gaussian distribution was found in Wang and Jivkov Wang and Jivkov (2015) for peak stress which indicates the stress distribution is not completely random, however, substantially more samples are necessary to yield conclusive results for softening behaviour due to the relatively high coefficient of variance of stress (standard deviation/mean value).Monte Carlo simulations with 50 realisations per set of parameters were carried out in this study using the High Performance Computing (HPC) cluster at Computational Shared Facility (CSF), University of Manchester.A typical Monte Carlo analysis of 50 realisations using the mesh size of 1 mm and loading time of 0.005 s took approximately 6 hours by ABAQUS default parallel processing with 32 CPUs.

Comparison with 2D models
As mentioned in the introduction, most of existing work on mesoscale modelling of damage and failure of concrete is performed with 2D models.It is therefore interesting to assess the difference between the responses of 2D and 3D meso-structures.To this end, 50 2D models of concrete with elliptical aggregates and circular voids were generated using the algorithms from Wang et al. (2015b).All meso-scale parameters were fixed to the values used in the 3D model described in Section 3.4.Afterwards, Monte Carlo simulations are carried out to estimate the mechanical behaviour of these 2D concrete samples under uniaxial tension.
For comparison, the mean responses from 2D and 3D simulations are given together in Fig. 22.The mean peak stresses predicted by 2D and 3D modelling are 2.65 MPa and 3.49 MPa, respectively, i.e. the more realistic 3D meso-structure predicts 24.1% higher peak stress.The standard deviations of peak stress are 0.11 MPa and 0.05 MPa, respectively, i.e. the 3D model provide 54.5% lower standard deviation.In the softening phase of the responses, the 3D models show generally larger standard deviation of stress than the 2D models.The increase of peak stress in 3D is attributed to the constraint effect in the thickness direction, not captured by the 2D model.The decrease of standard deviation of peak stress is attributed to the larger ratio between surfaces available for cracking and volume in 3D than line segments available for cracking and area in 2D.This leads to more uniformly distributed pre-peak micro-cracks in the 3D volume.The outcome is that the spatial distribution of features has substantially weaker effect on the pre-peak response in the 3D models than in the 2D models.In the softening branch, the larger standard deviation of stress predicted by the 3D models is related to the process of microcrack coalescences into macroscopic cracks.Here the more uniformly distributed micro-cracks from the pre-peak phase in 3D have significantly more options to coalesce depending on the spatial distribution of the features.As a result, the spatial distribution of features guiding micro-crack coalescence has substantially stronger effect on the post-peak response in the 3D models than in the 2D models.
The results presented suggest that the use of 2D approximation can bring significant underestimation of strength and toughness of concrete, where toughness is understood as the area under the stress-displacement curve.Clearly, the third dimension cannot be neglected due to the nature of fracture process.The importance of the thickness effect has been pointed out experimentally by Van Mier and Schlangen (1989).

Crack patterns
Two typical macro-crack patterns are observed from the Monte Carlo simulations under uniaxial tension: Type I cracking with only one dominant crack illustrated in Fig. 23; and Type II cracking with two or more dominant cracks illustrated in Fig. 24.This finding is consistent with that observed by Roubin (2014)  The damage index is represented in red colour.In the early stages of loading, a large number of micro-cracks initiate at mortaraggregate interfaces.As loading increases, some cracks coalescence into the dominant cracks, while the other cracks arrest due to stress redistribution.Type I cracking arises when micro-cracks coalesce into a single dominant crack.Type II cracking occurs when two or more independent cracks form during the coalescence process.
Fig. 26 shows the stress-displacement curves corresponding to two different types of cracking, red dashed lines for Type I cracking and blue solid lines for Type II cracking.Among 50 Monte Carlo samples investigated in this particular study, 38 samples behave as Type I cracking and 12 samples behave as Type II cracking.It can be seen that, there is an almost identical response of monotonous increase in the pre-peak range.For both Type I and Type II, the average peak stress is very close to each other, about 3.50 MPa.However, a clear difference is observed between the post-peak responses.
If a single crack develops (Fig. 23), the stress-displacement curve presents a rapid softening (see red curves in Fig. 26).If two or more macro-cracks form concurrently (Fig. 24), the softening is more graceful (see blue curves in Fig. 26).This indicates that the post-peak softening is related to cracking type.The emergence of multiple concurrently growing macro-cracks leads to larger apparent toughness of the representative volume analysed.

Effect of aggregate and void content
Based on Monte Carlo simulations, the effect of aggregate content is investigated.All the parameters except for aggregate volume fraction are fixed to the values used in Section 3.4.Fig. 27(a) shows the mean stress-displacement curves for samples with different aggregate volume fractions.The corresponding mean peak stress against aggregate volume fraction is plotted in Fig. 27(b).It can be observed that the mean load capacity decreases 14% in strength, from 3.99 MPa to 3.49 MPa, with the increase of aggregate volume fraction from 0% to 30%.This is because the increase of weak  aggregate-mortar interfaces has a greater influence than the increase in strong aggregates as the cohesive strength of interfacial CIEs is assumed to be only half of mortar CIEs in this study.Fig. 28 shows the effect of void content on the concrete response.The mean tensile strength decreases from 3.59 MPa to 3.39 MPa, as the void content increases from 0% to 4%.This is attributed to the fact that the voids provide easier pathways for crack propagation.Although the differences observed for void contents rising from 0% to 4% are not significant, the results indicate that the voids existing in concrete should not be neglected in modelling of mechanical properties and fracture of concrete, in particular for concrete with high void content, such as porous concrete.

Effect of aggregate and void shape
The effect of aggregate and void shape is investigated by undertaking 50 Monte Carlo simulations with varying aggregate and void shape while all other parameters are fixed to those used in Section 3.4.Fig. 29 illustrates two typical crack types observed in concrete with spherical/polyhedral aggregates and ellipsoidal voids.Both types of cracking can be found in concrete samples with different shape of aggregates and voids.This suggests that the type of cracking is independent of aggregate and void shape.
Fig. 30 shows the simulated mean stress-displacement curves for concrete samples with different aggregate and void shape under uniaxial tension.The results indicate that the load capacity of concrete containing spherical and ellipsoidal aggregates is greater than that including polyhedral aggregates.The difference of tensile strength between samples with spherical aggregates and elliptical aggregates is about 1% while its difference between samples with spherical aggregates and polyhedral aggregates is about 3%.This may be explained by two factors.On the one hand, concrete with unsmooth aggregates has more mortar-aggregate interface elements if aggregate and void content are constant.This is presented in Table 3.On the other hand, the local stresses are enhanced by the higher stress concentration at the corners of polyhedral aggregates, while the smooth edges of spherical and ellipsoidal aggregates have a more benign stress distribution which delays the fracture process and increases the tensile strength.

Analysis of cohesive interface effects
In this sub-section, the effect of principal material parameters, critical stress and cohesive energy, used for aggregate-aggregate, mortar-mortar and mortar-aggregate interface elements on the simulation results are investigated.All other geometric parameters are fixed to those used in Section 3.4.It can be seen from Fig. 31 that the load capacity is greatly affected by the critical tensile stress of mortar-aggregate interfaces.The decrease of tensile stress of aggregate-mortar interfaces leads to a rapid loss of load capacity.This is because the aggregate-mortar interfaces are weak relative to aggregates and mortar, and become preferred crack paths.different critical tensile stresses of mortar-aggregate interfaces.At the lowest stress of 1 MPa of mortar-aggregate interface, the crack path appears to be along mortar-aggregate interfaces.For higher mortar-aggregate interface critical stress, the crack paths tend to be independent of the critical stress of mortar-aggregate interfaces.
(a) t n =2 MPa (b) t n =4 MPa   Fig. 33 shows that the load capacity is hardly affected by the critical stress of the aggregate-aggregate interfaces until it becomes equal to the strength of mortar-aggregate interfaces.This is because aggregates are stronger than mortar and interfaces, and the load capacity is not sensitive to the tensile strength of aggregates.The crack paths for different tensile strength of aggregate-aggregate interfaces are shown in Fig. 34.It can be seen that, as the aggregate strength is reduced, cracks seem to a greater tendency to propagate across the aggregates leading to straighter paths with less bridging and branching.strength of mortar-mortar interfaces is the same as that of mortaraggregate interfaces (see Fig. 36(a)).This is because the mortaraggregate interfaces are no longer the weakest part.However, the crack paths seem to be insensitive to the tensile strength of mortarmortar interfaces when the mortar is strong.

Effect of cohesive energy
Figs. 37, 39, and 41 depict the simulation results obtained by varying the cohesive energy of interface elements.For each curve in the figures, the legend includes three values, which denote the cohesive energy assigned to the mortar-aggregate, aggregate-aggregate and mortar-mortar interfaces, respectively.As in the previous subsection, a missing value for the aggregate-aggregate interface indicates that cracking through aggregates is not allowed.
It can be seen from Fig. 37, that the cohesive energy of mortaraggregate interfaces has a significant effect on the softening response of the concrete.The results shown in Fig. 41 indicate that the cohesive energy of mortar-mortar interfaces has a pronounced effect on shape of softening branch and a slight effect on peak stress.The stress in the postpeak range diminishes at a lower rate when the cohesive energy increases, which leads to a larger critical displacement.This suggests that the increase of cohesive energy for mortar-mortar interfaces result in a softer response of concrete samples under tension.The effects of meso-scale material parameters, critical stress and cohesive energy of meso-structure constituents, are logical and explain the variability in tensile strength and failure energy density (toughness) which is observed experimentally using variable mortar and aggregates.

Conclusions
Effective and efficient algorithms for generating threedimensional meso-structures and corresponding finite element models with cohesive elements are developed.Meso-structures of concrete contain spherical, ellipsoidal or polyhedral aggregates, and spherical or ellipsoidal voids, dispersed in mortar, with intra-phase (aggregate/aggregate and mortar/mortar) and inter-phase (aggregate/mortar) cohesive elements.Monte Carlo approach is proposed for investigating the effects of 3D scaling, aggregate and void content, aggregate and void shape, and cohesive properties on longer-scale mechanical response.Core of the approach is statistical analysis of responses of 50 model realisations with given meso-structure parameters.3. Compared to 2D modelling, the 3D modelling demonstrates: a larger mean peak stress and a smaller standard deviation in the pre-peak response, attributed to more uniformly distributed micro-cracks within the volume; a larger standard deviation in the post-peak response, attributed to the larger number of possibilities for micro-crack coalescence under the constraint of randomly distributed features.4. The cohesive properties of concrete constituents affect strongly the macroscopic mechanical response and the cracking patterns: critical tensile stress of constituents has direct effect on the peak stress and differences between constituents' strengths change the patterns of micro-crack coalescence and macro-crack propagation; cohesive energy has little effect on the pre-peak behaviour, but affects strongly the softening behaviour and crack paths of concrete under tension.
It is acknowledged that all conclusions of this work are derived from meso-scale models of concrete subjected to uniaxial tension.The selection of simple zero-thickness cohesive elements with uncoupled tensile and shear responses has been dictated by this simple loading scenario.Development of advanced cohesive elements, specifically for coupling compressive and shear responses for simulations under more complex loading conditions, is a subject of ongoing research.

A3. Polyhedron
The polyhedrons can be generated by randomly picking the nodes on the spheres using spherical coordinate system as follows, xi = r cos(θ ) sin(ϕ ) + x 0 ȳi = r sin(θ ) cos(ϕ ) + y 0 zi = r cos(ϕ ) + z 0 (A.11)where xi , ȳi , zi are the coordinates of the ith vertex of the polyhedron.x 0 , y 0 , z 0 are the coordinates of centre of the sphere, r' is the radius of the sphere.
where θ is the polar angle measured from a fixed zenith direction, ϕ is the azimuth angle of its orthogonal projection on a reference plane that passes through the origin and is orthogonal to the zenith, measured from a fixed reference direction on that plane.r 0 and r 1 are the minimum and maximum radius in the size distribution, η, ζ , ξ are independent random numbers uniformly distributed between 0 and 1.
In order to make sure the randomness of generated polyhedrons, the number of vertexes of polyhedron is determined as uniformly distributed from N min to N max as where η is a random number uniformly distributed between 0 and 1.
The polyhedrons are generated by linking the nearest three vertexes to an area, and then make all the areas connected to a volume.To get a better quality mesh for finite element modelling, a minimum distance limitation between every two vertexes is set to 0.5r'.Fig. A2 shows a few examples of numerical polyhedrons with different random numbers.

Fig. 20 .
Fig. 20.Effect of the sample number on the peak stress.

Fig. 21 .
Fig. 21.Effect of the sample number on stress at displacement = 0.03.
via image-based modelling, as shown in Fig. 25.Figs.23 and 24 present the crack patterns at four loading stages.To clearly visualise the fracture surfaces, three different ways are used: models with damaged cohesive elements (left column, displacement magnification factor of 20), model without cohesive elements (middle column, displacement magnification factor of 200) and failed interfaces only (right column, displacement magnification factor of 20).

Fig. 23 .
Fig. 23.Typical Type I cracking evolutions at various loading stages shown in models with CIEs (left column), models without CIEs (middle column) and cracked interfaces (right column).

Fig. 24 .
Fig. 24.Typical Type II cracking evolutions at various loading stages shown in models with CIEs (left column), models without CIEs (middle column) and cracked interfaces (right column).

Fig. 26 .
Fig. 26.Comparison of stress-displacement curves for different types of cracking.(For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

Fig. 29 .
Fig. 29.Typical type I and type II cracking paths for samples with different aggregate and void shape (displacement magnification factor of 200).

Fig. 39 .
Fig.33shows that the load capacity is hardly affected by the critical stress of the aggregate-aggregate interfaces until it becomes equal to the strength of mortar-aggregate interfaces.This is because aggregates are stronger than mortar and interfaces, and the load capacity is not sensitive to the tensile strength of aggregates.The crack paths for different tensile strength of aggregate-aggregate interfaces are shown in Fig.34.It can be seen that, as the aggregate strength is reduced, cracks seem to a greater tendency to propagate across the aggregates leading to straighter paths with less bridging and branching.Fig.35illustrates that the tensile strength of mortar-mortar interfaces has a significant influence on the load capacity of concrete samples.The decrease of tensile strength of mortar-mortar interfaces results in a rapidly falling load capacity.Fig.36shows the crack paths corresponding to different tensile strength of mortar-mortar interfaces.The cracks show a tendency to bypass the aggregates and propagate in a straight path through mortar and voids when the tensile
Fig. 42 also shows that the crack paths are dependent on the cohesive energy of mortar-mortar interfaces.
1. Mechanical response and crack propagation in concrete are affected strongly by the meso-structure: load capacity decreases with increasing aggregate or void content and increasing aggregate roughness; pre-peak response and peak stress are relatively insensitive to the aggregate and void spatial distribution; postpeak softening response is sensitive to the spatial distribution of aggregates and voids.2. Two macro-crack (failure) patterns are observed in concrete under uniaxial tension regardless of the shape of aggregate and void: Type I -micro-cracks coalesce into one macro-crack propagating to failure.The result is rapid softening response and smaller failure energy density (apparent toughness of representative volume); Type II -micro-cracks coalesce into two or more macrocracks propagating concurrently to failure.The result is slower softening response and larger failure energy density (apparent toughness).
MCS are carried out to investigate the effects of main parameters, such as, shape, size, spatial distribution and volume fractions of aggregate and void, and cohesive laws, on mechanical behaviour of concrete in a quantitative manner.The outcomes are presented comprehensively and provide valuable information for improved understanding of damage and failure mechanisms of concrete.

Table 2
Material properties.

Table 3
Typical mortar-aggregate interface element number of specimens with different aggregate and void shape.Figs.31,33 and 35present the simulation results obtained by varying the critical stress of interface elements.For each curve in the figures, the legend includes three values, which correspond to the Fig. 30.Effect of aggregate and void shape.4.4.1.Effect of critical stressFig.31.Effect of tensile critical stress of mortar-aggregate interfaces.criticalstress assigned to the mortar-aggregate, aggregate-aggregate and mortar-mortar interfaces, respectively.A missing value for the aggregate-aggregate critical stress indicates that cracking of the aggregate is not allowed.