Steel fiber reinforced concrete (SFRC) is a multi-phase composite material in which steel fibers are incorporated into conventional concrete. It has been observed that the inclusion of steel fibers effectively inhibits the formation and propagation of microcracks in concrete, thereby improving the tensile, bending, fatigue and impact resistance of concrete. Therefore, SFRC is widely used in various fields such as construction engineering, bridge engineering, water conservancy engineering and tunnel engineering (Meng et al. 2016). In tunnel and underground engineering construction, SFRC is widely employed as a supporting structure through high-speed injection molding of compressed air (Wang et al. 2021). The lining of underground engineering not only bears ground stresses but also experiences dynamic loading events such as earthquake, explosion or excavation, leading to varying degrees of dynamic response (Tsinidis 2018; Feng et al. 2022). Therefore, it is of great significance to study the dynamic response and mechanical behavior of SFRC under static and dynamic loading for underground engineering support.

In the past few decades, numerous scholars have studied the mechanical behavior of steel fiber reinforced concrete under static loading conditions, resulting in a wealth of research outcomes (Luo et al. 2022; Nataraja et al. 1999; Lu and Hsu 2006). Recently, many scholars have commenced conducting diverse experimental studies on SFRC to explore its dynamic mechanical behavior. Marar et al. (2001) used a drop hammer impact tester to conduct impact tests on SFRC and found a logarithmic relationship between toughness energy and impact force. Lok and Zhao 2003; 2004) and Zhao and Lok (2005) studied the effect of strain rate on dynamic compressive strength and toughness of SFRC with different fiber volume fraction using a split Hopkinson pressure bar (SHPB). These studies explored the application of SHPB technology to SFRC and improved the experimental techniques to overcome challenges encountered in conventional SHPB facility. Wang et al. (2011) conducted static and dynamic compression tests on steel fiber reinforced concrete with different steel fiber contents. Their research shows that as the steel fiber content increases, the failure mode of concrete transitions from brittle to ductile, with the toughness energy being proportional to the fiber content under dynamic compression. Sun et al. (2018) employed a SHPB device with a diameter of 75 mm to perform dynamic compression tests on SFRC with six different steel fiber volume fractions. They investigated the relationships between the dynamic strength increase factor, fiber enhancement factor, critical strain and energy absorption with respect to strain rate and fiber volume fraction. Liao et al. (2020) utilized the SHPB device to perform dynamic compression tests on SFRC specimens with three steel fiber aspect ratios, investigating the influence of aspect ratio on the dynamic mechanical behavior of steel fiber reinforced concrete. Zhao et al. (2023) examined the effect of steel fiber on the dynamic compressive properties of SFRC at different strain rates and employed CT scanning technology to elucidate its enhancement mechanism at both macroscopic and microscopic levels.

Owing to the limitation of experimental technology and cost associated with dynamic test and blasting experiments on SFRC, it becomes challenging and impractical to investigate the dynamic characteristics and failure mechanism of SFRC solely rely on a large number of experiments. Therefore, numerical simulation is considered to be an effective method to study the dynamic response and damage mechanism of SFRC. With the aid of computing technology and commercial software, such as LS-DYNA (Yoo and Yoon 2014), AUTODYN (Nyström and Gylltoft 2011), ABAQUS (Swaddiwudhipong and Seow 2006), etc., a large number of models based on finite element method have been established to study the static mechanical characteristics of steel fiber reinforced concrete. In recent years, researchers have also begun employing mesoscale models based on the finite element method to study the mechanical behavior of steel fiber reinforced concrete under dynamic loading conditions. Xu et al. (2012) proposed a 2D axisymmetric SFRC model consisting of circular aggregates, fibers, and mortar to modeling the dynamic compressive properties of SFRC, and studied the effect of steel fiber on dynamic strength increase factor. Fang and Zhang (2013) considered the random position and discrete orientation of straight circular steel fibers in the model, developed a three-dimensional model consisting of steel fibers and concrete matrix, and investigated the dynamic growth factors of SFRC materials with different fiber volume percentages under high strain rate loading as well as the failure modes under blast loading. Su et al. (2017) established a two-phase model composed of concrete matrix and steel fiber to study the tensile mechanical behavior of SFRC under static and dynamic loads. Wu et al. (2020) employed a 3D mesoscale modeling approach to simulate the SHPB test for SFRC under high-rate loadings, observing that the cracking failure primarily occurred at the specimen edge, and the dispersed steel fibers bridged and interconnected the concrete fragments.

In general, the mesoscale model of steel fiber reinforced concrete based on finite element method has been developed and extensively studied. Zhang et al. (2021) stated that, as a continuum-based method, the traditional finite element method does not accommodate discontinuity within the element, nor does it allow the separation of adjacent elements sharing boundary nodes. These limitations present challenges in simulating crack formation within or along the element boundary, rendering it intricate to faithfully replicate the fracture process of steel fiber reinforced concrete from a meso-fracture standpoint. Moreover, the finite element method barely accounts for the interaction between steel fibers, and between steel fibers and coarse aggregates. Consequently, it fails to capture the collision, rebound, and stability processes between steel fibers and aggregates during mixing, resulting in an inaccurate representation of the actual distribution of steel fibers and aggregates within the matrix.

The discrete element method (DEM) is a numerical simulation technique based on particle dynamics, making it suitable for simulating the nonlinear mechanics and fracture behavior of complex geotechnical materials. Over the past few decades, many researchers have employed DEM to investigate the strength characteristics and cracking behavior of concrete under static loading (Nitka and Tejchman 2015; Tran et al. 2011), and have analyzed its influence on the strength and failure characteristics of concrete from the mesoscopic level of aggregate shape (Wang et al. 2020) and ITZ strength (Lei et al. 2022). Moreover, DEM has proven to be an effective approach for studying the dynamic response of materials. Several scholars (Kou et al. 2019; Cai et al. 2023; Zhou et al. 2022, 2023) have initiated the construction of concrete models based on DEM to study the dynamic mechanical behavior of concrete.

