A 3D FEM Mesoscale Numerical Analysis of Concrete Tensile Strength Behaviour

A three-dimensional (3D) finite element method (FEM) based on an inserted cohesive element numerical analysis procedure was developed for concrete mesoscale systems on the ABAQUS platform with python scripts. Aggregates were generated based on dividing the existing geometrical element algorithms to randomize arbitrary spheres. Simultaneously, randomizations of the maximum aggregate size and uniform distributions of aggregate particles were also considered. An FEM for the mortar phase in concrete mesoscale systems was generated along with the interfacial transition zone (ITZ) by inserting a cohesive element. Numerical parameter analyses were performed for nine different concrete systems by varying the coarse aggregate volume fraction (α) and the ITZ tension strength (ITZ-S). The mechanical performance of concrete systems with the coupling effects of α and ITZ-S was evaluated for simulated tensile loading. The results of the numerical simulations for mechanical properties, such as the simulated tensile strengths and tension damage behaviour of concrete systems, were verified with experimental results. The proposed aggregate and ITZ generation approach and numerical simulation procedure can be used by researchers to better understand how aggregate volume fraction and ITZ strength affect the tensile behaviour of concrete mesoscale systems.


Introduction
At mesoscale, concrete is composed of three or four phases, namely, mortar with cementation, aggregate with a skeleton or filling effect, interfacial transition zone (ITZ) connecting mortar and aggregate, and random defects [1]. e ITZ forms a honeycomb wall between aggregate and mortar. Due to adhesion, the mechanical properties of the ITZ have a substantial effect on the elasticity and strength of concrete [2,3]. However, at the initial stage of pouring, the fluidity of mortar near the surface of the aggregate is weakened, and the chemical potential energy caused by different ion concentrations results in higher porosity and lower strength of the ITZ than bulk cement paste [4,5]. is makes the ITZ the cause of dry shrinkage of concrete, creep, and even damage due to cracking [6][7][8].
Also, with the rapid development of high-performance, lightweight, and recycled aggregate concrete (RAC), higher requirements are placed on the research of the role of the ITZ in concrete [9][10][11]. In particular, coupling new and old ITZs in RAC results in complex mechanical properties of concrete [12][13][14][15]. To enhance the mechanical properties of RAC, some researchers have added polypropylene, basalt, steel fiber, and other materials, which further makes concrete a multiphase composite material [16][17][18][19].
Several researchers have tested the physical and mechanical properties of the ITZ with a large number of laboratory micro tests. Among these, tests based on nanoindentation technology can directly obtain the elastic modulus of the ITZ on a mesoscale [20][21][22][23]. Scanning electron microscopy (SEM) can identify the two-dimensional boundary of the ITZ [24], and X-ray diffraction (XRD) can examine the chemical composition of hardened cement paste which helps to explain the macroscopic behavior of concrete [25]. Digital image correlation (DIC) can track and reproduce ITZ failures and cracking processes in real time [26]. ese methods provide a solid physical foundation for the theoretical analysis and numerical simulation of the ITZ.
A substantial amount of theoretical research on multiscale mechanics and micromechanics that satisfy the law of conservation in thermodynamics has been conducted for concrete containing an ITZ. At the initial stage, the volume fraction of the ITZ in concrete was calculated by a representative volume element (RVE) based on the thickness, after which the RVE was converted to the corresponding volume mortar [27]. ereafter, in-depth theoretical threephase analyses of concrete based on the Mori-Tanaka method for the Eshelby inclusion problem were conducted [28]. However, only simple boundary conditions have been used in theoretical research; besides, the interactions between multiple phases have been ignored. erefore, a finite element method (FEM) based on an inserted cohesive element numerical analysis procedure is now being used to conduct further research.
Over the past several decades, several 2D and 3D numerical studies have focused on the performance of the ITZ on mesoscale. Wang et al. [29] used a stochastic finite element model (SFEM) to evaluate the effect of creep on the behavior of concrete. Zhou and Hao [30] conducted indoor experiments and also used an SFEM to evaluate the elastic modulus of the ITZ for concrete, including recycled concrete. Du et al. analyzed the diffusion characteristics of chloride ions in the ITZ by SEM and SFEM [31,32]. Lu et al. investigated the effects of acid rain solutions and mineral admixtures on the compressive strength and elastic modulus of RAC [33]. Bernard and Kamali-Bernard et al. [2] examined the mesoscale behaviors of mortar and concrete. Gong et al. [34] developed a mesoscale model based on the rigid-body spring method to simulate the nonlinear mechanical behavior of concrete to investigate damages due to freeze-thaw cycles. e spatial effect of the ITZ is ignored in the RVE which equates the ITZ to mortar. Although 2D SFEM partially reflects the spatial effect, its essence is to consider the ITZ as a cylinder, which increases the effect of the ITZ on concrete, resulting in a significant difference between the calculated results of the model and test results. Although 3D SFEM accurately reflects the spatial effect of the ITZ [35], most existing numerical analyses are based on 2D SFEM. While some 3D multiscale computational models of concrete platforms have used CEMHYD3D and ABAQUS to conduct 3D FEM analysis of the ITZ and mortar, the element volumes were too large due to the limitations of the algorithms [3,5,8]. For typical thin-shell structures of the ITZ, if 3D solid elements are used and the number of elements is not increased substantially, then ill-conditioned meshes are inevitable. e thickness of the ITZ is usually 0.01 mm to 0.1 mm, and the coarse aggregate diameter is about 10 mm [36][37][38]. If solid elements are used to simulate all phases within the same model, then the number of unilateral elements is at least 10 mm/0.01 mm � 10 4 , and the number of three-dimensional elements is above 10 12 , which is unacceptable for FEM analysis. Although Kamali-Bernard [4] used connecting elements to simulate the ITZ, the elements still need to be equal to the mortar in size, which results in the model analysis being unable to fully consider the role of the aggregate.
To address this issue, the present study first used the Ran_gen_agg_3D_ABAQUS program developed using python in an ABAQUS 3D 100 mm × 100 mm × 100 mm characteristic FEM model, according to the Fuller grading curve randomly picked aggregate. Subsequently, with the COH_3D_ABAQUS program developed using python, a cohesive element without thickness was inserted between the randomly generated aggregate and the mortar to represent the ITZ. e traction-separation cracking criterion was applied to the cohesive element so that the three-phase interaction between the aggregate, mortar, and the ITZ in concrete could be expressed scientifically and reasonably in a multiscale form. On this basis, a parametric analysis of the effects of the aggregate volume fraction and the ITZ strength behavior on the elastic and tensile strength of concrete was conducted to verify the reliability of the FEM based on an inserted cohesive element numerical analysis. Figure 1 shows the 3D FEM mesoscale principle on a representative concrete cube unit. e side length of each cube unit was 100 mm and was divided into 40 equidistant units, each of which was 2.5 mm in length. Each cube unit consisted of three phases as indicated in Figure 1: (a) aggregate, (b) mortar, and (c) ITZ. e aggregate played a skeletal role in the strength of the unit, and its gradation and particle size determined the primary mechanical characteristics of the multiscale system. Mortar was used as a binder for the aggregate, restrained the freedom of aggregate, and improved the skeletal function of the aggregate on the overall strength of the unit. Being the transitional layer between the aggregate and the mortar, the ITZ facilitated the interaction between the mortar and the aggregate.

