Failure process of steel – polypropylene hybrid ﬁ ber - reinforced concrete based on numerical simulations

: In this work, we studied the failure mechanism of steel – polypropylene hybrid ﬁ ber reinforced concrete ( HFRC ) at the mesolevel. The uniaxial tensile test of HFRC was simulated using ABAQUS ﬁ nite element ana lysis software. Then, the relationship between the meso scale failure process and the mechanical properties was analyzed based on the simulation results. The results showed that the cracks ﬁ rst appeared in the interfacial transition zone and then gradually propagated into the mortar elements and intersected with adjacent cracks, forming major macroscopic cracks. According to the crack evolution process, the incorporation of steel ﬁ bers and polypropylene ﬁ bers changed the concrete crack expansion paths and served to inhibit crack expansion. Furthermore, the increase in the hybrid ﬁ ber volume had a positive e ﬀ ect on the mechanical properties, and the steel ﬁ bers dominated in providing reinforcement compared to the polypropylene ﬁ bers.


Introduction
Concrete is currently the most widely used building material because it offers suitable compressive strength and durability along with convenient material acquisition and easy preparation, giving it certain advantages over other traditional building materials. However, concrete also has disadvantages such as low tensile strength and a propensity for cracking, which greatly limit its engineering applications. Over the past few decades, fiberreinforced concrete (FRC) has been extensively studied through model tests, numerical simulations, and theoretical studies [1][2][3][4][5][6]. Several studies have shown that the incorporation of fibers in concrete can significantly improve the mechanical properties of concrete, such as tensile strength, bending strength, impact strength, toughness, and ductility [7][8][9][10]. Moreover, during crack initiation and expansion, incorporated steel fibers and polymer fibers can change the crack formation mechanism, thereby providing better control of the size and development of crack defects [11][12][13]. In addition, steel-polypropylene hybrid fiber-reinforced concrete (HFRC), which incorporates steel fibers and polypropylene fibers, has been widely used in practical engineering due to its excellent mechanical properties, and it has broad application prospects [14].
Numerous studies have shown that the contribution of FRC is more significant during tension than during compression on the fibers [15,16]. The tensile properties of composites are substantially better because the incorporated fibers provide bridging stresses, inhibit crack opening, and limit crack expansion in the concrete matrix. In addition, to evaluate how fibers enhance the mechanical properties of concrete, numerous experimental studies in recent years have been conducted with the goal of determining the effects of different fiber types and the parameters of various loading scenarios [17][18][19][20]. However, large numbers of laboratory experiments are inevitably constrained and are sometimes challenging to conduct. Thankfully, due to the development of supercomputers and commercial software such as ANSYS, ABAQUS, and LS-DYNA, numerical simulation techniques have become essential tools for studying concrete, and researchers have started to study the failure processes of FRC at the mesolevel with better results. Wittmann et al. [21] first used a multi-scale method to propose randomly distributed polygons to simulate a mesoscopic mechanical model of concrete aggregates. However, the polygons used in this model were relatively simple, making it difficult to accurately describe the mesostructure of the concrete. Xu et al. [22] established a two-dimensional (2D) mesoscopic numerical model of hook-shaped and spiral-shaped steel fiber circular aggregate concrete and carried out compressive impact loading simulation and dynamic tensile test simulations. Zeranka [23] developed a constitutive material model and performed numerical calculations with user material in the commercial finite element software ABAQUS. Tran et al. [24] used a numerical model to highlight the failure mode with the transversal cracks on the specimen surface utilizing the crack damage model for the concrete matrix. Sun et al. [25] simulated the effect of the basalt fiber length and content on the fundamental mechanical properties of concrete and developed a damage ontology model based on Mori-Tanaka homogenization theory and progressive damage theory to predict the properties of basalt FRC composites. Du et al. [26] used the damaged plasticity theory combined with the strainrate effect to describe the dynamic mechanical behavior of a mortar matrix, and the aggregate phase was assumed to be elastic. Rezakhani et al. [27] performed the numerical modeling of the presented experiments using the Lattice discrete particle model enriched with the fibers effect. Leite et al. [28] performed experimental simulations of uniaxial compression and splitting damage to 2D and three-dimensional steel FRC. Considering the above literature, a research basis has been established for the mesoscale aspects of HFRC. Because the destruction process of FRC involves complex damage evolution mechanisms, if the crack initiation, propagation, and instability processes were to be revealed, then these processes would provide a theoretical basis for the design of high-performance concrete.
The material properties and shape differences of the constituent substances of concrete are difficult to model. Accurate representation of the natural features of the structure with a reasonable computational workload is difficult to achieve, and adding fibers to the concrete fine view structural model becomes more challenging [29][30][31]. So far, most models have been simplified in different ways without taking into account the key geometrical features of the fine view structure, and fewer models deal directly with the irregular fine view structure of fibrous concrete. For example, homogenization methods were used to reduce the degree of inhomogeneity of the original structure, and aggregates were modeled with simple geometries such as spheres and ellipsoids [32][33][34][35], while the irregular shape of aggregates has a large influence on macroscopic and microscopic cracking processes. There are also some scholars who performed numerical simulations of mixed-fiber concrete, generating only two-phase models consisting of fibers and concrete matrix [36,37]. The purpose of this study was to propose a reasonable HFRC model at the mesoscale level and investigate the connection between the mechanical properties and macroscopic mechanical behavior of the model under uniaxial tension. We also analyzed the crack opening expansion process in the model in terms of damage as well as the macroscopic mechanical properties of the mesoscale unit, and we investigated the crack expansion law of the FRC, along with the fracture mechanism problem. In practical engineering applications, artificially crushed coarse aggregate is mostly in the shape of irregular convex polyhedra, and on the 2D scale, convex polygons can be used to approximate simulated crushed stone coarse aggregate. When programming a convex polygonal aggregate algorithm, the prevailing generation method places round aggregates in the model that match the conditions based on the number of aggregates obtained from the grading calculations. Then, an inner jointed triangle is randomly generated in each circle as a substitute for the circular aggregate. Finally, the triangular aggregate is extended according to the topology, and the numbers and shapes of the extensions are controlled by the area until an eligible convex polygonal aggregate is generated [38]. Accordingly, in this work, we generated convex polygons directly in the model region to approximate the crushed stone aggregate by controlling the deformation of the diamond-shaped parametric equation in accordance with the expansion factor method proposed by Laizhong et al. [39]. The specific steps were as follows.
(1) First, we examined the 2D rhombic parametric equation: x y cos sin sin cos cos cos sin sin cos sin where O x y , T 0 0 ( ) is the central coordinate of the rhombus and ∈ α π 0, 2 [ ] is the angle between the vector − − x x y y , 0 0 ( ) and the x-axis. The diameter of the rhombus could be defined as = R r r 2 max , 1 2 {| | | |}, as shown in Figure 1. The hypothesis was The expansion transformation of the rhombic parameter equation (equation (1)) yielded the new convex polygon parameter equation, as follows As shown in Figure 2, a polygon was obtained utilizing the rhombus parametric equation under the control of multiple superimposed deformations.

