Monte Carlo simulations of mesoscale fracture modelling of concrete with random aggregates and pores

(cid:2) A procedure is devised to generate mesoscale concrete samples with random multi-phases. (cid:2) Complex crack initiation and propagation is realised using cohesive interface elements. (cid:2) Samples in tension fail with 1 or 2 cracks, regardless of aggregates’ and pores’ shape and fraction. (cid:2) The effects of aggregate shape and porosity should not be neglected in meso-modelling. A procedure for generating two-dimensional heterogeneous meso-scale concrete samples is developed, in which the multi-phasic features including the shape, size, volume fraction and spatial distribution of aggregates and pores are randomised. Zero-thickness cohesive interface elements with softening traction–separation relations are pre-inserted within solid element meshes to simulate complex crack initiation and propagation. Extensive Monte Carlo simulations (MCS) of uniaxial tension tests were carried out to investigate the effects of key multi-phasic features on the fracture patterns and load-carrying capacities. It is found that the fracture behaviour and stress-displacement responses of the numerical specimens are highly dependent on the random mesostructures, especially the post-peak softening responses. The specimens fail with either one or two macro-cracks, regardless of the shapes and volume fractions of aggregates and pores. Assuming that the aggregate–mortar interface is weaker than the mortar, using polygonal rather than circular or elliptical aggregates, or increasing the aggregate volume fraction will reduce the tensile strength of specimens. The porosity is found to have severely adverse effects on the specimen strength and cannot be neglected in mesoscale fracture modelling of concrete. (cid:2) 2014 Elsevier Ltd. All reserved.


a b s t r a c t
A procedure for generating two-dimensional heterogeneous meso-scale concrete samples is developed, in which the multi-phasic features including the shape, size, volume fraction and spatial distribution of aggregates and pores are randomised. Zero-thickness cohesive interface elements with softening traction-separation relations are pre-inserted within solid element meshes to simulate complex crack initiation and propagation. Extensive Monte Carlo simulations (MCS) of uniaxial tension tests were carried out to investigate the effects of key multi-phasic features on the fracture patterns and loadcarrying capacities. It is found that the fracture behaviour and stress-displacement responses of the numerical specimens are highly dependent on the random mesostructures, especially the post-peak softening responses. The specimens fail with either one or two macro-cracks, regardless of the shapes and volume fractions of aggregates and pores. Assuming that the aggregate-mortar interface is weaker than the mortar, using polygonal rather than circular or elliptical aggregates, or increasing the aggregate volume fraction will reduce the tensile strength of specimens. The porosity is found to have severely adverse effects on the specimen strength and cannot be neglected in mesoscale fracture modelling of concrete.

