1 INTRODUCTION

To meet the needs of economic development, an increasing number of tunnel projects have been successively constructed under large burial depths and high crustal stresses. Under the condition of high crustal stress, rock is usually weak and displays obvious creep characteristics. Under this geological condition, it is difficult to predict the deformation of the surrounding rock during tunnel construction. Under the action of long-term loading, the deformation of the surrounding rock easily exceeds the engineering design value, resulting in initial tunnel support failure and even secondary lining cracking. Therefore, it is of great theoretical value and guiding importance for the construction of deep tunnel engineering and the deformation control of surrounding rock to study the aging mechanical properties of rock under high crustal stress and develop a reasonable numerical model of surrounding rock failure.

In recent years, many experts and scholars have conducted much research on the creep characteristics of rocks under long-term loads. In terms of research on the creep damage characteristics of rock, Ma et al. (2021) studied the effects of creep and fatigue on the hysteresis behavior, deformation ability, progressive damage, cycle life, and fracture mechanism of rock salt through laboratory tests. X. Li et al. (2022) explored the damage characteristics of high-temperature sandstone through a triaxial creep test. X. Yang and Jiang (2022) observed the creep characteristics of slate under different bedding angles and different numbers of freeze–thaw cycles. W. Li et al. (2007) studied the creep characteristics of composite hard and soft rock through laboratory tests and modified the Burgers model. Vladimir et al. (2022) constructed a porous elastic rheological damage model based on the nonlinear anisotropic damage rheological model. Wang et al. (2021) focused on the creep characteristics of salt rock under long-term cyclic loading. Zhang et al. (2021) studied the creep damage characteristics of granite under the action of wetting and drying cycles. Zhao et al. (2021) proposed the creep damage accumulation model of rock by the nonlinear summation method. Bai et al. (2021) analyzed the damage characteristics of rock under the action of freeze–thaw cycles.

In terms of applying the creep model, Ma et al. (2017) introduced nonlinear elements into the Burgers model and applied the resulting model to the stability analysis of a storage cavern. Siddiquee et al. (2015) used a complex nonlinear finite element method to analyze the creep characteristics of geosynthetic-reinforced soil retaining walls. Wang et al. (2022) constructed the nonlinear Maxwell model according to the creep test results of salt rock and applied the creep model to an engineering analysis of salt rock. Xu et al. (2019) studied the creep characteristics of melite through tests and applied the creep constitutive equation to numerical calculations. Cheng et al. (2021) applied the established fractional nonlinear creep constitutive equation to an engineering stability analysis of salt rock. Wu et al. (2022) improved the accuracy of creep parameters through acoustic emission tests and applied a creep equation to predict rock failure. Gutiérrez-Ch et al. (2022) combined the rate process theory with the discrete element method to analyze the surrounding rock stability of a deep tunnel. Zhao et al. (2021) proposed a new damage accumulation model and applied it to the long-term stability analysis of a compressed air energy storage salt cavern. Sun et al. (2022) constructed the nonlinear creep equation of granite and applied it to the stability analysis of a gold mine roadway. Thus far, rich achievements have been made in the research of creep damage characteristics and the application of constitutive models. However, most creep damage models are based on creep loading tests without creep unloading analysis, which cannot fully analyze the deformation properties of creep curves at various stages.

In view of this, with the red sandstone in the Yuezhishan Tunnel of China as the research object, this study performed triaxial graded loading–unloading creep tests under different stress levels. According to the test results, the relation function between the steady-state creep rate and stress level was constructed, and the long-term strength of the red sandstone was determined; meanwhile, the NCD constitutive equation was established. On this basis, the numerical NCD constitutive model and the BL model were developed by FLAC3D and verified by laboratory tests. Finally, the two models were jointly applied to a creep deformation and stability analysis of the surrounding rock of the Yuezhishan Tunnel.

2 THE OVERVIEW OF THE YUEZHISHAN TUNNEL PROJECT

2.1 Geological conditions

