Abstract
To study the long-term stability of the large-scale freezing project, the longest subway connecting channel constructed by the artificial freezing method in China is taken as the background. The uniaxial compressive strength test, triaxial shear test, and triaxial creep test of artificial frozen soil under different temperatures were carried out by using the MTS 370.25 mechanical testing machine. A series of triaxial creep tests for artificial frozen soil under different stresses and temperatures were carried out. Based on the creep test results, a fractional-order constitutive model of frozen silt with sand was proposed. Comparing the calculated results with the tested ones, it is found that the proposed fractional-order model can simulate the properties of frozen silt with sand very well. Compared with the other creep models of frozen soil, the proposed model here has higher accuracy and strong stress sensitivity. Through the simulation of the connecting channel construction and comparison with the observed data in site, the effectiveness of the model is verified, which further verifies that the model has important guiding significance for the freezing project.
1. Introduction
Frozen soils are compound materials consisting of solid mineral particles, ice crystals, liquid water, and gaseous inclusions. Compared with unfrozen soils, the mechanical property of frozen soils is more complex because of their components and sensibility to temperature. It is now well recognized that frozen soils show nonlinear viscoelastic-plastic behavior. Therefore, in some projects, for example, the superlong freezing connecting passage project studied in this paper, linear approximations cannot be used. The strength and creep characteristics of frozen soils are the most important mechanical properties for the engineering construction. In order to meet the needs of the application of engineering activities, a series of experimental studies on the mechanical properties of frozen soil had been conducted [1–9]. Jessberger [10] believed that frozen soil satisfies the Mohr–Coulomb strength criterion under low confining pressure. Yang et al. [11] used an improved Hoek–Brown criterion to describe the nonlinear strength characteristics of frozen soil under low confining pressure. Lai et al. [12, 13] carried out triaxial strength tests of frozen sand and proposed the elastoplastic damage constitutive model based on the continuous damage theory. Yang et al. [14] proposed an elastic-plastic constitutive model of frozen soil under high confining pressure. Lai et al. [15, 16] carried out static triaxial compression tests of frozen loess and introduced the plastic theory to propose an elastoplastic incremental frozen loess constitutive model under the triaxial state. Jia et al. and Zuo et al. [17, 18] used FLAC3D to study the frozen soil constitutive model based on the Mohr–Coulomb yield criterion.
For the research on the creep mechanical properties of frozen soil, Andersland and Akili [19] studied the effects of stress, temperature, and soil structure on the creep rate of frozen soil by differential creep test. Taking sand as an example, Vyalov et al. [20] pointed out that creep is caused by the generation and development of microcracks and other defects in soil. Based on the element model, Sun et al. [21] established a simple and applicable creep constitutive relation of frozen soil. Based on a large number of uniaxial and triaxial creep tests of artificial frozen soil, Li et al. [22–24] deduced a nonlinear creep damage coupled constitutive model of artificially frozen clay and frozen sand under unloading state and complex stress path conditions. Bai et al. [25–27] studied the creep effects caused by temperature, particle rearrangement, and particle contact surface characteristics. Chen and Qiao [28] studied the creep process and uniaxial compression test of artificial frozen soil and established a corresponding mechanical model. Yao et al. [29] proposed a creep model of frozen soil with temperature as an independent variable by analyzing the creep test results of frozen soil at different temperatures. Arenson and Springman [30] obtained the creep model of rich ice-frozen soil with a temperature close to 0°C through the frozen soil creep test. Wang et al. [31] used the non-Newtonian viscous element instead of the Newtonian viscous element to describe the one-dimensional accelerated creep state of frozen soil in the Nishihara model. Many studies on the creep of frozen soil show that the stress level is an important factor affecting the strength and deformation of frozen soil [32]. The application of the geotechnical constitutive model is very important. Based on the construction of the constitutive model, Yuan et al. [33, 34] put forward some works that should be paid attention to in engineering applications. Mengke et al. [35] conducted a series of triaxial creep tests on warm frozen silts, and they established a fractional-order rheological element model based on the creep test results.
Most of these studies show that the creep curves only include primary and secondary creep stages under different stress levels and ignore the influence of temperature effect. From the test results in this paper, the creep curves can be distinguished by three stages clearly when the stress level exceeds the stress threshold, and the influence of temperature effect will be involved in the characteristics of frozen soil and constitutive model. In this paper, the triaxial creep tests for artificial frozen soil under different stresses and temperatures were carried out. Based on creep test results, a fractional-order constitutive model of frozen silt with sand was proposed. Compared with the other creep models of frozen soil, the proposed model here has higher accuracy and strong stress sensitivity. Finally, the test data of frozen soils at different temperatures and test conditions were used to verify the proposed model, and it is found that the proposed fractional-order model can simulate the properties of frozen silt with sand very well and has a great engineering value and realistic significance.
2. Engineering Background and Test Plan
2.1. Engineering Background
The connecting passage of Fuzhou Zi-Wu Interval is the longest subway connecting passage constructed by the freezing method in China. The center distance of this connecting passage is 66 m, and the centerline is about 6.4 m away from the nearest pile foundation of the overpass pier. The artificial freezing method is used to reinforce the strata, and the underground excavation method is used for construction. The section of the freezing curtain of the connecting passage is shown in Figure 1. The effective thickness of the frozen curtain is 2.0 m, and the average temperature of frozen soil is less than −10°C.

