IMPROVED RANDOM AGGREGATE MODEL FOR NUMERICAL SIMULATIONS OF CONCRETE ENGINEERING SIMULATIONS OF CONCRETE ENGINEERING

Abstract In numerical simulations, concrete is usually considered as a three-phase material consisting of an aggregate, a cement matrix, and an interfacial transition zone (ITZ). Three-dimensional modeling of concrete usually requires extremely large computational requirements. In this study, an improved random aggregate model for numerical simulations of concrete is developed, which can minimize the number of elements, optimize the ITZ thickness, and create internal cracks and holes. Numerical investigations on the cracks form as well as deflection and tensile strength are also conducted based on three-point bending tests. The simulation results agree well with the experimental results.


Introduction
Concrete is a construction material that has been extensively studied in the field of civil engineering (Kudzys, Kliukas 2010). Various kinds of models have been developed by different researchers to investigate the fracture process and mechanical properties of concrete. These models can be classified into two groups: two-dimensional (2D) and three-dimensional (3D) models. The 2D models include a simple finiteelement 2D model (Berthelot, Fatmi 2004), three-phase composite lattice models (Prado, Mier 2003;Lilliu, van Mier 2007), a generalized beam (GB) lattice model (Liu et al. 2007), a 2D random aggregate notched concrete beam model (Skarżyń ski, Tejchman 2010), a high aggregate content concrete model (Leite et al. 2004). The 3D models include a 3D beam lattice model (Man, van Mier 2008), random aggregate structure models (Zhou, Hao 2008;Ma et al. 2007;Tian et al. 2010), and an irregular lattice model (Kim, Lim 2011). Berthelot and Fatmi (2004) investigated the generation of damage in concrete using a simple finite element 2D concrete model. Skarżyń ski and Tejchman (2010) built a random 2D aggregate notched concrete beam model and utilized this model to investigate the fracture process zones (FPZ). They also found that this model could reduce the number of finite elements and the corresponding calculation time. Liu et al. (2007) and Liu and Liang (2009) developed a new GB lattice model that could significantly reduce the computational time. The effects of particle overlay on the fracture process have also been discussed by him. In general, the 2D models require smaller data set and less computation time, but they cannot perfectly reflect the complexity of the internal structure of concrete.
Although the simulation results from 2D models are comparable to those from 3D models in many cases (Vořechovský, Sadílek 2008), 2D models cannot replace 3D models for a number of reasons. First, 2D models do not offer sectional views or capture the development of space fractures of the entire concrete model. Second, in 2D models, the initial defects in concrete are not taken into consideration.
A three-phase composite 3D lattice model was used to investigate the fracture process as well as the effects of thickness and size of the interfacial transition zone (ITZ) on the strength of concrete (Man, van Mier 2008;Lilliu, van Mier 2003). Leite et al. (2004) built a concrete model with high aggregate content and realistic distribution, and utilized this model to simulate direct tension and wedge-splitting. Berthelot et al. (2008) investigated the damage evolution in heterogeneous materials based on a 3D model. The advantages of using 3D models include: (1) more accurate reflection of the internal structure of concrete; (2) capability of simulating various states (such as different boundary conditions, loading types, etc.); and (3) better visualization of the numerical results than 2D models. However, using 3D models involve huge amounts of data, resulting in slower computation speed and longer computation time than using 2D models. Some researchers used parallel computers to increase the computation speed, but this would cause high costs and is not practical in common engineering practice.
It is desirable to develop a new model which can simulate the inner structure concrete correctly and completely while improving the calculation precision and speed. In this study, one such model is developed by improving an existing model. it is found that the improved model has high calculation precision, and requires less computation time and less computer memory. The detailed information of this model is listed as follows: a) A traditional random aggregate structure concrete model (Tian et al. 2010) is chosen as a base to build the improved model. The traditional model consists of three phases, namely, the aggregate, cement matrix, and interfacial transition zone (ITZ). The thickness of ITZ is usually very small in experiments, even few millimeters (Lilliu, van Mier 2003), but in the numerical simulation the thickness of ITZ always more larger than that in the experiment, and in numerical simulation, in order to reduce the number of elements, it is very hard to simulate such tiny element. The ITZ is a medium located between the aggregate and cement matrix. The strength of ITZ is lower than those of the aggregate or cement matrix. The thickness of ITZ affects the accuracy of the numerical simulation. Thus, the thickness of ITZ has to be optimized in the improved model. Meanwhile, less data should be used in the improved model to reduce the computational requirements; b) Concrete is a man-made material that contains randomly distributed internal cracks and holes (Ganguli et al. 2012). In most finite element numerical simulations, people ignored this feature of concrete material because they assumed that the elements in concrete were continuous. In this study, it is attempted to find a method of including these internal cracks and holes in the simulation. In addition, the effects of these cracks and holes on the fracture processes and mechanical properties of concrete are investigated. In this study, the simulation results are compared with the experimental results in terms of the deflection, cracks form and computational accuracy. It is found that the improved model developed in this study is viable, practical, and better than traditional models in many aspects.

