# Effects of Voids on Concrete Tensile Fracturing: A Mesoscale Study

Authors: Lei Xu; Yefei Huang
Journal: Advances in Materials Science and Engineering
 (2017)

---

## Abstract

A two-dimensional mesoscale modeling framework, which considers concrete as a four-phase material including voids, is developed for studying the effects of voids on concrete tensile fracturing under the plane stress condition. Aggregate is assumed to behave elastically, while a continuum damaged plasticity model is employed to describe the mechanical behaviors of mortar and ITZ. The effects of voids on the fracture mechanism of concrete under uniaxial tension are first detailed, followed by an extensive investigation of the effects of void volume fraction on concrete tensile fracturing. It is found that both the prepeak and postpeak mesoscale cracking in concrete are highly affected by voids, and there is not a straightforward relation between void volume fraction and the postpeak behavior due to the randomness of void distribution. The fracture pattern of concrete specimen with voids is controlled by both the aggregate arrangement and the distribution of voids, and two types of failure modes are identified for concrete specimens under uniaxial tension. It is suggested that voids should be explicitly modeled for the accurate fracturing simulation of concrete on the mesoscale.


---

## Body

## 1. Introduction

Concrete is widely used as a construction material and is traditionally treated as a homogeneous continuum on the structural scale (macroscale). This homogenization assumption can hold well as long as the mechanical response of concrete remains in the elastic regime [1, 2]. However, when fracturing occurs, the macroscale mechanical behavior of concrete is greatly controlled by its components and their interactions taking place on a finer scale (mesoscale) [3, 4], which means accurate modeling of concrete fracturing calls for the consideration of its mesostructure.Up to date, several mesoscale models have been developed to provide tools for a better understanding of concrete fracturing. From the simulation strategy point of view, most of the existing concrete mesoscale models can be broadly grouped into two types: the continuum model and the lattice model. In the continuum model, concrete is usually characterized by a continuum composite material with each component discretized by finite elements, while, for the lattice model, a discrete system composed of lattice elements is used to represent concrete. Moreover, the discrete element method (DEM) has been recently used to perform the mesoscale simulation of concrete [5], and it is shown that the discrete model requires a huge numerical effort that is necessary for this approach to obtain a reasonable representation of concrete mesostructure.Several researchers studied the concrete fracturing by employing the continuum modeling strategy, and representative contributions can be found in [6–12]. The most recent investigations following this strategy were carried out by Du et al. [13] who studied the dynamic tensile fracturing of concrete by assuming concrete to be composed of aggregate and mortar matrix, by Huang et al. [14] who performed a 3D mesoscale fracturing simulation based on the actual concrete mesostructure, and by Wang et al. [15, 16] who developed a computational technology using the interface element with a cohesive law to perform Monte Carlo simulations of concrete fracturing and to study the 3D mesostructure effects on concrete damage and failure. Overall, the principal merit of the continuum model lies in the detailed representation of concrete mesostructure, which ensures the ability to realize reasonable simulations of cracking initiation on the mesoscale and coalescence of multiple distributed cracks into localized macroscale cracks and fracture propagation. However, it tends to be computationally intensive even for laboratory-scale specimens, especially for three-dimensional cases.With respect to the lattice modeling strategy, representative studies were carried out in [17–20], and the most recent improvements are performed by Cusatis et al. [21, 22] who proposed a novel model named the lattice discrete particle model (LDPM) by exploiting the merits of both the lattice model and the discrete particle model. In contrast to the continuum model, the lattice model is considered computationally less demanding as concrete mesostructure is roughly represented by a discrete system with relatively less degrees of freedom and meanwhile can still possess the ability to capture the most important aspects of concrete fracturing. However, it is hard to investigate the interactions of concrete components in a real sense since the actual concrete mesostructure is not fully taken into account in the lattice model.Voids (or pores) with different sizes always exist in concrete and typically take up 2–6% of the total volume, and the use of entrained air void system is a common approach in concrete technology to resist cyclic freezing and thawing degradation [23]. However, the effects of voids on concrete fracturing on the mesoscale are still not well understood. Wang et al. [15] built numerical concrete samples with pores using interface elements and studied the effects of porosity on concrete loading-carrying capability under uniaxial tension, but the fracturing mechanism on the mesoscale was not detailed. Huang et al. [14] reported the distribution of voids greatly influences the tensile strength and crack patterns based on the simulation results of a single 3D specimen. On the whole, it has been recognized that the existence of voids affects the concrete mechanical behavior to a large extent, but further research is needed to reveal the effects of voids on concrete fracturing.With this in mind, a 2D finite element (FE) mesoscale modeling framework for concrete is proposed in this study in which concrete is considered as a four-phase material composed of aggregate, mortar, interfacial transitional zone (ITZ), and void, and the effects of voids on concrete tensile fracturing under the plane stress condition are detailed by performing several simulations. The rest of this paper is organized as follows: Section2 presents the generation procedures of concrete mesostructure; the FE modeling methodology including mesh discretization, insertion of ITZ elements, and constitutive modeling of mortar and ITZ is described in Section 3; in Section 4, the effects of voids on concrete tensile fracturing are discussed in detail based on the simulation results of several concrete specimens with different mesostructures; finally, the study is summarized with conclusions in Section 5.


## 2. Generation of Concrete Mesostructure

In this study, concrete is treated as a four-phase composite material, that is, coarse aggregate, mortar composed of cement matrix and fine aggregate, interfacial transitional zone (ITZ), and void randomly distributed in the mortar. Regarding aggregate generation, gravel is idealized as circle, while crushed aggregate is considered as polygon. Mortar is assumed as a homogenous continuum, and the interface with a specified thickness between coarse aggregate (hereinafter referred to as aggregate) and mortar is used to represent ITZ. Moreover, void is viewed as circle for simplicity.


### 2.1. Size Distribution of Aggregates and Voids