Fiber-generation algorithm
In the 2D plane, three parameters were used to determine the fiber position: the coordinates for the midpoint of the fiber (x, y), the length of the fiber, and the angle of inclination α (expressed in radians) between the fiber direction and the positive x-axis direction. We then followed the steps below to complete the fiber placement process during model generation.
(1) The quantities of the steel fibers and polypropylene fibers were calculated based on the volumetric fiber ratio and equation (5), respectively. (2) According to the rand(x) function in the rectangular region of the model and based on the horizontal coordinates (x min , x max ) and vertical coordinates (y min , y max ), two random numbers (x i , y i ) and the inclination angle α i between the fiber length direction and the positive direction of the x-axis were randomly generated to determine the planar position of the fiber according to the following random array: where x min and x max are the minimum and maximum values of the model in the x direction, and y min and y max are the minimum and maximum values of the model in the y direction, respectively. (3) The fiber parameters had to meet the boundary conditions and could not intersect with the aggregate. For the boundary conditions, the generated fiber endpoints had to be within the model, according to the expression below In this work, the vector fork multiplication method was used to determine the intersections of the fibers with convex polygonal aggregates. The flow chart for the random generation of aggregates and fibers is shown in Figure 3.

Random generation of aggregates
In the early twentieth century, Fuller, an engineer, proposed the Fuller grading curve [40] to determine the volume fraction for each range of particle sizes. In this work, the fuller grading curve was selected as the ideal maximum density grading curve according to the following equation: where D 0 is the diameter of the sieve hole, D max is the maximum aggregate particle size, and P is the mass fraction of aggregate through the sieve. Walraven and Reinhardt [41] transformed Fuller's formula into a 2D aggregate gradation curve, which is known as Walraven's formula: where P k is the percentage of the volume of coarse aggregate and fine aggregate in the total volume of the concrete model, D 0 is the diameter of the sieve hole, P c (D < D 0 ) is the probability that the particle diameter D in the section will be smaller than the sieve hole range D 0 , and D max is the maximum aggregate particle size. According to Walraven's formula, the number of aggregates in each particle size zone can be determined. In this work, the selected aggregate had a maximum particle size of 20 mm and a minimum particle size of 5 mm. Aggregate particle sizes ranging from 5 to 10 mm and 10 to 20 mm were generated with an area ratio of 1:1 and the thickness of the interfacial transition zone (ITZ) was 0.2 mm. The C45 concrete matrix mix proportion was obtained from the literature [42], as shown in Table 1.

