Mesoscopic Finite Element Modeling of Concrete Considering Geometric Boundaries of Actual Aggregates
Concrete is nonhomogeneous and comprises aggregate, mortar, and interfacial transition zones at the mesoscopic scale. The aggregate shapes significantly affect the development of microcracks. To deal with the problem of imprecise description of actual aggregate, an innovative method of modeling concrete is proposed in this study considering geometric boundaries of actual aggregate. First, the geometric feature points of the actual gravel aggregates, that is, the shape of the actual aggregate, are obtained by laser scanning. The geometric feature points are then moved randomly in the plane. Using this method, an aggregate library is established based on the actual aggregates. Finally, the front polygons-rear circumcircle conflict and overlap criteria are proposed, which can achieve a rapid placing process of the multicontrol point aggregate. Using this method, numerical uniaxial tensile and three-point bending beam tests are conducted and the results are compared with the round aggregate model. The results indicate that the geometric properties of aggregates have both blocking and guiding effects on crack development. Therefore, the proposed modeling method is better suited for analyzing crack development.
Traditional macroscopic models assume that concrete is isotropic based on continuum mechanics. However, the models cannot explain the crack propagation routes under external loads. Furthermore, it is also difficult to describe the concrete damage caused by stress concentrations . Therefore, the macroscopic model is not able to determine the relationship between the internal structure of concrete and the macroscopic mechanical properties.
With the advances in computing power, different modeling methods of concrete in the mesoscopic scale were proposed [2–9]. Among these methods, the random aggregate model assumes that concrete comprises aggregate, mortar, and interfacial transition zones [10–12]. It highlights the random distribution of the aggregate and gives an explicit characterization of the aggregate shape .
Pebble aggregate was typically simulated by the circle and the ellipse in two dimensions [14–17] or the sphere and the oval in three dimensions [18, 19]. For the gravel aggregate, the scholars typically reduced it to polygons or polyhedrons. For instance, Gao and Liu  used triangles and quadrilaterals as an aggregate base, and new aggregates were generated by adding new vertices and controlling the side lengths. Sun et al.  used random circles as the aggregate size to generate triangular aggregate bases and then grouped all the bases simultaneously and extended them randomly, thereby generating high aggregate content models. Ma et al.  used the Wallavin formula to determine the distribution of two-dimensional (2D) concrete specimens and then extended the internal polygon base until its area was the same as the circle aggregate. In the abovementioned studies, all the aggregate models were based on regular geometry, and specific mathematical algorithms were applied to structure the shape. However, this does not reflect the characteristics of actual aggregates.
Simulations based on actual aggregate have become increasingly popular . The primary method is to build an aggregate model using computed tomography (CT) slice images. These can better characterize the actual aggregate geometry. Ren et al.  used X-ray CT images to obtain the shape and location information of aggregates, and from this, they proposed the mesoscopic finite element fracture model. Qin and Du  used the bottom-hat transform methods to convert CT slices into three-dimensional (3D) images. He reconstructed concrete 3D models based on the volume data method; the model took porosity into consideration. Pan et al.  analyzed the size distribution of 2D aggregate based on the 3D sections and proposed a numerical method. Applying this method, the 2D mass-distribution function with arbitrary gradation can be determined. Fu et al.  used the polygons to approximately describe clear outlines on the CT slices. The aggregate sizes were normalized, and the double criterion was used to judge the aggregate intrusion. Li and Wang  established a relatively complete system for concrete mesoscopic mechanics analysis to simulate the process of early-age shrinkage cracking in high-performance concrete based on CT image reconstruction. Although the CT reconstruction technology can characterize the shape of the actual aggregate satisfactorily, the aggregate shape and gradation are limited by the specific slice. A significant number of experiments are needed to satisfy varying gradations; therefore, this method has a low simulation efficiency.
In order to improve the aggregate modeling, we proposed an innovative method. This study uses 2D modeling as an example and is organized as follows. Section 2 introduces the method of aggregate generation, including scanning of actual aggregate and the establishment of an aggregate library. In Section 3, the aggregate placement basic principle and overlap criteria are discussed. Two numerical mechanical tests are conducted and discussed in Section 4, and the advantages of the proposed model are analyzed. Finally, the conclusions are presented in Section 5.
2. Aggregate Generation
2.1. Scanning of Actual Aggregate
For actual aggregate, the outline shape of the single aggregate is complex and concave-convex. A Creaform handheld laser scanner was used to scan aggregate particles for the 3D shapes. The reconstructed aggregate images were obtained by scanning the actual aggregate, as shown in Figure 1.
The numerous contour points of the reconstructed aggregate can describe the actual aggregate accurately, as shown in Figure 2(a). However, this makes the meshing and calculation in finite element models highly complex. Therefore, it is necessary to deal with the contour points of the actual aggregate.
To simplify the contour points, the key control points should be extracted from the original contour points. These points can describe the basic shape of the actual aggregate using as few as possible contour points. The aggregate shape is simplified as shown in Figure 2(b) and is known as the aggregate base of the actual aggregate.
In this study, the control of the curvature and the minimum distance between two points are used to reduce the number of contour points, and the simplified aggregate shape can be generated by geometric control points and straight lines. Before determining the minimum distance between two points, the average particle size of each aggregate is calculated as follows:where is the average particle size of the aggregate, is the distance from the contour control point to the geometric center, and is the number of contour control points.
Assuming any point on the contour as the starting point, the distance between two adjacent control points satisfieswhere is the distance control coefficient.
The control conditions for curvature control using the three-point method are as follows: find the circumscribed circle of the triangle formed by three adjacent points, and the curvature radius of the circumscribed circle can be obtained. If the curvature radius is greater than a given value, the middle point is removed. The radius of curvature for any three consecutive points is calculated as follows:(1)Calculate the length of each side of the triangle according to the coordinates of three points: , , and .(2)Calculate the angle of the triangle vertex opposite to the circumscribed circle, , according to the cosine theorem.(3)Use the sine theorem to find the radius of curvature .(4)The curvature control radius can be taken as 0.15–0.2 times the average particle size.
2.2. Establishment of Aggregate Library
Because of the statistical similarities among different aggregates, the key control points of the actual aggregates can be used as an aggregate base to generate new aggregates. The method of generating new aggregates randomly is as follows:(1)Two circles are drawn with radii and , and the centroid of the aggregate is taken as the center of the circle. As shown in Figure 3(a), the control points move between the annular band comprising the two circles, and and can be determined according to the shape of the new aggregate.(2)Select a control point in all control points randomly, and allow it to move randomly along the double arrow line in Figure 3(a). The distance moved is controlled as , and the number of random moves “” is selected according to the deviation from the base shape of the original aggregate.
Figures 3(b)–3(d) show an example of the new aggregates generated after random movement, and the contours of the aggregates had been determined. However, the particle size and the location of the aggregate should be taken into consideration when it is placed into the designated area. Therefore, we need ways to determine the local coordinates of the aggregate particles and to adjust the particle size in accordance with the particle size distribution to treat the aggregate generated above.
The local coordinates of the aggregate particles are determined by subtracting the geometric center coordinates from the outline coordinates of the aggregate so that the geometric center point is at the origin. Therefore, this point is easily located when the aggregate is placed.
Adjusting the particle size in accordance with the particle size distribution allows the determination of the circumscribed circle diameter of each aggregate in the aggregate library, with aggregate coordinates as follows:
The aggregate sizes have been treated after the application of (3), and when placing aggregates, the aggregate coordinates are enlarged and located based on the gradation determined by the Walraven formula.
After scanning the different shapes of actual aggregate, the number of aggregate bases is determined. On this basis, the above methods are used to generate new aggregates and a statistically significant actual aggregate library is formed.
3. Aggregate Placement
3.1. Basic Principle
Aggregate placement requires that aggregates do not overlap and achieve a certain aggregate content. Previously, the bulk of random aggregate models were simple convex polygons. In the process of placing, scholars established the criterion of aggregate intrusion in order to avoid the overlapping of aggregates, including the arc discrimination method , the area discrimination method , the angle sum test method , the state matrix method , and the aggregate determination method based on a background grid . In actual aggregates, the contour shape is more complex and has a greater number of control points than the simple convex polygons. In addition, earlier discriminant methods for the placement were limited by computational inefficiencies. Therefore, it is necessary to establish an aggregate placement algorithm that adapts to the actual aggregate multicontrol points.
In this study, we propose an aggregate placement algorithm based on the intrusion judgment of front polygons-rear circumcircle. The algorithm is proposed as follows.
Determine the circumcircles of each aggregate in the actual aggregate library, and then place the aggregates from large to small according to gradation. When placing a new aggregate, determine if the new aggregate circumcircle and the placed aggregate contours overlap: the new aggregate is correctly located when they do not overlap. The algorithm effectively solves the problem of low discriminating efficiency between polygons with significant numbers of control points while simultaneously satisfying the gradation variation of the aggregates without creating large gaps between the aggregates. This method is a novel way to place actual aggregates with complex geometric boundaries.
The Walraven formula is used to calculate the area occupied by each gradation. The Walraven formula  is as follows:
An aggregate is randomly selected from the aggregate library, and the circumcenter circle is determined. In order to ensure the randomness of aggregates, the aggregates were subjected to aggregate rotation and particle size gradation before placing.
Aggregate rotation is performed by randomly rotating a selected aggregate around its circumscribed circle center:where are the coordinates of the control point before rotation, are the coordinates of the control point after rotation, and are the random numbers from 0 to 1.
Particle size grading is performed by randomly expanding aggregates within a range according to the graded distribution range of the aggregates, generating aggregates meeting the specified gradation requirements:where are the coordinates of the control point after expansion; are the maximum and minimum particle sizes in the grading range, respectively; and is a random number from 0 to 1.
The aggregates are placed from large to small until the cumulative area satisfies the gradation requirements. Once this gradation placement is completed, the next gradation placement will be done.
3.2. Aggregate Overlap Criteria
An aggregate placement process is complete when the circumcircle of the aggregate is placed into the designated area. Therefore, it is necessary to judge whether the circumcircles of the existing aggregates and the new aggregate overlap when placing, and this can be regarded as the intrusion relationship between polygons and circles. The method simplifies the judgment required between traditional polygons. Simultaneously, the principle of placing aggregates from large to small circumvents the problem of low aggregate content because of gaps between the circumscribed circles and the polygons.
When judging the location of the relationship between circles and polygons, it can be divided into three cases: (1) control point inside circle; (2) polygonal edges intersect with circles; and (3) circle inside the polygon, as shown in Figure 4. For the first case, it can be immediately determined if the aggregate control point is inside the circumcircle; for the second case, it is required to judge whether the distance from the circumcircle center to the straight line is smaller than the circumcircle radius and whether the vertical point falls on the line segment formed by the adjacent control point; for the third case, calculate each angle between the circumcircle center and two adjacent vertices, and then determine if the sum of these angles is equal to 360°.
3.3. Aggregate Generation by Gradation
The coarse aggregates of concrete are divided into four gradations according to the particle size: small stone (5–20 mm), stone (20–40 mm), large stone (40–80 mm), and extra-large stone (80–150 mm) . Using the method proposed in this study to generate two, three, and four-gradation concrete specimens, as shown in Figure 5, the size of the two-gradation aggregate model is 150 × 150 mm, and the aggregate content is 45.42%. The size of the three-gradation aggregate model is 300 × 300 mm, and the aggregate content is 52.17%. The size of the four-gradation aggregate model is 450 × 450 mm, and the aggregate content is 54.00%.
4. Numerical Test and Discussion
4.1. Concrete Uniaxial Tension Test
The 2D aggregate model established in this study is based on the actual aggregate boundary features and has a high degree of similarity with the actual aggregate. In order to study the effect of aggregate shape on the macroscopic mechanical properties and failure characteristics of concrete, numerical simulations of the uniaxial tensile test and three-point bending test of circle aggregate model and actual aggregate model were conducted. The aim was to determine the aggregate shape influence on macroscopic mechanical behavior.
According to the above algorithms for generation and placement, an actual aggregate model and a circular aggregate model are generated; the geometric model is shown in Figure 6. The specimen size is 100 × 100 mm and the aggregate size is 2–20 mm. The aggregate contents of the actual aggregate model and the circular aggregate model are 49.09% and 49.11%, respectively, and geometric shrinkage reduces the interface thickness by 0.5 mm to form an interface between aggregate and mortar.
The geometric model of the aggregate is imported into the ABAQUS finite element software. The mesh is divided by a quadrilateral-based freeform grid. Taking the actual aggregate model as an example, the model grid is divided into a grid as shown in Figure 7, and the total number of elements is 50,717. The total number of nodes is 50,149.
The ABAQUS damage plasticity model is used to describe the fracture behavior of concrete materials. This model can be used to characterize the inelastic behavior of concrete combined with isotropic elastic damage and isotropic tensile and compressive plasticity theory. In order to simplify the calculation, the evolution model of the double-fold damage variable  is adopted as the constitutive of each phase material, as shown in Figure 8.
The damage variables can be expressed aswhere is the tensile strength, is the residual tensile strength, is the principal tensile strain corresponding to , is the residual strain corresponding to the tensile residual strength , is the ultimate tensile strain, and is the largest primary tensile strain value of the loading history. The relationship between parameters is , , and .
The concrete phase material mechanical parameters are presented in Table 1, where was taken as 0.1, was taken as 10, and the aggregate, mortar, and interface were taken 5, 4, and 4, respectively. The concrete uniaxial tensile specimen calculation diagram is shown in Figure 9, where the left boundary is the applied horizontal constraint, the left middle is the applied vertical constraint, and the right boundary is the applied displacement load.
The macroscopic stress-strain curve of the concrete tensile specimen is shown in Figure 10. In the elastic phase, the gradients of the stress-strain curve of the circular aggregate model and the actual aggregate model are essentially identical, which indicates that the stress characteristics of the macrospecimens in the elastic phase are the same. At the end of the elastic phase, the strain of the specimen approaches , numerous cracks are initiated at the interface in the concrete specimens, and microcracks are formed.
After entering the softening stage, the damage area expands as the strain increases. In order to observe the development of macroscopic cracks, the strain of the test specimens was increased to greater than , where macroscopic cracks, greater than 0.02 mm, became visible. Apart from the small cracks that developed on the left side of the circular aggregate model, macroscopic cracks continued to develop along the aggregate interface and the mortar in the middle left of the specimen and finally formed a macroscopic crack perpendicular to the tensile direction of the specimen, as shown in Figure 11(a). For the actual aggregate model, the crack developed on the right of the specimen perpendicular to the direction of stretching. Unlike for the circular aggregate, the macroscopic cracks in the actual aggregate model developed simultaneously at the two interfaces, a phenomenon known as crack bifurcation, and finally met at the bottom, as shown in Figure 11(b). From the stress-strain curve, as the aggregate shape and aggregate distribution of the circular aggregate model and the actual aggregate model are different, the stress peak and the corresponding strain of the actual aggregate model are marginally greater than that of the circular aggregate, which led to single and bifurcation cracks. This phenomenon reflects the significant influence of aggregate geometry on the macroscopic mechanical properties of the specimens.
In the descending section of the stress-strain curve, the aggregate cracks continue to develop. The stress-strain curves of the actual aggregate model and the circular aggregate model decrease at a similar rate, and the softened end of the curves are coincident. At this point, the crack in the specimen has approximately penetrated the specimen, indicating that the geometric characteristics of the aggregate will result in different curves in the softening section. The ultimate stresses of the specimen were essentially identical.
4.2. Three-Point Bending Concrete Beam Test
The calculation diagram of the concrete three-point bending beam specimen is shown in Figure 12. The material parameters of each phase are the same as before, and the specimen size is taken as 400 × 100 mm. To reduce the difficulty of meshing and calculation, the middle of the specimen is taken as a nonuniform area and both sides of the specimen are taken as a uniform area. The model mesh is shown in Figure 13, with a total number of 62,988 units and 62,707 nodes.
The reaction force and the deflection of the mid-bottom point of the beam are extracted from the simulation results. The load-deflection curve is obtained as shown in Figure 14, where it can be seen that the gradient of the circular aggregate model in the elastic phase is slightly greater than that of the actual aggregate model, and the circular aggregate has a greater peak stress than the actual aggregate. When compared with the uniaxial tensile stress-strain curve, the softening section of the load-deflection curve is relatively flat while the curve of the circular aggregate model decreases faster than the actual aggregate model. However, the eventual trends of the two models are similar.
From the final crack distribution and crack developing process, as shown in Figure 15, it can be seen that under the influence of the loading and the boundary conditions, the macroscopic crack initiation zone is formed in the right-bottom area and points to the loading point. The crack of the circular aggregate model begins developing in small aggregates and mortar and reaches aggregate 1 as shown in Figure 15(a), where the development of cracks is blocked and continues to develop upward along the boundary of aggregate 1. The same phenomenon, known as the aggregate blocking effect, is seen for aggregates 2, 3, and 4. For the actual aggregate model, the cracks also occur in the interface and mortar, and the crack development changes the direction along the boundary of the aggregate depending on the shape of the aggregate, as can be seen for aggregates 1, 2, 3, and 4 in Figure 15(b). This was especially noticeable in aggregate 1, where the crack clearly changed the direction, indicating that the shape of the aggregate directs the crack development.
Therefore, it can be concluded that the aggregate geometry has the effect of blocking and guiding the crack in the mesoscale model, while the blocking effect tends to occur in the circle aggregate model and the guiding effect in the actual aggregate model. However, the blocking can also be observed in the actual aggregate model, as shown in Figure 15(b), and for aggregate 5, the effect of blocking cracking is more obvious when the aggregate length direction is perpendicular to the fracture direction. It can be seen that the softening end of the actual aggregate model load-deflection curve gradually flattens.
A model of concrete micromechanics based on the geometric boundaries of actual aggregate was proposed. As opposed to the polygon model, it uses the boundaries of actual aggregate as a base to generate new aggregates. In addition, the model could place the aggregates according to their gradation. To solve the problem of low discriminating efficiency of polygons with a great number of vertices, the front polygons-rear circumcircle conflict and overlap criteria were proposed.
Based on these proposed solutions, numerical tests of uniaxial tension and three-point moment beams were conducted. The results indicated that the geometric properties of aggregates have both blocking and guiding effects on crack development. As opposed to traditional models, the simulation of crack development in the model proposed in this study approximates actual crack development.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
X. Du and L. Jin, “A review on meso-mechanical method for studying the static-mechanical properties of concrete,” Advances in Mechanics, vol. 41, no. 4, pp. 411–426, 2011.View at: Google Scholar
A. Cundall and R. D. Hart, “Numerical modeling of discontinua,” in Proceedings of the 1st US Conference on Discrete Element Methods, Golden, CO, USA, 1989.View at: Google Scholar
A. R. Mohamed and W. Hansen, “Micromechanical modeling of concrete response under static loading: part I: model development and validation,” ACI Materials Journal, vol. 96, no. 2, pp. 196–203, 1999.View at: Google Scholar
A. R. Mohamed and W. Hansen, “Micromechanical modeling of concrete response under static loading: part II: model prediction for shear and compressive loading,” ACI Materials Journal, vol. 96, no. 3, pp. 354–358, 1999.View at: Google Scholar
C. Tang and W. Zhu, Concrete Damage and Fracture–Numerical Test, Science Press, Beijing, China, 2003.
G. Liu and Z. Wang, “Numerical simulation study of fracture of concrete materials using random aggregate model,” Journal of Tsinghua University, vol. 36, no. 1, pp. 84–89, 1996.View at: Google Scholar
C. Zhang, X. Tang et al., “State-of-the-art literature review on concrete meso-scale mechanics,” Journal of Hydroelectric Engineering, vol. 34, no. 12, pp. 1–18, 2015.View at: Google Scholar
Z. P. Bažant, M. R. Tabbara, M. T. Kazemi et al., “Random particle model for fracture of aggregate or fiber composites,” Journal of Engineering Mechanics, vol. 116, no. 8, 1990.View at: Google Scholar
X. Ruan and Z. Pan, “Mesoscopic simulation method of concrete carbonation process,” Structure and Infrastructure Engineering, vol. 8, no. 2, pp. 99–110, 2013.View at: Google Scholar
Z. Gao and G. Liu, “Two-dimensional random aggregate structure for concrete,” Tsinghua Science and Technology, vol. 43, no. 5, pp. 710–714, 2003.View at: Google Scholar
L. Sun, C. Du, and C. Dai, “Numerical simulation of random aggregate model for mass concrete,” Journal of Hehai University, vol. 33, no. 3, Article ID 291295, 2005.View at: Google Scholar
H. Ma, S. Mi, and H. Chen, “A generating approach of random convex polygon aggregate model,” Journal of China Institute of Water Resources and Hydropower Research, vol. 4, no. 3, pp. 196–201, 2006.View at: Google Scholar
A. Chen, Z. Pan, R. Ma et al., “New development of mesoscopic research on durability performance of structural concrete in bridges,” China Journal of Highway and Transport, vol. 29, no. 11, pp. 42–48, 2016.View at: Google Scholar
W. Ren, Z. Yang, and Y. Huang, “Meso-scale fracture modelling of concrete based on X-ray computed tomography images,” Journal of Hydraulic Engineering, vol. 46, no. 4, pp. 452–459, 2015.View at: Google Scholar
W. Qin and C. Du, “Meso-level model of three-dimensional concrete based on the CT slices,” Engineering Mechanics, vol. 7, pp. 186–193, 2012.View at: Google Scholar
Z. Pan, X. Ruan, and A. Chen, “Simulation method of random aggregate in two dimension based on arbitrary gradation,” Journal of Tongji University, vol. 41, no. 5, pp. 759–764, 2013.View at: Google Scholar
B. Fu, J. Li, and G. Lin, “Mesoscopic numerical model of concrete based on data base of real aggregate shapes,” Journal of Architecture and Civil Engineering, vol. 27, no. 2, pp. 10–17, 2010.View at: Google Scholar
Z. Wang, “Crack growth, computer strength and deformation of nonhomogeneous materials (concrete),” Tsinghua University, Beijing, China, 1996, Ph.D. thesis.View at: Google Scholar
X. Tang and C. Zhang, “Layering disposition and FE coordinate generation for random aggregate arrangements,” Tsinghua Science and Technology, vol. 48, no. 12, pp. 2048–2052, 2008.View at: Google Scholar
C. Qin, C. Guo, and C. Zhang, “A pre-processing scheme based on background grid approach for meso-concrete mechanics,” Journal of Hydraulic Engineering, vol. 42, no. 8, pp. 941–948, 2011.View at: Google Scholar
J. C. Walraven and H. W. Reinhardt, Theory and Experiments on the Mechanical Behaviour of Cracks in Plain and Reinforced Concrete Subjected to Shear Loading, Delft University of Technology, Delft, Netherlands, 1981.
H. Ma and H. Chen, Study on the Mechanism of Dynamic Damage Damage of Concrete Dam and Its Microscopic Mechanics Analysis Method, China Water and Power Press, Beijing, China, 2008.
H. Ma, H. Chen, and B. Li, “Meso-structure numerical simulation of concrete specimens,” Journal of Hydraulic Engineering, vol. 35, no. 10, pp. 27–35, 2004.View at: Google Scholar