Therefore, in this study, based on the discrete element software PFC3D, we have developed a novel three-dimensional mesoscale model consisting of four phases: concrete matrix, aggregate, interfacial transition zone and steel fiber. This innovative model framework considers the random shape of aggregates and the initial distribution of steel fibers within the model region, fully accounting for the real interaction between the two during the formation process of the sample. The corresponding material parameters and contact models have been established, and the bond-slip behavior between steel fiber and matrix has been simulated using an improved soft-bond model. After randomly generating aggregates within the model area through template importation, their size distribution remains fixed throughout subsequent calculations, although they retain the ability to move or rotate within the model upon colliding with steel fibers. Additionally, a comprehensive calibration procedure has been devised to determine the corresponding mesoscopic contact parameters based on the macroscopic mechanical parameters of concrete, ensuring the rationality of the parameters. By comparing the simulation results with experimental findings from a macroscopic mechanical perspective, a good agreement is observed between the numerical simulation results and the relevant experimental results. Building upon this foundation, SFRC has been simulated through static compression tests and SHPB dynamic compression tests, enabling analysis of the mechanical behavior and fracture mechanisms of SFRC under static and dynamic compression loads from multiple perspectives.

2.1 Size and distribution of aggregates

In concrete engineering, it is widely acknowledged that the aggregates employed typically consist of gravel or crushed materials, which exhibit random, angular, and irregular polyhedral shapes. Therefore, this study adopts a methodology that constructs polyhedron based on ellipsoidal surfaces, generating five distinct templates of irregular polyhedral aggregates using MATLAB.

Given an elliptical matrix, the coordinates of any point on the ellipsoidal surface within the spherical coordinate system can be determined using five parameters (\(R_{1}\), \(R_{2}\), \(R_{3}\), \(\theta\), \(\varphi\)). In the construction of a polyhedron composed of N vertices from an ellipsoidal primitive, the vertices are segregated into two sets, selected randomly from the upper and lower portions of the ellipsoidal primitive surface. Assuming the count of vertices chosen selected from the upper half of the ellipsoidal primitive surface is denoted as \(N_{1}\), then \(\theta_{i}\) and \(\varphi_{i}\) of these random points are:

where, \(\eta_{1}\) and \(\eta_{2}\) are two independent random variables uniformly distributed in the range [0,1]. \(\delta\) is a variable whose value is usually 0.3.

The remaining polyhedron vertices are randomly selected from the lower half of the surface of the ellipsoid primitive, the \(\theta_{i}\) and \(\varphi_{i}\) can be determined similarly employing the aforementioned equation.

By importing the aggregate template, the random particle size range, rotation angle, distribution type and volume fraction of the aggregate can be specified in the model domain. To achieve an equivalent particle size distribution to that observed in the experiment, PFC3D software regulates the size distribution of the generated aggregate, ensuring its homogeneous dispersion throughout the model domain. Subsequently, a specific time step is computed to establish equilibrium of the aggregate within the model domain (Figs. 1 and 2).

a Ellipsoidal primitives b Irregular polyhedrons

Irregular polyhedron templates and aggregates

2.2 Generation of random fibers

2.2.1 Steel fiber model

In the conventional finite element model, the steel fiber is simulated by the same continuum as its actual shape (Zhang et al. 2020). Although the interaction and deformation between the steel fiber and the matrix can be taken into account, the interaction among the steel fibers is disregarded. Conversely, in the discrete element approach, the basic element being a spherical particle prohibits the direct generation of cylindrical steel fibers within the model. Therefore, a deformable block cluster (referred to as a flexible cluster), consisting of multiple spherical particles, is employed to simulate the steel fibers. The cluster represents a combination of numerous particles. This method enables the simulation of various steel fiber shapes, such as straight fiber and hooked-end fibers (Fig. 3a). In this study, a straight steel fiber with an aspect ratio of 60 is adopted for the simulation. Different mechanical parameters can be assigned within the cluster, as well as between the clusters, and between the cluster and the external matrix. Upon being subjected to specific forces, the particle assembly within the cluster may be destroyed. This method offers several advantages over the finite element method, including the ability to simulate the interaction between steel fibers and provide a more realistic representation of the stress experienced by the steel fibers (Figs. 4, 5, 6 and 7).

a Single steel fiber b Contact connection between internal and external steel fiber

a Randomly distributed aggregates b Randomly distributed steel fibers c A mixture of steel fiber and aggregate interaction

Mesoscale model building process

Normal force–displacement contact law in bonded state a Linear-parallel-bond model b Soft-bond model

Force–displacement contact law in unbonded state a normal direction b tangential direction

2.2.2 Generation of random fibers

Due to the limitations of the PFC3D software in directly generating clusters for the random distribution of steel fibers within the model domain, an innovative approach has been proposed. Initially, clump templates comprising individual steel fibers are constructed within the designated model region. Each clump represents a rigid body formed by aggregating multiple spherical particles. Subsequently, utilizing the random number algorithm embedded in the PFC software, steel fibers are generated with random positions and orientations within the spatial region, based on the established steel fiber templates. The generation process ceases once the model porosity reaches a predefined value. The generated clumps are initially equilibrated within the model domain. Subsequently, employing the fish function, the contact between the steel fibers and the aggregates is detected to identify and eliminate any overlapping steel fibers within the aggregates, thus producing a mixture of steel fibers and aggregates.

The linear contact model is employed to simulate the interaction of collision and friction between steel fiber and aggregate during the sample forming process. The simulation continues until the unbalanced force between steel fiber and aggregate diminishes below a threshold of 1e-4, ensuring the accurate representation of the distribution of steel fibers and aggregates within the model. Finally, by utilizing the fish function, the pebbles within the clumps are replaced with individual spheres, and different mechanical parameters are assigned to them within and between clusters, thereby accomplishing the objective of generating randomly simulated cluster fibers.

2.3 Generation of cement matrix

The cement matrix is formed by generating ball particles with specified particle gradation within designated the model region. The volume fraction of the cement matrix is controlled by the porosity during generation. However, the random generation of ball particles across the model area may result in their distribution within the coarse aggregate, leading to an unbalanced model. Therefore, it is necessary to obtain the geometric position parameters of the ball and the aggregate clump through the fish function. If the distance between the ball particles and any pebble in the clump is less than the clump size, the ball particles are removed. In our model, ITZ is not considered as a separate entity. It is represented as a zero-thickness interface, representing the interaction between the cement matrix and the aggregate in the model. The construction process of the entire model is as follows:

Step 1. Generate random aggregates in the model domain.

Step 2. Generate random steel fiber in the model domain.

