The concrete microplane damage model (Caner and Bažant, 2013a; Caner and Bazant, 2013b) consists of three parts: physical mapping of “macro to micro physical variables,” establishment of constitutive relationship at the micro scale, and homogenization of “micro to macro physical variables,” as shown in Figure 1.

The microplane model portrays the mechanical behavior of concrete materials at the microscopic level in terms of stresses and strains in vector form, so it is necessary to transform the stresses or strains at the macroscopic level into the stresses or strains at the microscopic level. According to the treatment of the relationship between the macroscopic stress tensor or macroscopic strain tensor and the stress or strain vector on the microplane, they are usually categorized into static and kinematic constraints. They can be understood as the projection of the macroscopic stress tensor on the microplane to obtain the corresponding stress vector and the projection of the macroscopic strain tensor on the microplane to obtain the corresponding strain vector, respectively. Due to the strain-softening behavior of quasi-brittle materials such as concrete, the kinematic constraints shown in Figure 2 are used in the concrete microplane damage model to ensure the stability of the model when analyzing strain softening.

As shown in Figure 2, the strain vector e on the microplane i (normal vector is denoted as n ) is expressed as the projection of the macroscopic strain vector ϵ by Equation 1. e = n ⋅ ϵ = e N n + e T (1) where e N = ϵ N = N : ϵ is the normal strain component on the microplane, e T = ϵ T = T ⋅ ϵ is the tangential shear strain vector on the microplane, e T = ϵ K 2 + ϵ M 2 ; N and T are the corresponding projection operators. Further, the tangential strain e T in the microplane is used to characterize plasticity and friction; the normal strain e N is distinguished into the tensile strain (i.e., the part where e N > 0 ) to characterize the tensile ability, and the compressive strain (i.e., the part where e N < 0 ) to characterize the compressive ability, while the compressive strain can be decomposed into the volume component e V and the deviatoric component e D by Equation 2. e N = e V + e D = V : ϵ + D : ϵ (2) where V = 1 3 I is the volume component, D = n ⊗ n − 1 3 I denotes the bias obtained by subtracting the volume component from the total strain, and I is the third-order unit tensor.

After defining the stress-strain relationship on the microplane, the principle of virtual work is applied to establish the equation between the microplane stress vector and the macroscopic stress tensor, from which the macroscopic stress tensor σ i j satisfies Equation 17. 2 π 3 σ i j δ ε i j = ∫ Ω σ N δ ε N + σ L δ ε L + σ M δ ε M d Ω (17) where Ω is the surface of the unit hemisphere and 2 π / 3 is the volume of the hemisphere, the integral is regarded as a homogenization of the microplane contributions in different directions within the material. Due to δ ε N = N i j δ ε i j , δ ε L = L i j δ ε i j and δ ε M = M i j δ ε i j , from the arbitrariness of the variation δ ε i j , the equilibrium relation can be obtained as Equations 18, 19. σ i j = 3 2 π ∫ Ω s i j d Ω ≈ 6 ∑ μ = 1 N m w μ s i j μ (18) s i j = σ N N i j + σ L L i j + σ M M i j , (19)

The integral is approximated by the optimal Gaussian integral formula for the sphere, denoting the weighted sum of microplane in the direction n μ , with the weight w μ normalized, so ∑ μ w μ = 1 / 2 . For better accuracy in the far-peak post-softening, 37 microplane are recommended to be preferred.

The simplest way to extend the plain concrete microplane model to make it applicable to SFRC containing materials such as steel fibers is to adjust the stress-strain boundary conditions on the microplane (Caner et al., 2013). First, fibers usually increase the tensile capacity of the concrete material (Jiang et al., 2023; Lakavath et al., 2024), so the effect of steel fiber reinforcement needs to be introduced on the normal tensile stress-strain boundary. The contribution of fiber reinforcement is given by a simplified form of the Kholmyansky equation (Kholmyansky, 2002), as Equation 20. σ N f = E p 1 k 1 ϵ N / k 1 e − p 2 ϵ N / k 1 , if  ϵ N / k 1 < 1 / p 2 + p 4 . E p 1 k 1 / p 2 e − 1 , if  1 / p 2 + p 4 ⩽ ϵ N / k 1 < p 3 . E p 1 k 1 ϵ N / k 1 − p 3 + 1 / p 2 e − p 2 ϵ N / k 1 − p 3 + 1 / p 2 , if  p 3 ⩽ ϵ N / k 1 . (20)