The aggregate size distribution of concrete is described by Talbot’s equation as(1)Fpda=dadmaxn,where da is the size of aggregate, dmax is the maximum size of aggregate, Fpda represents the ratio of aggregates by weight passing through a sieve of characteristic size equal to da, and n is the exponent of Talbot’s equation. For n=0.5, the corresponding curve is known as Fuller’s curve extensively employed in concrete grading design for optimal packing properties.For a concrete specimen with total volumeV, the volume of aggregates within a grading segment [di,di+1] can be calculated by(2)Vaggdi,di+1=Fpdi+1-FpdiFpdmax-Fpdmin×AF×V, where dmin is the minimum size of aggregate and AF represents the aggregate volume fraction.Currently, the size distribution of voids in concrete has not been detailed. In general, these voids can be broadly grouped into two types according to different formation ways and the resulting different sizes: the (smaller) entrained voids with typical sizes on the order of 0.1 mm and the (larger) entrapped voids with typical sizes commonly more than 1 mm. In this study void size is considered to be uniformly distributed, and the same assumption is also employed by other researchers [15, 16]. Thus, denoting the size range of void by [dminv,dmaxv], the void size can be calculated by dv=dminv+P×(dmaxv-dminv) (P is a uniformly distributed random number between 0 and 1).


### 2.2. Generation and Placement of Aggregates and Voids

In order to build numerical concrete specimens automatically, a mesostructure generator for concrete (MGC) is developed using MATLAB based on the take-and-place method [24, 25].In the take-process, aggregates and voids, which will be placed into the specimen volume in the place-process, are generated separately. For the aggregate generation, the aggregate volume for each grading segment is first calculated according to (2). Then, starting with the grading segment with the maximum average size, the aggregates are generated one by one for each grading segment. For a certain grading segment [di,di+1], the generation of aggregates takes the following procedures.Step 1.
Generate a random number representing the aggregate sizeda, which is assumed to follow a uniform distribution and therefore can be taken as da=di+P×(di+1-di).Step 2.
For gravel, a circle with radius ofda/2 is defined to represent the aggregate, while, for crushed aggregate, a polygon with the random number of sides ranging from 4 to 10 and with the smallest width equal to da is generated to represent the aggregate (see [24] for more details). Then, the volume of the current generated aggregate is calculated.Step 3.
Repeat the previous two steps until the remaining volume left is less thanπda2/4, namely, not enough to generate a new aggregate.Step 4.
Transfer the remaining volume to the next grading segment.Following the similar procedures for generating gravel aggregates, the generation of voids can be performed with ease provided by the given void volume fraction and size range, which is followed by the placement of aggregates and voids (the place-process).In the place-process, the generated aggregates and voids are first sorted according to their volume, respectively. Then, for the convenience of mesh discretization discussed in Section3, the size of each aggregate is increased by a specified value (the thickness of ITZ, TITZ) to consider the surrounding ITZ, which means the aggregate finally placed consists of two parts (i.e., aggregate piece and the surrounding ITZ with a specified thickness). After the modification of aggregate size, all aggregates are placed into the specified specimen one by one starting with the aggregate with the largest volume, followed by the placement of voids starting from the biggest one. The procedures of the placement of aggregates and voids are detailed as follows.Step 1.
Define the shape of concrete specimen usingX and Y coordinates of the boundary vertices numbered clockwise or anticlockwise, which will be used in Step  3 to check if an aggregate is inside the concrete specimen.Step 2.
Generate random numbers to define the position (and orientation if polygon is used to represent the crushed aggregate) of the aggregate usingXmin, Xmax, Ymin, and Ymax, which represent the minimum and maximum value of X and Y coordinates of all boundary vertices, respectively.Step 3.
Perform the aggregate placing. The placement is considered to be successful if the following four conditions are satisfied: (1) the whole aggregate should be within the concrete specimen; (2) no overlapping/intersection occurs between the aggregate to be placed and any existing aggregate; (3) a minimum distance (defined by r1) should exist between the aggregate and the specimen boundary; and (4) a minimum gap (defined by r2) should exist between any two aggregates. If any of the four conditions is violated, Step  2 is repeated to make a new attempt until the placement of the aggregate is completed.Step 4.
Repeat Steps  2-3 until all aggregates are placed inside the specimen.After the placement of aggregates, voids should also be placed into the specimen, which can be carried out by following the similar steps given above. It is worth noting that voids are considered to be embedded in mortar in this study.Using MGC, numerical concrete specimens can be built with ease. The specimens shown in this paper are 100 mm squares, and the 4-segment Fuller curve is used to describe the aggregate grading for all specimens. For aggregates,dmax and dmin are set to 10 mm and 4 mm, respectively, while, for voids, dmaxv = 2 mm and dminv = 1 mm are used. In addition, TITZ is approximately set to 100 μm according to the experimental observation [1], and both r1 and r2 are taken as 0.1 times of the size of aggregate or void to be placed.Figures1(a) and 1(b) sketch two numerical samples using circular aggregates with the same void volume fraction (Av=2%) and the aggregate volume fraction AF=40% and 50%, respectively, whereas two numerical samples using polygonal aggregates with Av=2% and AF=40% and 50%, respectively, are plotted in Figures 1(c) and 1(d).Figure 1
Numerical concrete specimens. (a) Circular aggregate (AF=40%, Av=2%). (b) Circular aggregate (AF=50%, Av=2%). (c) Polygonal aggregate (AF=40%, Av=2%). (d) Polygonal aggregate (AF=50%, Av=2%).
(a)
(b)
(c)
(d)


## 2.1. Size Distribution of Aggregates and Voids

The aggregate size distribution of concrete is described by Talbot’s equation as(1)Fpda=dadmaxn,where da is the size of aggregate, dmax is the maximum size of aggregate, Fpda represents the ratio of aggregates by weight passing through a sieve of characteristic size equal to da, and n is the exponent of Talbot’s equation. For n=0.5, the corresponding curve is known as Fuller’s curve extensively employed in concrete grading design for optimal packing properties.For a concrete specimen with total volumeV, the volume of aggregates within a grading segment [di,di+1] can be calculated by(2)Vaggdi,di+1=Fpdi+1-FpdiFpdmax-Fpdmin×AF×V, where dmin is the minimum size of aggregate and AF represents the aggregate volume fraction.Currently, the size distribution of voids in concrete has not been detailed. In general, these voids can be broadly grouped into two types according to different formation ways and the resulting different sizes: the (smaller) entrained voids with typical sizes on the order of 0.1 mm and the (larger) entrapped voids with typical sizes commonly more than 1 mm. In this study void size is considered to be uniformly distributed, and the same assumption is also employed by other researchers [15, 16]. Thus, denoting the size range of void by [dminv,dmaxv], the void size can be calculated by dv=dminv+P×(dmaxv-dminv) (P is a uniformly distributed random number between 0 and 1).