Generation of Aggregates.
e Ran_gen_agg_3D_A-BAQUS program in python was used to simulate randomly distributed aggregates in concrete units. e program is mainly based on the Fuller grading curve [31] but with modifications, as shown in Figure 2.
Referring to Jayasuriya's algorithm [12], the aggregate generation process was as follows: First, the geometry of a concrete system was generated for the aggregate characteristics obtained from the ABAQUS model. Next, a uniform 3D cube mesh was generated in the ABAQUS model. Finally, the element sets that satisfy the Fuller grading curve in the generated mesh were selected using coded scripts (Ran_-gen_agg_3D_ABAQUS) written in python.
As stated by Wang et al. [39,40] and Jayasuriya et al. [11], both fine and coarse aggregates occupy between 60 and 80% of the total volume of concrete, and 40-60% of the volume fraction is occupied by coarse aggregates. e aggregate gradation feature is one of the key values that determines the mechanical properties of concrete. For a reasonable particlesize distribution, an appropriate gradation curve needs to be generated. In this paper, the aggregate particle distribution in the mesoscale system was based on the Fuller curve, as shown in the following equation: where P(d s ) is the cumulative ratio for a sieve with an aperture d s ; d max is the maximum aggregate diameter in the aggregate sets; and n is the gradation exponent that varies from 0.33 to 0.50 depending on the coarseness of the gradation. For a perfect maximum aggregate density in concrete, a gradation exponent between 0.40 and 0.45 was achieved. erefore, 0.40 was used as the Fuller curve exponent in this study. e representative coarse aggregate gradations for three levels, 4.75 mm, 9.25 mm, and 16 mm, are shown in Figure 2. e 3D mesoscale concrete systems were randomly generated, and aggregate particles were considered as approximate spheres such that each aggregate had a specific volume within a concrete boundary. e geometry of the concrete 3D cube, that is, width (W), height (H), and length (L), can be defined by a python script for ABAQUS, and aggregates were randomly generated within this 3D space. Representative aggregate volumes that are passed through two consecutive sieve sizes are defined in the following equation: where V agg [d (s+1) , d s ] is the volume of the aggregate within d s+1 and d s ; d max and d min are the maximum and minimum aggregate sizes, respectively; α is the aggregate volume fraction of a concrete system being generated; and V c is the total concrete volume (W × H × L).
Of these values, α is a key parameter that controls the aggregate generation ratio. e gradation range was discretized into consecutively smaller sieve aperture sizes, such that a target aggregate volume was obtained for each pair of successive sieve sizes. e discretization of each sieve aperture size pair over the selected gradation is significant when the aggregate packing density of the concrete structure is considered.