Definition of inherent defects
The notched concrete beam models built by Skarżyński and Tejchman (2010), Wu et al. (2006), Sagar and Prasad (2011) and Roesler et al. (2007) considered exterior defects on a concrete beam. Internal defects (such as cracks and holes) were not considered in these models. However, such internal defects do exist in concrete which is mixed by machine or manually. In this paper, these internal cracks and holes are referred to as inherent defects in concrete.

Mesoscale modeling approach
The flow chart for the model development is shown in Fig. 1. The key points are listed as follows: a) Based on programming and stochastic mathematical methods, a program is developed that can randomly generate the coordinates of every aggregate; b) The geometric dimensions of the concrete used in experiments and the numerical model are shown in Fig. 2(a). The number of elements is compressed. The lattice model ( Fig. 2 which accounts for 84.27% and 86.22% of the total number of elements and nodes in the entire concrete model, respectively; c) Hexahedral elements are widely used in traditional modeling methods. In turn, the ITZ in those models is very thick, resulting in low computational accuracy. In this study, tetrahedral elements are used to build the new model. Afterward, the file (which is generated in step ''a'') is loaded into the new model. Then the integrity concrete beam model is built using the method illustrated in Fig. 3. If every vertex of the lattice is located within the outline of the aggregate, this element will be considered as a component of the aggregate. If every vertex of the lattice is located outside the outline of the aggregate, this element will be considered as a component of the cement matrix. Otherwise, this element will be considered as a component of the ITZ. The three-phase random aggregate structure is then built based on this method; d) To optimize the ITZ thickness, the ITZ is subdivided (the flow chart is shown in Fig. 4), and step ''c'' is repeated. This optimizing step is repeated for two times during the development of the new model. Thus, the thickness of the ITZ in the new model is significantly reduced (   shows the comparison between the improved model and traditional models). While this optimization can effectively reduce the ITZ thickness, it should not be repeated for too many times because the repeated optimization step would produce more elements require better computer performance. After optimizing the ITZ, the random aggregate model is built; e) The random selection process is accomplished by programming. The selected element is named as the ''defective element n'', where n represents the sequence number of the element. In this paper, it is assumed that the defective elements comprise 1% of the total number of elements within the computational region. Therefore the total number of defective elements is 2356, thus n 0(1, 2, . . ., 2356).
To ensure that the lattices are continuous, the principle of ''ekill,n'' determines the elastic modulus for each defective element. The value of the elastic modulus should be very small. Thus, these defective elements can be regarded as inherent defects in the concrete material. Here, ''ekill'' means ''kill element.'' As the final step, the random aggregate concrete model with inherent defects, hereafter called the improved model, is built, and its sectional view is shown in Fig. 5(a). As a comparison, the sectional view of a typical traditional model is shown in Fig. 5(b).

Parameters and states
The parameters for the numerical simulation are listed in Table 1 (Appendix A). The load is applied at the top center of the computational region ( Fig. 2(a)). The states are shown in Table 2 (Appendix A).