## 2.2. Generation and Placement of Aggregates and Voids

In order to build numerical concrete specimens automatically, a mesostructure generator for concrete (MGC) is developed using MATLAB based on the take-and-place method [24, 25].In the take-process, aggregates and voids, which will be placed into the specimen volume in the place-process, are generated separately. For the aggregate generation, the aggregate volume for each grading segment is first calculated according to (2). Then, starting with the grading segment with the maximum average size, the aggregates are generated one by one for each grading segment. For a certain grading segment [di,di+1], the generation of aggregates takes the following procedures.Step 1.
Generate a random number representing the aggregate sizeda, which is assumed to follow a uniform distribution and therefore can be taken as da=di+P×(di+1-di).Step 2.
For gravel, a circle with radius ofda/2 is defined to represent the aggregate, while, for crushed aggregate, a polygon with the random number of sides ranging from 4 to 10 and with the smallest width equal to da is generated to represent the aggregate (see [24] for more details). Then, the volume of the current generated aggregate is calculated.Step 3.
Repeat the previous two steps until the remaining volume left is less thanπda2/4, namely, not enough to generate a new aggregate.Step 4.
Transfer the remaining volume to the next grading segment.Following the similar procedures for generating gravel aggregates, the generation of voids can be performed with ease provided by the given void volume fraction and size range, which is followed by the placement of aggregates and voids (the place-process).In the place-process, the generated aggregates and voids are first sorted according to their volume, respectively. Then, for the convenience of mesh discretization discussed in Section3, the size of each aggregate is increased by a specified value (the thickness of ITZ, TITZ) to consider the surrounding ITZ, which means the aggregate finally placed consists of two parts (i.e., aggregate piece and the surrounding ITZ with a specified thickness). After the modification of aggregate size, all aggregates are placed into the specified specimen one by one starting with the aggregate with the largest volume, followed by the placement of voids starting from the biggest one. The procedures of the placement of aggregates and voids are detailed as follows.Step 1.
Define the shape of concrete specimen usingX and Y coordinates of the boundary vertices numbered clockwise or anticlockwise, which will be used in Step  3 to check if an aggregate is inside the concrete specimen.Step 2.
Generate random numbers to define the position (and orientation if polygon is used to represent the crushed aggregate) of the aggregate usingXmin, Xmax, Ymin, and Ymax, which represent the minimum and maximum value of X and Y coordinates of all boundary vertices, respectively.Step 3.
Perform the aggregate placing. The placement is considered to be successful if the following four conditions are satisfied: (1) the whole aggregate should be within the concrete specimen; (2) no overlapping/intersection occurs between the aggregate to be placed and any existing aggregate; (3) a minimum distance (defined by r1) should exist between the aggregate and the specimen boundary; and (4) a minimum gap (defined by r2) should exist between any two aggregates. If any of the four conditions is violated, Step  2 is repeated to make a new attempt until the placement of the aggregate is completed.Step 4.
Repeat Steps  2-3 until all aggregates are placed inside the specimen.After the placement of aggregates, voids should also be placed into the specimen, which can be carried out by following the similar steps given above. It is worth noting that voids are considered to be embedded in mortar in this study.Using MGC, numerical concrete specimens can be built with ease. The specimens shown in this paper are 100 mm squares, and the 4-segment Fuller curve is used to describe the aggregate grading for all specimens. For aggregates,dmax and dmin are set to 10 mm and 4 mm, respectively, while, for voids, dmaxv = 2 mm and dminv = 1 mm are used. In addition, TITZ is approximately set to 100 μm according to the experimental observation [1], and both r1 and r2 are taken as 0.1 times of the size of aggregate or void to be placed.Figures1(a) and 1(b) sketch two numerical samples using circular aggregates with the same void volume fraction (Av=2%) and the aggregate volume fraction AF=40% and 50%, respectively, whereas two numerical samples using polygonal aggregates with Av=2% and AF=40% and 50%, respectively, are plotted in Figures 1(c) and 1(d).Figure 1
Numerical concrete specimens. (a) Circular aggregate (AF=40%, Av=2%). (b) Circular aggregate (AF=50%, Av=2%). (c) Polygonal aggregate (AF=40%, Av=2%). (d) Polygonal aggregate (AF=50%, Av=2%).
(a)
(b)
(c)
(d)


## 3. Finite Element Modeling Methodology

Once the concrete mesostructure is obtained, a corresponding FE model is required for performing numerical simulations. The details of the FE modeling methodology developed in this study are presented as follows.


### 3.1. Mesh Discretization

In order to automatically carry out the mesh discretization of the concrete specimen with complex mesostructure, a mesh generator is developed using MATLAB by exploiting the powerful preprocessing modules provided by the commercial finite element software ABAQUS. For a numerical concrete specimen with pregenerated mesostructure, a two-part python file, which defines the boundary of the specimen together with the locations and shapes of aggregates and voids using the first part of the file and specifies the mesh discretization parameters using the second part of the file, is first generated using the mesh generator by taking concrete mesostructure as input data. Then, a FE mesh composed of linear triangular elements can be obtained by the mesh generator through calling ABAQUS/CAE kernel to execute the generated python file. An example of the FE mesh discretization with aggregate elements highlighted is shown in Figure2. It must be noted that ITZ elements are not included in the original generated FE mesh as the tiny thickness of ITZ makes the mesh discretization harder and leads to a poor mesh quality, and therefore a modification of the original FE mesh is needed for obtaining the final FE mesh, which is discussed in detail in Section 3.2.Figure 2
FE mesh discretization (polygonal aggregate).


### 3.2. Insertion of ITZ Elements