Step 3. Allow for full interaction and equilibrium between the aggregates and steel fibers within the model domain, while removing any overlapping steel fibers within the aggregates.

Step 4. Generate cement matrix particles based on the desired volume fraction, while eliminating particles that are distributed within the aggregates. Balance the sample to eliminate any unbalanced forces between the aggregates, fibers, and matrix..

Step 5. Transform the matrix composed of balls into a set of clumps composed of individual pebbles, which enables the subsequent application of steel fiber contact parameters.

Step 6. Utilize the fish function to convert the steel fibers from a rigid cluster clump into flexible cluster, assigning different mechanical parameters between and within the cluster.

Step 7. Assign contact groupings based on interface interactions between the three materials, employing the fish function, and assigning distinct contact parameters.

In this study, the cement matrix is considered as a bonding material, and the aggregate is considered as an incompressible rigid cluster. Since there is no elastic–plastic model for simulating metal materials in the DEM and the steel fiber exhibits negligible fracture behavior during the compression process, it can be simplified as a linear elastic material. For the interface between the three materials, existing research suggests that the interface between the cement matrix and other components can be treated as a bonded interface, while the interface between steel fibers and between steel fiber and aggregate regarded as non-bonded interfaces. Thus, different models and parameters are assigned to these interfaces, which will be discussed in detail below.

3.1 Contact model employed in the DEM

At present, the Linear-parallel-bond model (Potyondy and Cundall 2004) is widely used in the simulation of bonding materials such as concrete and has demonstrated satisfactory outcomes. A parallel bond can be conceptualized as an array of elastic springs with uniform normal and shear stiffness, distributed evenly across a rectangular cross-section, situated at the contact plane and centered at the contact point. These springs operate in parallel with those in the linear component. Following the establishment of the parallel bond, if relative motion ensues at the contact, a force and moment are induced in the bond material, exerting influence on the two contacting surfaces. These forces and moments can be associated with the maximum normal and shear stresses experienced within the bond material at the bond edges. Should either of these maximum stresses exceed the corresponding bond strength, the parallel bond fractures, resulting in the bond being broken. Consequently, a linear model is adopted, as shown in Fig. 8.

Contact used in this model a Linear-parallel-bond model b Soft-bond model c Linear model

Based on previous research (Isla et al. 2022; Georgiadi-Stefanidi et al. 2015), from the microscopic point of view, the simulation of bond-slip behavior of the interface between steel fiber and matrix has been implemented by means of an improved soft-bond model in this study. When soft-bond model is in a bonded state, its behavior resembles that of the Linear-parallel-bond model, incorporating a frictional strength parameter that restricts the shear force. If the bond strength is exceeded either in shear or tension, there is a possibility of bond failure. Unlike the Linear-parallel-bond model, however, the bond is not immediately eliminated upon failure. Instead, its normal force and shear force may enter into a softening regime, progressively diminishing until the bond stress reaches a threshold value at which point the bond is considered broken. In this model, the slope of the softening stage and the tensile strength at failure are controlled by two parameters, the softening factor \((\xi )\) and the softening tensile strength factor \((\gamma )\).

In the event that the bond is within the softening regime, the normal stress at the bond periphery is examined against the softening envelope in the subsequent manner.

The maximum stress is denoted by:

where \(l_{\text{c}}\) is the critical bond elongation (at peak strength) and \(\xi\) is the bond softening factor, \(\delta_{\text{c}}\) is the bond tensile strength.

where the current bond elongation \(l\) is given by:

where \(\delta l\) represents the degree of the bond elongation since the onset of softening and \(\left| {{\mathbf{\delta \theta }}_{{\mathbf{b}}} } \right|\) denotes the extent of the accumulated bending since the onset of softening.

Upon the breaking of bonding, the soft-bond model undergoes a transition into an unbonded state. In the unbonded state, it provides the ability to transmit both a normal force and a shear force at the contact point, as shown in Fig. 9. In the event of contact slips, a coulomb friction criterion is enforced on the shear force such that:

Parameter calibration flow chart

The slip state is updated as:

If the slip state is true, signifying sliding contact, the contact is deemed to be in motion. Consequently, in this study, the non-zero residual frictional strength, which persists post complete debonding at the interface, is modelled by accounting for the interface contact friction coefficient subsequent to interface debonding. The elastic-softening-debonding mechanical behavior of steel fiber-matrix interface were simulated with this model.

In addition, the linear contact model is also employed to simulate the linear elastic friction behavior of the unbonded interface between fibers and between fibers and aggregates.

3.2 Selection of the material contact model

For steel fiber, a typical elastic–plastic material, there is no constitutive model in the discrete element method to simulate the material properties. Generally, it is simplified as a linear elastic material and simulated by the parallel bond model. The cement matrix uses a parallel bond model that is currently widely used to simulate cementitious materials such as cement. For the ITZ interface and the fiber-matrix interface, the bond-slip behavior of the interface tends to be nonlinear. To simplify the model, the soft bond is used to simulate the elastic bonding-softening-debonding behavior of the interface. For the interaction between fibers and between fibers and aggregates, the friction coefficient of the interface applies a Coulomb limit to adapt to the sliding, which is simulated using a linear model. The contact model used in each part is shown in Table 1.

3.3 Development of the SHPB system using FDM-DEM coupling

The SHPB test system comprises an incident bar and a transmitted bar, with a specimen sandwiched between them. A specially shaped striker impacts the incident bar to generate a compressive stress wave. In this paper, the SHPB simulation system is established by the coupling method of finite element and discrete element. The incident bar and transmitted bar are established by the zones in FLAC3D, and the bar adopts the linear elastic model. The SFRC model based on DEM locates between the incident bar and the transmitted bar. At the left end of the incident bar, the rigid rock block (Rblock) in the PFC3D is used to generate a striker with a specific shape, simulating the impact of the bullet. Measurement points are placed in the middle of the incident bar and the transmitted bar to monitor the stress signal. The complete SHPB model is presented in Fig. 10.

a Coupled SHPB system; b Discrete specimen containing a single flaw; c Coupled wall elements; d Rblock striker of a variable cross-section shape

