Lee et al. (Lee et al., 2010) showed that the interface shear stress is related to the relative slip between the steel fiber and UHPC matrix. Therefore, a two-dimensional cohesive element was selected to simulate the interface between the steel fiber and UHPC matrix (Figure 1A). The interface element consisted of four nodes; the element was composed of a top surface, bottom surface, and intermediate layer. When failure occurred, the element cracked from the intermediate layer. The cohesive model has a concise form, and can be used to simulate the damage and failure of heterogeneous materials. Among the available constitutive relations of the cohesive model, the exponential constitutive relation proposed by Xu and Needleman (Xu and Needleman, 1994) represents the softening of the material as an exponential function; therefore, it is very suitable for simulating the damage of UHPC materials containing steel fibers. Thus, the exponential cohesive constitutive model was selected to characterize the interface characteristics (Figure 1B): when there was no debonding, the interface stress increased linearly with increasing relative displacement of the interface; when the stress reached t s 0 , the element cracked from the intermediate layer. Subsequently, the stress began to decrease, and the element stiffness degraded irreversibly owing to material damage. In ABAQUS, the coefficient D in the range of (0, 1) (D = 0 means no damage, D = 1 means complete damage) was used to describe the degree of damage. The exponential softening relationship is given by D = 1 − { δ s 0 δ s } { 1 − 1 − exp [ − μ ( δ s − δ s 0 δ s f − δ s 0 ) ] 1 − exp ( − μ ) } , (1) where μ is the softening coefficient, which is a dimensionless parameter describing the rate of damage evolution (μ was determined experimentally); δ s 0 is the displacement when initial damage occurs; and δ s f is the displacement when the shear stress is zero.

Then, the actual element stiffness k s could be expressed as k s = ( 1 − D ) k s 0 , (2) where k s 0 is the initial stiffness of the element.

The second nominal strain criterion was used as the damage initiation criterion: { 〈 ε n 〉 ε n 0 } 2 + { 〈 ε s 〉 ε s 0 } 2 + { 〈 ε t 〉 ε t 0 } 2 = 1 (3)