As stated in Section2, the size of an aggregate is increased by a specified value determined by the thickness of ITZ before placing this aggregate into the specimen volume. Hence, the aggregate elements in the FE mesh generated in Section 3.1 not only occupy the volume of aggregates but also take up the volume of ITZs. In order to explicitly represent the surrounding ITZs and the actual sizes of aggregates in the FE model, ITZ elements should be inserted between aggregate elements and the corresponding mortar elements and the coordinates of nodes on the boundaries of aggregate elements are required to be adjusted for accurately representing the actual aggregate sizes.To this end, a four-step procedure is proposed. Firstly, the original nodes on the boundaries of amplified aggregates are identified, followed by the definition of new nodes based on the coordinates of original nodes on the boundaries and the given thickness of ITZ. Then, the connectivities of the aggregate elements associated with these nodes are redefined by replacing the number of the original nodes by the number of the corresponding new nodes. Subsequently, 4-noded ITZ elements are formulated one by one using the original nodes and the corresponding new nodes. Finally, an updated input file for ABAQUS, which contains final mesh data, is generated. An in-house MATLAB program, which follows the above procedure, is developed, and part of the final FE mesh discretization corresponding to Figure2 is illustrated as an example in Figure 3.Figure 3
FE mesh discretization after inserting and adjusting.


### 3.3. Continuum Damaged Plasticity Model

Without considering voids, it is well recognized that ITZ is weaker than aggregate and mortar, and consequently mesoscale cracking in concrete under loading is commonly considered to first appear in ITZs. After that, mesoscale cracks propagate into mortar and additional cracks initiate within mortar with the further increase of loading [1]. In general, aggregates behave elastically during the process of concrete fracturing. Hence, the isotropic linear elastic model is employed to model the mechanical behavior of aggregates, while a continuum damaged plasticity (CDP) model [26, 27] implemented in ABAQUS is utilized to describe the mechanical behavior of both mortar and ITZ, which is briefly summarized as follows.In CDP model, two independent hardening variables, that is, equivalent compressive and tensile plastic strainsε~cp and ε~tp, are introduced for considering compressive crushing and tensile cracking, respectively. Then, two independent damage variables, dc(ε~cp) and dt(ε~tp), are defined to represent the compressive and tensile damage states. Furthermore, in order to describe the overall damage in an isotropic manner, a scale variable d is defined as(3)d=1-1-stdc1-scdt, where st and sc are functions of the stress state that are introduced to represent stiffness recovery effects associated with stress reversals (see [27] for more details).Then, the damaged elastic modulusE depending on different failure mechanisms under tension and compression can be described by(4)E=1-dE0 in which E0 represents the initial elastic modulus.Based on the concept of damage mechanics, the effective stressσ¯ can be calculated by(5)σ¯=σ1-d, where σ is the Cauchy stress.The yield function of CDP model is given in the effective stress space as(6)Fσ¯,ε~p=11-αq¯-3αp¯+βε~pσ¯^max-γ-σ¯^max-σ¯cε~cp,where ε~p=[ε~tpε~cp]T; p¯ and q¯ are the effective hydrostatic pressure and the effective Mises equivalent deviatoric stress, respectively; σ¯^max is the algebraically maximum eigenvalue of σ¯; the brackets 〈·〉 are used in Macaulay sense; σ¯c(ε~cp) is the uniaxial compressive effective strength; α and γ are dimensionless material constants, which can be determined by comparing the initial equibiaxial and uniaxial compressive yield stress and by comparing the yield conditions along the tensile and compressive meridians, respectively; and β(ε~p) can be calculated by(7)βε~p=σ¯cε~cpσ¯tε~tp1-α-1+α in which σ¯t(ε~tp) is the uniaxial tensile effective strength. Figure 4 illustrates the yield surface in the case of plane stress.Figure 4
Yield surface in plane stress.In order to describe the dilatancy reasonably, nonassociated flow rule is employed in CDP model, and the flow potential takes the form as(8)Gσ¯=ϵσt0tan⁡ψ2+q¯2-p¯tan⁡ψ in which ϵ is a parameter defining the rate at which the function approaches the asymptote; σt0 is the uniaxial tensile stress at failure; ψ is the dilation angle measured in p¯-q¯ plane at high confining pressure.As presented earlier, the material softening under tension is defined by the relationship between the uniaxial tensile effective strength and equivalent tensile plastic strain (see (7)), which means mesh sensitivity will be encountered when applying the CDP model in FE simulations. Therefore, a stress-displacement relation is used in this study to define the tensile softening behavior for alleviating the influence of mesh sensitivity on the simulation results.


### 3.4. Numerical Solution Algorithm

Due to the highly nonlinear and softening behavior of concrete in the process of fracturing, the ABAQUS/Explicit solver is employed in all simulations with the aim of capturing the entire fracturing process of concrete.As is well known, the dynamic effect inevitably exists in an explicit FE analysis, and its influence on the solution of a quasi-static problem should be small enough to be neglected. In order to minimise the dynamic effect, the loading time should be large enough, while, on the other hand, the computational effort increases proportionally with the increase of loading time. Hence, a balance has to be made between the computational efficiency and simulation accuracy, which can be achieved through comparing the results under different loading time (or loading rates).


## 3.1. Mesh Discretization

In order to automatically carry out the mesh discretization of the concrete specimen with complex mesostructure, a mesh generator is developed using MATLAB by exploiting the powerful preprocessing modules provided by the commercial finite element software ABAQUS. For a numerical concrete specimen with pregenerated mesostructure, a two-part python file, which defines the boundary of the specimen together with the locations and shapes of aggregates and voids using the first part of the file and specifies the mesh discretization parameters using the second part of the file, is first generated using the mesh generator by taking concrete mesostructure as input data. Then, a FE mesh composed of linear triangular elements can be obtained by the mesh generator through calling ABAQUS/CAE kernel to execute the generated python file. An example of the FE mesh discretization with aggregate elements highlighted is shown in Figure2. It must be noted that ITZ elements are not included in the original generated FE mesh as the tiny thickness of ITZ makes the mesh discretization harder and leads to a poor mesh quality, and therefore a modification of the original FE mesh is needed for obtaining the final FE mesh, which is discussed in detail in Section 3.2.Figure 2
FE mesh discretization (polygonal aggregate).


## 3.2. Insertion of ITZ Elements

