Equivalent Fatigue Constitutive Model Based on Fatigue Damage Evolution of Concrete

Abstract

Concrete structures such as bridge decks and road pavements are subjected to repetitive loading and are susceptible to fatigue failure. A simplified stress–strain analysis method that can simulate concrete behavior with a sound physical basis, acceptable prediction precision, and reasonable computation cost is urgently needed to address the critical issue of high-cycle fatigue in structural engineering. An equivalent fatigue constitutive model at discrete loading cycles incorporated into the concrete damaged plasticity model (CDPM) in Abaqus is proposed based on fatigue damage evolution. A damage variable is constructed from maximum fatigue strains, and fatigue damage evolution is described by a general equation whose parameters’ physical meaning and value range are identified. With the descending branch of the monotonic stress–strain curve as the envelope of fatigue residual strength and fatigue damage evolution equation as shape function, fatigue residual strength, residual stiffness, and residual strain are calculated. The equivalent fatigue constitutive model is validated through comparison with experimental data, where satisfactory simulation results were obtained for axial compression and flexural tension fatigue. The model’s novelty lies in integrating the fatigue damage evolution equation with CDPM, explicitly explaining performance degradation caused by fatigue damage. The proposed model could accommodate various forms of concrete constitution and fatigue stress states and has a broad application prospect for fatigue analysis of concrete structures.

1. Introduction

Cyclic loading, a significant factor in the performance degradation of structures like highway bridges, road pavements, railway track slabs, and wind turbine foundations, induces a progressive, permanent, and internal structural change known as fatigue damage [1]. In concrete, the initiation and propagation of microcracks in the cement matrix, mortar, and interfacial transition zones are crucial causes of performance changes due to fatigue action [2]. Structural failure of concrete could occur when microcracks converge into macrocracks of the cross-section. As stress levels in concrete structures increase, fatigue is vital in affecting concrete’s durability, crack development, and even failure. This phenomenon underscores the increasing importance of investigating concrete fatigue in structural engineering.

Fatigue testing has been the traditional method to investigate the fatigue performance of concrete structures. However, the cost of structural component experiments is exceptionally high, and fatigue test results are usually restricted by experimental requirements such as structural dimensions, loading capacity, and boundary conditions. In recent years, scholars have attempted to propose concrete constitutive models under cyclic loading and combine them with numerical simulation to describe concrete’s stress–strain evolution and failure behavior [3]. Models based on plasticity theory and damage mechanics have been developed to evaluate the evolution of concrete’s fatigue damage [4,5,6]. These theoretical models could accurately simulate concrete’s observed behavior; however, predicting engineering structures’ fatigue performance requires software redevelopment. The intermediate parameters required to derive theoretical models are difficult to measure through traditional experiments. Employing concrete constitutive models available in commercial software for more accessible applications to practical and complex components is desirable.

The constitutive models that incorporate both damage and plasticity have successfully captured the nonlinear behavior of concrete [7]. Among them, the concrete damaged plasticity model (CDPM) proposed by Lubliner et al. [8] has proven to be particularly effective in reproducing various failure modes, including tensile cracking and compressive crushing. CDPM has been seamlessly integrated into the widely used general-purpose finite element software Abaqus [9]. Its versatility and robustness in successfully simulating monotonic, cyclic, and dynamic behavior is impressive. However, it has been reported that CDPM underestimates the fatigue resistance of concrete structures in some cases [10]. The drawback of high computational cost in a cycle-by-cycle simulation further underscores the urgent need for a simplified analysis method.

It is recognized that the fatigue problem experienced by most bridge engineering structures is high-cycle fatigue, where the stress level is relatively low, and the number of cycles before a fatigue failure is longer than 10⁴ cycles. Given the minimal fatigue damage induced in a single loading cycle for high-cycle fatigue problems, a simplified fatigue analysis method that overcomes the high computational cost and complicated simulation process associated with cycle-by-cycle fatigue analysis is practically possible. The cycle-jump acceleration algorithm [11,12,13], a method that accelerates the simulation process by skipping some cycles based on the assumption of gradual damage accumulation, can significantly improve the efficiency of the finite element method in simulating the evolution of concrete fatigue damage. Alternatively, the fatigue process can be discretized and analyzed at prescribed cycle ratios to express fatigue damage evolution in a standardized form, making it more convenient for engineering applications.

Zanuy et al. [14] proposed a simplified concrete fatigue model for axial compression. Bai et al. [15] conducted numerical research on the fatigue performance of wind turbine concrete foundations. They simulated how performance degradation of concrete caused by fatigue loads led to stress redistribution along the cross-section by combining several actual cycles into one equivalent virtual cycle. Some scholars found that the stress–strain curve of concrete under fatigue loading followed a similar form as that under monotonic loading and derived the equivalent fatigue damage model by substituting parameters such as stiffness and strength after specific loading cycles into the monotonic stress–strain relationship of concrete [10,16,17,18]. The above models discretize the fatigue analysis process and conduct concrete fatigue performance analysis in a jump-in-cycle manner. From an empirical and phenomenal perspective, this method is more convenient for engineering applications. However, these models have limitations in providing a physical explanation for fatigue damage evolution, and further research is crucial to understand the definition of fatigue damage variable and its evolution better.

The progressive accumulation of internal damage causes concrete performance degradation at the macro level. Macroscopic mechanical properties of concrete evolve as fatigue loading progresses; the evolution of these concrete properties can then represent internal fatigue damage. Different parameters, such as maximum fatigue strain [19], residual strain [20,21], secant modulus of elasticity [22], and residual strength [23], have been defined as concrete damage variables. Maximum fatigue strain (corresponding strain when the upper limit of fatigue loading is reached) reflects the deformation under fatigue loading and is convenient to measure in experiments [24]. Thus, defining concrete fatigue damage variable through maximum fatigue strain is valuable in revealing the damage deterioration mechanism [19]. However, the maximum fatigue strain measured during the first loading cycle is more significant than that generated in subsequent fatigue cycles due to initial defects such as pores and microcracks within the material. Initial damage induced by initial defects and preloading also shows strong dispersion affected by the manufacturing process and loading conditions. Therefore, defining fatigue damage variable with damage exclusively generated by fatigue loading is desirable.

Fatigue strain evolution is generally accepted to exhibit a three-stage development pattern [24,25], consisting of an initial rapid but gradually stabilizing stage, a secondary stabilized stage, and a final rapid growth stage. Multi-segment linear [26], polynomial [27], or exponential equations [28] have been applied to fit the fatigue damage evolution. When obtaining such equations, it is not easy to guarantee clear physical meaning and satisfactory fitting accuracy simultaneously. Because each fitted parameter’s physical meaning and value range in these evolution equations have yet to be identified, the regression equation derived from a specific experiment might apply poorly to other experimental data with a large dispersion of fitted coefficients. The applicability of a regression equation to other experiments with similar loading conditions remains to be verified.

To address a significant need in the field, we aim to develop a simplified, practical, and physically sound method to simulate the fatigue response of concrete based on Abaqus-integrated CDPM. This method should involve a discrete simulation of loading cycles and an equivalent fatigue stress–strain relationship that can reflect fatigue damage evolution. A damage variable is first constructed based on maximum fatigue strain, and parametrical analysis is conducted on the damage evolution equation. The physical meaning of each parameter is clarified, thus overcoming the limitation that fitted equations could only reflect the trend of specific test data. The proposed value range of each parameter and a family of the evolution curves make the fitted fatigue damage evolution equation under various stress levels valuable. A novel definition for several vital parameters, such as residual modulus of elasticity, residual strength, and residual strain, is proposed. The equivalent fatigue stress–strain relationship is subsequently derived by explicitly considering fatigue damage evolution, and CDPM is practically extended to high-cycle fatigue analysis. The proposed equivalent fatigue constitutive model can be applied to any concrete constitution, either empirical formulas commonly recommended by codes or fitted equations directly based on specific experiments. The validity of the proposed fatigue constitutive model is verified by its application to fatigue experiments for concrete under axial compression and flexural tension.

2. Methodology

2.1. Concrete Damaged Plasticity Model (CDPM)

CDPM treats concrete as a continuum material and characterizes its crushing and cracking properties based on isotropic damage. When describing the nonlinear behavior of concrete in tension and compression, CDPM explicitly considers the loss in the material’s modulus of elasticity due to damage, in addition to plastic deformation [8]. Compared to other concrete constitutive models involving plastic deformation (termed plasticity model for simplicity), CDPM differs because the plasticity problem is formulated in the effective stress space. Moreover, the damage is associated with equivalent plastic strains, and tensile or compressive damage is calibrated from uniaxial tensile or compressive stress–strain curves. The following briefly describes the main ingredients of CDPM. Readers could refer to the Abaqus users’ manual, reference [9], for more details.

2.1.1. Stress–Strain Relations in CDPM

The stress–strain relationship for general loading is given by the scalar damaged elasticity equation [9]:

$$\sigma =(1-d)D_0^{el}:(\epsilon -{\epsilon }^{pl})=D^{el}:(\epsilon -{\epsilon }^{pl}),$$

(1)

where $\sigma$ and $\epsilon$ are the stress and strain tensors, respectively; $D_0^{el}$ is the initial (undamaged) elastic stiffness tensor; $D^{el}$ is the degraded elastic stiffness tensor; ${\epsilon }^{pl}$ is the plastic strain tensor. $d$ is the scalar stiffness degradation variable, which can take values in the range from zero (undamaged material) to one (fully damaged material). Damage associated with the failure mechanisms of the concrete (cracking and crushing), therefore, results in a reduction in the elastic stiffness. Following the usual notions of continuum damage mechanics, the effective stress is defined as [9]:

$$\overline{\sigma }\stackrel{\mathrm{def}}{=}{D}_0^{el}:(\epsilon -{\epsilon }^{pl}).$$

(2)

The Cauchy stress is related to the effective stress through the scalar degradation relation:

$$\sigma =(1-d)\overline{\sigma }.$$

(3)

