Monte Carlo simulations of mesoscale fracture of concrete with random aggregates and pores: a size effect study Building

(cid:2) Size effect of concrete studied by realistic mesoscale models and Monte Carlo simulations. (cid:2) Key meso-control parameters (aggregate fraction and porosity) considered. (cid:2) Weibull size effect laws of strength obtained considering the key parameters. (cid:2) Equations relating strength and the key parameters obtained. (cid:2) Both aggregate fraction and porosity play signiﬁcant role. Size effect in concrete under tension is studied by Monte Carlo simulations of mesoscale ﬁnite element models containing random inclusions (aggregates and pores) with prescribed volume fractions, shapes and size distributions (called meso-structure controls). For a given size and a set of controls, a number of realisations with different spatial distribution of inclusions are simulated to produce statistical data for macroscopic load/stress–strain curves. The complex meso-crack initiation and propagation is captured by pre-inserted cohesive interface elements. The effects of specimen size and meso-structure controls on macroscopic strength and toughness are analysed, and empirical size-effect laws for their dependences are proposed by data regression. It is also shown that the mesoscale porosity affects both strength and toughness and should not be ignored in size effect studies of concrete. (cid:2) 2015 The Authors. Published by Elsevier Ltd. ThisisanopenaccessarticleundertheCCBYlicense(http://


Introduction
Concrete is a composite material with random arrangements of mesoscale (mm) constituents of different mechanical properties. The effect of the resultant heterogeneous meso-structure on the mechanical behaviour is strong when the specimen size is comparable to the sizes of characteristic constituents, making homogenisation questionable. The effect weakens and the homogenisation becomes increasingly valid as the specimen size increases. This is a ''structure-property'' explanation of the measurable changes in mechanical behaviour with specimen size, known as the size effect [1,2]. The size effect of concrete has been widely investigated experimentally and numerically since 1990s [1,[3][4][5][6][7][8][9][10]. It has been shown that the nominal strength of concrete can decrease significantly with increasing specimen size, which is now regarded as an inherent concrete behaviour [11][12][13].
Different theories have been developed to explain the size effect, amongst which the best known are the works of Weibull [14], Bazant [15] and Carpinteri [16]. It appears that the size effect needs to be investigated by a bottom-up approach, with explicit account of concrete non-homogeneous structure, replicated in specimens of increasing sizes. This is an experimentally challenging task, but can be tackled by increasingly realistic numerical models. For example, Van Mier et al. [1,17] and Grassl et al. [18,19] used discrete lattice model to study the size effect of quasi-brittle materials. Karihaloo et al. [20] and Duan et al. [6] used the fictitious crack model to investigate the size effect on concrete. However, mesoscale numerical studies of the size effect of concrete with statistical analyses are still challenging [21][22][23], mainly due to the complex multi-phase composition, the complicated nonlinear multi-cracking behaviour, and substantial computational cost from repetitive simulations. In this paper, we study the size effect in concrete under tension by Monte Carlo simulations (MCS) of mesoscale specimens of varying sizes. The direct parameterisation modelling technique [24][25][26] is used to build mesoscale concrete models with randomly-distributed aggregates and pores generated according to prescribed size distributions, shapes and volume fractions. The four different phases, i.e., mortar, aggregates, interfaces and pores are explicitly modelled, allowing for examination of their relative influences on meso/macro mechanical responses. The complicate meso-crack initiation, propagation and coalesce into macro-cracks is captured by pre-inserted cohesive interface elements into solid finite element (FE) meshes [23,27]. The detailed procedure and initial results have been presented in [23].
The size effect of concrete nominal strength is mostly analysed by pre-notched beams under three-point bending [20,[28][29][30]. However, such loading conditions complicate the analysis and understanding of the size effect on the failure processes due to the strain and stress gradients introduced. Therefore, this study is focused on analyses of specimens without notches under displacement-controlled uniaxial tension. Such analyses replicate standard tensile tests, where stress concentrator is not present initially. This allows us to study a spatially distributed micro-cracking prior to localisation into a macroscopic crack, rather than to predefine localisation with the presence of a concentrator. The size effect caused by stochasticity, namely more defects are included as the size increases, and the deterministic size effect caused by nonlinear fracture process, were both taken into account. The macroscopic properties of interest in this study are strength and toughness. These are calculated by statistical analyses of a number of realisations. The effects of specimen size and the main meso-structure controls (volume fractions and size distributions of aggregates and pores) on these properties are analysed and empirical expressions for their dependences are proposed.
The outline of the paper is as follows. In Section 2 we present the techniques for generation of inhomogeneous meso-structures and the basics of finite element modelling with cohesive interfacial elements. In Section 3 we report and discuss the results of size effects in concrete with different prescribed parameters and propose several size effect laws based on them. The main conclusions are drawn in Section 4.

Model and method
The detailed procedure of generating numerical concrete samples and inserting cohesive interface elements into solid FE meshes has been given in [23] and [27]. Herein only an outline is presented for the convenience of discussion.

Generation of numerical concrete samples
The size distribution of aggregates in concrete is often described by the Fuller curve [31], which is discretised into a number of segments. The aggregate size distribution reported by Hirsch [32] and summarised in Table 1 is used in this study. The concrete contains coarse and fine aggregates and the cut-off size between coarse and fine aggregates is taken to be 2.36 mm. Here only coarse aggregates are explicitly modelled as meso-scale features. The large number of fine aggregates together with the cement matrix is regarded as mortar with homogenised uniform mechanical properties. The coarse aggregates are considered to have elliptical shape. In most concretes, the volume density of coarse aggregates is between 40% and 50% [33]. Elliptical pores are introduced with uniformly distributed sizes in the range from 2 to 4 mm. The prescribed volume fractions of aggregates and pores (i.e., porosity) are met by generating spatially-distributed random aggregates and pores in a repeated manner until a target area is filled [23].
2D square numerical samples are generated. Fig. 1 shows five typical models with different sizes or length L and Fig. 2 shows five samples with L = 50 mm, using the following parameters: aggregate volume fraction P agg = 40%, porosity P pore = 2%, aspect ratio for elliptical aggregates and pores R 1 = R 2 = [2, 2.5], minimum space between the edge of an aggregate and the specimen boundary c 1 = 0.5 mm and minimum gap between any two aggregates/ pores c 2 = 0.5 mm.

Finite element mesh generation and fracture modelling
The numerical concrete models are then meshed in ANSYS automatically 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. The 4-noded cohesive interface elements (CIEs) with zero in-plane thickness are then inserted into the generated solid FE mesh (triangular plane stress elements) by an augmented procedure devised for homogeneous materials in [27] to account for multi-phases and interfaces. CIEs with different traction-separation softening laws are inserted inside the aggregates, inside the mortar, and on the aggregate-mortar interfaces, respectively, to simulate the complicated nonlinear fracture behaviour [23,34].
The cohesive element COH2D4 with zero in-plane thickness in ABAQUS is used in this model. Its constitutive behaviour is described by a damage initiation criterion and a damage evolution law. A bilinear cohesive zone model, illustrated in Fig. 3, was used in this work.
The damage is assumed to initiate when the following condition is met: where h i is the Macaulay bracket, i.e., The damage evolution is characterised by a scalar parameter, D, representing the overall extension of the crack across the element caused by all physical mechanisms. It is defined in terms of effective relative displacement d m given by: The definition of the damage parameter D is: where d m,max is the maximum effective relative displacement attained during the loading history, and d m0 and d mf are effective relative displacements corresponding to d n0 and d s0, and d nf and d sf in Fig. 3, respectively. D evolves monotonically from 0 to 1 upon further loading after the initiation of damage. The damage initiation and evolution degrades the unloading and reloading stiffness coefficients, k n and k s , calculated by: The tractions are also affected by damage according to: where t n and t s are the traction components predicted by the elastic traction-displacement behaviour for the current separation without damage.
Material properties, such as density, Young's modulus, Poisson's ratio, tensile strength and fracture energy, are set for continuum elements of aggregates and mortar, and for the three sets of interface cohesive elements. The material heterogeneity is represented by different phases with corresponding material properties.

Definition of macroscopic properties and statistical analyses
The size effects on the tensile strength and toughness are investigated in this study. The strength is defined as the mean value of peak stresses (the peak load divided by the cross-sectional area) on the stress-strain curves, obtained from fracture modelling of a number of numerical samples for each size. The toughness is defined as the mean value of energies consumed per unit volume up to failure, i.e., the area under the stress-strain curve up to the ultimate failure from a number of numerical samples for each size, measured in J/m 3 .
The Monte Carlo simulation results are evaluated statistically. The standard deviation (SD) s of n number of results x i (i = 1, n) is defined as: where x the mean of x i . The coefficient of variation (CoV = s/ x) is also calculated, to represent the fluctuation of measured property relative to its average value.

FE modelling setups
Uniaxial tension tests of square numerical specimens were modelled in this study (see Fig. 4). The out-of-plane thickness  was unit and plane stress was assumed. All the models were fixed at the left boundary and subjected to a uniformly distributed displacement at the right boundary, i.e., a displacement-controlled loading scheme was used.
The solid elements for aggregates and mortar were assumed to behave linear elastically. The linear tension/shear softening laws were used to model CIEs with quadratic nominal stress initiation criterion and linear damage evolution criterion. Similar material properties as in [24] were used in this study. Young's modulus was 7 Â 10 4 MPa for aggregates and 2.5 Â 10 4 MPa for mortar. Poisson's ratio of both aggregates and mortar was 0.2. The cohesive elements inside mortar had elastic stiffness k n = k s = 10 6 MPa/mm, tensile strength t n = 6 MPa and fracture energy G f = 0.06 N/mm. The cohesive elements on aggregatemortar interfaces had k n = k s = 10 6 MPa/mm, t n = 3 MPa and G f = 0.03 N/mm. Elastic behaviour without damage was assigned to the cohesive elements inside aggregates so that no cracks were allowed to initiate inside aggregates as they are much stronger than other phases in normal concrete. Due to the lack of experimental data, the shear fracture properties were assumed to be the same as the normal ones.
Extensive   statistically converged. Based on our previous study of mesoscale fracture modelling [23], the element length of 1 mm, the loading time of 0.01 s used in the ABAQUS/Explicit solver, and the final displacement of 0.1 mm were used for all sizes of samples considering the balance between accuracy and efficiency. The presented model has been validated previously [23] by comparison of numerical results with available experimental data and with other numerical studies. Table 2 shows the average finite element numbers of numerical samples of the five sizes. It can be seen that a considerable number of nonlinear cohesive interface elements are inserted to simulate the fracture process, making computation time-consuming. A typical MCS for L = 50 mm took about 5 h by parallel computation using 12 Intel Xeon CPUs @ 2.8 GHz.

Results for samples with L = 50 mm
The failure evolution simulated for a typical sample with L = 50 mm is illustrated in Fig. 5. Fig. 5a-d shows the deformed meso-structure at elastic stage, at peak stress, in softening branch, and prior to sample disintegration, respectively. Displacement magnification factors of 800, 400, 200 and 30 were used, respectively, to illustrate clearly the meso-crack initiation and propagation process. This evolution is typical for all 50 realisations with variations reflected by statistical analysis of peak stress and toughness as follows.
The stress-strain curves of all 50 samples with the mean curve are plotted in Fig. 6a. The effects of the sample number on the mean and the CoV of peak stress can be seen in Fig. 6b and c, respectively. The results show that 50 samples are enough to obtain a convergent mean peak stress with a stable fluctuation below 4%.
The results with respect to the toughness are shown in Fig. 7a-c, respectively. The number of samples (50) is sufficient to reach a converged mean toughness with a stable fluctuation of about 19%.
The mechanical behaviour of concrete is often characterised by a single specimen or a small number of specimens [24,26,[35][36][37][38]. The results here show that the stress-strain responses vary with different realisations of specimens, especially in the post-peak stage. Therefore the conclusions based on a small number of analyses could be drawn for the elastic stage only, and substantially more realisations are necessary to yield conclusive results for softening behaviour.
The influence of the sample number on the CoV of peak stress and toughness are shown in Fig. 8a and b respectively for samples of all sizes. It can be seen that 50 samples are enough to reach stationary responses for all sizes.    Fig. 9 compares the mean stress-strain curves from 50 realisations each for the specimens of 5 sizes with elliptical inclusions (P agg = 40%, P pore = 2%). It is evident that the strength (mean peak stress) decreases monotonically as the size increases. It can also be seen that the slopes of the post-peak softening parts decrease with increasing size, indicating more brittle behaviour of larger specimens. This is consistent with experimental observations [9,13]. Table 3 summaries the mean, standard deviation and coefficient of variation of the peak stress and the toughness. The same data are also shown in Fig. 10. As the CoV of peak stress decreases when sample size increases, the effect of random aggregate and pore distribution on the peak stress decreases.

Size effects of reference samples
The results in Fig. 10 were curve-fitted using the least square method according to Weibull's size effect law [14] ln where r is the mean nominal strength, D the specimen size, and a and b are regression coefficients. a = 1.52 and b = À0.07 were found with correlation factor R 2 = 0.94, indicating a good agreement with the Weibull theory (Fig. 11). The results with the calculated coefficients match well with Weibull size effect law for Weibull modulus m = 24. This comparison provides a support for the soundness of the proposed meso-structure modelling and Monte Carlo simulation approach.
The size effect on toughness is shown in Fig. 12. It can be seen that the sample size has significant effects on the toughness. It is three times more for L = 12.5 mm than L = 62.5 mm. It should be noted that the CoV of toughness is relatively high (about 17.5%) even after sampling 50 different realisations.

Influence of aggregate volume fraction
Concrete samples with L = 25, 37.5, 50 and 62.5 mm containing elliptical aggregates and pores were studied with variable aggregate volume fraction P agg . The size L = 12.5 mm is too small to generate samples with large aggregates in Fuller curve when a low P agg is assumed and thus not considered. All other meso-structure parameters were the same as aforementioned in Section 3.1. Typical samples with P agg = 20-50% are shown in Fig. 13. Fig. 14 shows the size effect on strength for P agg = 20%, 30%, 40% and 50%, respectively.
To derive a functional relation between the strength and the aggregate volume fraction, additional samples of different sizes with pure mortar (2% pores) were simulated. Fig. 15 shows the relations between the peak stress and the aggregate volume fraction for different sizes. By curve-fitting, it is found that the following quadratic function can well describe the effect of P agg for all the sizes: 12.5mm×12.5mm 25mm×25mm 37.5mm×37.5mm 50mm×50mm 62.5mm×62.5mm Fig. 9. Comparison of stress-strain curves (P agg = 40%, P pore = 2%).
(a) Mean peak stress (b) CoV of peak stress CoV of peak stress (%) Sample size (mm) Fig. 10. Size effect on peak stress. where r agg 0 is the strength of pure mortar sample, and c is an empirical coefficient (c = 0.73 in this study). It should be noted that Eq.
The size effect laws (Eq. (2)) for P agg = 20-50% were computed by the least square regression and shown in Fig. 16. It can be found that the slope decreases as P agg increases, indicating a reduction of size effect for larger aggregate volume fractions. This may be due to the smaller variation in aggregate spatial distribution for higher P agg .
The size effect on toughness for different P agg is shown in Fig. 17. For a given size, an increase of the aggregate volume fraction results in a decrease of toughness.

Influence of porosity
Numerical samples with varying porosity P pore = 0-6% were then investigated with all other meso-structure parameters kept the same as in Section 3.1. Fig. 18 shows a few typical samples with L = 50 mm.
The size effect on strength for samples with P pore = 0-6% is shown in Fig. 19. It can be seen that the strength decreases as the sample size or porosity increases.
A power function, similar to the one for compression strength in cement paste suggested by Powers [39], is proposed here to relate the specimen strength and the porosity: where r pore 0 is the specimen strength without pores and n is an empirical value. n = 4 is found in this study resulting in good  (a) P agg =20% (b) P agg =30% (c) P agg =40% (d) P agg =50% Fig. 13. Typical samples of different aggregate volume fraction (L = 50 mm, P pore = 2%).  agreement between Eq. (11) and the numerical results, as shown in Fig. 20 for all the sizes. The size effect laws obtained according to Eq. (2) for different porosity are shown in Fig. 21. No definitive trend for the porosity influence on the size effect can be observed. The reason may be that the porosities investigated in this paper (0-6%, typical for standard concrete) are relatively small and therefore the spatial arrangement of the pores has a pronounced impact on the concrete non-linear responses. This means that the simulated number of random distributions might be insufficient to obtain a consistent trend for the effect of porosity on the size effect of peak stress. A consistent trend is expected from simulations with higher porosities, e.g., for porous and foamed concrete, but this is beyond the scope of the present paper.
The size effect on toughness for different porosity is plotted in Fig. 22. For a constant specimen size, it can be observed that a porosity increase results in decreased toughness. Peak Stress (MPa) Aggregate volume fraction (%)

Proposed law
Numerical results Fig. 15. Effect of aggregate volume fraction on peak stress for different sizes (P pore = 2%). (a) P pore =0%

Conclusions
In this paper, we have investigated the size effect on strength and toughness in concrete specimens under tension, through extensive Monte Carlo simulations of mesoscale finite element models containing randomly-distributed aggregates and pores with prescribed volume fractions, shapes and size distributions. The complex meso-crack initiation and propagation is successfully captured by pre-inserted cohesive interface elements. It is found that the numerical results closely follow Weibull's stochastic size effect law, although more research is needed to clarify the contribution of deterministic size effect. Both the size-effect equations for different porosities and aggregate volume fractions, and the relations between them and the strength, are given by curve-fitting the Monte Carlo simulations results. It is also shown that as the aggregate volume fraction, the porosity significantly affect both strength and toughness of the specimens and should not be ignored in size effect studies of concrete. Although the crack patterns and mean stress-strain curves obtained appear realistic, the 2D analysis remains limited to planar crack morphologies, and further verification of these findings via 3D models is required and ongoing.