As stated in Section2, the size of an aggregate is increased by a specified value determined by the thickness of ITZ before placing this aggregate into the specimen volume. Hence, the aggregate elements in the FE mesh generated in Section 3.1 not only occupy the volume of aggregates but also take up the volume of ITZs. In order to explicitly represent the surrounding ITZs and the actual sizes of aggregates in the FE model, ITZ elements should be inserted between aggregate elements and the corresponding mortar elements and the coordinates of nodes on the boundaries of aggregate elements are required to be adjusted for accurately representing the actual aggregate sizes.To this end, a four-step procedure is proposed. Firstly, the original nodes on the boundaries of amplified aggregates are identified, followed by the definition of new nodes based on the coordinates of original nodes on the boundaries and the given thickness of ITZ. Then, the connectivities of the aggregate elements associated with these nodes are redefined by replacing the number of the original nodes by the number of the corresponding new nodes. Subsequently, 4-noded ITZ elements are formulated one by one using the original nodes and the corresponding new nodes. Finally, an updated input file for ABAQUS, which contains final mesh data, is generated. An in-house MATLAB program, which follows the above procedure, is developed, and part of the final FE mesh discretization corresponding to Figure2 is illustrated as an example in Figure 3.Figure 3
FE mesh discretization after inserting and adjusting.


## 3.3. Continuum Damaged Plasticity Model

Without considering voids, it is well recognized that ITZ is weaker than aggregate and mortar, and consequently mesoscale cracking in concrete under loading is commonly considered to first appear in ITZs. After that, mesoscale cracks propagate into mortar and additional cracks initiate within mortar with the further increase of loading [1]. In general, aggregates behave elastically during the process of concrete fracturing. Hence, the isotropic linear elastic model is employed to model the mechanical behavior of aggregates, while a continuum damaged plasticity (CDP) model [26, 27] implemented in ABAQUS is utilized to describe the mechanical behavior of both mortar and ITZ, which is briefly summarized as follows.In CDP model, two independent hardening variables, that is, equivalent compressive and tensile plastic strainsε~cp and ε~tp, are introduced for considering compressive crushing and tensile cracking, respectively. Then, two independent damage variables, dc(ε~cp) and dt(ε~tp), are defined to represent the compressive and tensile damage states. Furthermore, in order to describe the overall damage in an isotropic manner, a scale variable d is defined as(3)d=1-1-stdc1-scdt, where st and sc are functions of the stress state that are introduced to represent stiffness recovery effects associated with stress reversals (see [27] for more details).Then, the damaged elastic modulusE depending on different failure mechanisms under tension and compression can be described by(4)E=1-dE0 in which E0 represents the initial elastic modulus.Based on the concept of damage mechanics, the effective stressσ¯ can be calculated by(5)σ¯=σ1-d, where σ is the Cauchy stress.The yield function of CDP model is given in the effective stress space as(6)Fσ¯,ε~p=11-αq¯-3αp¯+βε~pσ¯^max-γ-σ¯^max-σ¯cε~cp,where ε~p=[ε~tpε~cp]T; p¯ and q¯ are the effective hydrostatic pressure and the effective Mises equivalent deviatoric stress, respectively; σ¯^max is the algebraically maximum eigenvalue of σ¯; the brackets 〈·〉 are used in Macaulay sense; σ¯c(ε~cp) is the uniaxial compressive effective strength; α and γ are dimensionless material constants, which can be determined by comparing the initial equibiaxial and uniaxial compressive yield stress and by comparing the yield conditions along the tensile and compressive meridians, respectively; and β(ε~p) can be calculated by(7)βε~p=σ¯cε~cpσ¯tε~tp1-α-1+α in which σ¯t(ε~tp) is the uniaxial tensile effective strength. Figure 4 illustrates the yield surface in the case of plane stress.Figure 4
Yield surface in plane stress.In order to describe the dilatancy reasonably, nonassociated flow rule is employed in CDP model, and the flow potential takes the form as(8)Gσ¯=ϵσt0tan⁡ψ2+q¯2-p¯tan⁡ψ in which ϵ is a parameter defining the rate at which the function approaches the asymptote; σt0 is the uniaxial tensile stress at failure; ψ is the dilation angle measured in p¯-q¯ plane at high confining pressure.As presented earlier, the material softening under tension is defined by the relationship between the uniaxial tensile effective strength and equivalent tensile plastic strain (see (7)), which means mesh sensitivity will be encountered when applying the CDP model in FE simulations. Therefore, a stress-displacement relation is used in this study to define the tensile softening behavior for alleviating the influence of mesh sensitivity on the simulation results.


## 3.4. Numerical Solution Algorithm

Due to the highly nonlinear and softening behavior of concrete in the process of fracturing, the ABAQUS/Explicit solver is employed in all simulations with the aim of capturing the entire fracturing process of concrete.As is well known, the dynamic effect inevitably exists in an explicit FE analysis, and its influence on the solution of a quasi-static problem should be small enough to be neglected. In order to minimise the dynamic effect, the loading time should be large enough, while, on the other hand, the computational effort increases proportionally with the increase of loading time. Hence, a balance has to be made between the computational efficiency and simulation accuracy, which can be achieved through comparing the results under different loading time (or loading rates).


## 4. Results and Discussion

### 4.1. Numerical Specimens and Mechanical Properties

Aiming to investigate the effects of voids on the tensile fracturing of concrete with different aggregate volume fractions, three sets of numerical concrete specimens with dimensions of 100 mm × 100 mm using polygonal aggregates are generated, and each set contains four specimens with the same aggregate arrangement and differentAv (0%, 2%, 4%, and 6%, resp.). For Set I, AF is set to 30%, while AF=40% and 50% are chosen for Sets II and III, respectively. The typical element size is chosen as 0.4 times that of the minimum aggregate size; however the number of both elements and nodes of numerical specimens increases with the increase of Av. For example, the number of elements and nodes for the specimen with AF=40% and Av=2% is 15796 and 9136, respectively, and the number of elements and nodes for the specimen with AF=40% and Av=6% is 20528 and 12126, respectively.For each specimen, uniaxial tensile fracturing is simulated. In FE simulations, the left end of concrete specimen is fixed in the horizontal direction, while the opposite end is subjected to a uniformly distributed horizontal displacement up to 0.06 mm, namely, a displacement-controlled loading scheme is used. Following the strategy discussed in Section3.4, the loading time is set to 0.012 s, which corresponds to a loading rate 5 mm/s.The same mechanical properties of aggregate, mortar, and ITZ are adopted for all specimens, as listed in Table1. It is noted that the mechanical properties of mortar are directly obtained from Chinese code GB 50010-2002 (code for design of concrete structures), and the compression hardening curve and the tension softening curve are shown in Figures 5 and 6, respectively. Moreover, due to the lack of detailed experimental results for ITZ, the compressive and tensile strengths and elastic modulus of ITZ are assumed to be 75% of those of mortar for the sake of modeling the weaker ITZ, while the other mechanical properties of ITZ are taken as the same of mortar.Table 1
Mechanical properties of concrete components.
Material
E
0 (GPa)
Poisson’s ratio (—)
σ
t
0 (MPa)
ϵ (—)
ψ (°)
α (—)
γ (—)
Aggregate
60
0.2
—
—
—
—
—
Mortar
20
0.2
1.94
0.1
35
0.12
2.0
ITZ
15
0.2
1.46
0.1
35
0.12
2.0
Void
—
—
—
—
—
—
—Figure 5
Total stress-strain relation under uniaxial compression and compression hardening curve of mortar.Figure 6
Total stress-strain relation under uniaxial tension and tension softening curve of mortar.