This contribution is obtained by the gradual activation of the bridging of the fibers during the development of the crack, as shown in Figure 4. Assuming parallel coupling of the fibers and the matrix, the normal stress on the microplane is Equation 21. σ N b f = σ N b + σ N f . (21) where σ N b f is the total normal stress of the fiber-reinforced concrete, σ N b is the normal stress of the plain concrete matrix, and σ N f is the normal stress of steel fiber contribution given by Equation 20. Therefore, in the microplane model of fiber-reinforced concrete, Equation 12 is modified as Equation 22. σ N = max min σ N e , σ N b f , σ V b + σ D b . (22)

In addition, the tensile capacity that steel fiber-reinforced concrete can withstand changes (increases or decreases) before cracks develop. Therefore, the normal tensile stress-strain boundary Equation 8 of the plain concrete matrix needs to be adjusted to Equation 23. β 1 = − c 1 + c 17 e − c 19 ϵ e − c 18 + p 0 V f . (23)

Secondly, the addition of steel fibers changes the compressive capacity of the concrete material, especially the shear expansion deformation, so the compressive deviatoric stress-strain boundary condition Equation 9 needs to be adjusted to Equation 24. β 2 = c 5 γ 1 + c 7 V f . (24)

Finally, steel fibers change the mechanical behavior of concrete in triaxial compression, so the plastic-friction boundary Equation 11 needs to be adjusted to Equation 25. σ T b = F T − σ N = E T k 1 k 2 V f c 10 V f − σ N + σ N 0 E T k 1 k 2 V f + c 10 V f − σ N + σ N 0 , (25) where σ N 0 = E T k 1 c 11 V f − c 12 V f ϵ V .

The microplane model of steel fiber-reinforced concrete is based on the microplane model of plain concrete. Therefore to calibrate its parameters, it is necessary to calibrate the parameters of plain concrete first. Then from the uniaxial tensile data of SFRC, p 0 V f , p 1 V f , p 2 V f , p 3 V f , p 4 V f can be identified. From Equation 23, it can be found that the parameter p 0 V f determines whether the fibers in the concrete matrix are bonded or not. The parameter p 1 V f controls the proportion of steel fiber contribution, while p 2 V f controls the pullout of fibers connecting the open cracks. Parameters p 3 V f , p 4 V f determine the length of the plastic plateau in the stress-strain relationship (shown in Figure 4). c 7 V f , c 12 V f can be determined from the uniaxial compression data and k 2 V f , c 10 V f , c 11 V f can be determined from triaxial compression data. In many practical applications, fewer triaxial compression problems are involved. Therefore, it is mostly sufficient to study the uniaxial compression and uniaxial tension of fiber concrete. If uniaxial tension tests are difficult to perform, the parameters can be calibrated indirectly by means of the notched three-point loaded bending concrete beam test.

In the concrete microplane damage model, taking the stretch shown in Equation 5 as an example, the damage behavior of concrete under several cycles can be described by assuming that the elastic modulus E N in the μ -th microplane will undergo damage by Equation 26. E N μ E N 0 = exp − A ϵ N 0 + μ . (26) where E N 0 is the undamaged modulus, A is the material parameter which is called c 13 in Equation 5, and ϵ N 0 + μ is the maximum strain reached so far in the μ -th microplane. This approach can characterize the damage to the concrete for a few cycles or so, but it does not capture fatigue. It is worth noting that once damage begins to occur in the material, strain localization occurs, and strain calculations at this point are affected by the finite element mesh size in a way that is difficult to ignore. Therefore, it is necessary to analyze the microplane model in conjunction with nonlocal theory or the crack band model (Bažant and Oh, 1983) (used here). In order to effectively predict the fatigue response, it is necessary to introduce a variable that measures the fatigue damage history paths (Kirane and Bažant, 2015; Baktheer et al., 2021; Aguilar et al., 2022).