Aggregate Generation Algorithm.
In general, current aggregate generation algorithms first randomly divide multiunconnected subspaces in Euclidean space according to a certain probability distribution and then generate meshes based on these subspaces. ese kinds of aggregate generation Advances in Materials Science and Engineering algorithms limit the shape of the grid by the aggregate subspace and do not generate uniform cubical elements in 3D or generate uniform pixel elements. e regions need to be discretized very densely, which leads to an increase in the amounts of calculations in a geometric series, some of which are impossible to execute. e algorithm in this paper is entirely different from the above methods. First, a group of uniform pixel elements were generated in Euclidean subspaces, and then a group of virtual multi-unconnected subspaces was generated according to the set's probability distribution function. e generated element sets were used as random aggregate sets, to avoid the element singularity problem caused by the adaptation of the grid to the boundaries of the random subspaces. e detailed method is described below.
In this study, a 3D cube with 100 mm sides as the Euclidean space was divided into 40 grids on each side. e virtual multi-unconnected subspaces were all spherical with diameters satisfied by equation (1) and volume fractions satisfied by equation (2). In Euclidean space, the virtual subspace representing aggregates was mapped randomly in stages. First, the spherical subspace representing large aggregates was mapped, and then the spherical subspace with small aggregates was mapped, step by step among the large aggregates. e specific process is as follows: first, a spherical space with a diameter of 16 mm was mapped into the cube with a side length of 100 mm. Since the aggregates cannot overlap each other, the distance between the subspaces must be slightly larger than the sum of the radii of adjacent spherical spaces. In this study, a subspace size of 1.1 (d s+1 + d s )/2 has been considered. Next, the spherical subspaces with diameters of 9.5 mm and 4.75 mm were mapped step by step based on the above principles. Since this research simulated the concrete cube unit extracted from the Euclidean space instead of the self-made concrete experimental blocks, the aggregate could pass through the boundary of the cube, which resulted in a portion of the aggregate not being completely spherical, as shown in Figure 1. However, this mapping strategy was more representative and more conforming to the actual state of stress of the concrete, especially for tension, and the general failure modes of concrete include cracking at the outer edges of large aggregates. If the aggregates were not allowed to cross the boundary, cohesive elements could only be calculated inside the unit, which would have been a partial deviation from the actual simulation. In the process of mapping, a set of three uniformly distributed random number generations was adopted to locate the coordinates of the centres. e process that was used to generate the subspaces is explained in a flowchart shown in Figure 3.
e Ran_gen_agg_3D_ABAQUS program based on the above algorithm using a python script generated a seed size � 2.5 mm, with 40 units on each side of the concrete cube, in an ABAQUS CAE 100 mm × 100 mm × 100 mm model space, as shown in Figure 3(a).
Virtual multi-unconnected subspaces mapped according to equations (1) and (2) are represented by different colours in Figure 3