### 4.2. Effects of Voids on Concrete Tensile Fracturing Mechanism

In order to study the effects of voids on concrete tensile fracturing mechanism, the simulation results of two specimens in Set II, which include one without void and the other one with voids taking up2% of the specimen volume, are discussed in detail in the following.The macroscopic stress versus displacement curve of the concrete specimen without void (AF=40%) under uniaxial tension is shown in Figure 7, and the elements with nonzero equivalent tensile plastic strain at different loading levels are illustrated in Figure 8, in which the elements with equivalent tensile plastic strain bigger than 100 microstrains are highlighted in red color.Figure 7
Stress-displacement curve (AF=40%, without void).Figure 8
Mesoscale cracking development at different loading levels of the concrete specimen without void(AF=40%). (a), (b), (c), (d), (e), and (f) show the equivalent tensile plastic strain distributions at points a–f, respectively.
(a)
(b)
(c)
(d)
(e)
(f)It can be observed that although concrete specimen initially exhibits elastic response on the macroscale, mesoscale cracking still occurs in ITZs (see Figure8(a)) due to their lower fracture properties than mortar. With the increase of applied displacement, cracking develops in ITZs and subsequently propagates into mortar. The mesoscopic crack at the peak stress is shown in Figure 8(b), and it can be found that a macroscale crack close to the right end is formed as a result of the coalescence of mesoscopic cracks. In the softening stage, strain localization, which is accompanied by the decrease of macroscopic stress (see Figure 7) and finally leads to macroscopic tensile failure, can be clearly identified from Figures 8(c), 8(d), 8(e), and 8(f).The macroscopic stress versus displacement curve of the concrete specimen (AF=40% and Av=2%) under uniaxial tension is shown in Figure 9, and the elements with nonzero equivalent tensile plastic strain at different loading levels are illustrated in Figure 10. Compared to the case without void discussed above, mesoscale cracking appears in both ITZs and the mortar around voids in the macroscale elastic stage (see Figure 10(a)), and, additionally, cracked mortar generally suffers from bigger equivalent tensile plastic strain (see Figure 10(b)) at the peak stress as a result of the existence of voids. Similar to the specimen without void, a macroscale crack is formed at the peak stress point due to the coalescence of mesoscopic cracks. However, the fracture pattern is quite different from the one observed in the specimen without void even if the same aggregate arrangement is employed in both specimens, which indicates the existence of void dominates the fracturing behavior of concrete to a large extent. The phenomenon of strain localization, which is accompanied by the decrease of macroscopic stress (see Figure 9) and leads to the final tensile failure, is visualized in Figures 10(c), 10(d), 10(e), and 10(f).Figure 9
Stress-displacement curve (AF=40%, Av=2%).Figure 10
Mesoscale cracking development at different loading levels of concrete specimen withAF=40% and Av=2%. (a), (b), (c), (d), (e), and (f) show the equivalent tensile plastic strain distributions at points a–f, respectively.
(a)
(b)
(c)
(d)
(e)
(f)Overall, different fracturing mechanisms can be observed for the two concrete specimens by comparing the development processes of mesoscale cracking shown in Figures8 and 10. For the concrete specimen without void, mesoscale crack first appears in ITZ and then propagates into the neighboring mortar, while, for the concrete specimen with voids, mesoscale cracks are first found in both ITZ and mortar around voids and then coalesce in the fracturing process.


### 4.3. Effects of Void Volume Fraction

In this section, the effects of aggregate volume fraction (AF) on concrete tensile fracturing are first investigated, followed by the detailed discussion on the effects of void volume fraction on concrete with different aggregate volume fractions.The specimens without void in Sets I, II, and III are simulated, and the macroscopic stress versus displacement curves are depicted in Figure11. It is shown that the elastic responses of the three specimens are almost identical, which are in general linear and can be approximately characterized by the same elastic modulus. As AF increases, the reached peak tensile stress slightly decreases (similar to the results reported in [15]), which can be attributed to the fact that more ITZs exist in the case of higher AF. On the other hand, the postpeak stress in the case of AF=40% generally drops more quickly than in the case of AF=30%, as expected. However, in the case of AF=50%, the postpeak drops less quickly than in the other two cases, which may be due to the fact that this specimen finally fails with a different failure mode (two dominant macroscale cracks; see Figure 15(a)) from that of the other two specimens (one dominant macroscale crack; see Figures 13(a) and 14(a)), leading to higher fracture energy dissipated in tensile fracturing. From the above discussion, it can be concluded that the postpeak behavior of concrete highly depends on the failure mode, namely, the more macroscale cracks appear, the higher fracture energy should be expected.Figure 11
Stress-displacement curves of different aggregate volume fractions.Provided the simulation results of Sets I, II, and III, the effects of void volume fraction (Av) on concrete tensile fracturing are analyzed in the following. Figures 12(a), 12(b), and 12(c) present the stress-displacement relations of concrete specimens with AF = 30%, 40%, and 50%, respectively, and Figures 13, 14, and 15 illustrate the fracture patterns of concrete specimens with AF = 30%, 40%, and 50%, respectively. In Figures 12(a), 12(b), and 12(c), it is shown that the macroscale elastic modulus in prepeak stage can be considered to be independent of Av, while the peak stress decreases as Av increases. Moreover, the relation between Av and the postpeak behavior is not straightforward as depicted. As shown in Figures 13, 14, and 15, the fracture patterns of concrete specimens with the same aggregate arrangement (represented by highlighted red polygons) and different Av differ from each other. For concrete specimens without void, fracture pattern is mainly controlled by the distribution of aggregates (or ITZs), whereas, for concrete specimens with voids, the fracture pattern is dominated not only by the aggregate arrangement but also by the distribution of voids in mortar, leading to the differences of fracture patterns illustrated in Figures 13, 14, and 15. Consequently, accurate fracturing analysis of concrete calls for the explicit modeling of voids. Furthermore, both the concrete specimen without void and the concrete specimens with voids present two types of failure modes: one is with a single macroscale crack (Type I, typically shown in Figure 15(b)), and the other is with two macroscale cracks (Type II; see Figure 13(b), e.g.). Owing to the longer crack length, the higher fracture energy is dissipated in Type II, and therefore the postpeak stress-displacement curve of Type II drops less quickly than that of Type I, which can be found in Figure 12.Figure 12
Stress-displacement curves for (a)AF=30%; (b) AF=40%; (c) AF=50%.
(a)
(b)
(c)Figure 13
Fracture patterns of concrete specimens withAF=30% for (a) Av=0%; (b) Av=2%; (c) Av=4%; (d) Av=6%.
(a)
(b)
(c)
(d)Figure 14
Fracture patterns of concrete specimens withAF=40% for (a) Av=0%; (b) Av=2%; (c) Av=4%; (d) Av=6%.
(a)
(b)
(c)
(d)Figure 15
Fracture patterns of concrete specimens withAF=50% for (a) Av=0%; (b) Av=2%; (c) Av=4%; (d) Av=6%.
(a)
(b)
(c)
(d)


