Abstract
Generation and packing algorithms are developed to create models of mesoscale heterogeneous concrete with randomly distributed elliptical/polygonal aggregates and circular/elliptical voids in two dimensions (2D) or ellipsoidal/polyhedral aggregates and spherical/ellipsoidal voids in three dimensions (3D). The generation process is based on the Monte Carlo simulation method wherein the aggregates and voids are generated from prescribed distributions of their size, shape, and volume fraction. A combined numerical-statistical method is proposed to investigate damage and failure of mesoscale heterogeneous concrete: the geometrical models are first generated and meshed automatically, simulated by using cohesive zone model, and then results are statistically analysed. Zero-thickness cohesive elements with different traction-separation laws within the mortar, within the aggregates, and at the interfaces between these phases are preinserted inside solid element meshes to represent potential cracks. The proposed methodology provides an effective and efficient tool for damage and failure analysis of mesoscale heterogeneous concrete, and a comprehensive study was conducted for both 2D and 3D concrete in this paper.
1. Introduction
Concrete is the most widely used construction material in the world due to its good strength and durability relative to its cost. Numerical simulations coupled with theory and experiment are considered to be an important tool for examining the mechanical behaviour through computational materials science. Wittmann [1] proposed three levels of observation: microlevel (10−6 m), mesolevel (10−3 m), and macrolevel (100 m). The microlevel represents the most basic level, in which the internal structure of cement paste is considered and where physical and chemical forces dominate. In the mesolevel, concrete is usually represented as a material made up of coarse aggregates, mortar with fine aggregates and cement paste embedded inside, and interfaces between coarse aggregates and mortar. At the macroscale, concrete is treated as a homogenous and continuous material in which the internal structure is not considered [2]. As concrete is generally used in large-sized structures and its dependence of mechanical behaviours on mesostructures is significant, mesoscale modelling of concrete becomes essential and was investigated in this study.
An extensive literature review recently carried out by Wang et al. [3] shows that numerical image processing technique and parameterization modelling technique are two main approaches to generate mesostructure models of concrete. Although the mesh generated from images is closest to nature, it is presently very costly and time consuming to generate sufficient scanned images and convert them to mesh for meaningful statistical analyses (e.g., [4–6]). In the parameterization modelling technique, both direct (e.g., [3, 7–9]) and indirect approaches (e.g., [10–13]) to characterize material random heterogeneity have been used. The indirect approach is able to generate a large number of samples with ease but the key multiphasic parameters which could significantly influence the mechanical behaviour cannot be considered. The direct approach can take into account most of the key parameters such as shape, size, gradation and spatial distribution of voids and aggregates, phase volume fractions, and aggregate-mortar interfaces. So among all the approaches, it seems that the direct approach is particularly suitable for statistical analysis of concrete samples and was used in this study.
Identification and generation of unit cell geometry are a vital step in the mesoscale analysis of concrete. Both shape and size of aggregates have significant influences on the stress distribution, crack initiation, and damage accumulation up to the macroscopic failure within the concrete material [16]. Most of the existing algorithms for mesostructure generation of concrete only consider random aggregates but neglect voids [2, 7, 17–19]. However, the X-ray CT images [4, 14, 20] clearly show that voids exist in concrete at this scale and have severely adverse effects on the specimen strength [3, 9]. So automatic generation of morphological details of materials with randomly distributed different shape of aggregates and voids poses new challenges. In this paper, a computationally efficient and versatile in-house program of heterogeneous material generator (HMG) has been developed considering both 2D and 3D random aggregates and voids based on the prescribed parameters.
Several numerical models for crack propagation at mesoscopic level have been used to study the heterogeneity of concrete. Continuum based finite element models are the main approaches employed in the literature [21–24], primarily based on cohesive zone model [25]. Other alternative approaches have also been developed, such as discrete element model (DEM) [26] and lattice model [11, 12]. The biggest disadvantage of DEM and lattice model is the difficulty in determining the equivalent model parameters, which is relatively straight forward for continuum based finite element models [27]. The numerical model used in this work for crack propagation at mesoscopic level is based on the cohesive zone model [3]. It is becoming more and more popular in modelling crack propagation due to its simple formation and easy implementation in the form of cohesive interface elements (CIEs).
With the preprocessing functionality in ANSYS [29] and explicit solver in ABAQUS [30], nonlinear modelling of multiple crack propagation with a considerable number of concrete samples has been performed by high performance computing in this study. Statistically, analysis of mechanical behaviour of both 2D and 3D specimens has been conducted using the combined numerical-statistical method, leading to valuable statistical results that may help improve designs of concrete materials.
2. 2D and 3D Heterogeneous Concrete
2.1. Aggregate and Void Size Distribution
The size distribution of aggregates in concrete is often described by the Fuller curve [16], which is discretised into a certain number of segments. The aggregate size distributions found in Hirsch [28] and summarised in Table 1 are used in this study. In this study, concrete is treated as a material consisting of coarse aggregates, mortar incorporating fine aggregates dissolved in it, voids, and interfaces between aggregates and mortar. Here only coarse aggregates are explicitly modelled as mesoscale features. For normal strength concrete cast with a mould in the laboratory, coarse aggregates usually represent 30–50% of the concrete volume.
From X-ray CT images, it was found that gravel aggregates have circular/elliptical shapes in 2D and spherical/ellipsoidal shapes in 3D, while crushed aggregates have polygonal shapes in 2D and polyhedral shapes in 3D [14]. A number of voids can be clearly observed in concrete at mesoscale from CT images (see Figure 1). Therefore, voids with size ranging from 2 to 4 mm as presented in [31] are included in the mesostructure of concrete in this study. Coarse aggregates are considered to have elliptical/polygonal shapes in 2D and ellipsoidal/polyhedral shapes in 3D, and the shape of voids is assumed to be circle/ellipse in 2D or sphere/ellipsoid in 3D.
2.2. Specimen Generation
The basic idea is to create the aggregates and voids in the concrete in a repeated manner, until the target area/volume is filled with aggregates and voids. The generation starts with the information input process and is followed by taking and placing voids within size range; the aggregates within the grading segments are produced last. The “input” process reads the controlling parameters for generating a heterogeneous concrete; the “taking” process generates an individual void or aggregate in accordance with the random size and shape descriptions. The “placing” process positions the aggregates and voids into the predefined area/volume in a random manner, subject to the prescribed physical constraints. There are three physical conditions to be satisfied simultaneously: (1) the whole inclusion (void or aggregate) must be within the concrete, (2) there is no overlapping of inclusions, and (3) there is no intersection between any two inclusions. A three-level hierarchical method is proposed to reduce the computational cost and is outlined in Figure 2. The procedure was programmed using MATLAB [32] in this study, and the detailed algorithms and procedure in generating mesostructures can be found in our previous work [3, 33].
2.3. Mesostructure Models
Elliptical aggregates and circular/elliptical voids in 2D or ellipsoidal aggregates and spherical/ellipsoidal voids in 3D were used for concrete with gravel aggregates. A series of concrete samples with dimensions of 50 mm × 50 mm in 2D and 50 mm × 50 mm × 50 mm in 3D are generated. Figure 3 shows the numerically generated models with gravel aggregates (aggregates in green and voids in red). Aggregate gradation in Table 1, aggregate volume fraction , void content , aspect ratios for elliptical/ellipsoidal aggregates and voids (the ratio of the major axis to the minor axis) , minimum space between the edge of an aggregate and the boundary of the concrete specimen mm, and minimum gap width between any two aggregates mm are adopted here [31].
(a) Elliptical aggregates and circular voids
(b) Elliptical aggregates and elliptical voids
(c) Ellipsoidal aggregates and spherical voids
(d) Ellipsoidal aggregates and ellipsoidal voids
Polygonal aggregates and circular/elliptical voids in 2D or polyhedral aggregates and spherical/ellipsoidal voids in 3D are considered for concrete with crushed aggregates. A series of concrete samples with the same dimensions as the previous gravel aggregate samples are generated. Figure 4 shows the numerically generated models with crushed aggregates. The minimum and maximum edges of the polygons/polyhedrons are set to , . All the other parameters are the same as the values used for gravel aggregates.
(a) Polygonal aggregates and circular voids
(b) Polygonal aggregates and elliptical voids
(c) Polyhedral aggregates and spherical voids
(d) Polyhedral aggregates and ellipsoidal voids
3. Combined Numerical-Statistical Method
3.1. Description of the Method
A combined numerical-statistical method is proposed in this paper to study the material behaviour of concrete in a statistical sense. The detailed procedure is as follows:(1)Generate a model with prescribed variables, for example, sample size, aggregate volume fraction, void content, and aggregate shape.(2)Perform a finite element simulation of the sample for given boundary conditions.(3)Compute the mean value, standard deviation, and coefficient of variation (CoV) of effective property for the considered sample size.(4)Repeat steps (1) to (3) for sufficient number of random samples to meet the required precision, and conduct statistical analysis.
This procedure is automated by running a batch file in this study.
Results from all realisations are evaluated statistically. The standard deviation within a series of samples is given bywhere is the series average result and is the result from sample .
To compare results obtained with different sample sizes quantitatively, we use the coefficient of variation given byThis expresses the fluctuation of measured property relative to its average value.
3.2. Cohesive Zone Model in ABAQUS
The cohesive zone model developed by Barenblatt [34, 35] and Dugdale [36] enables the simulation of the energy dissipation process in the fracture process zone (FPZ) during fracture, as illustrated in Figure 5. A bilinear cohesive crack model is used here to predict discrete crack propagation in concrete specimens under tension loading. A linear ascending branch is added in the softening curve to model the initially uncracked material.
The separation displacement is difficult to derive from experiments, so cohesive fracture energy and cohesive strength are usually used. Among them, only two parameters are independent, and the fracture energy can be obtained aswhere is the cohesive fracture energy, is the cohesive strength, and is the separation displacement.
In the 3D cohesive zone model, it is assumed that there exist a normal traction and two tangential tractions (shear cohesion) and across the crack surfaces, through mechanisms such as material bonding, aggregate interlocking, and surface friction in the fracture process zone. Figures 6(a) and 6(b) illustrate typical linear softening curves for and , where and are the critical relative displacements when the tractions diminish for normal traction and tangential traction components, respectively. The unloading paths are also indicated.
(a) curve in normal direction
(b) curve in tangential direction
It is worth noting that the initial tensile stiffness should be high enough to ensure displacement continuity across the cohesive interfaces before the tensile strength is reached, but not too high to cause numerical instability due to ill-conditioned global stiffness matrix. Reasonable initial shear stiffness, and , is also needed before the shear strength, or , is reached. The effects of initial stiffness on computational results have been investigated previously [37, 38]. The following relation is suggested in [38] as a guideline for initial stiffness selection:where and are Young’s modulus and Poisson’s ratio, is the characteristic size of elements, and is taken as 10~100 from the experience in [38].
The areas under the curves in Figures 6(a) and 6(b) (calculated by (3)) represent the normal facture energy and twice the tangential fracture energy, and , respectively:
Cohesive elements in ABAQUS based on the cohesive zone model are used here. The damage is characterised by a scalar index representing the overall damage of the crack caused by all physical mechanisms. The effective relative displacements combining the effects of , , and can be obtained aswhere is the Macaulay bracket and
The damage evolution law is given bywhere is the maximum effective relative displacement attained during the loading history. and are effective relative displacements at damage initiation and complete failure, respectively. It is obvious that the damage variable evolves monotonically from 0 to 1 upon further loading after the initiation of damage.
The damage initiation and evolution will degrade the unloading and reloading stiffness coefficients , , and , which can be calculated asThe tractions are also affected by the damage according towhere and are the traction components predicted by the elastic traction-displacement behaviour for the current separation without damage.
In this study, damage in the cohesive zone model is assumed to initiate when a quadratic interaction function involving the nominal stress ratios reaches a value of one
The 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, three different interface cohesive elements. The material heterogeneity is investigated by considering different phases in the concrete specimen with corresponding material properties.
In the fracture process zone for a 2D case, tractions exist in the normal direction and shear direction across the crack interface, and the corresponding relative displacements are and . So the 2D constitutive law could be easily deduced from 3D characterization described above by taking one of the shear tractions and displacements out.
3.3. Cohesive Interface Elements Insertion
In this study, all finite element meshing is performed with the preprocessing functionality in commercial FE package ANSYS [29]. To mesh the mesoscopic structure of concrete, different material parts (mortar and aggregates) should maintain continuity on their surfaces. Hence, the finite element boundaries are coincident with material surfaces and there are no material discontinuities within the elements. The mortar and aggregates are meshed with triangular elements (plain stress) in 2D and hexahedral elements in 3D so that more realistic crack paths can be obtained. The specially developed ANSYS parametric design language (APDL) programs in combination with the ANSYS batch processing provide a powerful tool of automatic mesh generation for a large number of samples required for statistical analysis. Following the meshing procedure, the generated mesh will automatically have shared nodes at the interfaces between two elements. If the interfaces surrounding the elements are to be modelled explicitly as potential microcrack sources, a duplicate set of nodes will be required at the interface locations. Here 4-node or 6-node zero in-plane thickness CIEs are preinserted into the existing element interfaces in 2D or 3D. This is conducted by a purpose written and flexible in-house computer program cohesive interface elements insertion (CIEIN). Three sets of CIEs with different traction-separation softening laws are inserted (see Figure 7), namely, CIE-AGG inside the aggregates (grey in Figures 7(b) and 7(d)), CIE-MOR inside the mortar (green in Figures 7(b) and 7(d)), and CIE-INT on the aggregate-mortar interfaces (yellow in Figures 7(b) and 7(d)). The element and node numbers are denoted by and , respectively. The detailed numbering of elements and nodes in the initial mesh and the mesh with inserted cohesive elements shows the insertion procedure with the new nodes generated at the same positions and CIEs between the continuum elements.
(a) 2D initial mesh
(b) 2D mesh with zero-thickness CIEs
(c) 3D initial mesh
(d) 3D mesh with zero-thickness CIEs
4. Numerical Simulation and Results
4.1. Description of the Numerical Model
2D concrete specimens with elliptical aggregates and circular voids and 3D specimens with ellipsoidal aggregates and spherical voids (, , , 0.2 mm) were modelled, loaded under uniaxial tension. 25 mm × 25 mm concrete square in 2D and 25 mm × 25 mm × 25 mm concrete cubic in 3D are used here. All the other mesoscale parameters were fixed to be the same for 2D and 3D models. Figure 8 shows the geometry, boundary, and loading conditions of the numerical model. Horizontal displacements were prescribed to all nodes on the left and right surfaces of the specimen, with a value equal to zero for the left surface, and a uniformly distributed displacement on the right surface. Vertical displacements for the same nodes are left free, except for the nodes at the left lower corner of the specimen, which are fixed to prevent rigid body translation. A displacement controlled loading scheme was used and all the analyses were ended at a displacement mm and loading time s [3]. Considering the balance between accuracy and efficiency, the element length 1 mm was used for all the meshes in this study [3].
(a) 2D
(b) 3D
Generally, aggregates have much higher strength than mortar and interfaces in normal concrete. In this section, no cracks are allowed to initiate inside the aggregates by assuming elastic behaviour without damage in CIE-AGG. The linear tension/shear softening laws described above were used to model CIEs with quadratic nominal stress initiation criterion and linear damage evolution criterion. Similar material properties extracted from [39] were used in this study, which are listed in Table 2. The fracture properties related to the shear behaviour were assumed to be the same as the normal ones due to the lack of experimental data. Both initial tensile and shear stiffness for all the cohesive elements were set to 106 MPa/mm after trial and error.
ABAQUS/explicit with small time increments (typically s) was utilised to solve the highly nonlinear equation systems by the high performance computing (HPC) cluster at the University of Manchester computational shared facility (CSF). A typical run of simulations with 50 samples in this study takes about 6 hours by parallel computation with 32 CPUs.
4.2. Tensile Behaviour
Due to the statistical nature of mesostructure models, an extensive series of numerical simulations would be necessary to capture the range of behaviours. Fifty random samples were modelled in this study to ensure that the results are statistically converged.
In Figure 9, resulting mean stress-displacement and toughness-displacement curves, together with the variation range for both 2D and 3D, are depicted. It can be seen that the predicted stress-displacement curves are qualitatively similar with a clear peak and sharp initial postpeak drop followed by a long shallow tail. Numerical mean stress for every curve is obtained by summing the nodal reaction forces in the constraints and dividing by the specimen cross section area. The mean peak stress predicted by 2D and 3D modelling is 2.65 MPa and 3.49 MPa, respectively, representing an increase of 24.1%. The mean toughness (energy absorbed per unit volume up to the break point, i.e., the area under the stress-strain curve up to the break point, measured in J/m3) predicted in 2D and 3D modelling is 2.47 KJ/m3 and 3.15 KJ/m3, respectively, representing an increase of 21.6%; Figure 10 shows the standard deviation-displacement curves for both 2D and 3D modelling. The standard deviation of peak stress is 0.11 MPa and 0.05 MPa, respectively, representing a decrease of 54.5%. These increases in peak stress or decreases in standard deviation are due to the constraint effects in the thickness direction in 3D and the less smooth cracking surfaces from 3D modelling. However, the standard deviation of 3D modelling at softening range is larger than that in 2D modelling which indicates that 3D modelling is more complicated and dependent on different samples. This is due to the 3D heterogonous distribution of aggregates and voids. It is clear that the third dimension cannot be neglected due to the nature of fracture process. The importance of the thickness effect has also been pointed out experimentally by Van Mier and Schlangen [40].
(a) Stress-displacement curves
(b) Toughness-displacement curves
The mean curve, mean value, and standard deviation of the stress and toughness shown in Figures 9 and 10 are extracted from the stress-displacement curves for 50 samples with different aggregate and void distribution. In order to evaluate the effect of sample number on simulation results as required for the proposed combined numerical-statistical method, two special loading points, namely, point A at peak stress and point B at displacement 0.03mm where largest standard deviation is observed (denoted in Figures 9(a) and 10), are investigated. The influence of sample number on the coefficient of variation (CoV) of stress at points A and B is illustrated in Figure 11. It can be seen that 50 samples are enough to get convergent values of CoV of stress with a stable fluctuation. This is an indication that the combined numerical-statistical method with 50 random realisations offers a good balance between general applicability and testing of different samples. It can be observed that COV of peak stress for different void and aggregate distribution is low (4.2% for 2D modelling and 1.5% for 3D modelling), which demonstrates that the peak stress is relatively insensitive to the void and aggregate distribution. The high CoV of stress at point B further demonstrates the necessity to conduct a statistical analysis as the results vary greatly for different aggregate and void distribution in softening range.
Figure 12 shows the statistical distribution of peak stress, together with the probability density and the best fit Gaussian Probability Density Function (PDF). The probability density curve can be used to calculate structural reliability or failure reliability against given external loadings and material properties for structural design.
(a) 2D
(b) 3D
4.3. Crack Patterns
The complex mesoscale crack propagation is realistically simulated using the proposed method, and typical crack patterns for both 2D and 3D at failure are shown in Figure 13. To clearly visualise the fracture surfaces, models with damaged cohesive elements, models without cohesive elements, and failed interfaces only are used for 2D models (see Figures 13(a)–13(c)), and models with damaged cohesive elements, models without cohesive elements, and morphology of failed surface are used for 3D models (see Figures 13(d)–13(f)). The cracks shown in Figures 13(a), 13(c), and 13(d) are represented by red CIEs with damage index ( means complete failure). A magnification factor of 20 and 200 was used for the deformed models with damaged CIEs and those without damaged CIEs, respectively. The noticeable difference of 2D and 3D modelling discussed above could also be summarized in that the 3D modeling predicts more realistic fracture surfaces in the thickness direction whereas 2D modeling only assumes plane fracture problems. This result confirms the importance of including the third dimension into the analysis.
(a) 2D model with damaged CIEs
(b) 2D models without CIEs
(c) 2D models with failed interfaces
(d) 3D model with damaged CIEs
(e) 3D model without CIEs
(f) 3D morphology of failed surface
5. Conclusions
Models of numerical concrete with random mesostructures comprising elliptical/polygonal aggregates and circular/elliptical voids in 2D and ellipsoidal/polyhedral aggregates and spherical/ellipsoidal voids in 3D have been generated in this study. Numerical-statistical analysis was carried out using a cohesive zone model to simulate damage and failure of concrete. Crack patterns are realistically simulated using the technique of preembedding cohesive interface elements. The main conclusions, based on the results obtained from the statistical analysis under a uniaxial tension loading, are as follows: (1) statistical analysis is necessary in both 2D and 3D due to the high dependence of material behaviour on different aggregate and void spatial distribution; (2) the third dimension is demonstrated to have a pronounced influence on both macroscopic mechanical properties and crack patterns in tension; (3) compared to 2D modelling, 3D modelling demonstrates a larger mean peak stress and a smaller standard deviation in the prepeak response, attributed to more uniformly distributed microcracks within the specimen; a larger standard deviation/CoV of stress in the postpeak response is attributed to the larger number of possibilities for microcrack coalescence under the constraint of randomly distributed features. It has to be pointed out that the conclusions obtained are based on mesoscale modelling with specific specimen and aggregate size. The variation on resulting mechanical behaviour is also associated with the size of specimen and aggregates, and this phenomenon may not exist when specimen is large enough, for example, at the length scale of engineering structures.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is supported by the EPS Faculty Ph.D. Studentship from the University of Manchester for Xiaofeng Wang. The authors would like to acknowledge the assistance given by IT Services and the use of the Computational Shared Facility (CSF) at The University of Manchester.