In order to meet the basic assumptions of the SHPB test, a dynamic local damping value of 0 is set for the incident bar and the transmitted bar. The right boundary of the transmitted bar is designated as a quiet boundary to absorb the stress wave. Additionally, a coupled wall element is implemented as the medium for stress wave propagation between the continuous bar and the discrete sample, allowing for equivalently transmission of the stress wave into the specimen, as shown in Fig. 11. This SHPB model based on FDM-DEM coupling has gradually found application in simulating dynamic impact tests on concrete (Zhou et al. 2022, 2023; Wang et al. 2022).

Propagation of stress wave in bar

For the discrete element model, since the macroscopic mechanical parameters cannot be directly used for the mesoscopic contact parameters, it is necessary to calibrate the parameters (Shang et al. 2018). This paper presents a developed and comprehensive calibration procedure to validate the microscopic contact parameters required in our model. Specifically, the calibration process primarily relies on the experimental work conducted by Zhao et al. (2023), as depicted in Fig. 9.

4.1 Calibration of parameters using direct tensile test for steel fiber

Initially, the macroscopic elastic modulus and tensile strength are derived by fitting the relevant experimental macroscopic mechanical parameters, such as the elastic modulus and tensile strength in the experiments (Zhao et al. 2023). To achieve this, a steel fiber consisting of 30 particles is constructed and subsequently simulated the direct tensile test of the steel fiber. By applying opposing velocities to the particles at both ends, the fundamental contact parameters of the steel fiber were obtained (Fig. 12).

a Direct tensile test of steel fiber b Comparison of the stress–strain curve from laboratory experiments and DEM simulation

4.2 Conducting uniaxial compression tests on mortar

Cylindrical specimens of 50 mm × 100 mm were prepared according to the mass ratio of cement: water: sand of 1:0.5:1.8, as reported in the literature (Zhao et al. 2023). Subsequently, uniaxial compressive strength were conducted by a universal testing machine after 28 days of curing. To begin with, the normal shear stiffness ratio of linear contact and parallel bonding was calibrated by aligning the Poisson's ratio of mortar. Next, the cohesion, tensile strength, and normal tangential bond strength ratio of the parallel bond model were continuously adjusted to achieve a match with the average uniaxial compressive strength and macroscopic failure mode of the mortar. This calibration process is illustrated in Fig. 13.

Calibration of the matrix mesoscopic parameters. a Comparison of the stress–strain curves from laboratory experiments and DEM simulation; b Failure patterns of the simulated results and the laboratory samples

4.3 Execution of uniaxial compression tests on concrete

According to the test parameters in the literature, a cylindrical specimen of dimensions 37 mm × 75 mm was established in the PFC software, and the uniaxial compression test was simulated and calibrated with the test results, as shown in Fig. 14. According to the existing research, the elastic modulus of the ITZ has minimal influence on the overall elastic modulus of the model. Therefore, it is assumed that the elastic modulus of the ITZ is equal to that of the matrix. Initially, the effective modulus and normal tangential stiffness ratio of the matrix are calibrated according to the elastic modulus of the uniaxial compression curve. Subsequently, through continuous adjustments of the cohesion and tensile strength of the ITZ, the stiffness and strength micro-contact parameters of the ITZ interface were determined to match the uniaxial compressive strength of the concrete.

Calibration of the ITZ mesoscopic parameters. a Comparison of the stress–strain curves from laboratory experiment (Zhao et al. 2023) and DEM simulation; b Failure patterns of the simulated results and the laboratory samples (Zhao et al. 2023)

4.4 Implementation of pull-out tests on single steel fibers embedded in concrete

At present, the bond strength between steel fiber and matrix is generally determined by the pull-out test of a single steel fiber in the matrix. Consequently, this study adopts a simulation approach based on the representative pull-out test conducted by Isla et al. (2015) to determine the micro-mechanical parameters governing the interaction between steel fibers and cement matrix in the developed model. The same steel fiber drawing model (Figs. 15a and b) employed in the aforementioned literature was established, and a cylindrical specimen of 100 mm × 80 mm was generated. Subsequently, a 1mm diameter and 50mm length steel fiber is incorporated, and an initial velocity is assigned to initiate the fiber pull-out process. The stiffness parameters characterizing the fiber-matrix were determined based on the slope of the linear phase in the pull-out force–displacement curve, and the strength parameters were determined by the peak pull-out load. Additionally, the post-peak stage is utilized to assess the softening factor \((\xi )\) and softening tensile strength factor \((\gamma )\) of soft-bond model.

Calibration of the fiber-matrix interface mesoscopic parameters. a Diagram of experiment; b The established numerical model; c Comparison of the force–displacement curves from laboratory experiments (Shang et al. 2018) and DEM simulation

Despite employing the same cement matrix strength for testing purposes, variations in the mix ratios may yield discrepancies in the bond strength between the steel fiber and the cement matrix. Therefore, it is necessary to incorporate the obtained interface parameters into the model and compare them against the uniaxial compressive strength of steel fiber reinforced concrete with the initial volume fraction. Ultimately, the microscopic contact parameters obtained through the whole calibration program are shown in Tables 2 and 3.

5.1 Analysis of quasi-static compression testing

5.1.1 Stress–strain curve

Figure 16 compares the experimental and numerical quasi-static compressive stress–strain curves of SFRC with different fiber volume fractions. The different fiber volume fractions are represented by SFRC0, SFRC1, SFRC2, and SFRC3, corresponding to steel fiber volume contents of 0%, 1%, 2%, and 3% respectively. Both the experimental and numerical stress–strain curves of SFRC under static compression load exhibit three stages: elastic stage, elastic–plastic stage and descending stage. The numerical model accurately captures the changes in compressive strength of SFRC with varying steel fiber content and demonstrates the positive influence of steel fiber on improving the compressive strength and ductility of concrete. The numerical simulation curves are in good agreement with the experimental curve. From the perspective of uniaxial compressive strength, it is observed that the compressive strength of SFRC in the test initially increases and then decreases with the increase of steel fiber content. However, the compressive strength of SFRC in numerical simulation continues to increase. At high steel fiber contents, the test specimens exhibit a significant number of pores due to the evaporation of free water during the curing process. These pores reduce the bond strength between the steel fiber and the concrete matrix, thus diminishing the bridge effect provided by the steel fiber (Zhao et al. 2023).

Comparison of the stress–strain curves of SFRC with different steel fiber contents from experiment (Zhao et al. 2023) and DEM simulation

5.1.2 Failure mode