The Yuezhishan Tunnel is located along the Dadu River in the Jinkou District of Leshan city, Sichuan Province. The geomorphology of this section can be divided into two geomorphological units: tectonic denuded meso-alpine and valley terrace. The former unit runs through the Sichuan–Yunnan meridional tectonic system and is partially controlled by the combination of this system and the NW-trending Yunnan–Myanmar tectonic system on the Qinghai–Tibet Plateau. The western and eastern parts (bounded by the Anning River) belong to the snow mountain system and the Daliang Mountain system, with towering mountains. The river valley terraces are mainly Anning River valleys, with sharp river cutting and great differences in height. The geographical location of the Yuezhishan Tunnel is shown in Figure 1.

The Yuezhishan Tunnel has a total length of if 14 085 m and a maximum burial depth of 1810 m. The surrounding rock is red argillaceous sandstone with developed joints and fractures, and the local surrounding rock is broken. The stress of the surrounding rock is tested on-site by the hydraulic fracturing method, and the results show that the minimum and maximum principal stresses on the horizontal hole of the surrounding rock are 20 and 40.2 MPa, respectively. Because of the weak and broken surrounding rock and the high stress of the primary rock, the surrounding rock, which has been damaged many times, produces large rheological deformation and lining cracks. Accordingly, the local steel arch is severely distorted as shown in Figure 1.

2.3 Basic mechanical properties of the samples

To study the rock physical and mechanical characteristics of the Yuezhishan Tunnel, red sandstone pillars were drilled from deep depths. The rock pillars were processed and polished into Φ50 mm × 100 mm and Φ50 mm × 25 mm standard samples. Then uniaxial and triaxial compression tests were performed on the samples using a TAW2000 electrohydraulic servo rigidity test press as shown in Figure 2a. The test results show that the uniaxial compressive strength σc of the rock samples is 25.2 MPa, the tensile strength σt is 2.2 MPa, the cohesion c is 2.2 MPa, the internal friction angle ϕ is 35°, Young's modulus E is 6.2 GPa, Poisson's ratio υ is 0.21, and the compressive strength σ20 is 58.8 MPa.

3 CREEP TEST

3.1 Graded loading–unloading creep test procedure

To study the creep characteristics of the surrounding rock of the Yuezhishan Tunnel, a triaxial creep test of the red sandstone was performed. A global digital systems triaxial rheological tester was used to perform the graded loading‒unloading creep test on the red sandstone as shown in Figure 2b.

The confining pressure was set at 20 MPa. The minimum horizontal principal stress was measured in the field, and the axial deviatoric stress loading was divided into 8 levels (5, 10, 15, 20, 25, 27, 30, and 35 MPa). During the test, the confining pressure was first applied at a rate of 0.5 MPa/s. When the confining pressure was stable, the deviatoric stress was applied at a rate of 100 N/s. Once the target stress value was reached, the pressure was stabilized for 96 h, and then unloading was performed to maintain a contact stress of 100 Pa between the pressurized plate and the standard samples. The unloading deformation of the specimen within 24 h was recorded as shown in Figure 3.

According to Figure 3, the creep curve trend and particularly the creep characteristics of the red sandstone vary under different axial stresses. When the deviatoric stress is 5–20 MPa, the creep behavior of the red sandstone is viscoelastic, and the steady-state creep rate is 0. The instantaneous elastic and viscoelastic deformation recover after unloading, but residual deformation still persists, indicating that plastic deformation occurs during the instantaneous loading process. When the deviatoric stress is 25–30 MPa, the creep behavior of the red sandstone is viscoelastic‒plastic, and the steady-state creep rate increases with increasing stress level. The instantaneous elastic deformation recovers after unloading, and the instantaneous plastic deformation and the viscoplastic deformation generated during the creep process are residual. When the deviatoric stress is 35 MPa, the creep of the red sandstone also exhibits viscoelastic‒plastic behavior, and the creep acceleration stage appears. However, the sample is damaged if the loading process exceeds approximately 10 h, indicating that the stress level at this time far exceeds the starting value of the viscoplastic deformation.

3.2 Long-term strength of the red sandstone