This connecting passage is the longest one constructed by the freezing method in China, and its location environment is complex, so the construction risk level is grade I. Compared with conventional freezing engineering, this freezing engineering has the characteristics of long freezing time, large freezing volume, long excavation and construction time, complex surrounding environment, and high requirements for deformation control. The elastic model often used in previous engineering calculations cannot fully adapt to this project, so this paper adopts the fractional derivative viscoelastic-plastic constitutive model to improve the calculation.
The main soil layers where the connecting passage is located are mainly silty fine sand and silt with sand. In this paper, silt with sand is taken as the research object to carry out the conventional geotechnical tests and physical and mechanical properties’ tests of frozen soil, including uniaxial compressive strength test, triaxial shear test, and triaxial creep test.
2.2. Experimental Equipment and Samples
The test soil samples were taken from the strata of the connecting channel, and the sampling depth was 20 m–30 m below the underground. Routine geotechnical tests were performed on the soil samples, and the basic physical parameters required for the determination are shown in Table 1.
The mechanical property tests of artificial frozen soil are carried out by using the MTS 370.25 low-temperature material testing machine, and the test equipment is shown in Figure 2. The soil samples are placed into Ф50 mm × 100 mm cylindrical specimens and placed in a special incubator for curing.

2.3. Tests’ Plan
The loading and unloading process of the MTS testing machine is controlled by program. The compression strength tests of frozen soil adopt the displacement rate loading mode, and the displacement change rate is 1 mm/min. During the test, the axial strain reaches 20% or the peak stress decreases by 20%, and the test is terminated. The compression strength of frozen soil at −5°C, −7°C, −10°C, −15°C, and −20°C is measured by experiments, and the elastic modulus and Poisson’s ratio are calculated at the same time.
The triaxial shear tests adopt the strain loading method, and the strain rate is 1%/min. If the axial force has a peak value, after the peak point, shear 3%∼5% of the strain value, and then terminate the test. If the axial force continues to increase, then the test is terminated after shearing to 20% strain value. Four freezing temperatures of −5°C, −7°C, −10°C, and −15°C are adopted, and the temperature error is controlled within 0.1°C during the test. The test confining pressures are 0.5 MPa, 1.0 MPa, and 1.5 MPa, respectively. The data acquisition frequency is 2 Hz. In the creep test, the creep loading coefficient is 0.3, 0.5, and 0.7, and the test time is 12 h.
3. Test Results and Analysis
3.1. Uniaxial Compressive Strength Tests
3.1.1. Uniaxial Compressive Strength
The uniaxial compressive strength tests of frozen soil at five temperature levels of −5°C, −7°C, −10°C, −15°C, and −20°C were carried out. The test results are shown in Table 2.
Due to the existence of ice crystals in frozen soil, its structure, strength, and other characteristics are different from those of conventional soil. The change of temperature directly affects the compressive strength of frozen soil. With the decrease of temperature, the internal composition of frozen soil gradually changes, and unfrozen water changes into ice crystals, which makes the compressive strength of soil increase. The relationship curve between uniaxial compressive strength and temperature of frozen soil is shown in Figure 3.

In the temperature range from −5°C to −20°C, the relationship between uniaxial compressive strength and temperature of frozen soil can be fitted by a linear function, and the correlation coefficient R = 0.99, indicating that the compressive strength of frozen soil increases almost linearly with the decrease of temperature.
3.1.2. Elastic Modulus
The elastic modulus of frozen soil obtained from the test is shown in Table 3. The relationship curves between elastic modulus and temperature of frozen soil samples are shown in Figure 4.