For any given cross-section of the material, the factor $(1-d)$ represents the ratio of the effective load-carrying area (i.e., the overall area minus the damaged area) to the overall section area. In the absence of damage, the effective stress is equivalent to the Cauchy stress. When damage occurs, however, the effective stress is more representative than the Cauchy stress because the effective area resists the external loads. It is, therefore, convenient to formulate the plasticity problem in terms of the effective stress $\overline{\sigma }$ and hardening variables ${\tilde{\epsilon }}^{pl}$.

Damaged states in tension and compression are characterized independently by two hardening variables, ${\tilde{\epsilon }}_t^{pl}$ and ${\tilde{\epsilon }}_c^{pl}$, which are referred to as equivalent plastic strains in tension and compression, respectively. In general, the evolution of the hardening variables in rate format is given by an expression of the form [9]:

$${\dot{\tilde{\epsilon }}}^{pl}=\mathbf{h}(\overline{\sigma },{\tilde{\epsilon }}^{pl})\cdot {\dot{\epsilon }}^{pl}.$$

(4)

Please note that an additive strain rate decomposition is assumed for the rate-independent CDPM as:

$$\dot{\epsilon }={\dot{\epsilon }}^{el}+{\dot{\epsilon }}^{pl},$$

(5)

where $\dot{\epsilon }$ is the total strain rate, ${\dot{\epsilon }}^{el}$ is the elastic part of the strain rate, and ${\dot{\epsilon }}^{pl}$ is the plastic part of the strain rate.

Microcracking and crushing in the concrete are represented by increasing values of the hardening variables. These variables control the evolution of the yield surface (or failure surface) and the degradation of the elastic stiffness. Abaqus assumes that the equivalent plastic strain rates are evaluated according to the expressions [9]:

$${\dot{\tilde{\epsilon }}}_t^{pl}\stackrel{\mathrm{def}}{=}r(\widehat{\overline{\sigma }}){\widehat{\dot{\epsilon }}}_{\max }^{pl};{\dot{\tilde{\epsilon }}}_c^{pl}\stackrel{\mathrm{def}}{=}-\left(1-r(\widehat{\overline{\sigma }})\right){\widehat{\dot{\epsilon }}}_{\min }^{pl},$$

(6)

where ${\widehat{\dot{\epsilon }}}_{\max }^{pl}$ and ${\widehat{\dot{\epsilon }}}_{\min }^{pl}$ are the maximum and minimum eigenvalues of the plastic strain rate tensor, respectively.

$r(\widehat{\overline{\sigma }})\stackrel{\mathrm{def}}{=}{\displaystyle \frac{{\displaystyle {\sum }_{i=1}^3}〈{\widehat{\overline{\sigma }}}_i〉}{{\displaystyle {\sum }_{i=1}^3}\left|{\widehat{\overline{\sigma }}}_i\right|}},0\le r(\widehat{\overline{\sigma }})\le 1$ is a stress weight factor that is equal to one if all principal effective stresses ${\widehat{\overline{\sigma }}}_i(i=1,2,3)$ are positive and equal to zero if they are negative. The Macauley bracket $〈·〉$ is defined by $〈x〉={\displaystyle \frac{1}{2}}(\left|x\right|+x)$. If the eigenvalues of the plastic strain rate tensor ${\widehat{\dot{\epsilon }}}_i(i=1,2,3)$ are ordered such that ${\widehat{\dot{\epsilon }}}_{\max }^{pl}={\widehat{\dot{\epsilon }}}_1\ge {\widehat{\dot{\epsilon }}}_2\ge {\widehat{\dot{\epsilon }}}_3={\widehat{\dot{\epsilon }}}_{\min }^{pl}$, the evolution equation for general multi-axial stress conditions can be expressed in the following matrix form:

$${\dot{\tilde{\epsilon }}}^{pl}=\left[\begin{array}{c}{\dot{\tilde{\epsilon }}}_t^{pl}\\ {\dot{\tilde{\epsilon }}}_c^{pl}\end{array}\right]=\hat{\mathbf{h}}(\widehat{\overline{\sigma }},{\tilde{\epsilon }}^{pl})\cdot {\widehat{\dot{\epsilon }}}^{pl},\hat{\mathbf{h}}(\widehat{\overline{\sigma }},{\tilde{\epsilon }}^{pl})=\left[\begin{array}{ccc}r(\widehat{\overline{\sigma }})& 0& 0\\ 0& 0& -\left(1-r(\widehat{\overline{\sigma }})\right)\end{array}\right],{\widehat{\dot{\epsilon }}}^{pl}=\left[\begin{array}{c}{\widehat{\dot{\epsilon }}}_1\\ {\widehat{\dot{\epsilon }}}_2\\ {\widehat{\dot{\epsilon }}}_3\end{array}\right].$$

(7)

The yield function $F(\overline{\sigma },{\tilde{\epsilon }}^{pl})$ represents a surface in effective stress space, which determines the states of failure or damage. For the plastic-damage model,

$$F(\overline{\sigma },{\tilde{\epsilon }}^{pl})\le 0.$$

(8)

Plastic flow is governed by a flow potential G according to the flow rule [9]:

$${\dot{\epsilon }}^{pl}=\dot{\lambda }{\displaystyle \frac{\partial G(\overline{\sigma })}{\partial \overline{\sigma }}},$$

(9)

where $\dot{\lambda }$ is the nonnegative plastic multiplier. The plastic potential is also defined in the effective stress space, and the CDPM uses non-associated plasticity.

Within the context of the scalar-damage theory, the stiffness degradation is isotropic and characterized by a single degradation variable $d$. The evolution of the degradation variable is governed by hardening variables ${\tilde{\epsilon }}^{pl}$ and the effective stress $\overline{\sigma }$; that is, $d=d(\overline{\sigma },{\tilde{\epsilon }}^{pl})$. Similarly, by decomposing the degradation variable into tensile and compressive damage variables, Abaqus assumes:

$$1-d=(1-s_{\mathrm{t}}{d}_{\mathrm{c}})(1-s_{\mathrm{c}}{d}_{\mathrm{t}}),$$

(10)

where $s_{\mathrm{t}}$ and $s_{\mathrm{c}}$ are functions of the stress state and are used to simulate stiffness recovery effects due to stress reversals. They are defined as follows [9]:

$$s_{\mathrm{t}}\stackrel{\mathrm{def}}{=}1-{\omega }_{\mathrm{t}}r(\widehat{\overline{\sigma }})(0\le {\omega }_{\mathrm{t}}\le 1);s_c\stackrel{\mathrm{def}}{=}1-{\omega }_{\mathrm{c}}(1-r(\widehat{\overline{\sigma }}))(0\le {\omega }_{\mathrm{c}}\le 1),$$

(11)

where ${\omega }_{\mathrm{t}}$ and ${\omega }_{\mathrm{c}}$ are weight factors and control the recovery of the tensile and compressive stiffness upon load reversal, and the default values are ${\omega }_{\mathrm{t}}$ = 0 and ${\omega }_{\mathrm{c}}$ = 1 for concrete.

In summary, the CDPM’s elastic-plastic response is described in terms of the effective stress and the hardening variables by Equations (2), (4), (8) and (9), where $\dot{\lambda }$ and $F$ obey the Kuhn–Tucker conditions of $\dot{\lambda }F=0;\dot{\lambda }\ge 0;F\le 0$. The above constitutive relations for the elastic-plastic response are decoupled from the stiffness degradation response Equation (3), which makes the CDPM attractive for an effective numerical implementation.

2.1.2. Input of Plasticity Parameters

The plasticity parameters determine the yield function and the flow rule. In particular, the yield function of Abaqus CDPM is defined in Equation (12) to account for different tension and compression strength evolution.

$$F={\displaystyle \frac{1}{1-\alpha }}(\overline{q}-3\alpha \overline{p}+\beta 〈{\widehat{\overline{\sigma }}}_{\max }〉-\gamma 〈-{\widehat{\overline{\sigma }}}_{\max }〉)-{\overline{\sigma }}_{\mathrm{c}}\le 0.$$

(12)

$$\alpha ={\displaystyle \frac{{\sigma }_{b0}-{\sigma }_{c0}}{2{\sigma }_{b0}-{\sigma }_{c0}}};\beta ={\displaystyle \frac{{\overline{\sigma }}_{\mathrm{c}}}{{\overline{\sigma }}_{\mathrm{t}}}}(1-\alpha )-(1+\alpha );\gamma ={\displaystyle \frac{3(1-K_{\mathrm{c}})}{2K_{\mathrm{c}}-1}}.$$

(13)

In Equations (12) and (13), $\overline{q}$ represents the Mises equivalent effective stress and $\overline{p}$ is the effective hydrostatic pressure. As the typical experimental values of the ratio between biaxial and uniaxial compressive strength ${\sigma }_{b0}/{\sigma }_{c0}$ range between 1.10 and 1.16, the coefficient $\alpha$ can vary from 0.08 to 0.12. ${\widehat{\overline{\sigma }}}_{\max }$ stands for maximum principal effective stress, ${\overline{\sigma }}_{\mathrm{t}}$ and ${\overline{\sigma }}_{\mathrm{c}}$ are the effective tensile and compressive cohesion stresses (as explained in Section 2.1.3), respectively. The coefficient $\gamma$ enters the yield function only for stress states of triaxial compression, which can be determined by comparing the yield conditions along the tensile and compressive meridians. $K_c$ is the ratio of second stress invariant on tensile and compressive meridians that could be calculated by referring to the literature [7].

Regarding the flow rule, the Drucker–Prager hyperbolic function, as defined in Equation (14), is chosen for the flow potential

$$G=\sqrt{{(ϵf_{\mathrm{t}}\tan \psi )}^2+{\overline{q}}^2}-\overline{p}\tan \psi ,$$

(14)

where $\psi$ is the dilatancy angle measured in the p–q plane (deviatoric plane) at high confining pressure. f_t is the uniaxial tensile stress at failure. $ϵ$ is referred to as the eccentricity of the plastic potential surface.