Random generation of fibers
In the 2D planar model, if the shape and size of the fibers placed in the model remained unchanged, then the number of fibers could be determined according to the following equation: where N is the number of fibers, L is the fiber length (mm), d is the fiber diameter (mm), ρ f is the fiber volume fraction, W is the length of the model (mm), H is the width of the model (mm), and round(x) is the integral function.
During the model generation process in this work, to distinguish between the two types of fibers, the steel fibers were simulated with a rectangular shape and length of 32 mm, while the polypropylene fibers consisted of in-line segments with a length of 18 mm. The program was written to generate random positions and orientations for the fibers with the assumption that both the steel and polypropylene fibers were randomly and uniformly distributed in the model, and the generated fibers were discarded if they overlapped with the aggregate or were outside the boundary of the model. This process was repeated until the end, and in total, 16 steel fibers and 39 polypropylene fibers were generated.
Based on the above theory, the random placement of aggregates and fibers in the model was achieved and Figure 4 presents the generated 2D HFRC model with dimensions of 150 mm × 150 mm. The generated aggregate particle area was 10857.40 mm 2 , which accounted for 48.26% of the total area. The volume fraction of the steel fibers was 1.0%, the volume fraction of the polypropylene fibers was 0.1%, and all the random variables were assumed to follow a uniform distribution.

Determination of material mechanical parameters
In this work, the mechanical parameter values of each material were determined by constant inversion and fitting based on previous research and comprehensive experimental literature data, as shown in Table 2. The plastic parameters are shown in Table 3. The plasticity parameter values were used as recommended by the ABAQUS software.