From the failure mode comparison (shown in Fig. 17) of SFRC with different content, numerical simulation proves to be effective in predicting the failure mode of SRFC. The specimen with 0% content exhibited extensive damage and spalling. In terms of SFRC with 1% content, significant spalling occurred at the specimen's edge, while the central portion of SFRC with 2% and 3% content remained mostly intact, with only minor spalling observed at the sides. These comparisons demonstrate a satisfactory agreement between the numerical and experimental results (Zhao et al. 2023), indicating the reliability of the developed mesoscale discrete element model for simulating the static compression test of SFRC specimens.

Comparison of the failure mode of SFRC with different steel fiber contents from laboratory experiment (Zhao et al. 2023) and DEM simulation

5.1.3 Examination of interface failure processes

In this part, the crack generation algorithm has been modified by employing the fish language to visualize the expansion of micro-cracks and interface failures in different regions. The microcracks generated in the SFRC model are categorized into three regions: microcracks generated by ITZ failure, microcracks generated by matrix failure and microcracks generated by fiber-matrix interface failure (Fig. 18).

Comparison of interface failure process of SFRC0 and SFRC2 under quasi-static compression

Figure 15 displays the micro-cracks within the matrix, micro-cracks attributed to ITZ failure, and micro-cracks arising from fiber-matrix interface failure. Microcracks first initiated from the ITZ part and propagated into adjacent areas of the ITZ and the cement matrix. Only a small number of microcracks were generated at the fiber-matrix interface, and most of the cracks emerged from the external area of the specimen, indicating that the fracture of the specimen initiated from the external area. When the stress in the specimen reached the peak value (\(\varepsilon = 0.2\%\)), most of the ITZ has been destroyed and failed, leading to rapid expansion of micro-cracks within the matrix. Simultaneously, multiple fractures occurred at the fiber-matrix interface, initiating expansion along the interface. The bridging effect of steel fiber and matrix prevents the propagation and evolution of cracks in the center of the specimen. At a strain of 0.25%, the ITZ was essentially fully destroyed. Owing to the bridging effect of steel fibers, a ring-shaped pattern formed within the matrix. The fibers surrounding the specimen lost their bearing capacity due to the interface failure, while the steel fibers in the center of the specimen continue to bear the load, preventing the development of macroscopic cracks and thereby increasing the residual stress after reaching the peak.

5.2 Dynamic compression test

5.2.1 Dynamic compressive strength characteristics

Figure 19 compares the experimental and numerical dynamic compressive stress–strain curves of SFRC0 and SFRC1 at varying strain rates. It can be seen from the diagram that the peak strength and peak strain of SFRC in numerical simulation increase with the increase of strain rate, showing a good strain rate correlation. The numerical curve of plain concrete is in good agreement with the experimental curve. The numerical curve of SFRC1 closely resembles the shape of the experimental curve, but the peak strength at higher strain rates is slightly different. From the comparison of dynamic compressive strength, the dynamic compressive strength of SFRC in numerical simulation also increases with the increase of steel fiber content, which is not consistent with the experimental data. This inconsistency is attributed to the limitations of the experimental techniques employed. Specimens with higher steel fiber content tend to exhibit a greater number of initial pores, whereas our current numerical model is idealized and does not account for these initial pores and specimen damage. Consequently, the influence of steel fiber content on the dynamic compressive strength of the sample is readily apparent. The bridging effect of steel fibers serves to effectively fortify its dynamic compressive strength. The increase in steel fiber content results in a significant improvement in the dynamic compressive strength of the sample, regardless of the strain rate.

Comparison of the dynamic stress–strain curves of SFRC from experiments and DEM simulation a SFRC0 b SFRC1 c Comparison of SFRC dynamic compressive strength with different strain rates and steel fiber content

The dynamic influence factor (DIF) is used to quantify the effect of strain rate on the strength of concrete, which is determined by the ratio of compressive strength of concrete under dynamic and static loads. DIF represents a significant parameter to describe the dynamic mechanical properties of concrete materials. Therefore, the DIF values of different SFRCs at various strain rates were calculated and compared, as depicted in Fig. 10. In general, the DIF obtained through numerical simulation resembles the curve of SFRC0 in the test, with DIF increasing approximately linearly with the increase of strain rate. At the same strain rate, SFRC exhibits a slightly higher DIF than plain concrete, although the disparity remains insignificant, signifying that the inclusion of fibers in concrete does not exert a notable effect on compressive DIF. Therefore, the available DIF formula of plain concrete can be approximated to SFRC to simulate its high strain rate load. Moreover, the DIF resulting from numerical simulation outcomes is minimally influenced by variations in steel fiber content. The strain rate trend of SFRC is similar to that of plain concrete, with the effect of fiber volume fraction being negligible. This conclusion aligns with the previous research conclusion (Sun et al. 2017; Xu et al. 2012). In Fig. 20, a comparison is made between some previous research data (Lok et al. 2024; Tedescos and Ross 1998) with numerical simulation results. The DIF strain rate effect of numerical simulation is basically consistent with the experimental results of Lok and Zhao (2003) and the recommended value of CEB. At higher strain rates (above 50s⁻¹), DIF displays a linear increase with strain rate and remains within a reasonable range of both cases.

Comparison of DIF between numerical results and previous research. a Test data of Zhao et al. (2023); b Previous mathematical models

5.2.2 Assessment of failure characteristics

Figure 21 compares the macroscopic failure modes of the test (Zhao et al. 2023) and numerical simulation of SFRC with different strain rates and steel fiber contents. It is noteworthy that the numerical simulation results of SFRC shows the damage of the horizontal section. It can be seen from the diagram that the failed SFRC with the same fiber content exhibited poor integrity when subjected to higher strain rates. Under compressive load, SFRC0 experiences fragmentation when the strain rate ranged from 40 to 213s⁻¹, and the degree of fragmentation increases progressively. The macroscopic fragmentation diagram of numerical simulation also demonstrates a gradual escalation degree of fragmentation. For steel fiber reinforced concrete, specimens from SFRC2 (at a strain rate of 50s⁻¹) and SFRC3 (at a strain rate of 60s⁻¹), as shown in both the numerical simulation and experimental tests, indicated a high level of integrity with only minimal localized spalling. However, with the increase of strain rate, the central part of SFRC2 and SFRC3 gradually fragmented, but the whole specimen does not appear to be completely broken.