In terms of rock engineering, long-term strength is a key criterion for the long-term stability of the surrounding rock and a landmark parameter for rock creep characteristics (Okubo et al., 2010). At present, the commonly used methods to determine long-term strength include the transition creep method and the isochronal curve method (Innocente et al., 2021). However, the transition creep method has difficulty determining the specific value of the long-term strength and can only mark out the range of the long-term strength. The isochronal curve method requires subjective judgment when determining the inflection point of an isochronal curve, and accordingly the long-term strength obtained is not scientific.

According to the creep curve, in the process of creep loading, the red sandstone mainly exhibits viscoelasticity at low stress levels, only decay creep and steady-state creep occur, and the steady creep rate is 0. When the stress level is high, the red sandstone behaves as viscoelastic‒plastic material, and the steady-state creep rate is a constant value greater than 0. If the loading time is long enough, the creep curve will accelerate, causing the sample to be damaged. Therefore, the minimum stress level that makes the steady-state creep rate greater than 0 is the long-term strength, which is also the threshold for the viscoplastic deformation of rock. The long-term strength of the red sandstone can be solved by constructing the relationship between the stress levels and steady-state creep rates. The steady-state creep rate values at high stress levels are shown in Table 1. Since the structure of the sample is severely damaged during level 8 loading (axial stress, 55 MPa), the measured stress variables include not only the creep deformation of the specimen but also the staggered deformation of the rock block. As a result, this level is omitted, and only stress levels 4–7 are taken. The steady-state creep rate of each stress level was fitted, and the fitting curve is shown in Figure 4. The mathematical relationship between the deviatoric stress and steady-state rate is as follows:

(1)

where is the steady-state creep rate; σ denotes the axial deviatoric stress; and σ s refers to the long-term strength of the red sandstone under a confining pressure of 20 MPa.

Table 1. Steady-state creep rates.

Axial stress (MPa)	Deviatoric stress (MPa)	Steady-state creep variable in 96 h (10−4)	Steady-state creep rate (10−6∙h−1)
40	20	0.43	0
45	25	1.23	0.47
47	27	2.98	3.58
50	30	6.23	7.25

According to Equation (1), the deviatoric stress value corresponding to a steady-state creep rate equal to 0 can be calculated to be 23.76 MPa. Therefore, the long-term strength of the red sandstone is 43.76 MPa, which is approximately 74.2% of the compressive strength σ20, consistent with the empirical value of 60%–90% (Shen & Chen, 2011), and close to the maximum principal stress value measured in the field.

4 DEVELOPMENT OF THE NONLINEAR CREEP DAMAGE MODEL

4.1 Derivation of the constitutive equation

Numerous studies have shown that the creep process of rock is also a process of damage accumulation, and the viscoplastic deformation of rock exhibits nonlinear characteristics under long-term loading (X. Li et al., 2022; Ma et al., 2021; X. Yang & Jiang, 2022). However, the current classical element model is mostly used to describe the linear constitutive relationship. Although it can characterize decay creep and steady-state creep, it cannot reflect the accelerated creep process of rock. Therefore, in this paper, the damage factors in the creep process are considered, and the nonlinear damage viscous element with damage variable Ds is introduced to describe the nonlinear viscoplastic deformation of the red sandstone.

The element is essentially a nonlinear Newton body. According to Newton's law of viscosity, its constitutive equation is

(2)

where is the creep rate and η 1( D s) is the damage viscosity coefficient. The damage viscosity coefficient can be expressed as follows:

(3)

where η 1 is the initial viscosity coefficient at the viscoplastic stage and D s is the damage variable. Accordingly, the damage variable has the following relationship with time (W. Yang et al., 2014):

(4)

where α is the coefficient related to rock material properties; t denotes time; and 0 <  D s < 1. With Equations ( 3) and ( 4) being substituted into ( 2), the following Equation ( 5) can be obtained:

(5)

A physical model for describing the nonlinear creep of the red sandstone can be constructed by combining the results of the loading‒unloading creep test with the nonlinear damage viscous element. The NCD model includes an instantaneous elastic element, an instantaneous plastic element, a viscoelastic element, and a viscoplastic element, all of which are connected in series as shown in Figure 5.