From Figure 4, it can be concluded that the elastic modulus of frozen soil increases with the decrease of freezing temperature in the range of test temperature, and the elastic modulus increases linearly with the decrease of freezing temperature, and the correlation is good. Every time the temperature drops by one degree, the elastic modulus increases by 8.70 MPa on average. In the test temperature range, the elastic modulus of the corresponding soil layer at any temperature can be calculated by the interpolation method.
3.1.3. Poisson’s Ratio
The test results of Poisson’s ratio of frozen soil are shown in Table 4, and the curve of Poisson’s ratio and temperature is shown in Figure 5.

From Figure 5, it can be concluded that, within the test temperature range, Poisson’s ratio of the experimental frozen soil increases with the increase of freezing temperature, and it increases almost linearly. With one degree rise, Poisson’s ratio of frozen soil increases by 0.0053 on average.
3.2. Triaxial Shear Test
The stress-strain curves of frozen soil under different temperatures and confining pressures are shown in Figure 6.

(a)

(b)

(c)

(d)
The triaxial shear process of frozen soil can be roughly divided into three stages: elastic growth stage, plastic yield stage, and accelerated failure stage. At the initial stage of loading, the stress-strain curve approximately shows a linear growth trend. At this time, the deformation of the soil is mainly through the compression between the particles in the soil, which is the elastic growth stage. With the continuous increase of stress, the soil sample begins to produce irreversible plastic deformation, and the plastic deformation continues to increase. When the maximum stress is reached, fine cracks begin to appear on the surface of the soil, which is the plastic yield stage. After the sample reaches the maximum stress, the soil crack develops, extends, and penetrates, leading to the failure of the sample, which is the accelerated failure stage.
The triaxial shear strength parameters of frozen soil are shown in Table 5.
3.3. Triaxial Creep Tests
The relationship between triaxial shear strain and time of frozen soil under different temperatures is obtained by experiments. The creep curves are shown in Figure 7.

(a)

(b)

(c)

(d)
The creep of frozen soil is a typical nonlinear rheology; that is, the isochron curve of stress-strain is not a straight line or broken line. The change of stress level and time can cause nonlinearity, and the influence of stress and time on nonlinearity is coupled. Under the condition that the stress applied to the frozen soil remains constant, the strain will increase with time. In the beginning, the frozen soil shows instantaneous elastic and plastic deformation, followed by the creep stage, and enters the first stage, i.e., the unstable creep stage. The creep of frozen soil presents an attenuating deformation; that is, the deformation will gradually approach a stable value. After that, the strain rate remains constant approximately and enters the second stage, namely, the stable creep stage. When the stress applied to frozen soil is small, the second stage of creep lasts longer, and even the third stage (accelerated creep) does not appear. On the contrary, the second stage of creep is very short or even disappears.
4. Fractional Derivative Viscoelastic-Plastic Creep Model
Due to the advantages of fractional calculus in nonlinear dynamic systems [36, 37], fractional calculus is introduced to describe the nonlinear mechanical behavior of frozen soil creep. In this paper, an improved fractional-order viscoelastic-plastic creep constitutive model for frozen soil is proposed. On the one hand, the fractional viscous body is used to replace the classical viscous body in the ideal viscoplastic body to obtain the improved viscoplastic body; on the other hand, the influence of stress on the creep characteristics of frozen soil, especially the accelerated creep, is considered. A fractional viscoplastic body with an order greater than 1 is proposed, and its order can be changed according to the stress level, which makes it sensitive to the change of stress. Finally, the classical elastic body, classical viscous body, improved fractional-order elastoplastic body, and improved fractional-order viscoplastic body are connected in series, and the unsteady fractional-order differential integral creep model of frozen soil is obtained and verified.
4.1. Riemann–Liouville Fractional Differential
Firstly, Γ (γ) gamma function closely related to the fractional derivative is introduced:
For general continuous functions, the following Riemann–Liouville fractional differential is introduced. Let γ be any real number and n be an appropriate integer satisfying n − 1 < γ ≤ n; then, Riemann–Liouville fractional differential is defined as
When 0 <γ < 1 and N = 1,
4.2. Fractional Viscous Body
4.2.1. Steady Fractional Viscous Body
The stress-strain relationship of the ideal gluing pot element satisfies Newton’s law.where σ is the stress of the material; ε is the strain; η is the viscosity coefficient; and t is the time.
In essence, the rheological model theory using the fractional derivative is obtained by replacing Newton viscous fluid with Abel body in classical model theory or replacing Newton viscous fluid with the Abel viscosity pot. The fractional calculus viscous body is shown in Figure 8.