Comparison of failure modes of SFRC with different strain rates and steel fiber contents between experiment and numerical simulation

When considering the same strain rate, SFRC with higher fiber content shows better ductility and integrity in SHPB test. Under the second impact strain rate, the degree of fragmentation of specimens from SFRC1 to SFRC3 significantly reduced. Increasing the steel fiber content enhanced the ability of steel fiber reinforced concrete to withstand deformation under dynamic loading. The numerical simulation is in good agreement with the experimental failure mode, thereby confirming the reliability of the numerical model in predicting the failure mode of SFRC.

Figure 22 shows the deformation characteristics of steel fibers in SFRC with different strain rates and steel fiber contents. The deformation of steel fiber increased as the increase of strain rate. At lower strain rates (e.g., 60s⁻¹ for SFRC2 and 58s⁻¹ for SFRC3), the deformation of steel fibers was minimal, with little to no bending observed. Conversely, at higher strain rates (e.g., 206s⁻¹ for SFRC1 and 203s⁻¹ for SFRC2), the deformation of steel fibers increased significantly, with bending occurring in certain steel fibers surrounding the specimen. With the increase of steel fiber content, the degree of deformation in the steel fibers diminished noticeably, and the steel fibers within the central region of the SFRC3 specimen exhibited almost no deformation. These findings closely align with the deformation results observed in the study conducted by Zhao et al. (2023), indicating that the model is capable of accurately reflecting the deformation characteristics of steel fiber under dynamic compression.

Displacement deformation of steel fiber in numerical simulation

5.2.3 Microcrack propagation process

Figure 23 shows the three microcrack propagation processes of SFRC0 and SFRC2 at varying strains under a strain rate of 160s. Figure 24 shows the change of the number of cracks in SFRC0 and SFRC2 with increasing strain at a strain rate of 175 s. Analysis of the ITZ crack reveals that the initiation and propagation characteristics of ITZ cracks in SFRC0 and SFRC2 are fundamentally similar. However, the propagation speed of ITZ cracks in the central region of SFRC2 was marginally slower than that in SFRC0. This observation suggests that steel fibers possess the ability to partially inhibit the expansion of ITZ cracks towards the specimen's center. From the perspective of matrix crack, the cracks in both SFRC0 and SFRC2 first initiated near the periphery of the specimen and gradually propagated from the exterior towards the interior.

Comparison of interface failure process of SFRC0 and SFRC2 under strain rate of 160

a Crack-strain curve of SFRC2 under 160s⁻¹. b The relationship between crack statistics and fiber content under different strain rates

The development speed of the central crack in SFRC2 was obviously slower than that of SFRC0, resulting a ring-shaped crack morphology, which is similar to the static law. However, the number of matrix cracks in SFRC2 exceeded that in SFRC0, indicating that while the bridging effect of steel fibers somewhat hampers the speed of microcrack development from the periphery towards the center of the specimen, it intensifies the crushing effect around the matrix, enabling the matrix absorb more energy. During the dynamic loading process, the fiber-matrix crack eventually emerged, with the number of cracks initially increasing and then stabilizing. This pattern signifies that the fiber-matrix interface experiences the latest crack initiation, and at higher strain rates, the interface essentially undergoes complete debonding leading to failure.

Figure analysis reveals that the number of ITZ cracks remains consistent and unaffected by variations in steel fiber content, since the ITZ interface as the ITZ interface consistently experiences complete failure at higher strain rates. In contrast, the number of matrix cracks in steel fiber-reinforced concrete is significantly higher than in plain concrete under all three strain rates. However, increasing the steel fiber content does not result in an increase in the number of matrix cracks. Instead, higher steel fiber content primarily leads to more fiber-matrix interfaces, thereby amplifying the occurrence of fiber-matrix cracks due to damage.

From this point of view, the addition of steel fibers to concrete induces a greater presence of microcracks within the matrix vicinity, enhancing the matrix's capacity to absorb energy, as elaborated in detail in Sect. 6.2.5. For steel fiber reinforced concrete under different strain rates, it is obvious that the increase of strain rate deepened the damage of matrix and interface, resulting in a higher number of matrix cracks and fiber-matrix cracks.

5.2.4 The evolution of contact force fields

PFC enables the analysis of the contact force change and fracture development process of the contact force chain, so that the failure process of steel fiber reinforced concrete can be analyzed by examining the evolution process of the contact force field (Yan et al. 2021). Figure 25 shows the evolution process of typical contact force field for SFRC2 and SFRC0 at different strain times under a strain rate of 160s⁻¹. Figure 26 shows the stress of the axial and radial steel fibers of SFRC2 along the loading direction during loading. During the dynamic compression process, the SFRC0 specimen experiences evident compressive stress within a strain range from 0.01 to 0.018. The contact around the specimen gradually disappeared, and the specimen began to detach from the surrounding. When the strain reached 0.022, the contact force field broke, resulting in complete specimen failure. Conversely, the matrix in SFRC2 continued to bear compressive stress. Due to the steel fiber's deformation speed lagging behind that of the matrix, the matrix generated tensile stress on the steel fiber through the fiber-matrix interface.

The contact force field evolution process of SFRC2 under 160 strain rate

Stress state of steel fibers in SFRC2 under dynamic compression

Most of the steel fibers on the horizontal plane experience tensile stress, while the vertically distributed steel fibers, approximately parallel to the loading direction, undergo compression caused by dynamic compression waves. In addition, the compressive and tensile stresses experienced by steel fibers in the edge region surpass those in the central region. This observation suggests that the steel fiber bridging effect effectively inhibits concrete cracking under dynamic loads. From both radial and axial perspectives, when the strain ranges from 0.01 to 0.018, the residual contact of SFRC is significantly greater than that of plain concrete, indicating superior integrity. Some randomly distributed steel fibers withstood the compressive load, while others restrict the matrix's deformation rate by bridging it, delaying internal contact fracture and thereby enhancing concrete ductility. At a strain of 0.022, the central portion of SFRC2 retained a certain level of integrity, whereas SFRC0 has completely failed, aligning closely with existing test results.