Introduction
Concrete is a composite material with multiple phases including mortar, aggregates, interfaces and various defects such as pores and weak inclusions. Its mechanical behaviour is a result of multiple mechanisms occurring at different length scales, from macro, meso to microscale and even atomic or sub-atomic scales. Therefore, multi-scale experiments and simulations of concrete has become a hot research area offering more accurate insights to the mechanical behaviour of concrete material and structures [1][2][3][4]. A current focus has been stochastic heterogeneous modelling of concrete considering its multiple phases and interactions at micro and meso scales [5][6][7][8][9][10][11][12][13]. This is desirable, as such small-scale models take into account accurate internal structures so that more fundamental understanding of damage and failure mechanisms can be gained, and more representative homogenised material properties can be used in higher-scale models [13][14][15][16]. As concrete is generally used in large-sized structures, meso-scale rather than micro-scale modelling seems a more feasible choice and will be the topic of this study.
There are basically two approaches to generate meso-scale models of concrete: the digital image based approach and the http://dx.doi.org/10.1016/j.conbuildmat.2014.09.069 0950-0618/Ó 2014 Elsevier Ltd. All rights reserved. parameterization modelling approach. In the former, cameras, microscopes, or more advanced 3D techniques such as nuclear magnetic resonance spectroscopy and X-ray computed tomography (XCT) scanners are used to capture digital images of the mesostructures; the images are then converted into finite element (FE) meshes modelling multi-phases with real sizes, shapes and distributions [12,[17][18][19]. The image-based models are real reflections of the material mesostructures and thus offer tremendous potential. However, it is still costly and time-consuming to conduct 3D tests, and even more so when a large number of samples are usually needed to make meaningful statistical analyses. In addition, image processing and mesh generation of complicated mesostructures is not a trivial task.
In the parameterization modelling approach, indirect and direct algorithms can be used to generate mesostructures. In the indirect algorithms, different phases of concrete are not explicitly modelled; instead, the heterogeneous material properties are modelled as spatially-varying random fields assigned to conventional FE meshes (e.g. [6,7]), or by lattice elements randomly assigned as aggregates and mortar phases in lattice models (e.g. [19,20]). In the direct algorithms, aggregates of different shapes and sizes are generated, randomly packed into a space, and incorporated with mortar to form digital concrete specimens (e.g. [8][9][10][11]). The indirect algorithms are able to generate a large number of samples with ease and thus are very promising for global statistical analyses. The direct algorithms can take into account key multi-phasic parameters such as the shape, size, gradation and spatial distribution of pores and aggregates, phase volume fractions and aggregatemortar interfaces, and their effects on mechanical behaviour of specimens. This makes the direct algorithms attractive in mesoscale modelling, especially when accurate understanding of detailed failure mechanisms is required.
Recent XCT scanned images [19,21,25] show that a large number of pores exist in concrete at micro and meso scales. The pores may have significant effects on the damage and fracture behaviour of concrete specimens. However, they are rarely considered in existing meso-scale models, in which concrete is treated as either a two-phase (mortar and aggregates) or three-phase (mortar, aggregates and interfaces) material.
This paper presents a systematic study of statistical effects of key multi-phasic parameters on the complex fracture behaviour and load-carrying capacities of concrete specimens, through extensive Monte Carlo simulations (MCS) of direct meso-scale finite element models. A generate-and-place procedure similar to [10] is developed to generate 2D mesostructures with randomly packed aggregates and pores of different shapes and fractions. The mesostructures are then meshed automatically using the pre-processing functionalities in ANSYS. Zero-thickness cohesive interface elements governed by nonlinear normal and shear traction-separation laws are then inserted within aggregate, mortar and their interfaces to simulate complex crack initiation and propagation, using a procedure augmented from that in [6]. Finally, the nonlinear finite element models are solved using ABAQUS. The MCS results are critically analysed, leading to valuable statistical results that may help improve designs of concrete material and structures.

Size distribution of aggregates and pores
The optimal gradation or size distribution of aggregates is often determined by the Fuller curve as [13] PðdÞ ¼ 100 where P(d) is the cumulative percentage passing a sieve with aperture diameter d, d max is the maximum size of aggregates and n is a constant (n = 0.45-0.70), respectively. Eq. (1) can be discretised into a number of segments, and the area of aggregates within a grading segment [d i , d i+1 ] is [10] A agg ½d iþ1 À d i ¼ where d min is the maximum and minimum sizes of aggregates, P agg is the area fraction of all aggregates with respect to the total area of concrete sample A.
The aggregates generally occupy 60-80% in volume of the concrete. The aggregate size distribution in [22] is used in this study. Table 1 lists the sieve size ranges with passing and retaining percentages, with the corresponding four-segment gradation curve shown in Fig. 1. For simplicity, only coarse aggregates larger than 2.36 mm are modelled in this study, while the large number of smaller fine aggregates together with the cement matrix is treated as mortar, whose mechanical properties are assumed to be uniform. For normal strength concrete, coarse aggregates usually represent 40-50% of the concrete volume [13].
For simplicity, it is assumed that the pores are circular or elliptical, with a single size range 2-4 mm. The pore area fraction is specified by the porosity P pore .