Fracture criterion and constitutive model discussion
Concrete is a quasi-brittle material. The cracks caused by the loading usually leads to a non-linear stress-strain curve. Thus the concrete elastic damage model has been used in the simulation to describe the meso-mechanics properties of the concrete. This model describes the material degradation using a single damage variable D growing monotonically from 0 (undamaged material) to 1 (completely damaged material) (Marzec et al. 2007). The damage variable D is associated with material degradation due to the propagation and coalescence of microcracks and holes. The stress-strain relationship is represented by: r ¼Ẽe; (1) where: E 0 is the initial elastic modulus,Ẽ is the elastic modulus after the damage is produced. s and o are the nominal stress and nominal strain, respectively. The bilinear damage evolvement model is used in the simulation. The damage variable D is determined by Eq. (2) and plotted in Fig. 6: where: o 0 is the principal strain when the tensile strength is equal to the tensile strength of concrete f t ; f tr 0l × f t (0Bl51) is the residual tensile strength of concrete; l is the residual tensile strength factor; o r is the residual strain when the tensile strength is equal to f tr ; o r 0h × o 0 (1Bh 55); h is the residual strain factor; o u 0j × o 0 (j h) is the ultimate tensile strain; j is the ultimate tensile strain factor; and o max is the maximum value of principal strain in the course of the loading.
If o max Bo 0 , the element is undamaged. If o 0 Bo max 5o r , the element is in the first damage stage. If o r Bo max 5o u , the element is in the second damage stage. The second damage stage appears after the first The numerical damage results obtained using the improved model and traditional models are shown in Figs 7 and 8, respectively. According to Equation (2), there is no damaged element in the model when D 00. Therefore, elements retain their original color in the sectional view (aggregate 0light blue, cement matrix 0red, ITZ 0purple). When 0BD B1, there are damaged elements at the bottom of the model. These elements are not completely damaged. However, their mechanical properties have been altered because their elastic modulus has been significantly reduced. In the sectional view, these elements are shown in multicolor. When D 01, some elements are completely damaged and cracks are formed. The cracks appear in white color in the sectional view.
It can be seen from Figs 7 and 8 that, for both improve model and traditional models, damaged elements always appear initially at the bottom of the model in the early course of loading.
In addition, the cracks always appears initially at the bottom center of these damaged elements (such as in step 10). However, there are some differences between the two kinds of models: (1) In step 10, the crack has a tendency to move to inherent defects in the improved model, as shown in Fig. 9(a). The stress concentration area (red area) appears next to the inherent defects and cracks (white area), which will enlarge the area of the inherent defects and allow these defects to merge with the cracks in step 12 ( Fig. 9(b)). It has been proven that inherent defects can significantly affect the form of the cracks.
(2) The cracks are formed faster in the improved model than in the traditional model, because defects are found to be connected with each other when the improved model is used. The connection of defects is found in the entire course of loading in the improved model, but not in the traditional model.
In step 16, the concrete model has reached the tensile strength f t . Most inherent defects located near the bottom of the model are enlarged and connected to big cracks (such as A in Fig. 7(a)). However, the inherent defects located away from the bottom do not show significant changes (such as B in Fig. 8(b)). This indicates that a new trend of the formation of the cracks can be simulated in the improved model, but not by the traditional model.
The geometric dimensions of the concrete specimens and the boundary conditions in the acoustic emission (AE) experiment conducted by Wang et al. (2011) are shown in Fig. 2(a). The inherent defects do not change at the beginning of the AE experiment, indicating the absence of an acoustic emission. Then acoustic emissions occur at the bottom center of the specimen (Fig. 10(a)) due to the initial formation of cracks in this region. As the experiment continues, a and 10(c)), due to the increase in the size of inherent defects and the formation of connections between some of the inherent defects. Meanwhile, a small number of inherent defects connect to the big cracks which originate from the bottom center of the specimen.
The cracks and inherent defects of the concrete model are extracted in the numerical simulations, as shown in Fig. 11. Figs 11(a) and 11(b) clearly show the increase in the size of the inherent defects (such as A) in the right corner and at the middle center of the computational region. Some of the inherent defects are connected to each other (such as B in Fig. 11(b)), whereas some of the inherent defects are connected to the big cracks (such as C Fig. 11(b)). These results are consistent with those of the AE experiment. However, such results cannot be found by using the traditional model (Figs. 11(c) and 11(d)).
The numerical simulation results obtained by the improved model are more visualized and credible than those by the traditional model. Therefore, the improved model is a feasible replacement for the traditional model in cracks simulations.