Figure 26 shows the ratio of residual contact number to initial contact number after failure under different steel fiber content and strain rate, which can reflect the integrity of the specimen after failure. Regardless of the strain rate, the contact number ratio of plain concrete is always significantly lower than that of steel fiber reinforced concrete. With the increasing of strain rate, the ratio of the residual contact number to the initial contact number of SFRC gradually decreases, and the fragmentation degree of the specimen increases. At low strain rates, the contact number ratio of SFRC is very high, especially SFRC3 reaches nearly 90%, which indicates that steel fiber can effectively prevent the cracking of concrete matrix at low strain rates, so that the whole specimen shows a high degree of integrity. At the same strain rate, the ratio of residual contact number to initial contact number increases obviously with the increase of steel fiber content. The ratio of residual contact number to initial contact number of SFRC1, SFRC2 and SFRC3 increased by about 14%, 24% and 31% on average compared with SFRC0. The existence of steel fiber hinders the expansion of macroscopic cracks to the central area of the specimen, so that the integrity of the central part is higher.

5.2.5 Energy analysis

The PFC discrete element facilitates the tracking of alterations in energy storage and release during the simulation process, thereby significantly enhancing the understanding and analysis of experimental outcomes (Luo et al. 2020).

Strain energy (abbreviated as \(E_{\text{c}}\)) is defined as the energy stored in the springs. The total normal contact force results from the amalgamation of a linear spring contribution. The formula is as follows:

where \(k_{\text{s}}\) is the shear stiffness, \(k_{\text{n}}\) is the normal stiffness, \(N_{\text{c}}\) is the contact number, \(F_{\text{s}}^{\text{l}}\) is the linear partial shear force, \(F_{\text{n}}^{\text{l}}\) is the linear partial normal force.

Bond strain energy (abbreviated as \(E_{\text{b}}\)) is defined as the cumulative strain energy across all parallel bonds. The formula is calculated as follows:

In the formula: \(N_{\text{b}}\) represents the number of parallel bonds, \(A\) represents the cross-sectional area of parallel bond, \(\overline{I}\) is the moment of inertia of cross-section, \(\overline{J}\) is the polar moment of inertia of the parallel bond cross-section, \(\overline{F}_{\text{n}}\) is the parallel bond normal force, \(\overline{F}_{\text{s}}\) is the parallel bond tangential force, \(\overline{M}_{{\text{t}}}\) is the twisting moment, \(\overline{M}_{\text{b}}\) is the bending moment.

Slip energy (abbreviated as \(E_{\mu }\)) represents the energy expended due to sliding friction. The calculation formula is expressed as follows:

where \((F_{{\text{s}}}^{\text{l}} )_{\text{o}}\) is the linear shear force at the beginning of the timestep, and the adjusted relative shear-displacement increment has been decomposed into an elastic \(\Delta \delta_{\text{s}}^{\text{k}}\) component and a slip \(\Delta \delta_{\text{s}}^{\mu }\) component.

The dashpot energy (abbreviated as \(E_{\beta }\)) is defined as the total energy dissipated by the dashpots. The formula is as follows:

where \(\dot{\delta }\) is the relative translational velocity, \(F_{\text{d}}\) is the dashpot force.

Kinetic energy (abbreviated as \(E_{\text{k}}\)) is defined as the energy expended by both translational and rotational motions of all particles. The formula stands as follows:

In the formula:\(\varsigma_{i}\) is the generalized mass,\(v_{i}\) is the generalized velocity.

Including: \(E_{\text{c}}\), \(E_{\text{b}}\), \(E_{\mu }\) and \(E_{\beta }\) constitute the total energy stored in contact. And \(E_{\text{k}}\) is the energy of particles.

For this model, five types of absorbed energy are defined: total absorbed energy, ITZ absorbed energy, matrix absorbed energy, fiber-matrix interface absorbed energy and fiber strain energy. The total absorption energy is the sum of all the contact energy and the kinetic energy of the particles in the model. The energy absorbed by ITZ is the sum of the total energy in contact with the matrix. Matrix energy is the sum of the total energy of cement matrix contact and the kinetic energy of cement particles. Fiber-matrix energy is the sum of the total energy of the contact between steel fiber and matrix. Finally, the fiber strain energy includes the deformation energy of the internal contact of the steel fiber. For plain concrete, the total absorbed energy includes ITZ absorbed energy and matrix absorbed energy. For steel fiber reinforced concrete, the total absorbed energy includes ITZ energy, matrix energy, fiber-matrix energy and fiber strain energy.

Figure 27 shows the energy-time curves of SFRC0, SFRC1, SFRC2 and SFRC3 at strain rate of 170. Figure 28 depicts the peak curves of various energies with different steel fiber contents. It can be seen from the figure that the total energy absorbed by SFRC is significantly more than that of plain concrete, and the absorbed energy of SFRC is proportional to the content of steel fiber, which is consistent with the conclusion of Wang et al. (2022). The addition of steel fiber to concrete can increase the energy absorption threshold of the matrix by approximately 40%. However, increasing the steel fiber content does not substantially enhance the energy absorption capacity of the matrix. Irrespective of the steel fiber content, the ITZ interface consistently experiences complete failure, resulting in a similar absorption energy for the ITZ. The increase of energy absorption of steel fiber reinforced concrete with higher fiber content primarily attributes to more interface damage absorption energy and steel fiber deformation energy. From the microscopic point of view, a higher steel fiber content results in a greater number of fiber-matrix interfaces. During dynamic loading process, most of these interfaces debond and fail to absorb additional energy. Moreover, a higher fiber content enables a greater number of steel fibers to exert their deformation capabilities, resisting specimen failure and allowing for increased elastic–plastic deformation of the steel fibers themselves, thereby absorbing more energy. Therefore, within the allowable range of construction conditions, increasing the steel fiber content in shotcrete steel fiber reinforced concrete linings for underground engineering can enhance the energy absorption effectiveness of SRFC. It may support structures to a certain extent and reduce the occurrence of dynamic disasters in underground engineering.

The ratio of residual contact number to initial contact number under different steel fiber content and strain rate

Energy-time curve of SFRC under the same strain rate. a SFRC0 b SFRC1 c SFRC2 d SFRC3

Afterwards, the enhancement mechanism of SFRC on the dynamic compressive properties of concrete from the relationship between stress accumulation, energy absorption and crack propagation was analyzed. Additionally, the absorbed energy was defined for different energy types. The combined strain energy and bond strain energy are considered as the total strain energy. Broken energy encompasses dissipation energy resulting from crack generation, particle sliding friction, and damping dissipation. Broken energy corresponds to the summation slip energy and dashpot energy. The kinetic energy refers to the energy carried by the fragments.