In Figure 5, the instantaneous elastic element is used to describe the instantaneous elastic stage in the creep test curve. The instantaneous plastic element is used to describe the instantaneous plastic stage. The viscoelastic element is used to describe the decay creep stage. The viscoplastic element is used to describe the steady-state creep and accelerated creep stages. The constitutive relations of each element are described below.

• 1.

  Instantaneous elastic element

For the three-dimensional elastic element, its stress tensor σij can be decomposed into a spherical stress tensor δijσ m and a deviatoric stress tensor sij, which change the volume and shape of the object, respectively. The instantaneous elastic element strain is expressed as follows:

(6)

where K and G denote the bulk modulus and the shear modulus, respectively.

• 2.

  Instantaneous plastic element

According to the loading‒unloading creep test curve of the red sandstone, obvious instantaneous plastic deformation occurs in the sandstone sample at all stress levels. Therefore, by assumption, the threshold of the plastic slider is small, and the element will deform under any stress level. In this case, the constitutive equation of the instantaneous plastic element can be obtained by analogy with the constitutive equation of the elastic element. Thus, the instantaneous plastic element strain can be expressed as follows:

(7)

where and are the deviatoric and spherical stress tensors of the plastic slider, respectively; and is the mean normal stress of the slider. When i =  j = 1,

(8)

• 3.

  Viscoelastic element

The spherical stress tensor mainly causes the volume change of the sample. However, in the viscoelastic stage, the sample volume almost remains constant. Accordingly, the creep deformation caused by the spherical stress tensor can be neglected. The deformation of the viscoelastic element is mainly contributed by the deviatoric stress tensor. The viscoelastic element strain is expressed as follows:

(9)

• 4.

  Viscoplastic element

The plastic flow law should be considered to solve the constitutive equation of the viscoplastic element under three-dimensional stress. Therefore, the constitutive equation of the viscoplastic element can be expressed as follows:

(10)

where denotes the strain rate of the viscoplastic element; F refers to the yield function; Q represents the plastic potential function; F 0 denotes the initial value of the yield function; The P function has the following properties:

(11)

For rock materials, x is usually 1. When F ≥ 0, it can be known from relevant flow laws that F =  Q. The viscoplastic element strain can be obtained by integrating both sides of Equation ( 10).

(12)

The initial value of yield function F 0 is set to 1, and according to the Mises yield criterion

(13)

where J 2 is the second invariant of the deviatoric stress tensor, MPa.

Substituting Equation ( 13) into Equation ( 12) yields the following:

(14)

Since each element is connected in series, the total strain variable of the model should be the sum of the strain variables of each element. The constitutive equation of the NCD model can be obtained by combining Equations ( 6), ( 7), ( 9), and ( 14). The total strain εij is expressed as follows:

(15)

(16)

In the triaxial compression test, the three-dimensional normal stress is the three-dimensional principal stress, which is σ 1, σ 2, and σ 3, and σ 2 =  σ 3. Accordingly, it can be obtained when i =  j = 1:

(17)

By substituting Equations ( 8) and ( 17) into Equations ( 15) and ( 16), the axial constitutive equations of the NCD model can be obtained as follows:

(18)

(19)

4.2 Determination of the NCD model parameters through the inversion method

To further determine the unknown parameters in the NCD model, the model parameters must be inverted by using the triaxial creep test results. At present, for nonlinear function parameter inversion problems, the least-squares method is most commonly used. Since many unknown parameters are in the model, MATLAB is used to determine the NCD model parameters so as to avoid setting their initial values. The inversion results of the parameters at all stress levels are shown in Table 2, and the fitting curves are shown in Figure 6.

Table 2. Parameters of the NCD model.