Figures 3(e), 3(g)
, and 3(f ), respectively. Based on the above program, by the order of the diameter from large to small, the first of the large aggregate sets occupied positions in the main quadrants and ensured large gaps, while the secondary generated medium aggregate sets and fine aggregate sets gradually filled the gaps, causing the skeletal function of the aggregate to increase in strength.

Model of the ITZ
e method to simulate the ITZ that surrounds each coarse aggregate uses the particle sets previously generated. It is primarily based on the cohesive element technique of the ABAQUS [41,42] but with modifications. e cohesive element used in this research was an 8-node 3D cohesive element (COH3D8) as shown in Figure 4. Although ABAQUS CAE can insert cohesive elements between regular geometries, it is impossible to insert cohesive elements without thickness between aggregates and mortar, for a large number of irregular faces as shown in Figure 3. Hence, it was necessary to use the python script driver COH_3D_ABAQUS developed by us to re-edit the ABA-QUS.inp files.
Generate next "n" random points in mesh surplus spaces.
Map "n" subspaces in mesh space on the centres as "s" aggregates.

Traction-Separation Law of the ITZ.
Under tension, the cracking process for the cohesive elements of the ITZ is frequently in the mixed-mode I + II as shown in Figures 6  and 7.
δ 0 n , δ 0 s , and δ 0 t are the initial displacements in the normal and shear directions, respectively, when the ITZ is damaged; δ f n , δ f s , and δ f t are the displacements at eventual failure in the normal and shear directions, respectively; K is the stiffness of traction-separation; and G C n and G C s are the mode I or II interfacial fracture energies.

Damage Initiation.
e traction-separation constitutive equation of the ITZ represents the damage mechanism. e initial damage condition refers to the first degradation of the response of a material point.
e process of degradation starts as the combined stresses satisfy the quadratic nominal stress criterion, which can be defined as the composite stress ratio reaching a value of one: <> is the Macaulay bracket and guarantees that a purely compressive stress state does not initiate damage; t 0 n , t 0 s , and t 0 t are the peak stress values for deformation purely normal   to the interface, in the s shear direction and the t shear direction, respectively, as shown in Figure 6.

Damage Evolution.
e characteristics of damage evolution are represented by the rate at which the ITZ stiffness declines once the corresponding initiation criterion is reached. Figure 7 shows the scalar damage variable D, which represents the damage factor at the ITZ. Damage evolution is based on the energy that is dissipated as the result of the fracture. e fracture energy G C n /G C s is equal to the area under the traction-separation curve, as shown in Figure 7. In this study, the Benzeggagh fracture criterion [41] was used. e critical fracture energies during deformation purely along the s and t shear directions are the same for aggregates (i.e., G C s � G C t ). is is given by where G S � G S + G t , G T � G n + G S , and η is a material parameter. According to Nasiri and Liu [42], G C n � 40 N/m, G C s � 400 N/m, and η � 1.5. e linear softening damage variable D reduces to where δ f � 2G C /T 0 eff , T 0 eff is the effective traction at damage initiation, and δ max is the maximum effective displacement during loading. e bond strength, fracture energy, and corresponding strain of the ITZ were determined from experimental results in the literature [43,44].

Mechanical Properties and
Constitutive Models 4.1. Aggregate. e tensile properties of aggregates used in this analysis were based on Brazilian tensile tests [14,49]. e parameters of the aggregate are presented in Table 1. Table 1 presents elastic modulus E m , Poisson's ratio μ m , fracture energy G m , and tensile strength τ m of mortar element sets in a mesoscale concrete system under uniaxial tensile stress. e concrete damage plasticity (CDP) equation was adopted to calculate the response of the mortar phase to uniaxial loading with tension σ t and corresponding strain ε t . e relations between σ t and ε t are shown in Figure 4 and show the plastic behavior of the mortar in concrete under tension [50].