The constitutive relation of the fractional viscous body is as follows:where γ is the fractional differential order and ηγ is the viscosity coefficient.
When γ = 0 and γ = 1, the fractional viscous body can represent the ideal elastic body and the ideal fluid material, respectively; when 0 <γ < 1, it can describe the mechanical properties of materials with properties between them.
When the stress σ is constant, according to the Riemann–Liouville fractional calculus theory, the fractional integration of equation (5) is carried out.
4.2.2. Unsteady Fractional Viscous Body
Based on the theory of fractional calculus, a kind of fractional-order viscous element with unsteady differential order is proposed. When the stress of frozen soil exceeds its yield limit, the frozen soil enters into the accelerated creep stage. The relationship between fractional differential order γ and stress σ is as follows:where is the positive real number related to the nature of frozen soil; σs is the yield limit of frozen soil; and σ0 is the unit stress with the value of 1, and the dimension is the same as that of σ.
The improved fractional calculus viscous body is shown in Figure 9.

When the stress σ is constant, according to the Riemann–Liouville fractional calculus operator theory, the fractional integral is obtained:
Therein,
4.3. Fractional Derivative Viscoelastic-Plastic Creep Model
The fractional derivative viscoelastic-plastic creep model is shown in Figure 10.

The nonlinear model consists of the Maxwell model (part 1), a steady fractional viscoelastic body (part 2), and an unsteady fractional-order viscoplastic body (part 3).
The creep equation of the Maxwell model is as follows:where E1 is the elastic modulus and η1 is the viscosity coefficient.
The creep equation of the steady fractional viscoelastic body is as follows:where η2 is the viscosity coefficient.
The creep equation of the unsteady fractional viscoplastic body is as follows:
Therein,where η3 is the viscosity coefficient and σs is the yield stress.
When , the third part is a rigid body, which does not affect the constitutive model.
When , the switch of part 3 is turned on, which can describe the accelerated creep stage of frozen soil.
4.4. Parameters’ Inversion and Model Verification
According to the creep test results, the model parameters of frozen soil creep test curves under different temperatures are identified. Based on the particle swarm optimization algorithm and the least square fitting function parameters, the results are shown in Table 6.
Figure 11 shows the creep test results and the fitting curves of the creep model of artificial frozen soil under different temperature conditions. It can be seen that the theoretical model proposed in this paper can well describe the test results under different temperature conditions, including the complete creep process such as initial creep, steady-state creep, and accelerated creep, especially the accelerated creep stage. The high degree of fit indicates the correctness and applicability of the nonlinear viscoelastic-plastic model established in this paper.

(a)

(b)

(c)

(d)
5. Finite Element Calculation
5.1. Modeling
Using ANSYS finite element calculation software and the viscoelastic-plastic constitutive equation of frozen soil derived in this paper, the displacement field and stress field of the frozen wall are simulated. SOLID186 three-dimensional solid element is adopted, and PCG solver is selected. Due to the long distance of the connecting passage, the large volume of freezing, long excavation and construction time, the complex surrounding environment, and high requirements for deformation control, the plastic and creep deformation of frozen soil should be considered in the calculation. Therefore, the fractional derivative viscoelastic-plastic model is used to improve the calculation. The three-dimensional mechanical model is used for the mechanical analysis of the frozen curtain. The material parameters are shown in Table 7, and the mechanical characteristic parameters of frozen soil are taken as the mechanical properties of frozen soil at −10°C. According to the structural symmetry, 1/4 of the structure is taken as the calculation model, as shown in Figure 12.

5.2. Displacement and Creep
The total displacement and frozen curtain displacement caused by the connecting passage excavation are calculated, respectively. The calculation results are shown in Figures 13 and 14. From the total displacement nephogram, the maximum deformation of the ground up the middle of the connecting passage is 34.3 mm; from the displacement nephogram of the frozen curtain, the maximum deformation of the frozen curtain is 21.2 mm, and the frozen soil deformation accounts for 61.8% of the total ground deformation.


Taking the center of the tunnel as the origin and the distance from the center of the tunnel as the variable, the displacement curve of the upper frozen curtain of the connecting passage at different times is shown in Figure 15. Taking the frozen soil in the middle of the connecting passage as the research object, the creep law of frozen soil within 24 hours after the excavation is shown in Figure 16.