In the numerical implementation, it is necessary to incorporate the crack band model to reduce the mesh size sensitivity of the computed results (Bažant and Oh, 1983; Červenka et al., 2005). In ABAQUS commercial finite element software, VUMAT subroutine is written for numerical implementation. If the stress tried exceeds the stress-strain boundary, it is necessary to keep the strain and limit the stress to the stress-strain boundary. In the finite element method program, if the strain increment ∆ ε i j at the current step is known, as well as the strain ε i j 0 and stress σ i j 0 at the end of the previous step, the new stress σ i j is obtained at the end of the current step by the following steps:

Step 1

First, the strain and strain increment on the microplane are calculated by Equation 33 according to Equation 1, ∆ ε N = N i j ∆ ε i j , ∆ ε M = M i j ∆ ε i j , ∆ ε L = L i j ∆ ε i j . (33)

Step 2

Calculate the volumetric strain and its increment at the end of the previous and current steps, based on the given strain and its increment by Equation 34. ε V o = ε k k / 3 , ∆ ε V = ∆ ε k k / 3 , ε V = ε V o + ∆ ε V . (34)

Then, calculate ε e = − σ V 0 / E N 0 , while the corresponding volumetric stress-strain boundaries σ V b are later calculated by Equation 10. The deviatoric strain can be obtained from total strain and volumetric strain by Equation 35. ∆ ε D = ∆ ε N − ∆ ε V , ε D 0 = ε N − ε V 0 , ε D = ε D 0 + ∆ ε D , (35) and σ D b is calculated by Equation 9.

Step 3

Calculate the ε N = ε V + ε D retrieve the history variables ε N 0 + , ε N 0 − of the largest magnitude stored so far. Next, the current damage degree of the normal elastic modulus of the microplane is calculated according to E N 0 = E / 1 − 2 ν by Equations 5, 6. Calculating the normal elastic stress of the microplane by Equation 36. σ N e = σ N o + E N ∆ ε N (36)

Step 4

Calculate the normal tensile boundary σ N b of the microplane using Equation 8 when σ N b > 0 .

Step 5

Holding the strain constant and allowing the stress to drop vertically to the normal tensile stress-strain boundary by Equation 37. σ N = max min σ N e , σ N b , σ V b + σ D b (37)

Meanwhile, update the history maxima ε N 0 + and ε N 0 + .

Step 6

Calculate an approximation of the current volumetric stress by Equation 38. σ V ≈ 1 2 π ∑ μ = 1 N m w μ σ N (38)

Retrieve the originally stored microplane shear stresses σ L o , σ M o , followed by estimating E T according to Equation 4, and σ ^ N 0 = E T k 1 c 11 − c 12 ε V . If we want to simulate cell failure, the recommended deletion criterion of the adopted cell is ε I ≥ 0.5 % . The shear boundary is then calculated by Equation 11.

Step 7

Calculate the shear response upon return to the stress-strain boundary by Equations 12–16.

Step 8

The stress σ i j at the end of the current step is obtained by calculating the sum of stresses on all microplanes through Equations 18, 19, while recording the variables σ V , σ N , σ L , σ M , ε N 0 + , ε N 0 − at the end of the current step.

In the microplane model (M7), the shape of the response curve is determined by five free parameters and twenty-one fixed parameters. Tables 1, 2 provide a concise explanation of the meaning and default value of each parameter. The fixed parameters, as detailed in Table 1, are calibrated using the uniaxial compressive strength f c ′ = 36 MPa and the axial normal strain ε p = 0.0036 at peak stress. The calibration of the uniaxial compressive strength f c * , and the corresponding strain ε p * , of a particular concrete using a microplane model requires only that the reference values of the free parameter k 1 , and the elastic modulus E , be modified to, k 1 * = k 1 ε p * ε p , E * = E f c * f c ′ ε p ε p * . (39)