Mortar.
e CDP model assumes that the two main failure mechanisms are tensile cracking and compressive crushing of concrete. However, only the tensile cracking property was investigated in this study. e evolution of the failure surface is controlled by the hardening variable ε pl t (tensile equivalent plastic strain) linked to the failure mechanism under tensile loading. Under uniaxial tension, the stress-strain response follows a linear elastic relationship until the value of the failure stress, σ t0 , is reached. e failure stress corresponds to the onset of microcracking in the mortar phase. Beyond the failure stress, the formation of microcracks is represented macroscopically with a softening stress-strain response, which induces strain localization in the mortar phase. As shown in Figure 8, when the mortar phase is unloaded at any point on the strain-softening branch of the stress-strain curves, the unloading response is weakened and the elastic stiffness of the material appears to be reduced. e degradation of the elastic stiffness is characterized by the damage variable, d t , which is assumed to be a function of the plastic strain. If E 0 is the initial (undamaged) elastic stiffness of the material, the stress-strain relation under uniaxial tension is given by e "effective" tensile stress is e effective cohesion stresses determine the size of the failure surface. e plastic flow potential function and the yield surface make use of two stress invariants of the effective stress tensor, namely, the hydrostatic pressure stress, and the von Mises equivalent effective stress, where S is the effective stress deviator, defined as S � σ + PI. (10) e CDP model assumes nonassociated potential plastic flow. e flow potential G used for this model is where ψ is the dilation angle measured in the p-q plane at a high confining pressure and ς is a parameter referred to as the eccentricity, that defines the rate at which the function approaches the asymptote.

General Conditions.
e ABAQUS software was employed for FEM with inserted cohesive element numerical analysis. e aggregate, mortar, and ITZ were all simulated as 8-node hexahedral elements with three translational degrees of freedom (DoFs) at each node. For the mechanical analysis, the C3D8 element type was adopted for the aggregate and the mortar, and the COH3D8 cohesive element type was adopted for the ITZ.
In ABAQUS, the cohesive element (COH3D8) simulates the ITZ in detail, including the properties of tractionseparation, by using an adhesive layer with zero thickness in Euclidean space. e mortar and the ITZ around the aggregates are the weak locations of the concrete system. erefore, the tensile strengths of the ITZ and the mortar are important in determining the tensile strength of concrete. If the stress is greater than the tensile strength, the concrete cracks. e mesh design was optimized to determine the appropriate mesh density to obtain reliable data. e FEM based on an inserted cohesive element mesoscale algorithm used for this research involves large-scale parallel computation. Hence, ABAQUS/explicit multiple processors for the algorithm were adopted, and 48 core Intel® Xeon® Gold 624br CPU @ 3.00 GHz was used for the parallel computation. e geometry involved a 100 mm × 100 mm × 100 mm cube, and the model with 40 elements on each side required more than 72 hours to complete the analysis.

