Abstract
Freeze-thaw cycle is a type of fatigue loading, and rock stress relaxation under freeze-thaw cycles takes into account the influence of the freeze-thaw cycle damage and deterioration. Rock stress relaxation under freeze-thaw cycles is one of the paramount issues in tunnel and slope stability research. To accurately describe the mechanical behaviour of stress relaxation of rocks under freeze-thaw, the software element is constructed based on the theory of fractional calculus to replace the ideal viscous element in the traditional element model. The freeze-thaw damage degradation of viscosity coefficient is considered. A new three-element model is established to better reflect the nonlinear stress relaxation behavior of rocks under freeze-thaw. The freeze-thaw and stress relaxation of rock are simulated by ABAQUS, the relevant model parameters are determined, and the stress relaxation equation of rock under freeze-thaw cycle is obtained based on numerical simulation results. The research shows that the test results are consistent with the calculated results, indicating that the constitutive equation can better describe the stress relaxation characteristics of rocks under freeze-thaw and provide theoretical basis for surrounding rock support in cold region.
1. Introduction
In the construction of the infrastructure construction and operation, there is no lack of long-term stability of rock masses. Stress relaxation refers to the phenomenon that the stress of rock decays with time under constant deformation. It is common in underground caverns, roadways, slopes, and other rock mass projects. Many projects are often damaged due to stress relaxation, which is manifested in soft rock projects. Especially obvious, such as Dadu 15% of the lithology of the diversion tunnel of Hedanba Hydropower Station is quartz mica schist, which has significant creep and stress relaxation characteristics, which has a significant impact on the long-term strength and safety of surrounding rock. Rock stress relaxation is the basic standard to reflect the long-term stability of rock engineering, which is closely related to slope instability, bridge foundation settlement, and the surrounding rock deformation in tunnel [1–3]. The damage and stress relaxation of rocks in the freeze-thaw cycles are more serious in the construction of rock engineering in the extremely frigid zones. This behaviour results in negative effects on the long-term stability of rock engineering in the Alpine area. Therefore, establishing the constitutive model of rock stress relaxation in freeze-thaw cycles and studying the effect of freeze-thaw cycles, stress levels, and time on the evolution of rock stress relaxation is of immense significance to the study of tunnel and slope stability.
At present, some researchers have obtained some results and practical experiences in the research of rock freeze-thaw damage mechanism and freeze-thaw test [4–6]. Also, an increasing number of researchers have focused on the study of rock stress relaxation model. Tian et al. [7] carried out the stress relaxation test on argillaceous siltstone under the influence of confining pressure in the laboratory and obtained the curve influence characteristics of confining pressure on indoor rock stress relaxation from the test data; Tian et al. [8] also conducted shear relaxation laboratory tests on three structural planes with different angles made of cement mortar and analyzed the influence of normal stress of climbing angle on shear relaxation characteristics; and the triaxial postpeak stress relaxation tests of red sandstone were carried out by Liu et al. [9]. According to the characteristics of the curve obtained from the test data, the effects of rock degree, strain loading rate, and confining pressure on stress relaxation behaviour were discussed.
Apart from integral calculus, the order of fractional calculus can be fractional, irrational, or even complex. Therefore, there are significant advantages in describing memory in time and path dependence in space. Meanwhile, fractional calculus can preferably describe the nonlinear mechanical behaviour of rock and soil and provide a new dimension for the study of rock relaxation constitutive law due to the advantages of good global correlation and clear physical meaning [10–18]. Few studies have applied the theory of fractional calculus to describe the rheological nonlinearity of rock and soil mass with good results. For example, Zhang et al. [19] established the stress relaxation equation of expansive soil using fractional calculus; Liu and Xu [20] have computed the relaxation function for higher-order fractional rheological models (they contain three different fractional parameters, and the maximum order of equations is more than one) involving more parameters than in our analysis; Zhang et al. proposed a novel four-element fractional viscoplastic (FVP) model.
However, the existing stress relaxation studies are mainly focused on normal temperature conditions, and there are few reports on the study of rock stress relaxation features in freeze-thaw cycles. The stress relaxation process of rock engineering in cold regions is mostly affected by freeze-thaw cycles, which aggravates the damage and destruction of rocks and affects the stress relaxation features of rocks. Therefore, in this paper, the Newton viscous element in the viscoelastic model is replaced with soft element, and the damage variable is introduced to represent the mechanical properties of rock under freeze-thaw cycles based on fractional calculus operator and theory. Thus, a constitutive model is established to preferably reflect the stress relaxation of rock under freeze-thaw cycles.
2. Fractional Stress Relaxation Model
2.1. Fractional Calculus
Compared to traditional calculus, the order of factional calculus is not just an integer, but the entire complex number. When the value of the integer derivative of a point is calculated, the result merely depends on the locality of the point. But the fractional derivative is integrated in the whole range of values and has a strong dependence on the lower limit of the integral. For example, the stress change in rock depends not only on the current state but also on the mechanical state of rock in a fixed time. It is difficult to model and analyze using the classical differential equations for the memory system; however, the fractional derivative incorporates memory effects due to the nonlocality. Consequently, fractional calculus is quite a useful tool for analysing these systems.
Currently, there are many forms of the definition of fractional calculus, and this paper introduces the Riemann–Liouville [21] definition as follows:where β is the order of fractional calculus, n is order of fluxionary calculus, and Γ(β) is the gamma function defined as follows:where Re is the real part of a real number.
The constitutive equation for the model containing the fractional derivative with respect to the Riemann–Liouville definition can be given as follows:where σ is stress and ε is strain.
2.2. Fractional-Order Viscous Element with Freeze-Thaw Cycle Damage
It is assumed that the freeze-thaw damage of rock is only caused by frozen-heave force and frozen-heave force is a unidirectional tensile load, so the freeze-thaw cyclic damage to the rock can be equivalently replaced by the fatigue of the rock under the action of cyclic load. The freeze-thaw damage model of rock iswhere D is the damage variable, Nf is the number of freeze-thaw cycles, N is the number of freeze-thaw cycles at failure, and V is the material-related parameter.
It is assumed that the rock exhibits isotropic damage, the damage law of each parameter in the relaxation model is the same, the damage variable D is introduced to improve the H-K model, and the improved model is shown in Figure 1.