Figure 29 shows the energy-crack-stress curves of SFRC0 and SFRC2 specimens defined according to different absorbed energy, and Fig. 30 shows the energy-time curves defined according to different energy partitions. Combining the two curves reveals that energy absorption precedes crack formation within the specimen, predominantly stemming from the elastic strain energy stored in the specimen. When the stress reaches the peak, the number of cracks increases sharply, resulting in a significant rise in absorbed energy. Ultimately, as the stress in SFRC0 approaches zero, crack formation ceases, and energy absorption terminates. However, due to the load-bearing capacity of the steel fibers, SFRC2 retains residual stress, gradually diminishing the generation of new cracks and leading to a gradual cessation of energy absorption. SFRC2 exhibits significantly higher total strain energy and broken energy compared to SFRC0, while the kinetic energy carried by SFRC2 is significantly lower than that of SFRC0 (Fig. 31).

The relationship between absorbed energy and steel fiber content under the same level strain rates

Stress-crack-energy evolution curve. a SFRC0 b SFRC2

Energy-time curve of three types. a SFRC0 b SFRC2

Under the influence of external dynamic loads, deformation energy progressively accumulates within the concrete specimen, leading to a gradual increase in stress. Upon reaching the energy storage limit, characterized by the deformation limit, the strain energy stored within the specimen is rapidly released, triggering the generation of cracks within the matrix that absorb a substantial amount of energy. This process continues until the specimen is completely fractured, losing its load-bearing capacity and ceasing energy absorption and dissipation. In SFRC, the presence of steel fibers decelerates the deformation rate of the matrix, prolonging the time required for the matrix to reach its elastic–plastic deformation limit (the stress peak of SFRC2 shifts to the right compared to SFRC0). This, to some extent, improves the energy storage limit of the matrix. In addition, the steel fibers themselves store a portion of the deformation energy. These combined effects result in increased deformation energy stored in SFRC, manifesting as a higher stress peak.

When the matrix in SFRC reaches its energy storage limit, more strain energy stored in the matrix is released, leading to rapid crack initiation in the surrounding region. However, due to the bridging effect of steel fibers, macroscopic cracks struggle to propagate from the surrounding area to the central area. Instead, an increased number of microcracks are generated in the surrounding area of the matrix, resulting in greater fragmentation. Furthermore, the presence of more microcracks and friction between fragments leads to increased energy dissipated by the matrix. Moreover, because the steel fiber limits the deformation and failure of the central region of the matrix, and the interface between the steel fiber and the matrix absorbs a portion of the energy, the total volume of the fragments produced by SFRC is reduced, thereby reducing the kinetic energy carried by SFRC.

This paper constructs a novel three-dimensional mesoscale model of steel fiber reinforced concrete (SFRC) in tunnel employing the discrete element method. The model comprises a concrete matrix, aggregate, interfacial transition zone (ITZ), and steel fiber. It considers the stochastic shape of the coarse aggregate and the random distribution of steel fibers, effectively accounting for the microscopic interactions among the coarse aggregate, steel fibers, and matrix. The microscopic parameters necessary for the model are obtained through a comprehensive procedure for calibrating detailed parameters, and the model's reliability is verified by comparing it with experimental findings. Through the execution of static compression and SHPB dynamic compression tests, the mechanical properties of SFRC are simulated and investigated. The mechanical properties and the mechanism of crack arrest in SFRC under both static and dynamic compression loads are analyzed at the microscopic level. The primary conclusions are summarized as follows:

1. (1)

   In the simulation of quasi-static compression test, the compressive strength of SFRC exhibits an increase with higher steel fiber content. Compared to SFRC0, SFRC1, SFRC2 and SFRC3 display respective increases of 12.6%, 26.7% and 38%, respectively. At the same time, a higher steel fiber content leads to SFRC specimens retaining greater integrity after failure, accompanied by enhanced ductility. The bridging effect of steel fibers and matrix prevents the propagation of macroscopic cracks from the periphery towards the center of the specimen, thereby inhibits deformation and failure while improving ductility.

2. (2)

   In the simulation of dynamic compression test, the dynamic compressive strength of SFRC continues to rise with increasing steel fiber content. Furthermore, steel fiber reinforced concrete exhibits an evident strain rate effect, with the dynamic compressive strength increasing alongside the strain rate. The numerical simulation results also perform the calculation of the dynamic increase factor (DIF) of SFRC specimens. It is found that the compressive DIF of SFRC is only slightly larger than that of plain concrete, and the addition of fiber content within the range of 1–3% has negligible effect on the DIF of SFRC.

3. (3)

   Based on the micro-crack propagation and micro-force field evolution process, the dynamic compression experimental results are analyzed from the microscopic point of view. It is determined that ITZ cracks and matrix cracks are generated almost simultaneously, while crack initiation at the fiber-matrix interface occurs later. At higher strain rates, the ITZ interface consistently experiences complete failure. Although the bridging effect between the steel fiber and matrix consistently experiences complete failure speed of the micro-cracks from the surroundings towards the center of the specimen, it results in the generation of more micro-cracks in the matrix vicinity. The presence of steel fibers hinders the propagation of macrocracks towards the central region of the specimen, thus enhancing the integrity of the central part. The presence of steel fibers can be construed as effectively inhibiting deformation failure in concrete under dynamic loading.

4. (4)

   Considering the stress and deformation of steel fibers, it is found that vertical fibers aligned with the loading direction in SFRC primarily endure compressive stress, while the fibers located on the horizontal plane mainly experience tensile stress. The stress experienced by steel fibers in the edge region significantly exceeds that in the central area. At higher strain rates, the steel fibers surrounding the specimen may undergo bending deformation, and the overall deformation of steel fibers decreases with an increase in steel fiber content.

5. (5)

   In terms of energy perspective, the addition of steel fiber to concrete enhances its energy absorption capability, thereby improving toughness. The inclusion of steel fiber leads to approximately a 40% increase in the dynamic load energy absorbed by the matrix. However, the increase in steel fiber content has no significant impact on the energy absorption of the matrix. Nevertheless, a higher steel fiber content results in more damage to the fiber-matrix interface and increased steel fiber deformation, thereby raising the total energy absorbed by steel fiber reinforced concrete.