Aggregate and pore generation
A similar procedure presented by Wang et al. [10] is devised. The basic idea is to generate and place aggregates and pores in a repeated manner, until the target area is fully packed. The aggregate shape depends on its type, namely, circular or elliptical for gravels and polygonal for crushed ones. The detailed procedure is as follows.  (1) Input controlling parameters. The key parameters include the aggregate type and shape, the Fuller curve, the aggregate area/volume fraction P agg , the pore shape, the porosity P pore , the size range of pores (d 0 = 2 mm, d 1 = 4 mm), the minimum distance of aggregates/pores to the specimen boundary c 1 , and the minimum distance between aggregates/pores or film thickness c 2 . The minimum distance c 2 reflects that each aggregate is coated all around with a layer of mortar. The aspect ratio ranges for elliptical aggregates (R 1 ) and pores (R 2 ), and the number of sides (N) and angle range (h) for polygonal aggregates, are also specified. (2) Generate a pore. Assuming a uniform distribution between d 0 and d 1 , the pore size is calculated by d where g is a random number between 0 and 1. If the pore is elliptical, a random orientation angle is also generated.
(3) Place the pore. Random numbers are generated to define the position (and orientation if elliptical) of the newly generated pore, and the following conditions are checked: (1) the whole pore must be inside the concrete area and (2) there is no overlapping/intersection between this pore and all the existing ones. For circular pores, the condition (2) can be checked easily by comparing the distance between pore centres and the sum of two radii. For elliptical and polygonal pores, the search algorithm in [23] and the convex hull algorithm in MATLAB are used. When both conditions are satisfied, the pore is placed; or otherwise another set of random numbers are generated to place the pore in a new position. (4) Generate and place all pores. Steps 2 and 3 are repeated until the area of pore remaining to be generated is less than pd 0 2 /4, i.e. there is no space to place another pore.
where g is a random number between 0 and 1.
(6) Place the aggregate. Random numbers are generated to define the position (and orientation if elliptical) of the aggregate. Four conditions are checked before placement: (1) the whole aggregate must be within the concrete area; (2) there is no overlapping/intersection between this aggregate and any existing pores and aggregates; (3) there is a minimum distance between the edge of an aggregate (c 1 ) and the specimen boundary; and (4) there is a minimum gap between any two aggregates (c 2 ). The last two conditions are checked after the aggregate size is enlarged to (1 + c 1 ) and (1 + c 2 ) times its size, respectively. If either condition is violated, this aggregate is disregarded and another set of random numbers are generated to place it in a new position. When all the conditions are satisfied, the aggregate is placed.
(7) Generate and place all aggregates. Steps 5 and 6 are repeated until the remaining area to be generated is less than pd i 2 /4.
The remaining area is then transferred to the next smaller grading segment.
Step 5 to step 7 are repeated for all other size grading segments until the last aggregate of smallest size is generated and placed.
Using the above generation-placing procedure, numerical concrete specimens consisting of various shapes, sizes and distributions of aggregates/pores can be built with ease. All the numerical specimens shown in this paper are 50 mm squares, and the aspect ratio ranges R 1 = R 2 = [2, 2.5] and the film thickness c 1 = c 2 = 0.5 mm are used when elliptical aggregates and pores exist. For crushed polygonal aggregates, N = [3,7] and h = [15°, 180°] are used. The 4-segment aggregate gradation curve in Fig. 1 are used for all the samples. Fig. 2 shows a few examples of numerical samples with combinations of different shapes of gravel aggregates (P agg = 45% area fraction) and pores of varying porosity. Fig. 3 shows numerical samples with the same porosity 2% and aggregate area fraction 40%, 45% and 50%. Numerical samples with the same aggregate area fraction 40% and the porosity 1%, 3% and 6% are shown in Fig. 4. The effect of aspect ratios of elliptical aggregates can be seen in Fig. 5. Some numerical specimens with crushed aggregates are shown in Fig. 6, compared with an XCT image of a real concrete specimen [21].

Finite element mesh generation and fracture modelling
The pre-processing functionalities of ANSYS are utilized to mesh the numerical samples. As the regions of aggregates and pores are already known, the mortar regions are identified by applying the ''overlapping'' Boolean operations. The samples are then meshed using solid FE elements. The mesh generation process is automated by running a batch file of APDL programs, so that a large number of samples can be meshed quickly for Monte Carlo simulations and statistical analyses.
To simulate realistic fracture processes, 4-noded cohesive interface elements (CIEs) with zero in-plane thickness are then inserted 50 mm 50 mm   into the generated solid FE mesh. The detailed CIE insertion procedure devised for homogeneous materials in [6] is extended to account for multi-phases and interfaces. Three sets of CIEs with different traction-separation softening laws are inserted, namely, CIE_AGG inside the aggregates, CIE_MOR inside the mortar, and CIE_INT on the aggregate-mortar interfaces. As the aggregates have much higher strength than the cement and the interfaces in normal concrete, no cracks are allowed to initiate inside the aggregates by assuming elastic behaviour without damage in CIE_AGG. However, it is possible to model crack propagation through aggregates (e.g. in lightweight concrete) by assigning damage properties to CIEs in aggregates. A final FE mesh (element size = 1 mm) after the insertion of CIEs is shown in Fig. 7. The CIEs are highlighted as red lines. Triangular solid FE elements are used because they are flexible in modelling realistic, smooth crack paths.

