# Fractional Stress Relaxation Model of Rock Freeze-Thaw Damage

Authors: Q. Liu; W. Chen; J. K. Guo; R. F. Li; D. Ke; Y. Wu; W. Tian; X. Z. Li
Journal: Advances in Materials Science and Engineering
 (2021)

---

## 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.


---

## Body

## 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:(1)dβftdtβ=Dt0tβft=dndtn∫t0tt−τn−β−1Γn−βfτdτ,where β is the order of fractional calculus, n is order of fluxionary calculus, and Γ(β) is the gamma function defined as follows:(2)Γβ=∫0∞e−ttβ−1dt,Reβ≥0,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:(3)σt=ξdβdtβεt,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 is(4)D=1−1−NfN1/V+1,where 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 variableD is introduced to improve the H-K model, and the improved model is shown in Figure 1.Figure 1
Rock stress relaxation model.Figure1 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,(5)ε=ε1+ε2,σ=σ1+σ2,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:(6)σ=Gtε,where G(t) is relaxation modulus, which can measure stress relaxation deformability.(7)Gt=11/E1+1/E2+ξdβ/dtβ,where E1 and E2 are the elastic moduli of the two springs, respectively, and ξ is the coefficient of viscosity. Making Pt=1/Gt, we can obtain(8)Pt=1E1+1E2+ξdβ/dtβ,where Pt is a parameter to simplify the calculation through Laplace transformation:(9)Pt=L−11E1+1E2+ξdβ/dtβ,(10)L−11E1⋅s=1s.First, the second term can be expanded based on Taylor series and then subsequently perform the inverse Laplace transformation:(11)L−11E2+ξdβ/dtβ=∑n=0∞E2n−1−ξ1+nt1+nβ−1Γ1+nβ.Combining equations (6), (7), and (11), we obtain(12)Gt=1∑n=0∞E2n−1−ξ1+nt1+nβ−1/Γ1+nβ,(13)Ei=1−DE0,i=1,2,ξ=1−Dξ0,where 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).


## 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:(1)dβftdtβ=Dt0tβft=dndtn∫t0tt−τn−β−1Γn−βfτdτ,where β is the order of fractional calculus, n is order of fluxionary calculus, and Γ(β) is the gamma function defined as follows:(2)Γβ=∫0∞e−ttβ−1dt,Reβ≥0,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:(3)σt=ξdβdtβεt,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 is(4)D=1−1−NfN1/V+1,where 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 variableD is introduced to improve the H-K model, and the improved model is shown in Figure 1.Figure 1
Rock stress relaxation model.Figure1 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,(5)ε=ε1+ε2,σ=σ1+σ2,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:(6)σ=Gtε,where G(t) is relaxation modulus, which can measure stress relaxation deformability.(7)Gt=11/E1+1/E2+ξdβ/dtβ,where E1 and E2 are the elastic moduli of the two springs, respectively, and ξ is the coefficient of viscosity. Making Pt=1/Gt, we can obtain(8)Pt=1E1+1E2+ξdβ/dtβ,where Pt is a parameter to simplify the calculation through Laplace transformation:(9)Pt=L−11E1+1E2+ξdβ/dtβ,(10)L−11E1⋅s=1s.First, the second term can be expanded based on Taylor series and then subsequently perform the inverse Laplace transformation:(11)L−11E2+ξdβ/dtβ=∑n=0∞E2n−1−ξ1+nt1+nβ−1Γ1+nβ.Combining equations (6), (7), and (11), we obtain(12)Gt=1∑n=0∞E2n−1−ξ1+nt1+nβ−1/Γ1+nβ,(13)Ei=1−DE0,i=1,2,ξ=1−Dξ0,where 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 (Figure2).Figure 2
The cube test block of rock.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 Figure3. The physical parameters of the rock are shown in Table 1.Figure 3
Freeze-thaw cycles.
(a)(b)(c)Table 1
The physical parameters of the rock.
Density (kg/m3)Elastic modulus (MPa)Poisson’s ratioThermal conductivity (W/m°C)Specific heat capacity (J/kg°C)203843200.30.2291e4Figure3 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 as(14)σ=Etε,(15)Et=E01−∑i=1ng1−e−t/τii,where E(t) is the relaxation function of rock, which can be expressed in the form of Prony series, gi is the shear relaxation of rock, and τi is the Prony delay time constant.The temperature field and static stress are shown in Figures4 and 5:Figure 4
Temperature distribution map of rock freeze-thaw cycle: (a) 5 freeze-thaw cycles, (b) 10 freeze-thaw cycles, and (c) 20 freeze-thaw cycles.
(a)(b)(c)Figure 5
The stress distribution map of rock freeze-thaw cycles: (a) 5 freeze-thaw cycles, (b) 10 freeze-thaw cycles, and (c) 20 freeze-thaw cycles.
(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 Figure6.Figure 6
Rock stress relaxation curve by ABAQUS.


### 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 Figure6 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.Table 2
Stress parameters of the model.
NE1 (MPa)E2 (MPa)ξ (MPa ∗ h)β56.040.5470.2210.307105.680.5180.2080.358205.170.4680.1910.416Figure 7
Comparison of rock stress between the experimental data and predicted results.PuttingE1, E2, and ξ in Table 2 into equation (12), we obtain(16)D5=0.16,D10=0.21,D20=0.28.PutD5 = 0.16, D10 = 0.21, and D20 = 0.28 into equation (4) to obtain N and V.


## 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 (Figure2).Figure 2
The cube test block of rock.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 Figure3. The physical parameters of the rock are shown in Table 1.Figure 3
Freeze-thaw cycles.
(a)(b)(c)Table 1
The physical parameters of the rock.
Density (kg/m3)Elastic modulus (MPa)Poisson’s ratioThermal conductivity (W/m°C)Specific heat capacity (J/kg°C)203843200.30.2291e4Figure3 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 as(14)σ=Etε,(15)Et=E01−∑i=1ng1−e−t/τii,where E(t) is the relaxation function of rock, which can be expressed in the form of Prony series, gi is the shear relaxation of rock, and τi is the Prony delay time constant.The temperature field and static stress are shown in Figures4 and 5:Figure 4
Temperature distribution map of rock freeze-thaw cycle: (a) 5 freeze-thaw cycles, (b) 10 freeze-thaw cycles, and (c) 20 freeze-thaw cycles.
(a)(b)(c)Figure 5
The stress distribution map of rock freeze-thaw cycles: (a) 5 freeze-thaw cycles, (b) 10 freeze-thaw cycles, and (c) 20 freeze-thaw cycles.
(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 Figure6.Figure 6
Rock stress relaxation curve by ABAQUS.


## 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 Figure6 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.Table 2
Stress parameters of the model.
NE1 (MPa)E2 (MPa)ξ (MPa ∗ h)β56.040.5470.2210.307105.680.5180.2080.358205.170.4680.1910.416Figure 7
Comparison of rock stress between the experimental data and predicted results.PuttingE1, E2, and ξ in Table 2 into equation (12), we obtain(16)D5=0.16,D10=0.21,D20=0.28.PutD5 = 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.



---

*Source: 3936968-2021-02-13.xml*