According to Equations (12)–(14), the spatial stress behavior of concrete depends on four plasticity parameters $K_{\mathrm{c}}$, $\psi$, ${\sigma }_{b0}/{\sigma }_{c0}$, $ϵ$, whose value must be specified as input. The default values for the four parameters in Abaqus are 2/3, 35°, 1.16, and 0.1 [9], and users can adjust parameters according to specific loading conditions.

2.1.3. Input of Uniaxial Behavior

The scalar degradation variable $d$ is reduced to $d=d_{\mathrm{t}}$ for uniaxial tension and $d=d_{\mathrm{c}}$ for uniaxial compression. The stress–strain relationship response of concrete subjected to uniaxial tension or compression is shown in Figure 1, where the tension part is enlarged for a clearer view.

The tensile and compressive stresses in concrete in Figure 1 can be expressed as [9]

$${\sigma }_{\mathrm{t}}=(1-d_{\mathrm{t}})E_0({\epsilon }_{\mathrm{t}}-{\tilde{\epsilon }}_{\mathrm{t}}^{\mathrm{pl}});{\sigma }_{\mathrm{c}}=(1-d_{\mathrm{c}})E_0({\epsilon }_{\mathrm{c}}-{\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{pl}}),$$

(15)

and the effective uniaxial cohesion stresses, ${\overline{\sigma }}_{\mathrm{t}}$ and ${\overline{\sigma }}_{\mathrm{c}}$, are given as:

$${\overline{\sigma }}_{\mathrm{t}}={\displaystyle \frac{{\sigma }_{\mathrm{t}}}{(1-d_{\mathrm{t}})}}=E_0({\epsilon }_{\mathrm{t}}-{\tilde{\epsilon }}_{\mathrm{t}}^{\mathrm{pl}});{\overline{\sigma }}_{\mathrm{c}}={\displaystyle \frac{{\sigma }_{\mathrm{c}}}{(1-d_{\mathrm{c}})}}=E_0({\epsilon }_{\mathrm{c}}-{\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{pl}}).$$

(16)

d_t and d_c are damage variables of concrete under tension and compression, respectively, which could be derived as $d_{\mathrm{t}}=1-\sqrt{{\sigma }_{\mathrm{t}}/(E_0{\epsilon }_{\mathrm{t}})}$ and $d_{\mathrm{c}}=1-\sqrt{{\sigma }_{\mathrm{c}}/(E_0{\epsilon }_{\mathrm{c}})}$ based on the Sidoroff energy equivalence principle [29], E₀ is concrete’s initial modulus of elasticity; ε_t and ε_c are the tensile and compressive strains of concrete, respectively; ${\tilde{\epsilon }}_{\mathrm{t}}^{\mathrm{pl}}$ and ${\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{pl}}$ are the equivalent plastic strains of concrete.

The equivalent plastic strains are defined as a function of cracking strain ${\tilde{\epsilon }}_{\mathrm{t}}^{\mathrm{ck}}$ and inelastic strain ${\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{in}}$, as shown in Equation (17):

$${\tilde{\epsilon }}_{\mathrm{t}}^{\mathrm{pl}}={\tilde{\epsilon }}_{\mathrm{t}}^{\mathrm{ck}}-{\displaystyle \frac{d_{\mathrm{t}}}{1-d_{\mathrm{t}}}}{\displaystyle \frac{{\sigma }_{\mathrm{t}}}{E_0}};{\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{pl}}={\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{in}}-{\displaystyle \frac{d_c}{1-d_{\mathrm{c}}}}{\displaystyle \frac{{\sigma }_{\mathrm{c}}}{E_0}}.$$

(17)

The cracking strain ${\tilde{\epsilon }}_{\mathrm{t}}^{\mathrm{ck}}$ and inelastic strain ${\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{in}}$ are defined as the total strain minus the elastic strain in the undamaged state and can be expressed as:

$${\tilde{\epsilon }}_{\mathrm{t}}^{\mathrm{ck}}={\epsilon }_{\mathrm{t}}-{\displaystyle \frac{{\sigma }_{\mathrm{t}}}{E_0}};{\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{in}}={\epsilon }_{\mathrm{c}}-{\displaystyle \frac{{\sigma }_{\mathrm{c}}}{E_0}}.$$

(18)

Equations (15)–(18) determine the input of uniaxial mechanical behavior parameters. Data pairs of (${\sigma }_{\mathrm{t}}$,${\tilde{\epsilon }}_{\mathrm{t}}^{\mathrm{ck}}$) and ($d_{\mathrm{t}}$,${\tilde{\epsilon }}_{\mathrm{t}}^{\mathrm{ck}}$) define CDPM’s tension stiffening and damage; data pairs of (${\sigma }_{\mathrm{c}}$,${\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{in}}$) and ($d_{\mathrm{c}}$,${\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{in}}$) define stress hardening, followed by strain softening and damage in compression.

The constitutive theory in this section can capture the effects of irreversible damage associated with the failure mechanisms that occur in concrete as long as the confining pressures are less than four or five times the ultimate compressive strength in uniaxial compressive loading. Although CDPM can characterize the nonlinear mechanical behavior of concrete under cyclic or dynamic loading, the computational cost associated with cycle-by-cycle simulation in high-cycle fatigue is enormously high [10]. Therefore, an equivalent fatigue damage constitution of concrete that can directly describe concrete’s residual mechanical behavior after specific fatigue loading cycles is desirable. Such equivalent fatigue constitution can be obtained by introducing a fatigue damage variable that reflects the damage evolution law during fatigue loading into concrete’s monotonic constitution [30].

2.2. Fatigue Damage Evolution

2.2.1. Fatigue Damage Evolution Equation

The fatigue damage variable can be constructed based on the maximum fatigue strains of each cycle. Because of initial defects and the preloading effect, the strain obtained from the first loading cycle usually includes some initial strain of discrete nature besides that caused by cyclic loading. Therefore, to estimate the fatigue damage accumulated by cyclic loading more accurately, a fatigue damage variable in a normalized linear form is proposed, as shown in Equation (19):

$$D_{\mathrm{f}}={\displaystyle \frac{{\epsilon }_{\max }^{\mathrm{n}}-{\epsilon }_{\max }^1}{{\epsilon }_{\max }^{\mathrm{f}}-{\epsilon }_{\max }^1}}.$$

(19)

${\epsilon }_{\max }^1$ is the maximum fatigue strain during the first loading cycle, corresponding to the strain at maximum fatigue load; ${\epsilon }_{\max }^{\mathrm{n}}$ is the maximum fatigue strain during the nth loading cycle; ${\epsilon }_{\max }^{\mathrm{f}}$ is the maximum fatigue strain during the last loading cycle when fatigue failure occurs. These fatigue strains are shown schematically in Figure 2.

Since maximum fatigue strain typically follows three-stage development, the damage variable from Equation (19) also presents the three-stage development characteristic. In contrast to the piece-wise description widely used in the literature [26,27,28], a continuous nonlinear equation is proposed to describe the three-stage evolution of the fatigue damage variable, as expressed in Equation (20):

$$D_{\mathrm{f}}(n/N_{\mathrm{f}})=a{({\displaystyle \frac{b}{b-n/N_{\mathrm{f}}}}-1)}^{1/c}.$$

(20)

D_f(n/N_f) is the damage variable corresponding to the nth fatigue loading; a, b, and c are unitless equation parameters obtained by fitting the experimental data whose physical meaning and value range will be identified in the following sections. For normalizing loading cycles, the cycle ratio n/N_f is defined as the number of applied cycles n divided by fatigue life N_f, the total number of cycles a specimen can sustain before fatigue failure. Although any division of cycle ratio is acceptable, only values at given cycle ratios (i.e., n/N_f = 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 or n/N_f = 0, 0.2, 0.4, 0.6, 0.8, 1.0) are presented in this study for illustration purpose.

2.2.2. Analysis of Critical Parameters

By definition, there is no damage at the beginning of fatigue loading, and complete damage (i.e., D_f = 1) occurs when the cycle ratio reaches 1. The parameter b in Equation (20) can thus be interpreted as a function of parameters a and c by simple transformation, and the correlation is expressed as in Equation (21):

$$b=a^c+1.$$

(21)

Since Equation (21) is a power or exponential function of a and c, parameter b is called the composite factor. The analysis of critical parameters can be simplified to the impact analysis of a and c, whose effects on the damage variables (Figure 3) are explored by changing one parameter while keeping the other parameter constant.

Figure 3a shows the effect of various values of a on the damage evolution curve with a constant c = 4. As a increases, the slope of stage Ⅰ of the evolution curve grows more rapidly, and stage Ⅰ accounts for a more significant proportion of the entire curve. Stage Ⅲ of the evolution curve develops in the opposite direction; with an increase of a, stage III approaches unity more slowly, and the characteristics may gradually disappear. Since a exhibits significant influence on the development rate and proportion of stages Ⅰ and Ⅲ, a is therefore called the stability factor of the beginning and end stages.

Based on previous analysis, Figure 3b similarly explores the effect of c on the damage evolution curve while setting a = 0.4. The three-stage evolution trend becomes more apparent as c increases; unlike the effect of a, an increase in c results in an increased development rate of both stages Ⅰ and Ⅲ. However, compared with the influence of a on stages Ⅰ and Ⅲ, c shows a more significant impact on the damage stabilization development stage. Since increasing c decelerates the slope of stage Ⅱ, c can be regarded as a development rate factor.

With various combinations of a and c, a family of evolution curves that could reflect the three-stage fatigue damage development are generated and plotted in Figure 4.

The fatigue damage variable and its evolution equation in Equation (20) can be obtained from maximum fatigue strains of fatigue tests under different loading conditions. It should be noted that only a few references have reported the original strain data collected during the whole process of fatigue loading. Specimens from axial compression and flexural tension fatigue tests in the literature [31,32] are taken as examples to derive typical fatigue damage evolution equations.

2.2.3. Evolution Equation under Axial Compression