Monte Carlo simulations
Uniaxial tension tests of 50 mm Â 50 mm numerical specimens were modelled in this study (see Fig. 8). All the models were fixed at the left boundary and were subjected to a uniformly distributed displacement at the right boundary, i.e. a displacement-controlled loading scheme was used. All analyses were ended at a displacement d = 0.1 mm or strain e = 0.002.
The solid elements for aggregates and mortar were assumed to behave linear elastically. The linear tension/shear softening laws were used to model CIEs [6] with quadratic nominal stress initiation criterion and linear damage evolution criterion. For comparison of results, the same material properties as in [8] were used in this study. They are listed in Table 2. Due to the lack of experimental data, the shear fracture properties were assumed to be the same as the normal ones.
Extensive Monte Carlo simulations were conducted to investigate the effects of key multi-phasic parameters on the statistical responses of numerical specimens. The reference samples have a combination of elliptical aggregates and elliptical pores with P agg = 40% and P pore = 2% (Fig. 7). For each MCS, 100 samples were modelled to ensure that the results were statistically converged. A typical MSC using the ABAQUS/Explicit solver and a mesh in Fig. 7 took about 10 h by parallel computation using 12 Intel Xeon CPUs @ 2.8 GHz.

Mesh dependence
The method of pre-inserting CIEs may artificially weaken the sample although a very high elastic stiffness is used for CIEs  ( Table 2). The crack pattern may also be mesh-dependent as only CIEs can become cracks. Therefore, a mesh convergence study was conducted first. Three meshes with element length L e = 1 mm, 0.75 mm and 0.5 mm, respectively were modelled for a particular sample shown in Fig. 9. The mesh 1 has 20,346 nodes, 6782 solid elements and 10,004 cohesive elements. The numbers for the mesh 2 and 3 are 33,288, 11,096 and 16,423, and 71,652, 23,884 and 35,504, respectively. Fig. 10 shows that the final crack paths are very similar with little mesh-dependence. The cracks are represented by red cohesive interface elements with damage index D P 0.9 (D = 1 means complete failure). A magnification factor 10 was used for all deformed meshes unless specified otherwise.
For each element length, 100 numerical samples with random distribution of elliptical aggregates and pores were then generated and analysed. The mean stress-displacement (r-d) curves for each element length are compared in Fig. 11. The stress is the total horizontal reaction force of all the left boundary nodes divided by the specimen cross-section area. It is clear that there exists little mesh dependence as the curves are almost coincident. Considering the balance between accuracy and efficiency, the element length 1 mm (mesh 1) was used in all the meshes in the following simulations.

Effects of loading time
When the ABAQUS/Explicit solver is used to model quasi-static loading conditions, the loading/step time must be long enough to   Monte Carlo sample number Mean Peak Stress (MPa) Fig. 17. Effect of the sample number on the mean peak stress. minimise dynamic effects. However, a long loading time results in high computational cost. Fig. 12 shows the r-d curves from 100 random samples for three loading times. It can be seen that 0.001 s loading time results in dynamic oscillations in the curves, while smooth and identical curves are obtained for 0.01 s and 0.1 s. Therefore, a loading time 0.01 s was used in all the following simulations. Fig. 13 shows the r-d curves computed from 100 samples with random distribution of elliptical aggregates and pores. The mean curve, the mean and standard deviation of the peak stress (strength) are also shown. The scatter demonstrates the necessity of Monte Carlo simulations with random distributions of phases. Fig. 14 shows the probability density and the best fit Gaussian probability density function (PDF) of the peak stress, which can be used to calculate structural reliability or failure probability against given external loadings and material properties for structural design [6]. Fig. 15 shows the variation of standard deviation in stress with displacement. Before the peak stress is reached at about d = 0.01 mm (see Fig. 13), the standard deviation of the stress is lower than 0.1 MPa. It increases to over 0.5 MPa at d = 0.02 mm in the post-peak part and decreases afterwards. So the softening response is more sensitive to pore and aggregate distributions than the peak stress. Figs. 16 and 17 show the effect of the number of samples on the mean and standard deviation of peak stress, respectively. It is clear that 100 random samples are enough to achieve statistical convergence. This is consistent with the conclusion drawn from an MSC study using the indirect modelling approach based on Weibull random fields [6].

Stress-displacement curves and effect of sample number
In Fig. 18, the mean r-d curve from 100 samples is compared with the results from the experiment by Hordijk [24] and the FE simulation by Lopez et al. [8]. It can be seen that the peak loads and the post-peak softening stages of three curves are close. However, one should not intend to compare the curves directly, as they are responses of specimens with different sizes and different distributions and volume fractions of phases.