To optimize the fitting of a large amount of experimental data, it is not necessary to change all five free parameters simultaneously in the actual microplane model (Caner and Bažant, 2013a; Caner and Bazant, 2013b). The elastic modulus E and Poisson’s ratio v are determined through experimental measurement. First, all other parameters are assumed to be reference values, and subsequent parameter calibration is conducted. The parameter k 1 was calibrated based on the strain corresponding to the peak stress in the uniaxial compression experiment. If sufficient triaxial compression data are available, and the compression intensity is sufficiently strong to elicit an almost plastic response, the data are fitted by adjusting k 2 . The parameters k 3 and k 4 are calibrated based on hydrostatic pressure experimental data. If sufficient uniaxial, biaxial, and triaxial compression data at low hydrostatic pressure are available, k 5 should be determined by fitting. Otherwise, the default values are retained.

In this section, we consider a notched plain concrete beam of depth D = 107.8 mm , span S = 2.5 D , and a notch of length a 0 = D / 6 situated at the center of the beam. The width of the beam is 38.1 mm , as illustrated in Figure 9A. The parameters of plain concrete material: Young’s modulus E = 38.3 GPa , Poisson’s ratio ν = 0.18 , compressive strength f c = 90.3 MPa ,density ρ = 2400 kg / m 3 . The free parameters k 1 = 140 × 1 0 − 6 are calibrated by matching the peak loads of monotonic load, and the default values are adopted for all other parameters. Two finite element mesh sizes were considered, i.e., h size = 3 mm and h size = 5 mm , and the predicted peak loads were 5.32 MPa and 5.22 MPa , which were in good agreement with the experimental results 5.4 MPa (Bazant and Schell, 1993), as shown in Figure 9B. To reduce the computational cost, the finite element mesh size of the potential damage region was 5 mm in the subsequent analysis.

Next, the fatigue simulation of a three-point bending beam was performed. Cyclic loading was applied up to 84 % of the peak monotonic loading. The parameters of the fatigue damage law of the material are adjusted until the life prediction is satisfactorily close. The fatigue damage metric used was Equation 31, and with the adopted parameters s = 1800 and q = 2 , the model predicted failure after 225 cycles in excellent agreement with the tested fatigue life of 212 cycles. In the last cycle, the failure is characterized by a sudden increase in the overall deformation, as shown in Figure 9C.

Then, the fatigue damage behavior of steel fiber-reinforced concrete is analyzed. A three-point bending beam (Qing et al., 2023) with notched specimen size 440 × 100 × 100 m m 3 , a notch of width 4 mm is located at the bottom of the middle of the beam, the height of the notch is 40 mm , and the span at the bottom of the beam is 400 mm . Four groups of steel fiber-reinforced concrete (SFRC) with steel fiber volume admixture V f = 0 % , 0.5 % , 1.0 % , and 1.5 % were used. Fiber-reinforced concrete material parameters: elastic modulus E 0 = 27.738 GPa , Poisson’s ratio ν = 0.18 , tensile strength f t = 2.67 MPa , compressive strength f c = 34.39 MPa . The model’s free, fixed, and fiber parameters were calibrated, and the calibrated parameters are shown in Table 4. The peak monotonic loading predicted by the model is similar to the test, as shown in Table 5. The load capacity of the concrete beams was enhanced with the increase in the volume admixture of fiber-reinforced concrete.

Next, fatigue analysis of fiber-reinforced concrete three-point bending beams was performed. The maximum applied cyclic load was 85% of the monotonically loaded peak load. The fatigue damage law parameters were adjusted until the life prediction was close enough. The strains of SFRC three-point bending beam with 0.5% volume steel fiber admixture for the one cycle before fatigue loading failure are shown in Figure 10A, and the strains for the remaining three groups are similar. The fatigue damage metric used is Equation 31, and the calibrated fatigue damage law parameters are shown in Table 5 for various fiber-reinforced concretes. The corresponding material fatigue damage law is shown in Figure 10B. The fatigue resistance of concrete is gradually improved with the increase of steel fiber mixing, and the fatigue resistance improvement is gradually slowed down when the steel fiber admixture is more than 1%. Therefore, the present model can be well used to analyze the mechanical response of fiber-reinforced concrete under monotonic and fatigue loading.