## 4.1. Numerical Specimens and Mechanical Properties

Aiming to investigate the effects of voids on the tensile fracturing of concrete with different aggregate volume fractions, three sets of numerical concrete specimens with dimensions of 100 mm × 100 mm using polygonal aggregates are generated, and each set contains four specimens with the same aggregate arrangement and differentAv (0%, 2%, 4%, and 6%, resp.). For Set I, AF is set to 30%, while AF=40% and 50% are chosen for Sets II and III, respectively. The typical element size is chosen as 0.4 times that of the minimum aggregate size; however the number of both elements and nodes of numerical specimens increases with the increase of Av. For example, the number of elements and nodes for the specimen with AF=40% and Av=2% is 15796 and 9136, respectively, and the number of elements and nodes for the specimen with AF=40% and Av=6% is 20528 and 12126, respectively.For each specimen, uniaxial tensile fracturing is simulated. In FE simulations, the left end of concrete specimen is fixed in the horizontal direction, while the opposite end is subjected to a uniformly distributed horizontal displacement up to 0.06 mm, namely, a displacement-controlled loading scheme is used. Following the strategy discussed in Section3.4, the loading time is set to 0.012 s, which corresponds to a loading rate 5 mm/s.The same mechanical properties of aggregate, mortar, and ITZ are adopted for all specimens, as listed in Table1. It is noted that the mechanical properties of mortar are directly obtained from Chinese code GB 50010-2002 (code for design of concrete structures), and the compression hardening curve and the tension softening curve are shown in Figures 5 and 6, respectively. Moreover, due to the lack of detailed experimental results for ITZ, the compressive and tensile strengths and elastic modulus of ITZ are assumed to be 75% of those of mortar for the sake of modeling the weaker ITZ, while the other mechanical properties of ITZ are taken as the same of mortar.Table 1
Mechanical properties of concrete components.
Material
E
0 (GPa)
Poisson’s ratio (—)
σ
t
0 (MPa)
ϵ (—)
ψ (°)
α (—)
γ (—)
Aggregate
60
0.2
—
—
—
—
—
Mortar
20
0.2
1.94
0.1
35
0.12
2.0
ITZ
15
0.2
1.46
0.1
35
0.12
2.0
Void
—
—
—
—
—
—
—Figure 5
Total stress-strain relation under uniaxial compression and compression hardening curve of mortar.Figure 6
Total stress-strain relation under uniaxial tension and tension softening curve of mortar.


## 4.2. Effects of Voids on Concrete Tensile Fracturing Mechanism

In order to study the effects of voids on concrete tensile fracturing mechanism, the simulation results of two specimens in Set II, which include one without void and the other one with voids taking up2% of the specimen volume, are discussed in detail in the following.The macroscopic stress versus displacement curve of the concrete specimen without void (AF=40%) under uniaxial tension is shown in Figure 7, and the elements with nonzero equivalent tensile plastic strain at different loading levels are illustrated in Figure 8, in which the elements with equivalent tensile plastic strain bigger than 100 microstrains are highlighted in red color.Figure 7
Stress-displacement curve (AF=40%, without void).Figure 8
Mesoscale cracking development at different loading levels of the concrete specimen without void(AF=40%). (a), (b), (c), (d), (e), and (f) show the equivalent tensile plastic strain distributions at points a–f, respectively.
(a)
(b)
(c)
(d)
(e)
(f)It can be observed that although concrete specimen initially exhibits elastic response on the macroscale, mesoscale cracking still occurs in ITZs (see Figure8(a)) due to their lower fracture properties than mortar. With the increase of applied displacement, cracking develops in ITZs and subsequently propagates into mortar. The mesoscopic crack at the peak stress is shown in Figure 8(b), and it can be found that a macroscale crack close to the right end is formed as a result of the coalescence of mesoscopic cracks. In the softening stage, strain localization, which is accompanied by the decrease of macroscopic stress (see Figure 7) and finally leads to macroscopic tensile failure, can be clearly identified from Figures 8(c), 8(d), 8(e), and 8(f).The macroscopic stress versus displacement curve of the concrete specimen (AF=40% and Av=2%) under uniaxial tension is shown in Figure 9, and the elements with nonzero equivalent tensile plastic strain at different loading levels are illustrated in Figure 10. Compared to the case without void discussed above, mesoscale cracking appears in both ITZs and the mortar around voids in the macroscale elastic stage (see Figure 10(a)), and, additionally, cracked mortar generally suffers from bigger equivalent tensile plastic strain (see Figure 10(b)) at the peak stress as a result of the existence of voids. Similar to the specimen without void, a macroscale crack is formed at the peak stress point due to the coalescence of mesoscopic cracks. However, the fracture pattern is quite different from the one observed in the specimen without void even if the same aggregate arrangement is employed in both specimens, which indicates the existence of void dominates the fracturing behavior of concrete to a large extent. The phenomenon of strain localization, which is accompanied by the decrease of macroscopic stress (see Figure 9) and leads to the final tensile failure, is visualized in Figures 10(c), 10(d), 10(e), and 10(f).Figure 9
Stress-displacement curve (AF=40%, Av=2%).Figure 10
Mesoscale cracking development at different loading levels of concrete specimen withAF=40% and Av=2%. (a), (b), (c), (d), (e), and (f) show the equivalent tensile plastic strain distributions at points a–f, respectively.
(a)
(b)
(c)
(d)
(e)
(f)Overall, different fracturing mechanisms can be observed for the two concrete specimens by comparing the development processes of mesoscale cracking shown in Figures8 and 10. For the concrete specimen without void, mesoscale crack first appears in ITZ and then propagates into the neighboring mortar, while, for the concrete specimen with voids, mesoscale cracks are first found in both ITZ and mortar around voids and then coalesce in the fracturing process.