Tensile Stress-Strain Relationship for Concrete.
In the ABAQUS mesoscale calculation platform, the degrees of freedom of the upper and lower nodes were coupled to the reference point of the center of the surface, and all degrees of freedom of the lower surface reference point were constrained. Except for the vertical displacement of the upper surface, all other degrees of freedom were constrained. e loading was realized by forcing the upper surface coupling reference point displacement.
Based on mesoscale numerical simulations, the tensile strength-strain relationships of concrete systems under uniaxial tensile load are shown in Figure 9 for different volume fractions and ITZ strength. e tensile stress-strain relationships of concrete systems showed a linear material behavior up to the peak response; thereafter, they experienced a slight parabolic softening phase and failed at the corresponding ultimate strain capacities, for which the aggregate volume fraction α < 0.6, and the ITZ tensile strength is within the normal range. However, when the ITZ strength was particularly low, such as for the extreme case assumed in this study, which is 10% of the mortar strength value, the stress-strain linear relationship no longer existed, and strain softening commenced after the initial stage of the loading process.
Since this analysis focused on the effect of the ITZ tensile strength (ITZ-S) on the tensile strength of concrete, the maximum tensile strength of mortar for all specimens in the analysis was 2.95 MPa.
It can be seen from Figure 9(a) that when α � 0.4, and ITZ-S � 0.29 MPa, the stress-strain relationship of the concrete specimen weakened for an average strain greater than 0.2 × 10 −4 , and for a maximum average tensile strength of 1.97 MPa, which is 33% lower than the maximum tensile strength of 2.95 MPa of mortar, while the residual tensile strength was 0.66 MPa. e value of the strain corresponding to the residual tensile strength was 0.8 × 10 −4 , which is close to the yield stress of 0.58 MPa when the tensile strain was 0.8 × 10 −4 in the CDP constitutive model of mortar, and the decrease in strength was 1.31 MPa. It shows that when the ITZ-S is 10% of the mortar strength, the ITZ enters an open  Figure 8: Stress-strain curve. 8 Advances in Materials Science and Engineering state before the mortar cracks. When ITZ-S was 1.25 and 2.21 MPa, the maximum average tensile strength of the specimen was 2.18 MPa, which is 26% lower than that of mortar (2.95 MPa), and the residual strength was close to the above residual strength, which indicates that the ITZ-S is not significantly lower than the mortar strength. us, even if the ITZ-S is 50% of the mortar strength, it does not cause the tensile strength of the mesoscale concrete system to decrease significantly. When α � 0.5, as shown in Figure 9(b), the relationship of the stress-strain curve is consistent with that of α � 0.4, but for ITZ-S � 0.29 MPa, the maximum average tensile strength of the specimen was 1.92 MPa, which is 35% lower than that of mortar (2.95 MPa), and the residual strength of 0.905 MPa was 26% higher than that of 0.66 MPa for α � 0.4, which indicates that the effect of ITZ-S becomes increasingly significant with an increase in α.
When α � 0.6, as shown in Figure 9(c), the relationship of the stress-strain curve was different from that of α � 0.4 and 0.5. For the specimen with ITZ-S � 2.21 MPa, there was a large decrease in stress from 2.23 MPa to 0.0027 MPa, and the peak value was 1.4% lower than that of 2.26 MPa, for ITZ-S � 1.25 MPa. For the specimen with ITZ-S � 0.29 MPa, there were multiple yield steps, the peak value was 1.51 MPa, which is 48.8% lower than the tensile strength of mortar (2.95 MPa), while the residual strength was 0.849 MPa. is shows that the ITZ-S has a decisive effect on the strength of the entire mesoscale concrete system after α is increased to 0.6.
Since concrete is a heterogeneous material and the failure process inevitably involves geometrical uncertainties, Advances in Materials Science and Engineering three trials for each case were calculated to increase the reliability of the simulation results, and the peak values of strength with deviations under the above conditions are summarized in Figure 9(d). For different values of α, the characteristics of the inflection point of the peak of the tensile strength were consistent, and the effect of ITZ-S on the tensile strength of concrete was gradually enhanced with an increase in α.
For a single aggregate, the connection strength of mortar and aggregate showed a linear correlation with a decrease in ITZ-S. However, for the entire mesoscale concrete system, this linear correlation was not tenable. With an increase in α, the specific surface area of aggregate in mortar increased significantly in a nonlinear manner. is relationship was also significantly affected by the particle gradation of the aggregate. For finer aggregates, the surface area increased, based on the simulation results of this research. ese figures represent the DAMAGET and SDEG distributions at the corresponding peak response of each mesoscale concrete system that was simulated. e DAM-AGET and SDEG distribution contour limit values were from 0 to 1, considering the maximum tension damage capacity of concrete. Both DAMAGET and SDEG had a coupling effect on the mesoscale concrete system especially near the peak response because the concrete failure type is brittle under tension.