Damage analysis
The data from the numerical results are extracted to draw a histogram diagram describing the number of the damaged elements (Fig. 12). The histogram diagram shows that: (1) A small number of damaged elements in each state exist at the initial stage of the simulation [P(t)54 MPa], and this number increases with the loading rate. In addition, the number of the damaged elements is higher in the improved model than that in the traditional model under the same states. For example, the improved model has 1,689 damaged elements in state 3, when P(t)04 MPa, but the traditional model has only 1,502 damaged elements. This difference shows that, under the same states, using the improved model increases the number of damaged elements by 12.45%; (2) As the loading process continues, the number of damaged elements increases rapidly when P(t) $ 4.5 MPa. At this time, the specimen almost reaches its tensile strength. For example, the total number of damaged elements in stage 3 is increased from 2,871 to 4,649 when using the improved model and from 2,864 to 4,513 when using the traditional model. Using the improved model increases the number of damaged elements by 3.01%; (3) When P(t)06 MPa, the number of damaged elements in state 3 is 20,601 when using the improved model and 20,287 when using the traditional model, which shows a 1.55% difference in the number of damaged elements.
In general, inherent defects significantly affect the generation of the cracks and the total number of damaged elements during the early stages. However, this effect is gradually reduced afterwards because a large number of cracks are produced.

Strength analysis
In order to check the computational accuracy of the improved model, the load-deflection curve of a concrete specimen is obtained and shown in Fig. 13.
The x-axis is defined as the deflection of the concrete specimen and the y-axis as the loading value. The following conclusions can be drawn from the curve: (1) The characteristics of the load-deflection curves from both models (the improved model and the traditional model) are similar to those from experiments, indicating that the improved model has a great agreement with experiments; (2) During the entire loading process, the deflection of the improved model is always larger than those of the traditional model (OA in Fig. 13) under the same loading states (such as in state 1 during the early stage. And this trend is reversed during the later stage (AB in Fig. 13). The load-deflection curve from the improved model is more close to experimental curves  Fig. 13. LoadÁdeflection curves under different loading stages than that from the traditional model. This indicates that the inherent defects within the concrete reduce the tensile strength and increase the deflection of the concrete specimen in the early loading stage. In the later loading stage, however, the deflections of two numerical models are close to each other because the number of cracks is produced and the influence of inherent defects becomes negligible; (3) In order to compare the performance of the improved model with the traditional model, the results of the numerical simulations and those of the experiment conducted by Zhou et al. (2009Zhou et al. ( , 2010 are summarized, as shown in Table 3 (Appendix A). It is found that the numerical results of the deflection using the improved model are more accurate than those using the traditional model. This can be illustrated by the following example.
In state 1, the deflections corresponding to the maximum load on the concrete specimen under different loading states are compared. The accuracy of the deflection is 93.32% in the improved model and 89.09% in the traditional model, with a difference of 4.23%. In states 2 and 3, this difference is reduced to 1.22% and 1% respectively. The average accuracy (three states) of the deflection in the improved model is 89.24% as compared with 87.09% for the traditional model.
In summary, the numerical results of the improved model are closer to the experimental results than those of the traditional model. In addition, the computational accuracy of the improved model is higher than that of the traditional model. Results also show that the inherent defects cannot be ignored in numerical simulations. However, the inherent defects are ignored by traditional models.

Conclusions
In this study, the effects of inherent defects on concrete in terms of its cracks and deflection were investigated. Three loading states were simulated under three point bending experiments. The results from simulations matched the experimental results well, and the following conclusions can be drawn: (1) Inherent defects significantly affect the cracks and deflection of concrete, and the effects should be included in numerical simulations because they play a key role in determining the generation and development of cracks. The inherent defects significant affect the deflection values of concrete, especially in the early stage of loading. These phenomena have been observed in experiments but are not considered in the numerical simulations using the traditional model; (2) The computational accuracy of the improved model is higher than that of the traditional model. In addition, the numerical results from the improved model are more visualized than those from the traditional model. Therefore, the improved model is more practically applicable than the traditional model.