Deviatoric stress (MPa)	K1 (GPa)	G1 (GPa)	K2 (GPa)	G2 (GPa)	G3 (GPa)	η1 (GPa·h)	η2 (GPa·h)	α (10−3)	R2
5	11.14	33.26	11.67	34.39	39.32		122.3		0.9822
10	7.96	19.73	8.29	20.51	85.28		506.9		0.9583
15	6.87	6.08	6.84	3.91	17.41		225.2		0.9857
20	6.36	8.16	6.45	8.41	55.11		461.3		0.9755
25	7.83	5.80	7.83	6.21	65.74	300.0	300.0	5.420	0.9933
27	5.46	6.05	5.46	6.05	36.38	350.0	350.0	1.160	0.9856
30	5.23	5.23	5.23	5.22	18.60	337.9	193.1	0.001	0.9959
35	50.50	1.19	50.50	50.50	5.56	6867.0	42.5	676.200	0.9986

Although this model contains seven elements, compared with other models, it has fewer parameters, thus being more convenient to select, and better describing the properties of rock at various creep stages. According to Table 2 and Figure 6, all correlation coefficients of the creep test curves fitted by the NCD model are above 0.95, and the fitted curves are highly consistent with the test data, indicating that the NCD model can better reflect the creep characteristics of the red sandstone.

5 IMPLEMENTATION OF THE ALGORITHM FOR NUMERICAL SIMULATIONS

Based on the FLAC3D built-in Burgers creep constitutive model, the NCD model was embedded into the calculation program by using C++ to construct the numerical NCD constitutive model. At the same time, the combined use of the BL model was developed. The models were tested by experimental simulation and engineering calculation.

5.1 Secondary development of the NCD and BL models

5.1.1 Burgers model in FLAC3D

To solve the creep problem of geotechnical engineering, the Burgers (Kelvin model and Maxwell model in series) model is set in FLAC 3D. In the calculation process, the program iterates the strain increment continuously to solve the stress increment and combines the time increment (time step) Δ t to form the Burgers model algorithm. For the Kelvin model,

(20)

where u k, G k, and η k represent the strain, shear modulus, and the viscosity coefficient of the Kelvin model, respectively; F is the stress of the Burgers model; and the upper corner scripts ′ and o represent the old and new values before and after the change in Δ t, respectively.

For the Maxwell model,

(21)

where u m, G m, and η m represent the strain, shear modulus, and viscosity coefficient of the Maxwell model, respectively.

The strain increment Δ u of the Burgers model is shown below:

(22)

By substituting Equations ( 20) and ( 21) into Equation ( 22), the stress‒strain increment relation of the numerical Burgers model can be obtained as follows:

(23)

where

(24)

In the process of program calculation, the time step Δt and strain increment are automatically set; the initial value of uk is set to 0; and uk is continuously adjusted through Equation (20).

5.1.2 Numerical integration of the NCD constitutive model

According to the numerical Burgers model algorithm, although FLAC3D embeds the time parameter Δt into the calculation process through the constitutive relation of the Burgers model, the Burgers model has no plastic element, and the plastic state of the element cannot be characterized in the model calculation process. Moreover, the constitutive equation of the Burgers model is linear and cannot simulate the accelerated creep stage.

To solve these problems, based on the NCD constitutive model, the numerical Burgers model algorithm was improved through the Visual Studio development environment. The relationships among stress, strain, and time were adjusted, and the creep yield failure criterion was set. The improved principle and process are shown in Figure 7.

• 1.

  The new project was registered and named through Visual Studio, and the source and header files of the Burgers model provided by FLAC3D were imported into this project.

• 2.

  The NCD model parameters were declared, namely, the parameters shown in Table 2, and followed by parameter initialization. Accordingly, the user could assign values to the parameters in FLAC3D.

• 3.

  The NCD constitutive relation was written into the calculation program. At step i (i = 1, 2…), whether the set creep time had been reached was determined. If so, the calculation was finished; otherwise, it entered the constitutive regulation module. In the constitutive regulation model, whether F ≥ 0 was judged.

  
  

• a.

  If so, the constitutive relation containing the viscoplastic element was run, namely, Equation (16) (the formula needed to be converted into strain increment and stress increment format when it was introduced into the running program).

• b.

  Otherwise, the constitutive relation of the nonviscoplastic element was operated, which was Equation (15).