Element type selection
According to a series of HFRC mechanical tests conducted by our research group, because the aggregates often failed to reach the plastic stage after the tensile failure of the specimen model, the aggregate was set as an elastic material model according to the material properties. The mortar and ITZ was the potential cracking zone of the concrete, the mortar and ITZ were considered to be elastoplastic bodies, and the elastic-plastic damage model was used. In contrast, the steel fibers and polypropylene fibers both had high tensile strengths, and the fibers tended to be pulled out of the concrete matrix in the uniaxial tensile tests. Therefore, the steel fibers and polypropylene fibers were set as an elastic material model    regardless of their damage, simulating the constraint relationship between the fibers and the concrete matrix in the form of binding without considering the bond slip after the fibers lost their interfacial bonding. The thickness of the ITZ between the aggregate and mortar is usually between 20 and 50 µm in reality, and it is usually less than 100 µm. As an important part of the concrete material, the material properties of the ITZ will greatly affect the macro-mechanical properties of the concrete and the fracture failure process. However, an ITZ that is too small is difficult to mesh.
Consequently, in this study, the ITZ thickness was set to 0.2 mm, and then the aggregate, ITZ, steel fibers, polypropylene fibers, and matrix were meshed to generate a finite element model in which the aggregate, matrix, and ITZ were all CPS3 elements and the steel and polypropylene fibers were T2D2 elements. The meshed finite element model is shown in Figure 5.

Load application
The model was statistically loaded by applying both horizontal and vertical constraints on the left side of the model boundary and a horizontal displacement load on the right side of the model. The maximum tensile stress criterion was used as the basis for the damage determination, which was set to consider failure when the damage factor of the concrete matrix element reached 0.5 [43], and the failure of the element was regarded as the cracking of the model. The uniaxial tensile plane model of the HFRC is shown in Figure 6.

Tensile simulation test analysis
As shown by the stress-strain curve in Figure 7, both the peak stress and peak strain of the HFRC were significantly higher than those of the plain concrete (PC), with a higher peak strength exhibited during the residual stress phase of the curve. Before loading to the 94% peak stress in the curve, there were no significant differences in the rising trends of the curves between the PC and HFRC when the model was in the elastic phase. After exceeding this point (94% peak stress), we observed that the ITZ of some of the coarse aggregates started to accumulate varying degrees of damage, which was depicted macroscopically as the initiation of microcracks. Loading  continued until the stress approached the peak stress, the damage between the adjacent coarse aggregates converged, and the fibers started to exert their reinforcing effects. When the stress reached the peak stress, the peak stress of the HFRC was 4.08 MPa, the peak stress of the PC was 3.27 MPa, and the microcracks started to intersect with each other, forming the fracture surface. After the peak stress point was surpassed, the curves started to enter the softening phase in which the slopes of the curves decreased accordingly, the bearing capacity of the model continued to decrease, and a macroscopic crack formed. Finally, the curves entered the residual stress phase, by which point the main fracture surface of the model had formed. With the increasing displacement load, the fracture surface constantly expanded. Additionally, the curve trend gradually subsided with the increasing strain, and eventually stabilized at smaller values.

Model rationality verification
Two additional HFRC models with the same aggregate particle diameters, aggregate content, and fiber content (1.0% steel fiber content and 0.1% polypropylene fiber content) as the previous model were constructed to consider the effects of a random distribution of aggregates and fibers on the model. However, the aggregate and fiber distribution conditions were mutually different. The peak stress and peak strain obtained from the numerical calculations are shown in Table 4. These values were compared with the experimental values in the literature [33].
As shown in Table 4, by comparing the calculated peak stress and strain from the numerical simulations with the experimental values, we found that all the data errors were within 6%, and the dispersion between the numerical simulation results was small. The three stress-strain curve relationships obtained with the calculations are shown in Figure 8. As shown in the figure, the curves for the three HFRC models basically overlapped in the elastic phase, and when the peak stress was reached, the three models successively entered the softening phase. Although the amplitude of each curve varied, the same downward trend was maintained. Therefore, using the Monte Carlo random distribution algorithm, the random distribution of aggregates and fibers had less of an influence on the tensile softening behavior of the HFRC model, which also reflected the stability of the proposed model.