Here, ε n , ε s , and ε t are the normal and two tangential strain values of the cohesive element, respectively, and ε n 0 , ε s 0 , and ε t 0 are the normal and two tangential maximum allowable strain values of the cohesive element, respectively. The operator in the formula is the MacAuley operator: 〈 x 〉 = { 0    x < 0 x    x > 0 (4)

The plastic strengthening constitutive model was used to characterize the constitutive relationship of steel fiber: σ s = { E s ε s                          ε s ≤ ε y f y + k ( ε s − ε y )    ε y < ε s ≤ ε u 0                                ε s > ε u (5) where σ s , ε s , E s , f y , k, ε y , and ε u denote the stress, strain, elastic modulus, yield strength, slope of the hardened curve, yield strain, and the ultimate strain, respectively, of the steel fiber.

A concrete plastic damage model was used to simulate the damage to the UHPC matrix. To facilitate convergence in the calculations, the GFI in ABAQUS was used to define the tension behavior of the UHPC matrix, which means that the tensile softening curve was determined by the tensile strength and fracture energy. The following equation (Zhou, 2017) was used as the constitutive relationship of the UHPC matrix under compression: y = { a x + ( 5 − 4 a ) x 4 + ( 3 a − 4 ) 0 ≤ x ≤ 1 x b ( x − 1 ) 2 + x x > 1 (6)

Here, y is the ratio of the UHPC stress to peak stress, and x is the ratio of the UHPC strain to peak strain. a and b are material parameters obtained from the test. When the volume content of steel fiber is 0%, 1%, 2%, and 3%, the material parameters are a = 1.01, 1.05, 1.1, and 1.2, and b = 57.47, 11.56, 3.76, and 1.45 (Zhou, 2017). When simulating the pull-out of steel fiber from the UHPC matrix, the coefficients a and b in the stress–strain relationship of the UHPC matrix should be selected as the parameters when the steel fiber content is 0%, i.e., a = 1.01 and b = 57.47.

Based on the plastic damage model, after determining the uniaxial tension and compression behavior of the material, the uniaxial constitutive relation needs to be extended to the triaxial stress state. Table 1 shows the plastic damage material parameters required for extending the uniaxial constitutive relation to the triaxial stress state. The parameters commonly used in ordinary concrete (Lubliner et al., 1989) were used here; these parameters have been shown to simulate the plastic damage of UHPC well (Kueres et al., 2015).

The model shown in Figure 2 was used to simulate the pull-out of a fiber from the UHPC matrix. In this model, the UHPC matrix is a cylinder with a diameter of 40 mm and height of 30 mm. The diameter of the fibers was 0.3 mm, fiber length was 13 mm, and the fiber embedded length was 6.5 mm. The matrix and fiber element type was CAX4R (a four-node bilinear axisymmetric quadrilateral element), and the cohesive element described above with no thickness was inserted at the fiber–matrix interface (the element type was COH2D4). The cohesive element shared a node with the matrix and fiber elements (Figure 2). The shear stress distribution along the fiber at a certain loading time could be obtained by extracting the shear stress at the middle layer of the cohesive element. Displacement constraints were applied to the UHPC matrix, and a pulling load p was applied to the center of the fiber. The ABAQUS implicit solver was used to solve this problem.

The material parameters in the model were selected based on the work of Lee et al. (Lee et al., 2010) (Table 2). Lee et al. (Lee et al., 2010) conducted very comprehensive steel fiber pull-out tests, and the material characteristics of the UHPC matrix and fibers used in the tests were consistent with those of the UHPC material commonly used in engineering. Therefore, Lee et al.'s tests were used to verify our model.

In the exponential softening relationship, the value of the coefficient μ affects the shape of the softening curve. When μ = 0, the softening curve is a straight line, and when μ < 0, the softening curve is slightly convex, which closely reflects the material characteristics of UHPC. Therefore, an exponential softening curve with μ = 0, –1, and –2 was separately selected to represent the interface characteristics for the simulation. The simulated fiber pull-out P–δ curve is shown in Figure 3, and was compared with the results of Lee et al. (Lee et al., 2010). As shown in Figure 3, in the ascending section, the slope and peak value of the simulation curve agree well with the test results. The values of the peak load obtained from the test and simulation are respectively 27.75 N and 27.79 N—a difference of only 0.14%. In the descending section, the simulation result is closest to the test result when μ = –1. In this case, the simulated pull-out energy consumption is 105.1 N mm and the test pull-out energy consumption is 107.4 N mm (Lee et al., 2010)—a difference of only 2.2%. These results indicate that the pull-out energy consumption of the steel fiber in UHPC matrix can be simulated more appropriately when μ = –1. Therefore, it is recommended to choose the exponential cohesive model with μ = –1 as the constitutive model of the interface between the UHPC matrix and steel fiber.

Figure 4 shows a partially enlarged P–δ curve for pulling out the steel fiber. To determine the change in shear stress in the steel fiber during the pull-out process and analyze the debonding mechanism of the steel fiber, point A (non-debonding), point B (partial debonding), point C (two-way debonding), and point D (full debonding) were selected (Figure 4), and the shear stress distribution on the steel fiber at these four points was extracted. The pull-out displacements corresponding to points A, B, C, and D are 0.0003, 0.0036, 0.029, and 0.036 mm, respectively. Figure 5 shows the distribution along the buried part of the fiber and the Mises stress distribution in the matrix, at the four moments corresponding to points A, B, C, and D.

It can be seen that when the interface stress is in the elastic stage (point A), the fiber does not debond, and the shear stress at the interface is exponentially distributed along the embedded fiber. As the load increases, debonding begins (point B), the interfacial shear stress at the pull-out end of the fiber reaches the interface strength, and the fiber starts to debond from the matrix from the pull-out end. At this time, the interfacial shear stress at the embedded end of the fiber is still very small, and the interface that is not debonded is still in the elastic stress stage. With increasing pulling force, the shear stress is gradually transmitted to the embedded end, the fibers are debonded from the pull-out end to the embedded end, and the shear stresses on most of the interface and the embedded end reach the interface strength (point C). At this time, the fiber exhibits a two-way debonding phenomenon in which both the pull-out and embedded ends are debonded. Finally, the fibers are all debonded (point D), and the P–δ curve reaches a peak. At this time, the shear stress at almost the entire interface reaches the interface strength, and the fibers are gradually pulled out.

When establishing the UHPC model with a random distribution of steel fibers, it was considered that UHPC is composed of fibers and a matrix, and the randomly distributed steel fibers were modeled in the matrix based on the Monte Carlo method. In this model, the steel fibers were simulated as truss elements with a certain length, the positions and directions of the fibers were random, and the fiber distribution function was uniform. The number of fibers N in the model was calculated from the fiber length l f , fiber diameter d f , fiber volume content V f , and matrix volume V using the following equation: N = 4 × V × V f π × d f 2 × l f (7)

The positions of the fibers in the matrix were determined using the coordinates of the ends of the fibers (Figure 6). First, the coordinates ( x 1 ,   y 1 , z 1 ) of point A were randomly generated in the matrix space using the Monte Carlo method. Assuming that the angle between the fiber and the z-axis was α , and the angle between A ₁ B ₁ (which is the projection of the fiber on the xy plane) and the x-axis was β (both α and β were randomly generated within 0–360°), the coordinates of point B ( x 2 ,   y 2 , z 2 ) could be determined using the following equation: { x 2 = x 1 + l f ⁡ sin ⁡ α ⁡ cos ⁡ β y 2 = y 1 + l f ⁡ sin ⁡ α ⁡ sin ⁡ β z 2 = z 1 + l f ⁡ cos ⁡ α (8)

When point B was generated, it was additionally necessary to determine whether point B exceeded the boundary described by the boundary function; if it exceeds the boundary, it must be generated again. Finally, it was necessary to determine whether there was an overlap conflict between the two fibers in space. For this purpose, the minimum distance between the two fibers in the space was calculated and compared with the fiber diameter. If the former was greater, there were no conflicts. Figure 7 shows the flowchart of the algorithm applied to generate randomly distributed fibers.

In the UHPC model with randomly distributed fibers, the slip characteristics of the steel fiber–matrix interface were represented by the material properties of the steel fibers. By using the following equation, the P–δ relationship obtained in Figure 4 can be converted to the stress–strain relationship of the steel fiber: { σ = 4 P π d 2 ε = δ l f (9) where σ and ε are the equivalent stress and equivalent strain, respectively; P is the pull-out load; δ is the interface displacement; and l f is the fiber length.

The stress–strain relationship of steel fiber obtained by this transformation includes the interface characteristics, which can reflect the characteristics of the steel fiber pull-out and debonding from the UHPC matrix during the loading process.

The GFI provided by ABAQUS was used to set the tensile constitutive of the UHPC matrix, the softening behavior of the matrix was characterized by defining the fracture energy G _Fm, and the interface strength was characterized by the equivalent interface strength τ a v e . Therefore, the values of G _Fm and τ a v e determined the form of the softening curve of the matrix after crack formation. According to the method described in this section, a model was established to simulate a UHPC uniaxial tensile test (Yuan et al., 2009) conducted by our research group (Figure 8). Fixed constraints were applied to the left end of the model, and horizontal displacement was applied to the right end of the model. The material parameters were the same as those listed in Table 2.

The parameters G F m = 0.06 ,   0.08 , a n d   0.1   N / m m and τ a v e = 2.4 f t , 2.8 f t , and 3.2 f t were selected for the simulation. Because the random distribution of fibers influences the tensile behavior significantly, five simulations were performed with each parameter combination. The results show that the softening section of the curve agrees best with the test results when   τ a v e = 2.8 f t . In this case, the stress descending stage (after the peak) agrees most closely with the test results when G F m = 0.08 N/mm. Therefore, the value of G F m in the model is recommended to be 0.08 N/mm, and the value of τ a v e was selected as 2.8 f t .

Using this model, the UHPC uniaxial compression test reported in the literature (Zhou, 2017) was simulated, and the simulation parameters were selected according to the test results. Figure 9 compares the simulated and tested stress–strain curves, demonstrating that they agree well. The elastic moduli of the ascending phases of the two curves are very similar, and the trends of the descending phases of the curves are the same. The difference in compressive strength between the simulation results and the test results is 0.5%–4.4%, and the difference in peak strain is 0.03%–2.7%. These results indicate that the proposed mesoscopic model can simulate the mechanical behavior and failure process of UHPC well.

For UHPC with a fiber volume content of 2.5%, a tensile simulation of the notched specimens was conducted. Figure 11, which shows the crack distributions of the five UHPC specimens with different random fiber distributions, indicates that cracks in UHPC tend to pass by the fibers as the cracks propagate. When the crack enters the fiber-dense area, it tends to bypass this area and expands to the fiber-sparse area (the area inside the white circle in Figure 11). Therefore, if the fiber distribution is relatively even, the crack is relatively straight, such as in specimen one; if the fiber distribution is uneven, the crack is rugged, such as in specimen 4. It can also be seen that at 2.5% fiber volume content, multiple cracks tend to occur on the specimen. Because the fiber density is relatively large in this case, the fiber in front of the crack has a strong crack resistance effect, and it is difficult for the crack to propagate directly forward until specimen failure. New cracks are created, which propagate at other interfaces, causing multi-crack failure mode.

The stress–strain curves for each specimen are shown in Figure 12. The specimens exhibit obvious strain-hardening behavior, and the average tensile strength (8.88 MPa) and peak strain (1,560 × 10⁻⁶) are significantly higher than the average initial crack strength (5.87 MPa) and initial crack strain (150 × 10⁻⁶), respectively. Figure 12 further shows that the random distribution of fibers has little effect on the initial cracking strength and initial cracking strain of UHPC, but has greater effects on the tensile strength and peak strain. This finding indicates that the bridging effect of the fibers occurs only after the appearance of macroscopic cracks.

Figure 13 shows the stress cloud of the fibers of specimen 2 after the matrix cracks. Not only the fibers between the main cracks but also those in the tension zone are subjected to large stress. Thus, most of the steel fibers in the tension zone are involved in resisting tensile stress, which is the reason UHPC exhibits strain-hardening properties and high toughness. Figure 13 further shows that the fibers whose length direction is close to the tensile stress direction bear significantly more stress (the red fibers are shown in the figure), indicating that the steel fiber orientation affects the mechanical properties of UHPC significantly.

To study the effect of fiber orientation on the tensile behavior of UHPC, the tensions of UHPC specimens with different fiber inclination angles (the other parameters were the same) were simulated. The fiber inclination angle was defined as the angle between the fiber and tensile direction of the specimen. In the simulation, four fiber inclination angle regions were considered: 0°–15°, 15°–30°, 30°–45°, and 45°–60°.

The crack propagation of the specimens under different fiber inclination angles is shown in Figure 14A. It can be seen that the smaller the fiber inclination angle, the more difficult it is for the crack to bypass the fiber, and the more effective the crack resistance effect of the fiber, resulting in sinuous cracks and more fibers on the fracture surface. When the fiber inclination angle is 0°–15°, there are multiple cracks in the UHPC; when the fiber inclination angle is 15°–30°, the crack propagation in the UHPC is still sinuous and complicated, but is relatively gentle compared with the cracking of specimens with fiber inclination angles of 0°–15°. Further, when the fiber inclination angle is 30°–45°, the crack propagation is relatively gentle, and the number of fibers on the fracture surface decreases. When the fiber inclination angle is 45°–60°, the cracks develop almost perpendicular to the direction of the tensile stress.

Figure 15 shows the tensile stress–strain curves of UHPC with different parameters, and Table 3 summarizes the tensile mechanical properties of UHPC with different parameters. Figure 15A shows that when the fiber inclination angle is smaller than 45°, the stress–strain curves have an obvious plastic stage, and the smaller the fiber inclination angle, the greater the slope of the plastic stage. As shown in Table 3, with increasing fiber inclination angle, the initial crack strength and elastic modulus decrease gradually, and when the fiber inclination angle exceeds 30°, the tensile strength of the UHPC decreases significantly.

In general, upon comparing the specimens with fiber inclination angles of 15°–30° with the specimens with fiber inclination angles of 0°–15°, although the tensile and initial crack strengths of the former are slightly lower than those of the latter, the toughness and ductility of the former are obviously larger than those of the latter. Therefore, it is believed that UHPC can obtain the best tensile properties when the fiber inclination is in the range of 15°–30°.

Figure 14B shows the tensile crack modes of UHPC specimens with five different fiber volume contents (1%, 1.5%, 2%, 2.5%, and 3%). It can be seen that the crack propagation is relatively smooth in the specimen with a fiber content of 1%. When the fiber volume content exceeds 2%, the cracks in the specimen become sinuous. When the fiber content is greater than 2.5%, the specimens begin to exhibit multiple obvious cracks. When the fiber volume content is 3%, many cracks are discretely distributed in the weak interface area.

Figure 15B and Table 3 show that when the fiber volume content is greater than 2%, especially when the fiber content is greater than 2.5%, the stress–strain curve of the specimen exhibits strain-hardening behavior. In addition, the tensile strength and peak strain are obviously improved, the descending stages of the curves become gentle, and the toughness of the matrix is enhanced. Table 3 shows that the increase in fiber volume content has obvious effects on the tensile strength and peak strain of UHPC, but has little effect on the initial crack strength. The initial crack strength only increased by a small degree (5.33–6.04 MPa) with increasing fiber volume content.

The fiber volume content was set to 2%, and the fiber length was set to 8, 10.5, 13, 15.5, and 18 mm to study the effects of the fiber length on the tensile properties of UHPC. As shown in Figure 14C, when the fiber length is 8 mm, the anti-cracking effect of the fiber is not sufficiently strong, and most specimens crack along intermediate cracks. The peak load is reached soon after the initial crack occurs, and the bearing capacity begins to decrease (Figure 15C). As the fiber length increases, the crack propagation becomes sinuous, and correspondingly, the tensile strength and peak strain of the specimens increase significantly (Table 3). When the fiber length is greater than 15.5 mm, the stress–strain curve shows obvious strain-hardening characteristics, and multiple cracks develop in the specimens simultaneously (Figure 14C). This finding indicates that as the fibers in the matrix become longer, it becomes more difficult for the cracks to bypass the fibers, increasing the tensile energy consumption. When the fiber is sufficiently long (l _f > 15.5 mm in this study), it is difficult for the cracks to bypass the fiber. Therefore, new cracks must be generated at other weak interfaces of the matrix; thus, the specimen failure mode is multi-crack failure. Nevertheless, Table 3 shows that with the same fiber volume content, the change in fiber length has little effect on the initial crack strength and elastic modulus of UHPC.

With the same fiber distribution and other parameters, the equivalent interface strength τ a v e of the specimens was set to 10, 12, 14, 16, and 18 MPa in the tensile simulation to study the effects of the interface strength. Figure 14D shows that the change in the interface strength has little effect on the crack propagation of UHPC, which mainly depends on the fiber distribution. Figure 15D shows that with increasing interface strength, the descending stage of the stress–strain curve becomes smoother and the matrix toughness increases. The results in Table 3 indicate that the tensile strength and peak strain of the UHPC also increase with increasing interface strength.

To study the effect of matrix strength on the tensile behavior of UHPC, the matrix strength f t of the specimens was set to 3, 5, 7, and 9 MPa to conduct tensile simulations, and the fiber distribution and other parameters of the specimens were the same.

As shown in Figure 14E, when the strength of the matrix is low (3 MPa), once the crack development is resisted by the fibers, the cracks begin developing from other weak parts of the matrix instead of splitting the fiber–matrix interface. At this time, it is difficult for the fiber to exert the crack resistance effect. Therefore, although the specimens exhibit strain-hardening characteristics in this case, the tensile strength of the UHPC specimen is very low to exploit its material advantages. With increasing matrix strength, the initial crack strength and tensile strength of UHPC increase gradually (Table 3), but the crack propagation path is hardly affected by the matrix strength. When the matrix strength and interface strength are well matched (matrix strength ≥7 MPa), the crack resistance effect of the fibers can be fully exerted (Figure 15E). To further illustrate the effects of matrix strength on the crack resistance of fibers, the stress difference σ − σ 0 curve was drawn (Figure 15F), where σ 0 represents the tensile stress of the pure UHPC matrix (without fibers) specimens with various matrix strengths. It can be seen that the stress difference σ − σ 0 increases with increasing matrix strength, indicating that for the same fiber parameters, the higher the matrix strength, the more effective the fiber crack resistance.