• 4.

  This step visualized the unit failure state. When the strain rate exceeded 10 (10−6·h−1), the unit would yield failure. At this time, the characteristic state of the unit changed, which could be directly reflected in the generated graph. Notably, different rock materials have different failure creep rates, which should be set according to the rock creep test results.

• 5.

  The prepared constitutive model was encapsulated as a.dll file and installed into the exe64/cmodel folder in FLAC3D.

The core part of the program is composed of the constitutive regulation module and the unit state representation module. The former module is used to adjust the relationships among stress, strain, and time, and the latter module makes the calculation result more intuitive. After the program is installed, the NCD model can be directly invoked in FLAC3D, which is quick to use and convenient for engineering analysis.

5.1.3 Development of the BL model

The liner model is commonly used in FLAC 3D to simulate the tunnel lining structure. However, liner units are rigidly connected and cannot fracture, so the simulated lining is an elastic structure and cannot yield failure. To reflect the failure characteristics of the lining structure under creep action more accurately, the liner model was developed by using the Fish language to construct the BL model, with the specific process as follows (Sun et al., 2021).

• 1.

  The zone model was used to build the surrounding rock; the coordinates of the zone nodes were obtained by the Fish language and then stored in the Table database.

• 2.

  The independent liner units were built on the surface of the rock based on the coordinate information in the Table database. The liner units had different IDs (1, 2…, n), as shown in Figure 8a. At this time, the nodes of the liner units and the grid points of the zone units were connected through the link model, and the nodes of two adjacent liner units overlapped each other but were not connected.

• 3.

  The links between the liner units and the zone units were deleted; and the new links between the nodes of two adjacent liner units and between units and grid points were constructed by the Fish language, as shown in Figure 8b. This method realized the linkage between the zone model and liner model and established the stress transfer relationship between them.

• 4.

  Since the coordinate system of nodes of the liner units changes with the coordinate system of the surrounding rock, the coordinate system of nodes must be unified to facilitate the subsequent assignment of the link model. The node coordinate system was adjusted by the Fish language, as shown in Figure 8c. As was specified, the tangential direction of the liner nodes was the X direction, the axial direction was the Y direction, and the radial direction of the tunnel was the Z direction.

• 5.

  The links between nodes were assigned according to the node coordinate system, which included compressive strength, tensile strength, shear strength, flexural strength, and tensile strength.

The BL model automatically monitors the stress state of the links in the calculation process. When the stress of a link reaches the failure strength, the link is deleted by the Fish language to simulate the lining rupture.

5.2 Numerical simulation results and model validations

5.2.1 Validation of the numerical NCD model

To further test the simulation quality of the numerical NCD model, FLAC3D was used to perform the triaxial creep simulation test and compare it with the laboratory test. The model size and boundary conditions are shown in Figure 9. The model has a total of 20 736 zone units, with displacement constraint applied to the bottom, axial load σ1 applied to the top, and confining pressure σ3 applied to the circumference.

In the calculation process, the Burgers and NCD models were assigned to the sample model, and the deformation laws of the monitoring point on the top of the sample under different deviatoric stresses were monitored and recorded. After the calculation, the calculated data were extracted and compared with the laboratory test results as shown in Figure 10.

The comparison of the calculation results with the test results shows that the numerical NCD model is highly consistent with the laboratory test results. Therefore, it can accurately simulate the creep characteristics of the rock. However, the curve calculated by the Burgers model in the accelerated creep stage deviates considerably from the laboratory test curve. Accordingly, the numerical NCD model is more suitable for rock engineering creep failure analysis than the Burgers model.

5.2.2 Validation of the BL model

Concrete cube samples of 150 mm × 150 mm × 150 mm and cuboid samples of 150 mm × 150 mm × 550 mm were made indoors, and then concrete compression and flexural tests were performed on them. To test its rationality, the BL model was used for experimental simulation and compared with the laboratory test as shown in Figure 11.