## 4.3. Effects of Void Volume Fraction

In this section, the effects of aggregate volume fraction (AF) on concrete tensile fracturing are first investigated, followed by the detailed discussion on the effects of void volume fraction on concrete with different aggregate volume fractions.The specimens without void in Sets I, II, and III are simulated, and the macroscopic stress versus displacement curves are depicted in Figure11. It is shown that the elastic responses of the three specimens are almost identical, which are in general linear and can be approximately characterized by the same elastic modulus. As AF increases, the reached peak tensile stress slightly decreases (similar to the results reported in [15]), which can be attributed to the fact that more ITZs exist in the case of higher AF. On the other hand, the postpeak stress in the case of AF=40% generally drops more quickly than in the case of AF=30%, as expected. However, in the case of AF=50%, the postpeak drops less quickly than in the other two cases, which may be due to the fact that this specimen finally fails with a different failure mode (two dominant macroscale cracks; see Figure 15(a)) from that of the other two specimens (one dominant macroscale crack; see Figures 13(a) and 14(a)), leading to higher fracture energy dissipated in tensile fracturing. From the above discussion, it can be concluded that the postpeak behavior of concrete highly depends on the failure mode, namely, the more macroscale cracks appear, the higher fracture energy should be expected.Figure 11
Stress-displacement curves of different aggregate volume fractions.Provided the simulation results of Sets I, II, and III, the effects of void volume fraction (Av) on concrete tensile fracturing are analyzed in the following. Figures 12(a), 12(b), and 12(c) present the stress-displacement relations of concrete specimens with AF = 30%, 40%, and 50%, respectively, and Figures 13, 14, and 15 illustrate the fracture patterns of concrete specimens with AF = 30%, 40%, and 50%, respectively. In Figures 12(a), 12(b), and 12(c), it is shown that the macroscale elastic modulus in prepeak stage can be considered to be independent of Av, while the peak stress decreases as Av increases. Moreover, the relation between Av and the postpeak behavior is not straightforward as depicted. As shown in Figures 13, 14, and 15, the fracture patterns of concrete specimens with the same aggregate arrangement (represented by highlighted red polygons) and different Av differ from each other. For concrete specimens without void, fracture pattern is mainly controlled by the distribution of aggregates (or ITZs), whereas, for concrete specimens with voids, the fracture pattern is dominated not only by the aggregate arrangement but also by the distribution of voids in mortar, leading to the differences of fracture patterns illustrated in Figures 13, 14, and 15. Consequently, accurate fracturing analysis of concrete calls for the explicit modeling of voids. Furthermore, both the concrete specimen without void and the concrete specimens with voids present two types of failure modes: one is with a single macroscale crack (Type I, typically shown in Figure 15(b)), and the other is with two macroscale cracks (Type II; see Figure 13(b), e.g.). Owing to the longer crack length, the higher fracture energy is dissipated in Type II, and therefore the postpeak stress-displacement curve of Type II drops less quickly than that of Type I, which can be found in Figure 12.Figure 12
Stress-displacement curves for (a)AF=30%; (b) AF=40%; (c) AF=50%.
(a)
(b)
(c)Figure 13
Fracture patterns of concrete specimens withAF=30% for (a) Av=0%; (b) Av=2%; (c) Av=4%; (d) Av=6%.
(a)
(b)
(c)
(d)Figure 14
Fracture patterns of concrete specimens withAF=40% for (a) Av=0%; (b) Av=2%; (c) Av=4%; (d) Av=6%.
(a)
(b)
(c)
(d)Figure 15
Fracture patterns of concrete specimens withAF=50% for (a) Av=0%; (b) Av=2%; (c) Av=4%; (d) Av=6%.
(a)
(b)
(c)
(d)


## 5. Conclusions

A finite element modeling strategy of concrete with random mesostructure explicitly taking void into consideration has been proposed in the present work. The tensile fracturing mechanism of concrete with voids is detailed on the mesoscale by comparing the simulation results of two specimens consisting of one without void and the other one with voids with the same aggregate arrangement. Then, several simulations are carried out with prime attention placed on the effects of void volume fraction on concrete tensile fracturing. The main conclusions are as follows:(i)
Different fracturing mechanisms are observed for the two concrete specimens with the same aggregate arrangement including one without void and the other one with voids, and the fracture pattern of concrete specimen with voids is controlled by both the aggregate arrangement and the distribution of voids.(ii)
Compared to aggregate volume fraction, void volume fraction has a larger influence on concrete tension strength. The elastic modulus of concrete in the prepeak stage can be considered to be independent of both aggregate volume fraction and void volume fraction.(iii)
The relation between concrete postpeak behavior and void volume fraction is not straightforward due to the randomness of void distribution and the resulting fracture pattern.(iv)
Two types of failure modes are identified for concrete specimens under uniaxial tension, in which Type I is characterized by a single macroscale crack and Type II by two. Due to the longer crack length, the postpeak stress of Type II drops less quickly than that of Type I.(v)
It is necessary to model void explicitly for the accurate fracturing simulation of concrete on the mesoscale.



---

*Source: 7989346-2017-04-24.xml*