Tensile Damage
As is shown in Figure 10, for uniaxial tension, mesoscale concrete systems perform quite satisfactorily since the DAMAGET and SDEG have already reached the peak value of 1.0. erefore, it is confirmed that the ITZ and the mortar are the two material phases subjected to higher strain levels.
As shown in Figures 10(a), 10(c), and 10(e), the corresponding uniaxial tensile stresses of mortar reached ultimate strength, and the ITZ fracturing was initiated even before ultimate crushing occurred. is was identified as the material having undergone a softening phase where it can no longer resist the applied strain deformation due to the drastic change of the overall material stiffness. In the mesoscale system, a main cracking seam occurred in the mortar phase along the large aggregate, with a group of small cracks around it. Moreover, with an increase in the value of α, as seen by comparing Figures 10(a) and 10(e), the number of small cracks increased. e reason is that with an increase in α, the effect of ITZ-S on the tensile strength of the entire concrete system is more significant than that of mortar, and the number of ITZ cracks also increases. After the ITZ cracks, the tensile stress is rapidly transferred to the adjacent mortar, resulting in the formation of more small cracks. is phenomenon also confirms the reason for the several steps shown in Figure 9(c).
During this period, based on the mesoscopic simulations performed, it is inferred that microcracks were initiated at the ITZ and eventually led to failure due to interfacial cracking, while the material response continued to proceed beyond the peak level as shown in Figures 10(b), 10(d), and 10(f ). A similar material behavior was observed in other experimental findings that include X-ray, Brazilian, and Hopkinson pressure bar techniques [51][52][53][54][55], which further confirms that the results of the current numerical simulation are realistic.
In addition to the effect of α and ITZ-S on the tensile strength of concrete, the sequential order of the cracking pattern is dependent on the geometric size of an aggregate. It varies from large-to small-diameter aggregates that pertain to an ITZ along the loading direction, which was verified using a frame-by-frame hybrid analysis. Alternations in the strain transfer mechanism were attributed to strain localization near the intersection of the ITZ phases. Rapid tensile strain localization was observed at the locations where the aggregates connecting the plane were aligned in the direction of the loading plane; the ITZ formed a weak plane, causing accelerated damage progression. us, more tensile strain deformations accumulated in these areas, and potential interfacial cracking failures were aggravated. e three mesoscale concrete systems at α � 0.4, 0.5, and 0.6 showed that larger values of α tend to induce higher tensile strain deformations (i.e., cracking potential) as depicted by the simulations in Figures 10(e) and 10(f ). As shown in Figures 9(a)-9(c), the higher the ITZ-S, the more significant the change in tensile strain. is was due to the higher effective curvature of the ITZ cracking path associated with larger aggregate sizes. is result is similar to [56]. It leads to poor performance of the concrete system which showed low tensile strength. Further, it can be seen from the simulations in Figures 10(b), 10(d), and 10(f ) that when the concrete system reached its peak strength, the principal tensile strain created a clear cracking surface on the aggregate due to the smooth formation and interconnectedness of tensile strain contours around round aggregates.
is explains the result that the mechanical performance of concrete systems with larger, round aggregate particles and low adhered mortar content levels is highly vulnerable to getting damaged under applied tension loading.

Conclusions
A 3D FEM mesoscale model developed on the ABAQUS platform to generate random aggregates and insert ITZ cohesive elements for concrete systems has been presented in this paper. A python-based ABAQUS script was developed to divide the aggregate regions in the Euclidean space and insert cohesive elements as the ITZ around the boundary of the aggregate regions. Corresponding tensile stress-strain relationships were studied for each concrete system. e following conclusions were drawn from the numerical analyses performed with the mesoscale concrete systems: (1) Both the aggregate volume fraction α and the ITZ strength have important effects on the tensile strength of concrete, of which the aggregate volume fraction α determines the ITZ area. Hence, α is the main factor determining the tensile strength of concrete.
(2) Higher ITZ strength levels contribute to more stable strength around the coarse aggregates during tension loading processes, which facilitate concrete systems to reach higher strain capacities before softening.  (3) Extremely low ITZ strength, such as being 10% of the mortar strength, leads to the entire concrete mesoscale system strength becoming unstable; strain softening occurs at the initial stage of loading, and with an increase in the volume fraction of coarse aggregates, this phenomenon is significantly enhanced. (4) e general failure path due to uniaxial tensile deformation of concrete mesoscale systems initiates from ITZ traction-separation around large-diameter aggregates to small ones before tension damage occurs to the mortar. (5) e failure mode of a mesoscale concrete system is usually through the main crack at the edges of large aggregates, and there are a few small nonpenetrating cracks around the main crack. ese small cracks increase with volume fraction. When the ITZ strength is low, there are several steps in the stressstrain relationship curve.
e FEM based on an inserted cohesive element mesoscale method developed in this paper provides a reasonable approach for researchers to study the gradation and mechanical factors that affect the response of concrete mesoscale systems.

Data Availability
All data, models, and code included in this study are available upon request by contact with the corresponding author.

Conflicts of Interest
e authors declare that they have no conflicts of interest regarding the publication of this paper.