Fatigue damage modeling ultimately aims to fine-tune the analysis of the whole process of fatigue damage in concrete structures. This section simulates the fatigue damage of reinforced concrete beams using the proposed model as an illustrative study. The cross-section height of the concrete beam is h = 400 mm , the cross-section width is b = 300 mm , the protective layer thickness of the steel reinforcement is 40 mm , and the cross-section reinforcement ratio of the longitudinal reinforcement ρ ≈ 1.26 % . The geometry and loading of the concrete beam are shown in Figure 11. To show the potential of the model for capturing the fatigue damage behaviors of reinforced concrete structures, the steel reinforcement is assumed to be ideally elastic-plastic with a yield strength of f y = 400 MPa , and the material parameters of the SFRC are the same as in the previous section.

First, a numerical simulation of the four-point bending of reinforced concrete beams ( V f = 0 % , 0.5 % , 1.0 % , 1.5 % ) under monotonic loading is carried out. The load-displacement curves are shown in Figure 12. The ultimate load capacity of reinforced concrete increases gradually with the increase of fiber volume admixture, and the ultimate load capacity of reinforced concrete beams with a fiber volume admixture of 1% is the largest. Unlike the damage mechanics intrinsic model, the fatigue damage of a structure is usually accounted for by the development of strain in the microplane theory. It is worth noting that the bottom strain of plain reinforced concrete beams reaches up to 0.4361% at a displacement loading of 3 mm, and the stressed portion of the tensile zone gradually changes from concrete to steel reinforcement, as shown in Figure 13. With the increase in fiber volume admixture V f = 0.5 % , 1.0 % , 1.5 % , the maximum strains of the resinforced concrete beams were reduced to 0.2315%, 0.0927%, and 0.065% with displacement loading of 3 m. At a displacement loading of 7.5 mm, the steel reinforcement basically stressed the tensile zone of the concrete beams, and there were some areas where the strains were close to or more than 1%. However, the concrete structure would not be subjected to such large deformations under service conditions. Afterward, we will pay attention to the response of the concrete when cracked or just cracked. In this paper, the load capacity corresponding to a plain reinforced concrete beam subjected to a displacement loading of 3 mm is defined as F u = 124.72 kN , and the maximum value of the fatigue load magnitude is 0.9 F u .

Next, fatigue damage simulations of reinforced concrete structures ( V f = 0 % , 0.5 % , 1 % , 1.5 % ) under 200 cyclic loadings were performed. During cyclic loading, the fatigue history variables are gradually accumulated, resulting in gradual degradation of the material stiffness and gradual development of the strain in the concrete beams. The strain distribution after 200 loading cycles is shown in Figure 14. For a reinforced concrete beam composed of plain concrete with a fatigue life of 176 cyclic loadings, the strains at one time before structural failure are shown in Figure 14A. The strains at the bottom of the beam are more than 1% in most regions. With the increase in fiber dosage, the strain of the concrete structure improved during fatigue loading due to the improvement in tensile capacity, and the cracking pattern was transformed. The beam strain distribution after 200 loading cycles is shown in Figures 14B–D. The maximum strains of reinforced concrete beams with steel fiber volume admixture V f = 0.5 % , 1 % , 1.5 % are 0.1313%,0.1231%, and 0.1194%, respectively. The displacements at the loading points also show a significant reduction trend, as shown in Figure 15. Therefore, the microplane model considering the material stiffness degradation can be used to reflect the fatigue damage behavior of plain and fiber-reinforced concrete materials.

It is worth noting that for the fatigue of reinforced concrete structures, to capture the various mechanical behaviors in the tests, the slip between the concrete and the reinforcement also needs to be considered, and the degradation of the fatigue properties of materials such as reinforcement and fibers needs to be considered. However, this is not the focus of this study and will be illustrated in future studies. This model is expected to be used for the whole-process analysis of fatigue damage of plain concrete and steel fiber-reinforced concrete structures under complex loading conditions and structural forms, which will facilitate engineering design, evaluation, and optimization.