Crack damage evolution process
The damage to the HFRC specimens was essentially due to the process of crack emergence, expansion, and destabilization. Combined with the damage diagram, this fact can be used to more clearly and intuitively explain the whole process from the beginning of the loading to the emergence of the microcrack and the expansion until the final damage. Accordingly, in this research, when the damage factor of the matrix unit reached 0.5, the unit could be considered to have failed, and the failure of the unit was considered to have been the cracking of the specimen. With the application of the displacement load, the number of units reaching the damage threshold (0.5) gradually increased, which was manifested in the macroscopic process of the gradual expansion of the cracks.  In this work, crack expansion simulation tests were carried out on PC and HFRC. In combination with Figures 9 and 10, the entire process was determined from the beginning of loading until crack initiation, expansion, and, finally, catastrophic failure.
During the early stages of the loading of the HFRC model (incremental step 25), the mortar section was stressed first, and initial damage occurred at the ITZ, which formed microcracks around the aggregate. Moreover, we found that areas with a high density of smaller aggregates were more likely to be damaged first. As the loading continued to increase (incremental step 100), the failed units started to release stress, resulting in the expansion of microcracks from the ITZ to the mortar, while new microcracks continued to initiate elsewhere in the model. This caused adjacent microcracks to evolve, develop, and converge, forming short cracks. Additionally, the polypropylene fibers along the crack propagation path broke as they reached the tensile limit. During the final loading stage (incremental step 200), short cracks started to develop rapidly and intersect with the adjacent cracks. The model started to fracture along the main macroscopic cracks, eventually converging into a penetrating crack that was roughly perpendicular to the loading direction of the model.  Through the crack opening process in the PC, we observed that the initial damage locations were consistent with the HFRC, and these locations were all in areas with dense concentrations of small-sized particle aggregates. However, the PC showed varying degrees of damage on the left side of the initial damage positions and in the lower end of the model, while the cracks in the middle region of the model continued to converge and expand until penetration. In addition, the initial cracks on the right side gradually stopped developing. In contrast to the HFRC, the crack opening process was more abrupt and uncertain, resulting in more stable crack development and a delay in the crack development process.

Mechanism of fiber crack resistance
The path of the crack expansion was further complicated by the presence of fibers, and the path avoided various fibers when the extended crack encountered steel fibers in three specific situations, listed below.
(1) The crack bypassed the fibers. When the crack encountered a fiber obstruction, the expansion path of the crack deflected at the end of the fiber, turned to the side farther away from the fiber, and then continued to expand around the fiber, as shown in Figure 11(a). (2) The cracks stopped at the fibers. As shown in Figure  11(b), there were two cracks in the diagram, where the expansion path of the left crack was deflected after encountering the fibers and then bypassed the fibers, while the right crack could not continue to expand after encountering the fibers due to the fiber barrier, resulting in the termination of expansion in the presence of the fibers. (3) The cracks ran across the fibers. When cracks were generated in the weak areas on both sides, and the fibers were in the expansion path of the crack, the steel fibers continued to have a bridging effect due to the high tensile strength of the steel fibers. This occurred even if the yield strength was reached but the fibers did not break, as shown in Figure 11(c). When the cracks on both sides of the fibers intersected, the fibers were pulled out of the concrete.
Steel fibers and polypropylene fibers were randomly distributed within the concrete matrix, and due to the random distribution of fibers, the development of cracks and fiber obstruction was random. Additionally, the expansion of the crack path was slow and uncertain, resulting in HFRC tensile damage that exhibited more substantial plasticity than ordinary concrete. Furthermore, from the perspective of fracture energy theory, the inclusion of hybrid fibers caused the model to indirectly fracture several times, accompanied by the continuous transfer of energy, which resulted in an overall higher fracture strain value and more substantial energy consumption. As a result, HFRC presented better overall performance than ordinary concrete.