Crack propagation and fracture patterns
Two typical fracture patterns were observed in the 100 samples under uniaxial tension. In all the samples, many microcracks initiate in the early stages of loading on the aggregate-mortar interfaces due to their lower fracture properties than the mortar. (a) d=0.009mm (peak stress) (b) d=0.016mm (softening) (c) d=0.035mm (later stage) (a) d=0.009mm (peak stress) The microcracks at peak loads are shown in Fig. 19(a) and 20(a), respectively. It can be seen that at the peak load, there is still no evident macrocracks. The predicted post-peak macrocrack propagation process is illustrated in Figs. 19(b-c) and 20(b-c), in which the different phases are also shown. As the displacement increases, some aggregate-mortar interfacial cracks continue to propagate and are gradually coalesced with newly formed cracks in the mortar. The specimen fails either with one dominant crack (type I in Fig. 19(c)) or with two (type II in Fig. 20(c)). It should be noted that the microcracks in Figs. 19(a) and 20(a) still exist, but they are not shown in Figs. 19c-d and 20c-d because their width is much smaller than that of the macrocracks. A deformation magnification factor = 50 is used in Figs. 19 and 20 to clearly illustrate the detailed microcrack initiation and propagation process. Among the 100 numerical samples investigated, 58 samples behave as Type I cracking and 42 as Type II cracking. Fig. 21 shows the r-d curves separately for both types of cracking. It can be seen that the pre-peak responses are almost identical, and there is little difference in the peak stress. However, the post-peak stress drops more quickly in Type I than Type II, leading to lower dissipated energy in Type I cracking. This behaviour may be attributed to smaller fracture area in the single crack in the former than the two cracks in the latter. Fig. 22 shows the mean r-d curves and corresponding peak stress for four aggregate volume fraction P agg ranging from 20% to 50% with 100 samples each. It can be seen that the initial elastic stiffness of the specimens increases as more aggregates are added, which is expected. It is surprising to see that the computed mean tensile strength decreases. This is because the tensile strength is not determined by the number or strength of aggregates, but by the cohesive strength of CIEs that form the microcracks (see Figs. 19a and 20a). A higher P agg means more aggregates and thus more aggregate-mortar interfacial CIEs in the sample. Therefore, more aggregate-mortar interfacial CIEs and fewer mortar CIEs   become microcracks to resist the external loading. Because it is assumed that the cohesive strength of the interfacial CIEs is only half that of the mortar CIEs (Table 2), the specimen strength is reduced as P agg increases. This conclusion only holds if such an assumption is valid. The effects of porosity on the load-carrying capacity are shown in Fig. 23. A higher porosity results in lower strength, because more pores reduce the effective area to resist the tensile force. An increase of porosity from 0% to 6% reduces the strength by 26% (3.80-2.80 MPa).

Effects of aggregate and pore shapes
Apart from the reference samples using elliptical aggregates and pores, additional five combinations of different shapes of aggregates and pores were investigated, with 100 samples generated and modelled for each combination. Some typical crack patterns are shown in Fig. 24. It can be seen that both types of cracking occur regardless of the shapes of aggregates and pores.     Table 3) that tend to reduce the tensile strength as the aggregate volume fraction does, as discussed in Section 3.5. It may also be caused by the higher stress concentration at the corners of the polygonal aggregates, while the smooth edges of the circular and elliptical aggregates have a more benign stress distribution which delays the fracture event and increases the tensile strength.

Conclusions
Numerical models of concrete with random mesostructures comprising circular, elliptical, or polygonal aggregates and circular or elliptical pores have been developed in this study. The complex mesoscale crack initiation and propagation was realistically simulated using the technique of pre-embedding cohesive interface elements. Extensive Monte Carlo simulations of uniaxial tension were carried out to investigate the statistical effects of shapes, volume fractions and random distributions of aggregates and pores. The main conclusions are: (1) The fracture behaviour and stress-displacement responses of the numerical specimens are highly dependent on the random mesostructures. The computed pre-peak responses and the peak stress (strength) are insensitive to the distribution of aggregates and pores, but the post-peak softening responses are much more sensitive. Therefore, Monte Carlo simulations of a sufficient number of samples should be simulated for accurate understanding of the post-peak responses; (2) two cracking types were noticed in the Monte Carlo simulations of 50 mm square specimens under the assumed uniaxial tensile boundary conditions, regardless of the shapes and volume fractions of aggregates and pores, namely, Type I with one dominant macrocrack and Type II with two. The latter shows more progressive post-peak softening responses than the former; (3) assuming that the aggregate-mortar interface strength is lower than the mortar strength, using polygonal rather than circular or elliptical aggregates, or increasing the aggregate volume fraction will reduce the tensile strength of specimens; and (4) the porosity has severely adverse effects on the specimen strength, and the pores cannot be neglected in mesoscale fracture modelling of concrete.