Figure 1 shows that the model consists of two parts connected in series, one part is a spring (H1), the other part is a spring (H2), and a sticky pot (N) in parallel (H|N). According to the principle of series connection of component models,where ε1 is the strain of spring (H1), ε2 is the strain of element (H|N2), σ1 is the stress of spring (H1), and σ2 is the stress of element (H|N2).
Combining equations (3) and (5) yields equation (6), as follows:where G(t) is relaxation modulus, which can measure stress relaxation deformability.where E1 and E2 are the elastic moduli of the two springs, respectively, and ξ is the coefficient of viscosity. Making , we can obtainwhere is a parameter to simplify the calculation through Laplace transformation:
First, the second term can be expanded based on Taylor series and then subsequently perform the inverse Laplace transformation:
Combining equations (6), (7), and (11), we obtainwhere E0 is the initial elastic modulus and ξ0 is the initial coefficient of viscosity. The relaxation model can be obtained by substituting equation (12) into equation (6).
3. Model Verification and Parameter Identification
3.1. Numerical Simulation
In this paper, the stress relaxation behaviour of rocks under freeze-thaw cycles was simulated by ABAQUS. The sequential thermal stress-coupled transient heat conduction analysis is used as the simulation analysis method for the freeze-thaw test.
The model uses a cube test block of 100 mm × 100 mm × 100 mm (Figure 2).

According to the regulation of the medium and slow freezing method in Chinese current standards, the temperature change range is −20 °C∼20 °C, the cycle is 24 hours, and the number of cycles is 5 times, 10 times, and 20 times [22].
In sequential-coupled thermal stress analysis, temperature directly affects the stress and strain fields, but is not limited by displacement. For this type of analysis, firstly, we analyze the heat transfer problem and then analyze the thermal stress based on the temperature field to obtain the calculation result. The transient heat conduction of the model structure was analyzed by ABAQUS. The curve of the node temperature with time during the freeze-thaw process of the model test block can be obtained, as shown in Figure 3. The physical parameters of the rock are shown in Table 1.

(a)

(b)

(c)
Figure 3 reveals the feasibility of stimulating. It is feasible to simulate the freezing and thawing process of rock by heat conduction analysis: The node temperature alternates between −20°C and 20°C, which is in line with the temperature change setting of the slow freeze-thaw machine. The length of time represented by the horizontal axis of the temperature change graph can reflect the number of freeze-thaw cycles of rock, which is consistent with the setting of the quick freeze-thaw machine.
After the heat conduction analysis of different freeze-thaw cycles, the temperature field analysis results are imported into a model that defines the plastic damage parameters of concrete for stress relaxation analysis. When the strain of the rock remains unchanged, the relationship between the stress and strain of the rock under the influence of stress relaxation can be expressed aswhere E(t) is the relaxation function of rock, which can be expressed in the form of Prony series, is the shear relaxation of rock, and τi is the Prony delay time constant.
The temperature field and static stress are shown in Figures 4 and 5:

(a)

(b)

(c)

(a)

(b)

(c)
The stress relaxation data in ABAQUS were imported into the numerical analysis, and the stress relaxation diagram of the rock is shown in Figure 6.

3.2. Determination of Parameters of Rock Stress Relaxation Model
In this paper, the L-M algorithm was used to determine the parameter values of the model. This approach eliminates the problem of improper selection of the initial value of the least square method. The L-M algorithm is used to nonlinearly simulate the experimental data results in Figure 6 based on equations (12) and (13). The model parameter results are shown in Table 2, and the comparison of rock stress between the experimental data and predicted results is shown in Figure 7.

Putting E1, E2, and ξ in Table 2 into equation (12), we obtain
Put D5 = 0.16, D10 = 0.21, and D20 = 0.28 into equation (4) to obtain N and V.
4. Results and Discussion
(1)The improved nonlinear model can effectively simulate the whole process of rock stress relaxation. The calculated results of the model are consistent with the numerically simulated results, which verify the reliability of the model.(2)Compared with the integral order stress relaxation model, the fractional-order model is more suitable for describing the historical memory of materials and more in line with the mechanical properties of rocks.(3)The rock stress relaxation model can be used for the deformation trend of surrounding rocks with time after tunnel excavation in cold area, thus providing theoretical basis for surrounding rock support and lining construction in extremely frigid zones.(4)This model was based on the macroscopic stress relaxation characteristics of rocks under freeze-thaw cycles; hence, it is recommended that future studies are directed towards its application at the microscale.
Data Availability
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also form part of an ongoing study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was supported by the National Natural Science Fund Project (Grant no. 41877230) and Open Fund of State Key Laboratory of Mechanical Behavior and System Safety of Traffic Engineering Structures (KF020-23) and the Fundamental Research Funds for the Central Universities CHD.