Research on hybrid fiber effect
To explore the influence of the hybrid effect on the mechanical properties of the fibers, we designed orthogonal experiments and carried out numerical simulations with 0.5, 1.0, and 1.5% steel fibers and 0.05, 0.1, and 0.15% polypropylene fibers. The equivalent elastic modulus, peak stress, and peak strain were obtained for different fiber contents, as shown in Figures 12-14, where the term "Vpf" in the figure represents the polypropylene fiber content.
According to the mechanics of the composite materials, the mechanical properties of composite materials can be expressed by where σ is the real-time force on the composite, ε is the strain on the material, E is the modulus of elasticity of the composite, E i is the modulus of elasticity of the material component, and ν i is the volume fraction in the material. Numerous experimental studies have shown the existence of either positive or negative hybrid effects [44][45][46][47][48][49] in HFRC. The positive effects enhance a certain performance condition, while negative effects attenuate the corresponding performance. Figure 11 shows the variation trend of the model's equivalent elastic modulus. When the steel fiber volume fractions were fixed at 1.0 and 1.5%, the equivalent elastic moduli of the model decreased by 16.5 and 20.7% as the polypropylene fiber volume fraction increased from 0.05 to 0.15%. An increase in the volume fraction of polypropylene fibers might have decreased the equivalent elastic modulus of the model, with a negative hybrid effect in terms of the elastic modulus. Additionally, after combining the peak stress variation graph in Figure 12 with the peak strain variation graph in Figure 13, we observed that the peak strain of the model increased with increase in both the steel and polypropylene fiber volume fractions. The best reinforcement was obtained when the steel fiber volume was 1.5% and the polypropylene fiber volume was 0.15%. Therefore, a positive hybrid effect was shown in both the peak stress and strain. Furthermore, the steel fibers were far superior to the polypropylene fibers in terms of enhancing the peak tensile strength and peak tensile strain.
Considering the synergic effect produced by the steel and polypropylene fibers, if the hybrid fibers were regarded as micro-reinforcement bars, the volume difference between the steel and polypropylene fibers was more than two orders of magnitude, forming two levels of micro-reinforcement in the matrix in which the steel fibers acted as the primary reinforcement with a high tensile capacity and the polypropylene fiber acted as secondary reinforcement with a high binding capacity. During the loading process of the HFRC, stress concentrations occurred in some areas, producing high-stress and low-stress areas. In the low-stress areas, the polypropylene fibers dominated and suppressed crack initiation, while in the high-stress areas, the steel fibers dominated, controlling crack development and expansion. Consequently, the steel and polypropylene fibers complemented each other and had superior strengths, ensuring that they could effectively improve the peak tensile strength and peak strain.

Conclusion
In this work, the crack damage evolution process and macroscopic mechanical properties of HFRC were assessed based on numerical simulations. The following main conclusions could be drawn.
(1) Based on the shape characteristics of crushed stone aggregate under realistic conditions, in this work, we conceived a convex polygon approximation to simulate the crushed stone aggregates with rhombic parametric equations. We then assumed that the ITZs were wrapped around the aggregates based on a fiber-aggregate intersection determination algorithm that was added to generate the HFRC model with randomly distributed aggregates and fibers according to the Monte Carlo method. (2) According to the numerical simulation analysis of the crack opening process, we found that initial damage generally occurred in the ITZ, and microcracks mostly formed in areas with dense distributions of small-sized particle aggregates, which confirmed that the ITZ was the weakest region in the entire composite system. (3) The addition of hybrid fibers had varying types of barrier effects on the cracks, such as cracks bypassing the fibers, cracks stopping at the fibers, and cracks crossing the fibers. Therefore, the hybrid fibers allowed the HFRC to form better transmission paths, which limited the generation and development of cracks well and allowed the model to exhibit greater plasticity when subjected to tensile damage. (4) The fiber volume fraction had a greater impact on the macroscopic mechanical properties of the model, as illustrated by the significant increases in both the peak stress and peak strain as the volume fraction of steel and polypropylene fibers increased. Moreover, the polypropylene fibers had a negative hybrid effect on the modulus of elasticity, and the addition of excessive polypropylene fiber decreased the modulus of elasticity, while the steel fiber dominated in enhancing the mechanical properties of the model, and their enhancement effects were superior to those of polypropylene fibers.