This section presents a comprehensive analysis of the fatigue behavior of concrete under axial compression from the literature [31]. The fatigue strains have been measured in the middle portion of the cylindrical specimens along the longitudinal direction, ensuring a thorough understanding of the behavior. Maximum stress level S_max or concrete strength f will be treated as a controlling factor in fitting the equation. The compressive specimens in the literature [31] were Φ100 mm × H200 mm concrete cylinders with three strength grades, 26 MPa, 52 Mpa, and 84 Mpa, respectively. These specimens were loaded under four maximum stress levels, 0.75, 0.80, 0.85, and 0.95 (75%, 80%, 85%, and 95% of the static strength), respectively, while the minimum stress level was kept constant at 0.25. The preparation of materials included a single source of Portland cement, crushed stone of maximum size 13 mm (coarse aggregates), and natural river sand (fine aggregates); the fineness moduli of fine and coarse aggregates were 2.89 and 6.08, respectively. In addition, an admixture (superplasticizer with retarder) was used in the mix for the strength of 52 Mpa and 84 Mpa to give good workability, and silica fume (powder) was used for the high-strength concrete of 84 Mpa. The mixture proportion parameters, including the amount of cement, water, aggregates, and any additives used, are shown in

Table 1. Mix proportions of the concrete used in fatigue tests under axial compression from the literature [31].

Compressive Strength (Mpa)	Cement (kg/m3)	Water (kg/m3)	Fine Aggregates (kg/m3)	Coarse Aggregates (kg/m3)	Admixture (%)	Slump (mm)
26	320	218	717	931	/	61
52	450	156	698	1047	3	81
84	442	156	633	1032	2	169

.

The damage variables were calculated by Equation (19), and the damage development trend of various specimens is shown in Figure 5. Specimens are designated as S-XX-C-XX, where S-XX indicates the maximum stress level, and C-XX specifies the strength grade. For example, S-75-C-26 means a specimen with a compressive strength of 26 Mpa and a maximum stress level S_max of 0.75.

Damage variables of various specimens are meticulously grouped according to stress levels or material strength and are fitted to fatigue damage evolution equations following Equation (20). The equation parameters are listed in

Table 2. Fitted parameters a, b, and c in the fatigue damage evolution model (Equation (20)) for concrete under compression from the literature [31], where R² is the coefficient of determination from regression analysis, indicating the model’s goodness of fit to the experimental data.

Controlling Factor		a	b	c	R2
Maximum stress level Smax	0.75	0.442	1.020	4.786	0.98
	0.80	0.491	1.038	4.475	0.97
	0.85	0.550	1.086	4.000	0.96
	0.95	0.662	1.182	3.800	0.96
Strength f (Mpa)	26	0.537	1.098	3.739	0.89
	52	0.563	1.060	4.954	0.95
	84	0.476	1.037	4.458	0.94

.

The column of R² in Table 2 shows that the determination coefficient of the damage evolution equation based on maximum stress level is above 0.95, whose overall correlation effect is better than that based on concrete strength. Figure 6 and Figure 7 show the distribution of damage variables and fitted damage evolution curves under different concrete strength grades and maximum stress levels, respectively.

Comparing Figure 6 and Figure 7, the scatter between damage variables of various specimens is larger when grouped under material strength, particularly for fatigue specimens made from 26 MPa concrete. When specimens are grouped under maximum stress levels, Figure 7 shows a more satisfactory correlation between the damage evolution curve and the data points of damage variables. Therefore, with the stress level as the controlling factor, Equation (20) can ensure acceptable fitting goodness and accurately reflect the evolution of fatigue damage. The finding underscores the importance of stress levels in predicting the evolution of fatigue damage. This conclusion is consistent with the observation that stress level is the most significant factor affecting concrete’s fatigue life [33].

Based on theoretical evolution curves reflecting three-stage fatigue damage development in Figure 4 and experiment data in Figure 7, a value range for applying the parameters in Equation (20) is suggested. The stability factor a can be taken as (0.3, 0.7), the development rate factor c ranges (3, 7), and the corresponding composite factor b can be taken as (1.0, 1.3).

2.2.4. Evolution Equation under Flexural Tension

The flexural specimens, 150 mm × 150 mm × 550 mm concrete prisms with a strength of 50 Mpa, were fatigue tested under four-point bending at maximum stress levels of 0.75 and 0.80 and a minimum stress level of 0.10 [32]. The tests were performed according to the Chinese standard GB/T 50081-2019 [34]. Fatigue strains were measured on the bottom surface of the specimen in the pure bending region along the longitudinal direction. The damage development of flexural specimens is shown in Figure 8, where the last digit in the legend represents the serial number of specimens.

Similar to the work presented in Section 2.2.3, the damage evolution equations of flexural tension fatigue at different maximum stress levels are derived and presented in Figure 9. It can be observed that the damage evolution curves of flexural tension fatigue could approximate the damage variables well, with R² all above 0.95. Comparison between Figure 7 and Figure 9 shows that the stability coefficient a of flexural tension fatigue damage evolution curves is generally smaller than that in Table 2. This observation indicates that the proportion of accumulated damage in the first development stage from flexural tension fatigue is less than that from axial compression fatigue. While stability coefficient a exhibits a proportional relationship to maximum stress level under axial compression, the relationship between them is opposite under flexural tension fatigue. Therefore, the lower limit of stability coefficient a can be extended to 0.2 when proposing the range of factors in flexural tension fatigue damage evolution models. In addition, the development rate factor c is inversely proportional to the maximum stress level in both stress states, which is consistent with the fact that the rate of damage development accelerates with increasing applied fatigue stress.

2.3. Equivalent Fatigue Damage Constitution

A typical stress–strain curve under fatigue loading shows that as the cycle elapses, the residual strain increases and the elasticity modulus decreases due to fatigue damage. For CDPM to consider fatigue damage, four parameters should be determined to establish the equivalent fatigue damage constitution after a certain number of fatigue loadings. They are fatigue residual strength and associated peak strain, fatigue residual elasticity modulus, and fatigue residual strain after n loading cycles. Calculation assumptions and analysis methods to determine these parameters are described as follows.

2.3.1. Calculation Assumptions

Experimental investigation demonstrates that strength deterioration of concrete in fatigue tests develops in a similar pattern with the descending branch of the monotonic stress–strain curve [14,15,16,17,18]. A typical stress–strain response under fatigue loading is depicted with a monotonic stress–strain curve, as shown in Figure 10. f_r and ε_p,r are the ultimate strength and associated strain in the monotonic curve where r = c, t represents concrete’s compressive and tensile properties. x(n) is the ratio between fatigue peak strain ε_p,fat associated with fatigue residual strength after n loading cycles, and the monotonic peak strain ε_p,r. The accumulation of deformation and damage increases with the cyclic loading process, which induces the decrease of concrete’s fatigue residual strength σ_r(n) and stiffness E_r(n). Fatigue failure (n = N_f) occurs when the reloading branch under cyclic loading encounters the descending branch of the monotonic loading curve, where the fatigue residual strength σ_r(n) degrades to the maximum fatigue stress σ_max and the maximum fatigue strain at fatigue failure ${\epsilon }_{\max }^{\mathrm{f}}$ is approximately equivalent to the strain corresponding to the maximum fatigue stress in the descending branch of monotonic stress–strain curve [10,24].

Based on this assumption of fatigue failure criterion, scholars introduced the concept of envelope curve to indicate the degradation process of fatigue residual strength. The descending branch of the monotonic stress–strain curve, termed the envelope curve, therefore describes the relationship between fatigue residual strength and the number of loading cycles.

2.3.2. Fatigue Residual Strength

If there is no direct experimental constitution describing the material behavior, the constitutive equations from various codes can be employed to describe the envelope curve of fatigue residual strength. According to the Chinese code for the design of concrete structures [35], the descending branch of the monotonic constitutive equation can be expressed as:

$${\displaystyle \frac{{\sigma }_{\mathrm{r}}}{f_{\mathrm{r}}}}={\displaystyle \frac{x}{{\alpha }_{\mathrm{r}}{[x-1]}^{\beta }+x}},x={\displaystyle \frac{\epsilon }{{\epsilon }_{\mathrm{p},\mathrm{r}}}}\ge 1,$$

(22)

where x is defined as the ratio of the strain in the descending branch to the peak strain in the constitutive curve; σ_r is stress in the descending branch; f_r is the axial compressive or tensile strength; α_r is the parameter of the descending branch, which can be taken with reference to the code [35]; β is the equation parameter, which is taken as 2 for compression and 1.7 for tension, respectively.

Transforming from Equation (22), the envelope curve of fatigue residual strength is derived by replacing σ_r with the residual strength σ_r(n) after n loading cycles and representing x with x(n). Based on the physical meaning of the envelope model (Figure 10), x(n) is defined as the following function of fatigue loading cycles:

$$x(n)=\left\{\begin{array}{cc}1& (n=1)\\ f(n,N_{\mathrm{f}})[x(N_{\mathrm{f}})-x(1)]+x(1)& (1\le n\le N_{\mathrm{f}})\end{array}\right.,$$

(23)

where f(n,N_f) is the shape function of x(n) during fatigue loading; x(N_f), as the maximum fatigue strain at fatigue failure ${\epsilon }_{\max }^{\mathrm{f}}$ divided by the monotonic peak strain ε_p,r, is a constant for a given fatigue specimen.

Due to the irreversibility of energy dissipation, the fatigue residual strength envelope is a monotonically decreasing function. The envelope curve must also satisfy the following boundary conditions: ${\sigma }_{\mathrm{r}}(x(1))=f_{\mathrm{r}}$; ${\sigma }_{\mathrm{r}}(x(N_{\mathrm{f}}))={\sigma }_{\max }$; $\sigma {\prime }_{\mathrm{r}}(n)\le 0$. Combining these boundary conditions and Equations (22) and (23), the fatigue residual strength of concrete can be derived and is expressed in Equation (24) as:

$${\sigma }_r(n)=f_{\mathrm{r}}{\displaystyle \frac{f(n,N_{\mathrm{f}})[x(N_{\mathrm{f}})-x(1)]+x(1)}{{\alpha }_{\mathrm{r}}{\left(f(n,N_{\mathrm{f}})[x(N_{\mathrm{f}})-x(1)]\right)}^a+f(n,N_{\mathrm{f}})[x(N_{\mathrm{f}})-x(1)]+x(1)}}.$$

(24)

It is worth noting that the shape function in the form of cycle ratio $f(n,N_{\mathrm{f}})={\displaystyle \frac{x(n)-x(1)}{x(N_{\mathrm{f}})-x(1)}}={\displaystyle \frac{n}{N_{\mathrm{f}}}}$ or the relative ratio of logarithmic life $f(n,N_{\mathrm{f}})={\displaystyle \frac{x(n)-x(1)}{x(N_{\mathrm{f}})-x(1)}}={\displaystyle \frac{\lg n}{\lg N_{\mathrm{f}}}}$ has been used in the literature [10,16,17,18]. The shape function of the cycle ratio indicates linear accumulation of fatigue damage. For the shape function of relative logarithmic life, the second-order derivative of f(n,N_f) with respect to n yields $f(n,N_{\mathrm{f}})″=-{\displaystyle \frac{1}{n^2\ln 10\lg N_{\mathrm{f}}}}<0$, indicating that f(n,N_f) is a convex function, as illustrated by the dash-dot line in Figure 11. The dash-dot line shows that the rate of fatigue strain increase gradually decreases as loading cycles elapse. This phenomenon might be the case when maximum fatigue stress is less than the fatigue threshold. However, for specimens that failed under fatigue, residual strain gradually accumulated in the cyclic loading process, and fatigue strain grew rapidly in stage Ⅲ [19], as shown by the solid line in Figure 11. The cycle ratio or relative logarithmic life as shape function, therefore, deviates from experimental observations and fails to describe the fatigue residual strength correctly.

A revisit of Equation (23) yields the defining equation of the shape function f(n,N_f), as expressed in Equation (25):

$$f(n,N_{\mathrm{f}})={\displaystyle \frac{x(n)-x(1)}{x(N_{\mathrm{f}})-x(1)}}.$$

(25)

Comparing Equation (25) with Equation (19), it is found that the physical meaning between the shape function and the damage variable is very similar. As a result, we propose the shape function f(n,N_f) adopting the fatigue damage evolution model in Equation (20). This proposal marks a significant step forward in our understanding of fatigue damage evolution in concrete specimens. The method introduces the concept of damage variables when deriving fatigue residual strength and mathematically explains the strength degradation due to fatigue damage. Adopting the fatigue damage evolution equation as the shape function f(n,N_f) in contrast to the commonly used cycle ratio n/N_f or relative logarithmic life lgn/lgN_f is the most significant novelty in this study.

2.3.3. Fatigue Residual Elasticity Modulus

For cyclic loading involving inelastic deformation, two types of elasticity modulus, i.e., tangent stiffness and secant stiffness, should be differentiated. Fatigue residual secant stiffness, the slope of the straight line connecting the two points representing the maximum and minimum fatigue loading in a given loading cycle, is smaller than fatigue residual tangent stiffness, the tangent slope of the initial straight part of the loading branch. The degradation process of fatigue residual stiffness is exhibited in Figure 12 [30].

In Figure 12, E₀ is the initial stiffness; E_n is the fatigue residual stiffness after n cycles, considering plastic flow and microcrack damage. The maximum fatigue strain consists of three parts: ε_e is the elastic part calculated as σ_max/E₀; ε_d is the inelastic strain induced by microcrack damage evaluated by subtracting ε_e from σ_max/E_n; and ε_pl is the inelastic strain induced by plasticity, which equals ${\epsilon }_{\max }^{\mathrm{n}}-{\epsilon }_{\mathrm{e}}-{\epsilon }_{\mathrm{d}}$.

The fatigue residual stiffness could be defined as in Equation (26):

$$E_{\mathrm{n}}=E_0/(1+({\epsilon }_{\mathrm{d}}/{\epsilon }_{\mathrm{e}})).$$

(26)

Assuming that ε_pl is coaxial with ε_d and δ = ε_pl/ε_d [30,36], Equation (26) could be rewritten as Equation (27):

$$E_{\mathrm{n}}={\displaystyle \frac{E_0(1+\delta )}{1+\delta +({\epsilon }_{\mathrm{d}}+{\epsilon }_{\mathrm{pl}})/{\epsilon }_{\mathrm{e}}}}.$$

(27)

Since ε_pl and ε_d cannot be directly measured, substituting ${\epsilon }_{\max }^{\mathrm{n}}={\epsilon }_{\mathrm{e}}+{\epsilon }_{\mathrm{pl}}+{\epsilon }_{\mathrm{d}}$ and ε_e = σ_max/E₀ = f_rS_max/E₀ into Equation (27), it could be rewritten as Equation (28):

$$E_{\mathrm{n}}={\displaystyle \frac{E_0(1+\delta )}{\delta +E_0{\epsilon }_{\max }^{\mathrm{n}}/f_{\mathrm{r}}{S}_{\max }}}.$$

(28)

Now, think about the last loading cycle just before fatigue failure. Please note that the maximum strain of concrete at fatigue failure ${\epsilon }_{\max }^{\mathrm{f}}$ could be obtained by substituting σ_max into Equation (22), the descending branch of the monotonic constitutive equation. Based on Equations (26)–(28), and assuming the tangent or secant modulus follows a similar proportion, the elasticity modulus at fatigue failure can be calculated as in Equation (29):

$$E_{\mathrm{f}}={\displaystyle \frac{E_0{f}_{\mathrm{r}}{S}_{\max }(1+\delta )}{\delta f_{\mathrm{r}}{S}_{\max }+E_0{\epsilon }_{\max }^{\mathrm{f}}}}.$$

(29)

E₀ is the initial stiffness; f_r is the static strength; S_max is the maximum stress level; ${\epsilon }_{\max }^{\mathrm{f}}$ is the maximum fatigue strain at fatigue failure. Δ is the ratio of plasticity-induced inelastic strain to inelastic strain induced by damages such as microcracks and microvoids, which ranges [1, 1.5] according to experimental results, and 1.2 is taken in this study [30,36]. Please note that two types of stiffness exist; their difference grows as the fatigue loading cycles increase. Therefore, when calculating the fatigue residual tangent stiffness at fatigue failure $E_{\mathrm{t}}^{\mathrm{f}}$, the initial stiffness E₀ should be taken as its tangent value $E_{\mathrm{t}}^0$. In contrast, the initial secant stiffness $E_{\mathrm{s}}^0$ should be used to evaluate fatigue residual secant stiffness at fatigue failure $E_{\mathrm{s}}^{\mathrm{f}}$.

The development of Internal cumulative damage essentially causes stiffness degradation in concrete at the macro level. Therefore, fatigue residual tangent modulus after n cyclic loading can be deduced considering the fatigue damage evolution equation defined by Equation (20), as expressed in Equation (30):

$$E_{\mathrm{r}}\left(n\right)=E_{\mathrm{t}}^{\mathrm{f}}+(1-D_{\mathrm{f}}(n/N_{\mathrm{f}}))(E_{\mathrm{t}}^0-E_{\mathrm{t}}^{\mathrm{f}}),$$

(30)

where the fatigue residual tangent modulus just before fatigue failure E_r(N_f) equals $E_{\mathrm{t}}^{\mathrm{f}}$.

2.3.4. Fatigue Residual Strain

Under fatigue loading, concrete’s fatigue residual strength and elastic modulus continuously decrease, while the residual strain accumulates with the increase of loading cycles. Fatigue residual strain after n cyclic loading is usually fitted from fatigue experiment data [37], which consists of complicated equations with many uncertain parameters. Therefore, a simplified method is developed to estimate fatigue residual strain.

At fatigue failure, the nonlinear loading branch of the equivalent stress–strain curve is shown as a thin solid line (a) in Figure 13, and the fatigue residual secant modulus calculated from Equation (29) is shown as a dashed line (b) in Figure 13, in which fatigue residual strength envelope curve is shown as the thick solid line.

From Figure 13, the fatigue residual strain at fatigue failure is calculated from the limit state parameters at fatigue failure as:

$${\epsilon }_{\mathrm{res}}^{\mathrm{f}}={\epsilon }_{\max }^{\mathrm{f}}-{\displaystyle \frac{{\sigma }_{\mathrm{r}}(N_{\mathrm{f}})}{E_{\mathrm{s}}^{\mathrm{f}}}},$$

(31)

where σ_r(N_f) = σ_max and ${\epsilon }_{\max }^{\mathrm{f}}$ are the residual strength and maximum fatigue strain at fatigue failure, which can be obtained from the fatigue residual strength envelope curve, and $E_{\mathrm{s}}^{\mathrm{f}}$ is the fatigue residual secant modulus at fatigue failure.

Assuming fatigue residual strain evolves similarly to the maximum fatigue strain [30], the damage variable can also be described by residual strain as in Equation (32):

$$D_{\mathrm{f}}={\displaystyle \frac{{\epsilon }_{\mathrm{res}}^{\mathrm{n}}-{\epsilon }_{\mathrm{res}}^1}{{\epsilon }_{\mathrm{res}}^{\mathrm{f}}-{\epsilon }_{\mathrm{res}}^1}},$$

(32)

where ${\epsilon }_{\mathrm{res}}^{\mathrm{f}}$ is the residual strain at fatigue failure; ${\epsilon }_{\mathrm{res}}^{\mathrm{n}}$ is the fatigue residual strain at nth cyclic loading; ${\epsilon }_{\mathrm{res}}^1$ is the residual strain at the first cyclic loading, which could be taken as 0 for simplicity.

By transforming Equation (32), fatigue residual strain at nth cyclic loading is derived and expressed in Equation (33) as

$${\epsilon }_{\mathrm{res}}^{\mathrm{n}}={\epsilon }_{\mathrm{res}}^1+D_{\mathrm{f}}({\epsilon }_{\mathrm{res}}^{\mathrm{f}}-{\epsilon }_{\mathrm{res}}^1).$$

(33)

Substituting Equation (20) into Equation (33), we can rewrite fatigue residual strain at any cycle ratio as

$${\epsilon }_{\mathrm{res}}^{\mathrm{n}}(n/N_{\mathrm{f}})={\epsilon }_{\mathrm{res}}^1+a{({\displaystyle \frac{b}{b-n/N_{\mathrm{f}}}}-1)}^{1/c}({\epsilon }_{\mathrm{res}}^{\mathrm{f}}-{\epsilon }_{\mathrm{res}}^1).$$

(34)

2.3.5. Equivalent Fatigue Constitution Based on Chinese Code [35]

The above analysis methods can be integrated with various concrete stress–strain curves to derive the corresponding equivalent fatigue damage constitution. Based on the compressive and tensile constitutive equations recommended by Chinese code [35], the equivalent fatigue damage constitution of the concrete after n cycles of fatigue loading is defined as expressed in Equation (35) for compression and Equation (36) for tension:

$$\sigma ={\displaystyle \frac{{\rho }_{\mathrm{c}}m}{m-1+x^m}}{E}_{\mathrm{r}}\left(n\right)(\epsilon -{\epsilon }_{\mathrm{res}}^{\mathrm{n}})\hspace{1em}\hspace{1em}{\epsilon }_{\mathrm{res}}^{\mathrm{n}}\le \epsilon \le {\epsilon }_{\mathrm{p},\mathrm{fat}},$$

(35)

$$\sigma ={\rho }_{\mathrm{t}}[1.2-0.2x^5]E_{\mathrm{r}}\left(n\right)(\epsilon -{\epsilon }_{\mathrm{res}}^{\mathrm{n}})\hspace{1em}\hspace{1em}{\epsilon }_{\mathrm{res}}^{\mathrm{n}}\le \epsilon \le {\epsilon }_{\mathrm{p},\mathrm{fat}},$$

(36)

where $x={\displaystyle \frac{\epsilon -{\epsilon }_{\mathrm{res}}^{\mathrm{n}}}{{\epsilon }_{\mathrm{p},\mathrm{fat}}-{\epsilon }_{\mathrm{res}}^{\mathrm{n}}}}$, ${\rho }_{\mathrm{c},\mathrm{t}}={\displaystyle \frac{{\sigma }_{\mathrm{r}}(n)}{E_{\mathrm{r}}(n)({\epsilon }_{\mathrm{p},\mathrm{fat}}-{\epsilon }_{\mathrm{res}}^{\mathrm{n}})}}$, $m={\displaystyle \frac{E_{\mathrm{r}}\left(n\right)({\epsilon }_{\mathrm{p},\mathrm{fat}}-{\epsilon }_{\mathrm{res}}^{\mathrm{n}})}{E_{\mathrm{r}}\left(n\right)({\epsilon }_{\mathrm{p},\mathrm{fat}}-{\epsilon }_{\mathrm{res}}^{\mathrm{n}})-{\sigma }_{\mathrm{r}}(n)}}$, and ε_p,fat is the fatigue peak strain corresponding to fatigue residual strength σ_r(n) on the fatigue residual strength envelope.

In summary, with the damage variable defined as in Equation (19), fatigue damage evolution curves (Equation (20), shown in Figure 7 for axial compression and Figure 9 for flexural tension) were obtained from fatigue test results. One aspect of novelties is that the parameters’ physical meaning in the fatigue damage evolution equation is clearly described. The evolution equation was then innovatively introduced to calculate critical parameters in CDPM for equivalent fatigue damage constitution, explicitly explaining performance degradation caused by fatigue damage. These innovations significantly contribute to the field by providing a simplified and physically sound model for fatigue response analysis of concrete structures.

3. Verifications

3.1. Feasibility of the Proposed Method

To verify the proposed method, we derived the equivalent fatigue damage stress–strain curves at specified cycle ratios for fatigue specimens tested by Kim et al. [31].

The two specimens tested at a maximum stress level of 0.75 for concrete compressive strength grades of 26 Mpa and 84 Mpa were selected for comparison. The monotonic loading constitution (abbreviated as the Kim model) and maximum fatigue strain data were reported in the literature [31]. The theoretical behavior of these specimens was predicted by the equivalent fatigue damage constitutive model (with appropriate fatigue damage evolution equation parameters from Table 2) and is shown in Figure 14 as various dashed lines at given cycle ratios.

The monotonic stress–strain relationships of the concrete are shown as solid lines in Figure 14. It is observed that the strains corresponding to the maximum fatigue stress in the ascending and descending branches of the Kim model are approximate to the measured maximum fatigue strains during the first and last loading cycles. Since the Kim model can describe the mechanical behavior of these two specimens reasonably well, the equivalent fatigue damage model was derived using the descending branch of the Kim model as the fatigue residual strength envelope curve.

Figure 14 shows the specimens’ predicted fatigue damage constitutive curves after specific loading cycles. These equivalent stress–strain curves indicate that the fatigue damage evolution equation constructed based on maximum fatigue strains could effectively describe the degradation behavior of stiffness, strength, and residual strain accumulation. When the ascending branch of the equivalent fatigue stress–strain curve reaches the maximum fatigue stress, a theoretical strain is obtained at a given cycle ratio. These theoretical strains are then compared to the measured maximum fatigue strains, shown as red dots in Figure 14. The comparison is depicted in Figure 15 for more details.

As demonstrated in Figure 15, the maximum relative error between the theoretical and experimental values is 12%, and more than 95% of the relative errors do not exceed 5%. Regression analysis of the results shows that the root mean square error (RMSE) between experimental and theoretical values in Figure 15a is 0.00011, and the RMSE in Figure 15b is 0.00003. This comparison verifies that the proposed equivalent fatigue damage constitutive model can effectively explain the stress–strain relationship of concrete after considering fatigue damage.

3.2. Discussion on Different Forms of Monotonic Constitution

Based on the above analysis, the equivalent fatigue constitutive relationship derived from the Kim model can satisfactorily predict the fatigue strain evolution process of two specimens. To eliminate the specificity caused by the single constitutive form and verify the universality of the method for constructing equivalent fatigue constitutive models, a comparative analysis will be conducted on the results of predicting fatigue deformation behavior in the literature [31] based on fatigue constitutive relationships derived from different monotonic constitutive forms.

Choosing the experimental data of S-75-C-26 in the literature [31] for analysis, and the compressive stress–strain relationship was calculated according to the critical parameters of peak strength = 26 MPa, peak strain = 0.002074, and initial modulus of elasticity = 21,000 Mpa based on the Chinese Code GB50010-2010 [35] and the CEB-FIP 2010 [38], respectively. Three constitutive curves are shown in Figure 16.

Figure 16 demonstrates that the trend of the three curves in the ascending branch is the same for the three forms of constitutive curves. While the Chinese code and Kim model present a gentler descending branch, the strength degradation in the descending branch of the CEB-FIP 2010 model appears more rapid. Again, assuming the descending branch forms the envelope curve and determines the theoretical value of the maximum fatigue strain at fatigue failure [14], the equivalent stress–strain relationships for each cyclic loading stage based on the Chinese code and CEB-FIP 2010 were further derived according to the method proposed in Section 2.3, as illustrated in Figure 17 and Figure 18.

From the perspective of prediction performance, the RMSE between the predicted values and measured values in Figure 17 is 0.000128, while that is 0.00008 in Figure 18, demonstrating satisfactory prediction performance. It is worth noting that the theoretical values calculated based on the Kim model and Chinese code both exhibit that they are lower than experimental values in stage Ⅰ but higher than experimental values at the end of stage Ⅱ. In contrast, theoretical values in Figure 18 are smaller than experimental values throughout the whole fatigue strain evolution stage. This predicted variance is primarily related to the assumed envelope relationship between the constitutive curve and maximum fatigue strain at failure.

As shown in Figure 16, the ascending branch of three constitutive curves is on the outer side of the strain point that represents initial loading to maximum stress level; that is, the assumed initial theoretical value is smaller than the experimental value, which might lead to smaller predicted values in the initial stage. However, the theoretical value of maximum strain at fatigue failure in the Chinese code and Kim model is greater than the experimental value. Therefore, the theoretical strain development rate will be greater than the experimental one, leading to gradually higher accumulated theoretical values at the end stage. On the other hand, the stress–strain relationship described in Figure 18 shows that the theoretical maximum fatigue strains at both initial loading and fatigue failure are smaller than the corresponding experimental one, leading to consistently smaller theoretical values of maximum fatigue strains during the fatigue loading process compared to the experimental values.

It should be noted that the concrete fatigue strain evolution results are discrete, and the discussion in this section only covers one specimen. However, the discrepancy induced by different assumed monotonic stress–strain relationships is considered. The analysis indicates that the proposed method for obtaining an equivalent fatigue constitution could be applied to more than one specific equation form. Nonetheless, the predicted results depend on selecting a monotonic constitutive curve that could reasonably encompass the development of maximum fatigue strains.

3.3. Discussion on Concrete Parameters Affecting Fatigue Degradation

Apart from differences in predicted results caused by different constitutive forms, critical parameters, such as aggregate size, loading frequency, specimen size, and stress variation, also impact the fatigue degradation of concrete. Restricted by reported data, only studies involving loading frequency and size effects are detailed in this section.

3.3.1. Loading Frequency

Taking the axial compression fatigue test as an example, we have analyzed the effect of frequency on fatigue damage evolution by referring to the fatigue deformation behavior of concrete specimens under different combinations of test conditions in pieces of literature [39,40]. The main experimental condition combinations are shown in

Table 3. Loading conditions of fatigue tests under axial compression from the literature [39,40]. The maximum stress level, stress ratio (ratio between minimum and maximum fatigue stress), and loading frequency vary.

Maximum Stress Level Smax	Stress Ratio	Frequency (Hz)	Number of Specimens
0.85	0.3	1/16	2
0.85	0.3	¼	2
0.85	0.3	1	2
0.85	0.3	4	2
0.80	0.1	2	4
0.80	0.2	2	4
0.80	0.2	3	4
0.75	0.1	2	4
0.75	0.1	3	4
0.70	0.1	3	4

, where the first four rows of S_max = 0.85 are from the literature [40], and the rest are from the literature [39]. Similar to the discussion in Section 2.2.3, the damage variables are then summarized according to frequency and maximum stress level, and the damage evolution equation is fitted, as shown in

Table 4. Fitted parameters a, b, and c in the fatigue damage evolution model (Equation (20)) for concrete under compression from the literature [39,40], where R² is the coefficient of determination from regression analysis.

Controlling Factor		a	b	c	R2
Frequency	2 Hz	0.550	1.015	7.000	0.94
	3 Hz	0.550	1.009	7.848	0.95
Maximum stress level Smax	0.85	0.551	1.046	5.161	0.96
	0.80	0.532	1.004	8.684	0.95
	0.75	0.524	1.002	9.515	0.97

.

As shown in Table 4, the damage evolution trends of different specimens are more consistent under the same stress level conditions, where the fitting goodness of damage evolution reaches 0.97 at a lower stress level of S_max = 0.75. It is worth noting that the damage evolution appears more discrete at S_max = 0.85 compared to the stress level S_max = 0.75. This phenomenon may be because, at higher stress levels of S_max ≥ 0.80 [24], the frequency dramatically affects the number of failure cycles, therefore affecting the evolution of fatigue deformation, whereas, at lower S_max, the frequency has little effect on the evolution of fatigue deformation [41].

For a detailed analysis of the effect of loading frequency on fatigue behavior, damage evolution curves are fitted for different loading frequencies at the same stress level, as presented in Figure 19. In the figure, (S_max, R, f)-x represents the test results of the x-th specimen under test conditions (S_max, R, f). For example, (0.80, 0.1, 2)-1 represents the test results of the first specimen under the test conditions of a maximum stress level of 0.80, stress ratio of 0.1, and loading frequency of 2 Hz.

It is illustrated that the differences in damage variable evolution at different frequencies are not significant for lower stress levels, which is consistent with the previous findings. As shown by the curves in Figure 19 for higher stress levels, the extent of damage development decreases with increasing loading frequency, leading to relatively high theoretical fatigue residual strength at high loading frequency.

Scholars have pointed out that for lower stress levels S_max, the influence of frequency can be ignored, and the impact will increase as the stress level increases. Therefore, fitting the damage evolution equation based on stress levels does not need to consider the influence of frequency at lower stress levels. However, applying the damage evolution equation to analyze strength attenuation at higher stress levels should also consider the influence of different loading frequencies on the results.

3.3.2. Size Effect

The static strength varies with the size of concrete specimens, and specimen size may also affect the dynamic behavior of the material. Scholars have found that fatigue life generally increases with increasing specimen size [24], while discussion on the influence of specimen size on fatigue residual strength is scarce.

This section focuses on the influence of specimen size on fatigue behavior. For the axial compression fatigue test, the size effect is analyzed by comparing specimens’ fatigue damage evolution behavior from the literature [31,40,42]; the information on selected fatigue deformation data is shown in

Table 5. Data information for fatigue tests under axial compression from the literature [31,40,42] for which the effect of specimen size on fatigue damage evolution is studied. The loading frequency varies slightly between 1 and 5 Hz, a relatively standard value during fatigue loading.

Data Source	Specimen Size	Loading Frequency
Literature [31]	Φ100 mm × H200 mm	1 Hz
Literature [42]	Φ150 mm × H250 mm	5 Hz
Literature [40]	100 mm × 100 mm × 100 mm	1 Hz

.

Since loading frequency in the literature [31] and [42] are different, these data are used to analyze the size effect at lower stress levels. As discussed in the previous section, loading frequency has a more significant impact at higher stress levels. Therefore, the data from reference [40] is supplemented for comparative verification under S_max = 0.85. The comparison results are shown in Figure 20. For a given stress level, the smaller the specimen size, the higher the accumulated damage under the same cycle ratio. The damage development of smaller-sized concrete is relatively higher than that of larger specimens, resulting in a faster degradation rate of fatigue residual strength and, thus, a lower theoretical fatigue residual strength. However, this difference induced by the size effect is not particularly significant. Comparing the damage evolution curves in Figure 20, it is evident that stress level is more critical for fatigue behavior degradation than the size effect, which is also consistent with the findings in the literature [24].

4. Applications

4.1. Concrete Fatigue in Axial Compression

While the proposed fatigue damage evolution equation was constructed based on fatigue experiments conducted in the literature [31], it could reasonably be applied to reflect the fatigue performance of concrete specimens tested under similar loading conditions. The cylindrical specimens F1 and F3, with a diameter of 100 mm and a height of 250 mm under axial compression from another literature [42], were selected as an application example. Numerical simulation was performed using the finite element software Abaqus 2019, incorporating the proposed equivalent fatigue damage constitution and fatigue damage evolution equation parameters from Table 2.

4.1.1. Conversion of Material Parameters

Literature [42] did not provide monotonic loading constitutive relationships for specimens F1 and F3, and the equivalent fatigue damage model suggested in Section 2.3.5 was used. It should be noted that the concrete parameters (such as strength, peak strain, and tangent modulus) in constitutive equations suggested by code [35] are defined based on standard-size 150 mm × 150 mm × 300 mm prisms. Therefore, these parameters must be converted from cylindrical specimens based on empirical formulas [43].

As the prismatic compressive strength f_c increases monotonically with the cubic compressive strength f_cu, the ratio of f_c to f_cu grows with concrete strength grade. The proportional relationship between these two types of strength can be expressed by Equation (37):

$$f_{\mathrm{c}}=\left(0.70~0.92\right)f_{\mathrm{cu}}.$$

(37)

In the literature [44], the modification factor for converting between Φ100 mm × H250 mm cylindrical compressive strength $f_c\prime$ and f_cu was derived by an ultrasonic-rebound synthetic method as $f_{\mathrm{c}}\prime =\overline{{\eta }_1}{f}_{\mathrm{cu}}=0.85f_{\mathrm{cu}}$; and the relation between prismatic compressive strength f_ct of size 100 mm × 100 mm × 300 mm and f_cu expresses as $f_{\mathrm{ct}}=\overline{{\eta }_2}{f}_{\mathrm{cu}}=0.81f_{\mathrm{cu}}$. Since f_c decreases with the specimen’s depth-thickness ratio increase, the conversion between f_c and f_cu is determined as f_c = 0.90 f_cu for Φ100 mm × H250 mm cylinders.

The compressive peak strain of the prismatic specimen can be calculated by Equation (38), and the initial elasticity modulus is determined according to Equation (39) [45].

$${\epsilon }_{\mathrm{p},\mathrm{c}}=(700+172\sqrt{f_c})\times {10}^{-6},$$

(38)

$$E_0={\displaystyle \frac{{10}^5}{2.2+(33/f_{\mathrm{cu}})}}.$$

(39)

After combining Equations (37)–(39), the conversion from cylindrical parameters of specimens F1 and F3 into prismatic parameters that will be integrated into the constitutive equations suggested by code [35] is expressed in

Table 6. Converted material parameters in the monotonic constitutive curve (Equation (22)) for fatigue tests under axial compression from the literature [42]. f_c, ε_p,c, and E₀ are the prismatic compressive strength, peak strain, and initial modulus of elasticity, respectively.

Specimen	${\mathit{f}}_{\mathit{c}}\prime$ (MPa)	Converted fcu (MPa)	Converted fc (MPa)	Converted εp,c	Converted E0 (MPa)
F1	46.2	54	48.6	0.001899	35,573
F3	43.9	51	45.9	0.001865	35,124

.

4.1.2. Finite Element Simulation

The finite element model of the cylindrical fatigue specimen is shown in Figure 21. Rigid elements are provided at both ends of the cylinder to simulate the loading device. Only tangential and normal behavior contact properties are set when defining the interaction between the specimen and rigid elements. Reference points RP-1 and RP-2 are placed on the top and bottom surfaces of the model, respectively, and are connected to the rigid element surfaces by coupling constraint. The hinge constraint is applied on RP-2 to simulate the boundary conditions during the fatigue test, while the cyclic load is defined on RP-1 as a concentrated force. For specimen F1, the fatigue-loaded stress levels are S_max = 0.76 and S_min = 0.05; for specimen F3, S_max = 0.84 and S_min = 0.05.

In the finite element model, concrete is modeled as C3D8R (an 8-node linear hexahedral element with reduced integration and hourglass control) elements. A total of 1536 elements are used in the model with a standard mesh size of 20 mm along the periphery and height. The number of elements is chosen to strike a balance between computational efficiency and model accuracy, ensuring reliable results without excessive computational cost. The embedded CDPM is used to simulate material properties, where the plasticity parameters are defined as listed in

Table 7. Plasticity parameters of CDPM adopted in the analysis of fatigue tests under axial compression from the literature [42]. The definition of plasticity parameters and their default values in Abaqus are shown in Section 2.1.2.

Material	Dilatancy Angle (°)	Eccentricity	Ratio of Biaxial/Uniaxial Compressive Strength	Shape Factor of Yield Surface	Viscosity Parameter
Concrete	55	0.1	1.16	0.667	0.0005

[46]. A suggested value between 10⁻⁴ and 10⁻³ of the viscosity parameter is utilized in Abaqus to improve the convergence rate in strain softening and stiffness degradation [9]. The compressive behavior of the CDPM is calculated from converted parameters in Table 6 according to Equations (15)–(18).

The stress–strain simulation curves of specimens F1 and F3 under monotonic axial compression loading are shown in Figure 22.

Figure 22 shows that the compressive strength of F1 is 44.4 MPa, and the associated peak strain is 0.00174 from simulation; the simulated strength of F3 is 42.5 Mpa, and the peak strain is 0.001698. The simulated compressive strength is slightly smaller than the experiment value (shown in Table 6). However, the relative errors are within 5%. The literature does not specify the peak strain from monotonic loading [42].

The ratio of the maximum fatigue strains to the peak strain ${\epsilon }_{\max }^{\mathrm{n}}/{\epsilon }_{\mathrm{p},\mathrm{c}}$ for fatigue loading is reported in the literature [42]. Maximum fatigue strains under axial compression fatigue loading were obtained from simulation by incorporating the equivalent fatigue damage constitutive model into Abaqus. These maximum fatigue strains at specified cycle ratios are then divided by the simulated peak strain in Figure 22. The comparison between the numerical simulation results and experimental results of ${\epsilon }_{\max }^{\mathrm{n}}/{\epsilon }_{\mathrm{p},\mathrm{c}}$ is presented in Figure 23.

Figure 23 shows satisfactory simulation results for both specimens F1 and F3, and the relative errors between the simulation and experiment values of the two specimens do not exceed 7%. The RMSE of simulation and test values in Figure 23a is 0.04, and 0.02 in Figure 23b. This application example demonstrates that the fatigue damage evolution equation constructed based on the literature [31] could reflect the performance degradation of materials under similar conditions, and the proposed equivalent fatigue damage constitutive model could be applied to simulate fatigue loading under axial compression reasonably well.

4.2. Concrete Fatigue in Flexural Tension

The proposed equivalent fatigue damage constitutive model was applied to simulate the bending fatigue behavior of concrete in flexural tension for specimens available in the literature.

4.2.1. Equivalent Fatigue Damage Constitutive Curve

Literature [47] reported a four-point bending fatigue test of 100 mm × 100 mm × 400 mm specimens. Due to experimental condition limitations, only the stress–strain curve’s ascending branch under monotonic bending was measured [47]. The test confirmed that the maximum fatigue strain at failure is inversely proportional to the applied maximum fatigue stress. Thus, the descending branch, also the hypothetical fatigue residual strength envelope curve, was obtained by fitting Equation (22). Parameter α_r was determined based on the measured maximum fatigue strain at fatigue failure and the applied maximum fatigue stress. Specimens C39A5 from the literature [47] were selected for analysis. The specimens were made from a concrete batch with a mix ratio of 1:0.39:1.7:3.2, where C39 represents a water–cement ratio of 0.39, and A5 represents a 5% air content in the concrete. The monotonic constitutive model of specimen C39A5 is shown in Figure 24, which exhibits a good agreement with the experimental values at two critical positions, representing the initial and ultimate loading cycles, respectively.

The equivalent fatigue damage constitutive equations under two stress levels are derived based on the initial constitutive model and fatigue damage evolution equation from Figure 9. Figure 25 clearly shows the fatigue residual strength, fatigue residual strain, and stress–strain development curve at every tenth of the cycle ratio. The equivalent fatigue constitutive curves could reasonably reflect the fatigue strains at minimum and maximum stress levels. Compared to experimental values of fatigue strains, theoretical values of the maximum fatigue strain are higher than the experimental results during stage Ⅰ of fatigue damage development, gradually approaching the experimental observation by the end of stage Ⅲ. This discrepancy could be because the fatigue residual strength envelope curve based on the descending branch of the monotonic stress–strain curve is obtained from experiments under axial loading. As we know, a difference in stress distribution exists between these two tensile stress states. Because of the stress gradient, the probability of crack initiation in a beam specimen under bending would be lower than in axial tension. The stress reduction after peak strain under flexural tension is expected to be slower than under axial tension due to the redistribution of bending stress. The initial rapid decreasing rate of the envelope curve leads to a lower simulated fatigue residual strength calculated by Equation (24), which would cause faster strain development. The equivalent fatigue constitutive equations under two stress levels are then integrated into the CDPM in Abaqus to predict the residual fatigue behavior of concrete in flexural tension.

4.2.2. Finite Element Analysis

The finite element model of specimen C39A5 is shown in Figure 26, whose dimensions are identical to the 100 mm × 100 mm × 400 mm beam in the four-point bending fatigue test from the literature [47]. The loading and supporting device are rigid, and the concrete specimen is defined as deformable solid C3D8R elements. Tangential and normal behavior contact properties are set when defining the interaction between the supporting device and concrete specimen, and “Tie” constraints are defined between the loading device and the contacting nodes on the upper surface of the specimen.

The specimen is monotonically loaded under force control with a 12 kN/min loading rate in finite element analysis. The tensile bending stress at the bottom fiber in the constant moment region is then obtained and displayed in Figure 27. The bending stress increases linearly with the applied load initially. As tensile damage gradually develops due to microcracking, the stress increases more slowly. As shown in Figure 27, the ultimate bending load is simulated as 10.95 kN, beyond which the specimen is completely damaged.

With the obtained ultimate bending capacity, a numerical analysis of the fatigue behavior of specimens C39A5 was conducted at maximum stress levels of 0.80 and 0.75 and a stress ratio of 0.1, with a loading frequency of 5 Hz. Based on the equivalent fatigue damage constitutive model, the simulated strains corresponding to the maximum and minimum stress levels are obtained and compared with test values, as illustrated in Figure 28.

In Figure 28, the solid lines stand for the specimen’s simulated fatigue stress–strain responses at each predetermined cycle ratio when it is fatigue-loaded between S_min and S_max; the black dots are the projection of the ends of these space curves into the cycle ratio versus strain plane, representing the maximum and minimum fatigue strains, respectively. Regular and inverted triangles symbolize the measured values of maximum and minimum strains from the fatigue test.

Consistent with Figure 25, the simulated values of the maximum fatigue strain are slightly higher than the experimental values during the earlier stage of fatigue damage development, gradually approaching each other by the end of the fatigue test. Generally, the equivalent fatigue damage constitutive model presents satisfactory predictions. In Figure 28a, the RMSE between simulated and experimental values of maximum fatigue strain is 0.000017 and 0.000013 for minimum fatigue strain. In Figure 28b, the RMSE between simulated and experimental values of maximum and minimum fatigue strains are 0.000024 and 0.000013, respectively. This application example indicates that the proposed equivalent fatigue damage constitutive model could be applied to simulate fatigue behavior under flexural tension reasonably well.

4.3. Research Limitations and Future Work

While the proposed constitutive model has been validated with high-cycle fatigue experimental data in axial compression and flexural tension, it is important to note that several limitations should be addressed in future work. This section underscores the ongoing nature of our research and the need for continuous improvement.

Good agreement between predicted results and observed behavior depends on whether the selected monotonic stress–strain curve could reasonably predict the development of maximum fatigue strains: the loading branch considers the initial loading, and the unloading branch predicts the maximum fatigue strain just before failure. The method of determining fatigue residual strain could also be improved. Since the damage evolution equation is derived from uniaxial fatigue test data, only its application to particular cases of axial compression and flexural tension has been confirmed. Whether the proposed method can be applied to analyze more general cases of multi-axial fatigue and variable amplitude fatigue problems needs further investigation.

Furthermore, it is crucial to consider scatters in monotonic constitutive curves and fatigue damage evolution curves between specimens tested under the same conditions. Therefore, a reliability analysis that considers the probability of failure is desirable in the future. A similar constitutive relation should be constructed for reinforcement before the proposed model can be successfully applied to real-world scenarios for fatigue response analysis of reinforced concrete structures.

5. Conclusions

A novel simplified analysis method of equivalent fatigue stress–strain relationship is proposed by introducing fatigue damage evolution into the concrete damaged plasticity model (CDPM), extending the practical application of CDPM to high-cycle fatigue simulation. The following conclusions can be drawn:

(1) A damage variable was constructed based on the maximum fatigue strain, and a nonlinear fatigue damage evolution equation, Equation (20), was proposed to characterize the three-stage fatigue damage development. The effects of three parameters on the evolution equation were explored. The stability factor a influences the proportion of stages I and III, the development rate factor c affects the development rate of stage II, and the composite factor b is determined by (a^c + 1). These parameters could be fitted to experimental fatigue data, and the value range of each parameter was recommended as a ∈ (0.2, 0.7), c ∈ (3, 7), and the corresponding b ∈ (1.0, 1.3). Fatigue damage evolution equations are affected mainly by the maximum stress level and less affected by concrete strength grade, loading frequency, or specimen size.

(2) The fatigue damage evolution equation was incorporated when calculating critical CDPM parameters, and the equivalent fatigue constitutive model after multiple cycles was derived. Changing the shape function f(n,N_f) from cycle ratio n/N_f or relative logarithmic life lgn/lgN_f to fatigue damage evolution equation D_f(n/N_f), we obtained an improved estimate of fatigue residual strength. The ascending branch of the equivalent constitutive model follows the yield hardening rule, and the residual strain before fatigue failure can be calculated from the fatigue residual secant stiffness. The residual strain for any given number of fatigue cycles is then estimated based on the damage evolution equation. By modeling cyclic loading at discrete cycle ratios, the degradation process of fatigue residual strength, residual stiffness, and residual strain accumulation in high-cycle fatigue could be simulated with low computation cost.

(3) The proposed equivalent fatigue constitutive model was validated with experimental data, showing acceptable agreement with observed behavior for fatigue tests under axial compression and flexural tension stress conditions. Our work demonstrates for the first time that the fatigue residual strength envelope curve also applies to bending fatigue behavior analysis. The simplified analysis method is not restricted to any specific equation of concrete constitution; both empirical formulas in the codes and stress–strain relationship directly measured from experiments have derived adequate prediction results. However, the matching between assumed constitutive forms and the measured maximum fatigue strains affected the accuracy of predicted strains.

The proposed fatigue damage evolution equation and equivalent fatigue constitutive model laid the foundation for subsequent applications in analyzing the fatigue response of reinforced concrete structures. However, the proposed method belongs to the deterministic analysis method category and only applies to particular cases. Further research is needed to consider randomness and general cases of three-dimensional and variable amplitude loading fatigue.