In the simulation test, the load P was gradually increased until the model was damaged. The loading data were recorded and compared with the laboratory test data. As shown in Figure 11a, the model was damaged when the load P reached 32.2 kN, which was close to the value of 32.4 kN obtained in the laboratory flexural test, and the error was 0.62%. As shown in Figure 11b, in the compressive simulation test, the link in the middle of the concrete model was damaged when the load P reached 27.5 MPa, which was close to the value of 27.8 MPa obtained in the laboratory compressive test, and the error was 1.08%. As is shown by the comparison of the laboratory test results and those of numerical calculation, the BL model could simulate the failure characteristics of the actual concrete structure.

6 DISCUSSION

The creep deformation of the surrounding rock of the Yuezhishan Tunnel mainly occurs in the DK246 + 130–DK246 + 252 section, with a total length of approximately 120 m. According to the field stress distribution and construction scheme, a numerical model was constructed to analyze the creep deformation and the support failure characteristics of the surrounding rock. The model is shown in Figure 12, in which the numerical NCD model was selected as the constitutive model, the Cable model as the bolt support, the Beam model as the arch support, and the BL model as the lining support. Surrounding rock parameters are shown in Section 2.3. The support parameters are shown in Table 3.

To simulate the deformation and support failure characteristics of the surrounding rock in the creep section, the tunneling depth of the model was set as 120 m. The Mohr‒Coulomb model was used to calculate the instantaneous convergent deformation of the surrounding rock, and then the NCD model was used to calculate the creep deformation. The creep calculation time was 30 days. In the calculation process, four monitoring points were set on the monitoring section, which were located at the vault, tunnel shoulder, tunnel waist, and tunnel foot, respectively. After the calculation, the monitoring data were extracted to obtain the creep deformation curves of the tunnel section and compared with the on-site creep monitoring data within 30 days. The surrounding rock deformation is shown in Figure 13.

Figure 13a shows that without considering the creep deformation of the surrounding rock, the instantaneous convergent deformation of the surrounding rock of the vault is approximately 2.23 cm; the tunnel shoulder is approximately 1.21 cm; the tunnel waist is approximately 0.55 cm; and the tunnel foot is approximately 0.11 cm. When the creep deformation is considered, the converging deformation of the surrounding rock of the vault is approximately 29.76 cm on the 30th day, which exceeds the design value of the reserved deformation of the secondary lining. The deformation of the tunnel shoulder is approximately 18.58 cm; the tunnel waist is approximately 16.00 cm; and the tunnel foot is approximately 7.47 cm. In addition, the lining is damaged at the tunnel waist on the 20th day, and then the failure range gradually expands to both sides. According to Figure 13b, the instantaneous convergence of the surrounding rock of the vault measured on-site is approximately 2.27 cm (considering the constraint effect of the working face, the monitoring data when the working face exceeds the monitoring point by 10 m is taken); the tunnel shoulder is approximately 1.25 cm; the tunnel waist is approximately 0.53 cm; and the tunnel foot is approximately 0.11 cm. After 30 days of monitoring, the convergent deformation of the surrounding rock of the vault is approximately 29.83 cm; the tunnel shoulder is approximately 18.73 cm; the tunnel waist is approximately 16.48 cm; and the tunnel foot is approximately 7.66 cm. The first failure of the site lining also occurs at the arch waist, and it appears to expand to both sides. According to the calculation data and on-site monitoring data, the large deformation of the surrounding rock is mainly caused by creep, and the instantaneous deformation contributes little to the total deformation. By comparison, the instantaneous deformation and creep deformation of the numerical model are close to the actual monitoring value, and the maximum calculation error is 3.8%.

According to the calculation results, the NCD and BL models obtain basically the same results as on-site monitoring. Therefore, it can be summed up that the calculation method can effectively simulate the creep deformation characteristics of the surrounding rock and the failure law of the lining.

7 CONCLUSIONS

AUTHOR CONTRIBUTIONS

All authors contributed to the study, conception, and design. Material preparation and laboratory tests were performed by Dongxu Chen and Versaillot Pierre Darry. Data collection and theoretical analysis were performed by Laigui Wang and Chuang Sun. The first draft of the manuscript was written by Dongxu Chen, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflict of interest.