It can be seen from Figures 15 and 16 that, at the same time, the displacement curve of the frozen curtain at the upper part of the connecting passage is parabolic, and the maximum vertical displacement occurs in the middle of the connecting passage. The vertical displacement increases with time, and the closer to the middle, the greater the increase. After excavation, the calculated maximum vertical displacement is 22.5 mm after 8 hours, 23.2 mm after 20 hours, and 23.3 mm after 24 hours. The vertical displacement increases by 0.74 mm, 0.58 mm, 0.31 mm, 0.20 mm, 0.16 mm, and 0.13 mm, respectively, after every 4 hours, and the creep increase decreases gradually, showing an attenuation creep increase, and gradually tends to be stable. In the actual project, the tracking grouting was carried out according to the monitoring value. Therefore, the later measured values are smaller than the calculated values.
5.3. Stress and Strain
The von Mises stress distribution nephogram of frozen soil after the excavation is shown in Figure 17. The von Mises stress of the frozen soil in the middle-upper of the connecting passage within 24 h is shown in Figure 18.


It can be seen from the von Mises stress distribution nephogram of frozen soil that the stress distribution is generally uniform, most of which are in the range of 0.037 MPa–0.663 MPa. The maximum stress is 1.45 MPa, which is the stress concentration point and is located at the junction of the frozen soil and tunnel. It can be seen from Figure 18 that the von Mises stress of the frozen soil in the middle-upper part of the connecting passage presents a slow downward trend during the frozen soil creep, with a small decline, and finally tends to be stable.
According to the strain of the frozen curtain in the final stable state, the whole frozen curtain is divided into the elastic zone, plastic zone, and damage zone, as shown in Figure 19.

(a)

(b)

(c)
It can be seen from Figure 19 that the elastic zone, plastic zone, and damage area of the frozen curtain are cross distributed. Most of the frozen soil is in the elastic zone, and the stress-strain relationship of frozen soil in the elastic zone is linear elastic. The proportion of the plastic zone is relatively small, mainly concentrated in the inner curtain and bottom plate of the passage, and the stress-strain relationship of frozen soil in the plastic zone is nonlinear. The damage zone is very small, at the corner of the passage floor, and it is mainly caused by stress concentration.
During the excavation of the connecting passage, the deformation monitoring value of the ground and frozen soil excavation surface is the same as the model calculation value. Taking into account the requirements of deformation control, according to the monitoring and model calculation conditions, the deformation can meet the control requirements by grouting and temporary support. Therefore, this model has important engineering application value and can provide important guidance for freezing engineering.
6. Conclusions
Through the uniaxial compressive strength test, triaxial shear test, and triaxial creep test of frozen silt with the sand soil layer, many important physical and mechanical properties’ parameters of the soil layer are obtained, which provide basic conditions for the development of freezing engineering. By a series of triaxial creep tests for artificial frozen soil under different stresses and temperatures, a fractional-order constitutive model of frozen silt with sand is proposed. The following conclusions could be acquired:(1)The creep of frozen soil is a typical nonlinear rheology; that is, the isochron curve of stress-strain is not a straight line or broken line. The change of stress level and time can cause nonlinearity, and the influence of stress and time on nonlinearity is coupled. When the stress applied to frozen soil is small, the second stage of creep lasts longer, and even the third stage (accelerated creep) does not appear. On the contrary, the second stage of creep is very short or even disappears.(2)The proposed fractional derivative viscoelastic-plastic creep model can describe the three stages of artificial frozen soil creep, especially the accelerated creep stage. The creep curves under different temperature and stress conditions are fitted and analyzed by the particle swarm optimization algorithm and least square algorithm, and the model parameters are obtained.(3)Comparing the calculated results with the tested ones, it is found that the proposed fractional-order model can simulate the properties of frozen silt with sand very well. Compared with the other creep models of frozen soil, the proposed model here has higher accuracy and strong stress sensitivity.(4)Through the simulation of the connecting channel construction and comparison with the observed data in site, the effectiveness of the model was verified, which further verifies that the model has important guiding significance for the freezing project.
Data Availability
Some or all data, models, or codes that support the findings of this study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was supported by the National Natural Science Foundation of China (Grant nos. 42061011, 41977236, and 41672278), the Natural Science Foundation of Jiangxi Province (Grant no. 20192ACBL20002), the Doctoral Research Initiation Fund of East China University of Technology (Grant nos. DHBK2019229 and DHBK2019251), and the Science and Technology Project of XPCC (Grant no. 2020AB003). The authors gratefully acknowledge the support.