The Secondary Development and Application of the Improved Nishihara Creep Model in Soft Rock Tunnels

Abstract

Given the complexity and diversity of rock formations, existing constitutive models struggle to accurately portray their mechanical properties, leading to substantial discrepancies between numerical simulation outcomes and reality. This inadequacy fails to meet the demands of numerical analysis in practical engineering. This study first analyzes the physical and mechanical properties of thin-layered carbonaceous phyllite. Subsequently, an improved Nishihara rheological constitutive model is established based on these analyses. Utilizing the secondary development function offered by FLAC3D, the proposed model is further developed. The program’s correctness and reliability are confirmed through a numerical simulation using the triaxial creep test from existing research. Finally, the established constitutive model is applied in the numerical simulation of an actual soft rock tunnel engineering, obtaining results compared to real monitoring data. The results demonstrate that the improved Nishihara model is more effective at describing the creep deformation characteristics of soft rock. Moreover, the findings from this study can serve as a theoretical reference for predicting deformation in soft rock tunnel engineering.

1. Introduction

Nature is abundant with soft rock types such as mudstone, shale, sandstone, and phyllite. In practical engineering, these rock masses exhibit characteristics of softening and dilatancy as well as significant rheological properties. Engineering disasters, such as substantial non-linear deformation instability of the surrounding rock and the failure of support structures, often occur when constructing tunnels or underground structures in these soft rock strata [1,2]. This poses a risk to safety of their construction. Given the engineering construction risks in these soft rock strata, rheological research is of vital importance in geotechnical engineering. With a thorough understanding and study of the rheological properties of soft rocks like mudstone, shale, sandstone, and phyllite, their behavior can be better predicted and managed. This enables more effective planning and design in constructing tunnels or underground structures, thereby ensuring their long-term stability and safety [3,4,5].

Rheological properties of rocks relate to the time-dependent stress–strain relationship of the materials. The time effect in material deformation is referred to as the rheological phenomenon. Rheology, including creep, relaxation, and elastic aftereffect, is an important mechanical property of rock and soil. Numerous studies have focused on rock rheology research. For instance, Griggs [6] pioneered the introduction of rheology into rock mechanics, paving the way for rock rheology studies. Building upon this, Tan et al. [7] proposed a rheological theory and introduced it into rock mechanics. Li et al. [8] identified three usual stages of rock mass creep: deceleration, stable, and accelerated through uniaxial compression tests of four different types of soft rock. In fact, the primary focus of rock rheology research is to establish a suitable rheological model that can adequately reflect the rheological properties of rock mass [9,10]. This model can accurately simulate the stress–strain characteristics and the law of deformation displacement inherent in actual rock mass rheology in engineering contexts. Clearly, the rheological model serves a crucial role in predicting the mechanical deformation failure of rock masses [11,12,13].

The primary methodologies for constructing rheological constitutive models include empirical methods, component combination methods, and integral methods. Among them, the component combination method has garnered widespread application owing to its explicit interpretation and the intuitive understanding of its component parameters. Currently, many researchers have integrated the three fundamental components of elasticity, viscosity, and plasticity with specific tests, proposing a range of rheological constitutive models. These models are incorporated into FLAC3D via secondary development for numerical simulation purposes. Chen et al. [14] combined the damage element with the Poynting-Thomson creep model in FLAC3D to describe the accelerating stage of plastic creep. Wang et al. [15] proposed a nonlinear dashpot, replacing the linear dashpot in the Maxwell model and allowing the description of the entire process of salt rock creep. Chu et al. [16] examined the secondary development plug-in and the fundamental principles of FLAC3D 6.0 software operation, integrating these with the Nishihara rheological model. Using an actual tunnel project as the research object, Zou [17] integrated FLAC3D numerical simulations with the real project to verify the applicability of the Nishihara rheological model in engineering; the optimal timing of construction of the supporting structure is proposed. As scholars delve deeper into the study of rheological constitutive models, they have found that the Nishihara model can adequately reflect the elasticity, viscoelasticity, and viscoplasticity characteristics of rocks. Consequently, this model has become the most commonly used rheological model in rock engineering. However, with the increasingly complex characteristics of the surrounding rock exposed by practical engineering, the traditional Nishihara model occasionally fails to accurately depict the creep deformation of rock. Numerous efforts have been made in previous research to address the model’s lack of accuracy. Liu et al. [18] hypothesized that the viscosity coefficient in the creep model changes over time; based on this, they established an improved Nishihara rheological model connected to an enhanced dashpot. To compensate for the Nishihara model’s difficulties in accelerating creep, Yan et al. [19] connected a linear Hooke body, a fractional derivative Kelvin models, and a viscoplastic body composed of a nonlinear dashpot. Furthermore, some scholars employ damage theory to enhance the Nishihara model. Qi et al. [20] improved the model by incorporating a strain-triggered nonlinear dashpot into the Nishihara rheological model. Based on rheological experiments, Hou et al. [21] considered the impact of initial rock damage and proposed a nonlinear viscoplastic damage dashpot to replace the linear dashpot in the Nishihara model. Yang et al. [22] posited that the elastic modulus and viscosity coefficient of rock change with water content and time and also established a time-dependent creep damage model of rock. While these scholars’ research has compensated for some deficiencies in the Nishihara model, they have not further developed the improved Nishihara model or applied it in numerical simulation software. Therefore, the applicability of the improved Nishihara model is limited, making it difficult to promote its use.

Up to now, some research results in the secondary development of rock rheological constitutive models and numerical simulation software have been achieved. However, as construction reveals increasingly complex properties of the surrounding rock, existing constitutive models sometimes show significant deviations when analyzing the stability of tunnel surrounding rock. For example, when analyzing the stability of the surrounding rock in tunneling through the strata of thin-layered carbonaceous phyllite Rock, a type of soft rock formation, significant disparities are observed in the material parameters across different regions of the surrounding rock. Specifically, the cohesion within the plastic zone is markedly lower than that in the elastic zone [23,24]. This characteristic has not been sufficiently captured or reflected in the existing constitutive models. Current constitutive models predominantly focus on capturing the average properties of the surrounding rock. However, when dealing with rocks such as phyllite Rock, they may overlook the micro-structure and varied mechanical properties within different regions. Cohesion, as a critical mechanical parameter of rock, demonstrates significant variance between plastic and elastic zones and this discrepancy could have a substantial impact on the overall stability of the tunnel. Therefore, it is essential to establish corresponding constitutive models based on the characteristics of the surrounding rock in practical engineering situations [25,26,27].

In this study, thin-layered carbonaceous phyllite was chosen as the research subject and its physical and mechanical properties were initially analyzed. Based on this, a seven-element rheological model capable of describing the entire rheological process was established. This was then converted into a three-dimensional finite difference form and a custom constitutive model based on C++ was written using the secondary development function of FLAC3D software. The validity of this custom constitutive model was subsequently confirmed. Finally, the custom constitutive model was employed in the rheological analysis of tunnels. This study is significant in enhancing research on model secondary development.

2. Property Test of Thin-Layered Carbonaceous Phyllite

2.1. Longitudinal Wave Velocity Test

In fact, the longitudinal wave velocity can somewhat reflect the damage characteristics of the surrounding rock. A higher wave velocity implies lower internal damage and better integrity of the surrounding rock. Thus, the longitudinal wave velocity test can serve as an indirect measure of the surrounding rock’s degree of damage.

Wave velocity on the sidewall of the thin-layer carbonaceous phyllite tunnel was tested using the RSM-SY5 (T) non-metallic acoustic wave detector in conjunction with the ‘one-shot double-receive’ acoustic wave probe. During the test, the measurement hole was carefully maintained in a full water state, with water being used as the coupling medium. Three measurements were made for the same hole and the average value was subsequently taken. In this study, the longitudinal wave velocity of rock between the depths of 1~11 m was examined. The depth-longitudinal wave velocity diagram was then drawn, utilizing the longitudinal wave velocity values extracted from three measuring holes, as shown in Figure 1.

The test results in Figure 1 demonstrate that under construction disturbance, three sections can be classified in the longitudinal wave velocity of surrounding rock in thin-layered carbonaceous phyllite strata. A higher longitudinal wave velocity is exhibited by the rock mass outside 8 m of the hole depth with an average of 3.395 km/s, suggesting that the rock mass is relatively intact in this range, with lower levels of damage and minimal influence from construction disturbance. A decrease in the average longitudinal wave velocity of the surrounding rock from 3.395 km/s to 2.192 km/s is indicated within the 3.5 to 8.0 m range, indicating a significant influence from construction disturbance and a relatively higher degree of damage. A low average longitudinal wave velocity of only 1.552 km/s is observed within the 3.5 m hole depth range which can be attributed to severe damage to the rock mass in this range, leading to a low propagation velocity of acoustic waves. In conclusion, three levels of damage are primarily displayed by the thin-layered carbonaceous phyllite strata under construction disturbance.

2.2. Uniaxial Compression Test

Acquiring samples for testing is challenging due to the unique development of fractures and joints in thin-layered carbonaceous phyllite rock mass. Following tunnel excavation, rock samples larger in size were selected for further examination. The rock samples were then fashioned into cylindrical specimens with diameters and heights of 50 mm and 100 mm, respectively. The test rock sample is shown in Figure 2. Subsequently, uniaxial compression tests were conducted to simulate the deformation characteristics and principles of thin-layered carbonaceous phyllite. The test instrument is shown in Figure 3, the stress–strain curve can be plotted by using test data, as shown in Figure 4.

Figure 4 indicates that the deformation of thin-layered carbonaceous phyllite can be broken down into five stages, starting with Stage I: the compaction stage. This stage is characterized by low stress but substantial strain. This is primarily because thin-layered carbonaceous phyllite has numerous initial pores. The internal pores are closed under the action of vertical load. The next is Stage II: the elastic stage. This stage is marked by a linear relationship between stress and strain. Primarily, this is due to the compaction of initial pores within the test block. Vertical load induces elastic deformation of the internal micro-element, which can be considered an ideal elastomer. Stage III represents the elastic–plastic stage. During this stage, stress and strain deviate from their linear relationship with an increasing rate of strain change. This is because when external force reaches the yield limit, the internal micro-element begins to deteriorate, leading to plastic deformation. Stage IV denotes the softening stage. At this stage, the load-bearing capacity sharply declines, with an increased rate of softening. Primarily, this happens when stress reaches its peak, causing damage to the micro-element and subsequently reducing its load-bearing capacity. Stage V signifies the residual strength stage. Despite the destruction of the carbonaceous phyllite, it retains a certain load-bearing capacity.

The damage law of rock was examined from a macroscopic perspective using the longitudinal wave velocity test. Meanwhile, the deformation law of rock was investigated from a microscopic perspective using the uniaxial compression test. The degree of rock deformation increases proportionally with its damage. The results from the longitudinal wave velocity and uniaxial compression tests laid the groundwork for constructing the rheological constitutive model.

3. The Nonlinear Creep Equation of Thin-Layered Carbonaceous Phyllite

3.1. One-Dimensional Creep Constitutive Equation

The traditional Nishihara model is widely used in rheological analysis of rock mass which can describe the viscoelastic and viscoplastic creep stage of rock well. It consists of a Hooke body, Kelvin body, and viscoplastic body in series, as shown in Figure 5.

The one-dimensional creep constitutive equation of Nishihara model is as follows:

$$\left\{\begin{array}{l}\epsilon =\frac{\sigma }{E_{ie}}+\frac{\sigma }{E_{ve}}\left[1-\exp \left(-\frac{E_{ve}}{{\eta }_{ve}}t\right)\right],\hspace{0.33em}\hspace{0.33em}\sigma <{\sigma }_s,\\ \epsilon =\frac{\sigma }{E_{ie}}+\frac{\sigma }{E_{ve}}\left[1-\exp \left(-\frac{E_{ve}}{{\eta }_{ve}}t\right)\right]+\frac{\sigma -{\sigma }_s}{{\eta }_{vp}}t,\hspace{0.33em}\hspace{0.33em}\sigma \ge {\sigma }_s,\end{array}\right.$$

(1)

where $\sigma$ is the total stress of the model; $\epsilon$ is the total strain of the model; $E_{ie},E_{ve}$ are the elastic moduli of each elastomer in the element combination; ${\eta }_{ve},{\eta }_{vp}$ are the viscosity coefficients of each dashpot body in the element combination; and ${\sigma }_s$ is the long-term strength of rock.

However, the shortcomings of the traditional Nishihara model have also been pointed out by many scholars. For example, the model deformation develops rapidly at lower stress levels which leads to a shorter creep duration and makes it difficult to identify the creep mechanical behavior of soft rock. In order to describe the creep mechanical behavior of thin-layered carbonaceous phyllite, it needs to be improved.

Combined with the physical and mechanical properties of thin-layered carbonaceous phyllite test results in the second chapter and research results [23,28], after excavation disturbance the thin-layered carbonaceous phyllite tunnels have the following three features:

(1)

Under construction disturbance, the thin-layered carbonaceous phyllite stratum is mainly in three damage states;

(2)

Given the rock mass’s rheology, the partition mode of the surrounding rock is changed from the common ‘viscoelastic zone + viscoplastic zone ‘to‘ viscoelastic zone + viscoplastic zone + broken zone ‘, as shown in Figure 6. (At the same time, the analysis assumes the following: the tunnel’s excavation section is circular with an equivalent excavation radius of ‘a’; the surrounding rock is a homogeneous and isotropic viscoelastic–plastic body [29]; the tunnel is deeply buried in a hydrostatic stress field and the weight of the surrounding rock is neglected; and the initial support provides constant radial support resistance ‘pi’ and the support size remains constant.);

(3)

After tunnel excavation, the boundary stress of the viscoelastic zone and viscoplastic zone is the long-term strength of surrounding rock and the boundary stress of viscoplastic zone and broken zone is the temporary strength.

In view of them, if the stress after redistribution is greater than the temporary strength [28], the surrounding rock will form a broken zone, and the plasticity inside the rock will further increase. In practice, the surrounding rock of thin-layered carbonaceous phyllite formed a broken zone instead of a direct collapse after disturbance by construction. Therefore, this stage is defined as the broken creep stage. In order to describe the deformation of this stage, a modified viscoplastic model which is called the broken Bingham body was connected in series to improve the traditional Nishihara model and the stress threshold was temporary the strength; the improved model is named GNishihara model, as shown in Figure 7.

The one-dimensional creep constitutive equation of GNishihara model is as follows:

$$\left\{\begin{array}{l}\epsilon =\frac{\sigma }{E_{ie}}+\frac{\sigma }{E_{ve}}\left[1-\exp \left(-\frac{E_{ve}}{{\eta }_{ve}}t\right)\right],\hspace{0.33em}\hspace{0.33em}\sigma <{\sigma }_s,\\ \epsilon =\frac{\sigma }{E_{ie}}+\frac{\sigma }{E_{ve}}\left[1-\exp \left(-\frac{E_{ve}}{{\eta }_{ve}}t\right)\right]+\frac{\sigma -{\sigma }_s}{{\eta }_{vp}}t,\hspace{0.33em}\hspace{0.33em}{\sigma }_s\le \sigma <{\sigma }_b,\\ \epsilon =\frac{\sigma }{E_{ie}}+\frac{\sigma }{E_{ve}}\left[1-\exp \left(-\frac{E_{ve}}{{\eta }_{ve}}t\right)\right]+\frac{\sigma -{\sigma }_s}{{\eta }_{vp}}t+\frac{\sigma -{\sigma }_b}{{\eta }_b}t,\hspace{0.33em}\hspace{0.33em}{\sigma }_b<\sigma \end{array}\right.$$

(2)

where ${\eta }_b$ is the viscosity coefficient of the dashpot body in broken Bingham body and ${\sigma }_b$ is the temporary strength of rock.

3.2. D-P Strength Criterion Considering Softening Effect

In this paper, the Drucker–Prager criterion is selected as the yield criterion and expressed as follows:

$$\alpha I_1+\sqrt{J_2}-k=0$$

(3)

where I₁ is the first invariant of the stress tensor; J₂ is the second invariant of the stress tensor; α is a parameter related to cohesion; and k is a parameter related to cohesion and the internal friction angle, as shown in Equation (4).

$$\left\{\begin{array}{l}{I}_1={\sigma }_1+{\sigma }_2+{\sigma }_3\hfill \\ J_2=\frac{1}{6}[{({\sigma }_1-{\sigma }_2)}^2+{({\sigma }_2-{\sigma }_3)}^2+{({\sigma }_3-{\sigma }_1)}^2]\hfill \\ \alpha =\frac{2\sin \phi }{\sqrt{3}(3-\sin \phi )}\hfill \\ k=\frac{6c\cos \phi }{\sqrt{3}(3-\sin \phi )}\hfill \end{array}\right.$$

(4)

In rheological analysis of soft rock, such as thin-layered carbonaceous phyllite, it is vital not to overlook the significant impact of the intermediate principal stress. Hence, the intermediate principal stress coefficient was introduced, commonly used in coal mine roadway analysis to simplify the process, and expressed as follows:

$$n=\frac{{\sigma }_2-{\sigma }_3}{{\sigma }_1-{\sigma }_3}$$

(5)

The intermediate principal stress coefficient n is brought into Equation (3) and simplified to obtain Equation (6).

$${\sigma }_1+M{\sigma }_3-N=0$$

(6)

where

$$\left\{\begin{array}{c}M=\frac{2\alpha -n\alpha -m}{\alpha +n\alpha +m}\hfill \\ \begin{array}{l}N=\frac{k}{\alpha +n\alpha +m}\\ m=\sqrt{\frac{1}{3}(n^2-n+1)}\end{array}\end{array}\right.$$

(7)

According to Yu et al. [23] and Brady et al. [24], upon post-destruction of the rock there is a noticeable change in cohesion whereas the internal friction angle remains largely unaltered. Combined with a large number of experimental studies on carbonaceous phyllite and the engineering properties of on-site thin-layered carbonaceous phyllite, it assumes that the cohesion of thin-layered carbonaceous phyllite decreases with the increase in strain in the plastic stage and decreases to 0 after entering the crushing creep stage, as shown in Figure 8.

Figure 8 demonstrates the correlation between the damage inflicted by cohesion and strain. Prior to the viscoplastic stage, the surrounding rock remains undamaged and unsoftened. Upon entering the viscoplastic stage, the surrounding rock experiences damage and a gradual decrease in cohesion. When it transitions into the broken creep stage, the cohesion reduces to 0. Consequently, the damage factor ‘D’ of cohesion can be represented using the subsequent equation.

$$D=\frac{\epsilon -{\epsilon }^p}{{\epsilon }^b-{\epsilon }^p}$$

(8)

where ${\epsilon }^p$ is the strain at the boundary of the viscoplastic zone and ${\epsilon }^b$ is the strain at the boundary of the broken zone.

The softening factor S is introduced to simplify the calculation equation. The relationship between the softening factor and the damage factor is shown by the following equation:

$$S=1-D=1-\frac{\epsilon -{\epsilon }^p}{{\epsilon }^b-{\epsilon }^p}$$

(9)

The strain generated in different zones of surrounding rock is related to the location and the softening factor is further simplified and expressed as follows:

$$S=\frac{R_{p}r-R_p{R}_b}{\left(R_p-R_b\right)r}$$

(10)

where the value range of r is $R_p\le r\le R_b$; $R_p$ is the radius of the plastic zone; and $R_b$ is the radius of the broken zone.

The softening factor is a variable formulated to delineate the changes in cohesion. This variable, when multiplied with cohesion, can effectively portray the cohesion changes. In the D-P strength criterion, only parameters ‘k’ and ‘N’ have a relationship with cohesion ‘c’. Thus, the D-P strength criterion accounting for the softening factor can be established as follows:

$${\sigma }_1+M{\sigma }_3-NS=0$$

(11)

Equation (11) represents the D-P strength criterion, considering both the softening effect and intermediate principal stress. Within the viscoelastic zone, the softening factor equals 1 and Equation (6) exhibits the strength criterion expression. In the broken zone, the softening factor reduces to 0 and the strength criterion expression is shown in the succeeding equation:

$${\sigma }_1+M{\sigma }_3=0$$

(12)

4. Implementing the Nonlinear Creep Model in FLAC3D Software

4.1. Operating Principles of the FLAC3D Constitutive Model

FLAC3D6.0 (Fast Lagrange Analysis of Continua in 3D), developed by ITASCA, USA, is a finite difference program explicitly designed for the three-dimensional mechanical analysis of civil engineering issues. It effectively simulates the mechanical behavior of material damage and plastic flow. FLAC3D uses an explicit differencing scheme to calculate the stress state at time t + Δt which is based on the stress state at time t and the incremental strains in the time interval Δt. The specific computational procedure is illustrated in Figure 9. FLAC3D implements constitutive models as dynamic link library (.dll) files which primarily return updated stress tensors based on incremental strains. The main purpose of custom constitutive models is to compute the provided incremental strains and derive the corresponding updated stresses.

During the model’s secondary development, Microsoft Visual C++ should be used to compile the model into a.dll file which should be placed in the FLAC3D software’s installation directory. The custom constitutive model can be invoked in the main program using the “model load” command by specifying the name of the dynamically linked library file. This study used the FLAC3D development platform to develop a nonlinear creep model; the model’s accuracy was subsequently verified.

4.2. Finite Difference Form of the GNonlinear Creep Model

For the secondary development of the creep constitutive model introduced in this paper, the associated constitutive relation needs to be rewritten into a finite difference form for a three-dimensional state. This involves rewriting the increments of stress and strain during the solution process into a different form in terms of the creep time.

According to the series–parallel combination relationship of the components, the total strain in Figure 5 model can be expressed as the sum of four parts:

$${\epsilon }_{ij}={\epsilon }_{ij}^{ie}+{\epsilon }_{ij}^{ve}+{\epsilon }_{ij}^{vp}+{\epsilon }_{ij}^{vb}$$

(13)

Equation (13) is rewritten into the form of deviatoric strain rate, that is,

$${\dot{e}}_{ij}={\dot{e}}_{ij}^{ie}+{\dot{e}}_{ij}^{ve}+{\dot{e}}_{ij}^{vp}+{\dot{e}}_{ij}^{vb}$$

(14)

where ${\dot{e}}_{ij}$ is the total deviatoric strain rate; ${\dot{e}}_{ij}^{ie}$ is the deviatoric strain rate of the Hooke body; ${\dot{e}}_{ij}^{ve}$ is the deviatoric strain rate of the Kelvin body; ${\dot{e}}_{ij}^{vp}$ is the deviatoric strain rate of the Bingham body; and ${\dot{e}}_{ij}^{vb}$ is the deviatoric strain rate of breaking Bingham body.

Equation (14) should be written in increment form, that is,

$$\Delta e_{ij}=\Delta e_{ij}^{ie}+\Delta e_{ij}^{ve}+\Delta e_{ij}^{vp}+\Delta e_{ij}^{vb}$$

(15)

For the Hooke body in the first part, the constitutive relation in a one-dimensional state is

$$\sigma =E_{ie}{\epsilon }_{ie}$$

(16)

Equation (16) is extended to the three-dimensional case and the following can be obtained:

$$S_{ij}^{ie}=2G_1{\epsilon }_{ij}^{ie}$$

(17)

where $S_{ij}$ is deviatoric stress and $G_1$ is the shear modulus of the Hooke body.

Expressing the above equation as an incremental form, the following is obtained:

$$\Delta e_{ij}^{ie}=\frac{\Delta S_{ij}}{2G_1}$$

(18)

where $\Delta e_{ij}^{ie}=e_{ij}^{ie,N}-e_{ij}^{ie,O}$ $\Delta S_{ij}=e_{ij}^{ie,N}-e_{ij}^{ie,O}$ and the superscript N represents a new value in a time step increment while the superscript O represents an old value.

For the Kelvin body, the constitutive relation in a one-dimensional state is

$$\sigma =E_{ve}{\epsilon }_{ve}+{\eta }_{ve}{\dot{\epsilon }}_{ve}$$

(19)

Referring to Equation (17), the deviatoric stress–strain relationship of the Kelvin body can be obtained as follows:

$$S_{ij}=2G_2{\epsilon }_{ij}^{ve}+2{\eta }_{ve}{\dot{\epsilon }}_{ij}^{ve}$$

(20)

where $G_2$ is the shear modulus of the Kelvin body.

Equation (20) is written in incremental form:

$$\Delta e_{ij}^{ve}=e_{ij}^{ve,N}-e_{ij}^{ve,O}=\frac{1}{A}\left[B{e}_{ij}^{ve,O}+\frac{\Delta t}{4{\eta }_2}\left(S_{ij}^N+S_{ij}^O\right)\right]-e_{ij}^{ve,O}$$

(21)

where $A=1+\frac{G_2}{2{\eta }_{ve}}\Delta t$ and $B=1-\frac{G_2}{2{\eta }_{ve}}\Delta t$.

For the Bingham body and the broken Bingham body, the constitutive relation in the one-dimensional state is

$$\left\{\begin{array}{l}{\epsilon }_{vp}=\frac{\sigma -{\sigma }_s}{{\eta }_{vp}}t,\hspace{0.33em}\hspace{0.33em}\sigma -{\sigma }_s>0\\ {\epsilon }_{vb}=\frac{\sigma -{\sigma }_b}{{\eta }_{vb}}t,\hspace{0.33em}\hspace{0.33em}\sigma -{\sigma }_b>0\end{array}\right.$$

(22)

The Equation (22) is rewritten into a three-dimensional form, that is,

$$\left\{\begin{array}{l}{\epsilon }_{ij}^{vp}=\frac{\langle F_a\rangle }{2{\eta }_{vp}}\frac{\partial Q_a}{\partial {\sigma }_{ij}}t\\ {\epsilon }_{ij}^{vb}=\frac{\langle F_b\rangle }{2{\eta }_b}\frac{\partial Q_b}{\partial {\sigma }_{ij}}t\end{array}\right.$$

(23)

where $\langle F_a\rangle$ and $\langle F_b\rangle$ represent the switching function in two states, respectively, as shown in Equations (24) and (25).

$$\langle F_a\rangle \left\{\begin{array}{l}0,\hspace{0.33em}\left(F_a\le 0\right)\\ F_a,\hspace{0.33em}\left(F_a>0\right)\end{array}\right.,\hspace{0.33em}{F}_a=S_{ij}-S^{\infty }$$

(24)

$$\langle F_b\rangle \left\{\begin{array}{l}0,\hspace{0.33em}\left(F_b\le 0\right)\\ F_b,\hspace{0.33em}\left(F_b>0\right)\end{array}\right.,\hspace{0.33em}{F}_b=S_{ij}-S^b$$

(25)

where $S^{\infty }$ and $S^b$ are the long-term strength and temporary strength, respectively; $Q_a$ and $Q_b$ are the plastic potential functions in two states, respectively; and the associated flow criterion is adopted, allowing them to be yield criteria in two states, respectively, as shown in Equations (11) and (12).

Equation (23) can be expressed as follows:

$$\left\{\begin{array}{l}{\epsilon }_{ij}^{vp}=\frac{S_{ij}-S^{\infty }}{2{\eta }_{vp}}\left(1+M\right)t\\ {\epsilon }_{ij}^{vb}=\frac{S_{ij}-S^b}{2{\eta }_b}\left(1+M\right)t\end{array}\right.$$

(26)

Equation (26) is written in incremental form:

$$\left\{\begin{array}{l}\Delta {\epsilon }_{ij}^{vp}=\frac{\langle F_a\rangle }{2{\eta }_{vp}}\frac{\partial Q_a}{\partial {\sigma }_{ij}}\Delta t\\ \Delta {\epsilon }_{ij}^{vb}=\frac{\langle F_b\rangle }{2{\eta }_b}\frac{\partial Q_b}{\partial {\sigma }_{ij}}\Delta t\end{array}\right.$$

(27)

According to (18), (21), and (27), the three-dimensional difference form of the GNishihara model can be expressed as follows:

$$e_{ij}^N-e_{ij}^O=\frac{S_{ij}^N-S_{ij}^O}{2G_h}+\frac{1}{A}\left[B{e}_{ij}^{ve,O}+\frac{\Delta t}{4{\eta }_{ve}}\left(S_{ij}^N+S_{ij}^O\right)\right]-e_{ij}^{ve,O}+\frac{F_a}{2{\eta }_{vp}}\frac{\partial Q_a}{\partial {\sigma }_{ij}}\Delta t+\frac{F_b}{2{\eta }_b}\frac{\partial Q_b}{\partial {\sigma }_{ij}}\Delta t$$

(28)

Equation (28) is transformed into the form of deviatoric stress

$$S_{ij}^N=\frac{1}{g}\left\{e_{ij}^N-e_{ij}^O+f{S}_{ij}^O-\frac{B{e}_{ij}^{ve,O}}{A}+e_{ij}^{ve,O}-\frac{F_a}{2{\eta }_{vp}}\frac{\partial Q_a}{\partial {\sigma }_{ij}}\Delta t-\frac{F_b}{2{\eta }_b}\frac{\partial Q_b}{\partial {\sigma }_{ij}}\Delta t\right\}$$

(29)

where $g=\frac{1}{2G_h}+\frac{\Delta t}{4A{\eta }_{ve}}$ and $h=\frac{1}{2G_h}-\frac{\Delta t}{4A{\eta }_{ve}}$.

In fact, the spherical stress tensor has an effect on the surrounding rock before the plastic change which needs to be considered in the numerical simulation process. It only affects the elastomer and Kelvin body. Therefore, the spherical stress can be expressed as

$$S_m^N=S_m^O+K_1\Delta {\epsilon }_m^{ie}+K_2\Delta {\epsilon }_m^{ve}$$

(30)

where $K_1$ and $K_2$ are the bulk modulus of the Hooke and Kelvin bodies, respectively.

The differential form of the spherical stress of the Hooke and Kelvin bodies are as follows:

$$\left\{\begin{array}{l}{S}_m^N=S_m^O+3K_1\left({\epsilon }_m^{ie,N}-{\epsilon }_m^{ie,O}\right)\\ e_m^{ve,N}=\frac{1}{C}\left[D{e}_m^{ve,O}+\frac{\Delta t}{6{\eta }_2}\left(S_m^N+S_m^O\right)\right]\end{array}\right.$$

(31)

where $C=1+\frac{K_2}{2{\eta }_2}\Delta t$ and $D=1-\frac{K_2}{2{\eta }_2}\Delta t$.

At this point, the constitutive model relationship of the GNishihara model in the three-dimensional difference form is fully displayed and the dynamic link library can be compiled accordingly.

4.3. Program Verification

To affirm the precision of the GNishihara creep model’s secondary development used in this study, this section validates the simulation outcomes by creating a numerical model that aligns with the creep test per the research findings [30]. The simulation model is represented as a 50 mm × 100 mm cylinder partitioned into 1152 zones and 1350 nodes. In this setup, uniform vertical stress is imposed on the model’s top with constant circumferential stress on both sides (simulating confining pressure). Additionally, a displacement restriction is enacted at the bottom in the vertical direction. Figure 10 illustrates the specific numerical model and its corresponding boundary conditions. As illustrated in

Table 1. Calculating parameters obtained by inversion.

$\mathit{c}$ /kPa	$\mathit{\phi }$ /°	${\mathit{G}}_1$ /MPa	${\mathit{G}}_2$ /GPa	${\mathit{\eta }}_2$ /GPa·h−1	$\mathit{\eta }1$ /GPa·h−1	$\mathit{\eta }\mathit{b}$	ha	hb	$\mathit{n}$	${\mathit{p}}_0$ /MPa
69.4	8.274	841.8	1.085	1.531	3.790	4.469	1.26	1.51	0.688	7.765

, the rock material parameters in the numerical simulation are derived using the inversion method. In the simulation, an axial load is induced through a steady confining pressure and incremental loading. The confining pressures are assigned values of 3 MPa and 6 MPa. Each creep stage lasts 48 h with the axial load deviator stress for each stage set at 4.5 MPa, 6 MPa, 7.5 MPa, 9 MPa, 10.5 MPa, and 12 MPa, respectively.

The numerical simulation outcomes are compared with the laboratory experimental results [30]. Figure 11 and Figure 12 illustrate the variation of the specimen’s axial strain under diverse axial loading stresses over time. As shown in Figure 11 and Figure 12, under axial loads of 4.5 MPa, 6 MPa, 7.5 MPa, and 9 MPa, the GNishihara creep model’s calculations align closely with laboratory experiments, with the axial strain value discrepancy being less than 10%. Upon loading to 10.5 MPa, laboratory experiment results exceed GNishihara creep model calculations by approximately 15%. This discrepancy can be attributed to the reduction in actual parameters as the laboratory experiment rock sample nears failure under large loading. While the GNishihara creep model accounts for cohesion decay, it disregards the decay of the internal friction angle, thus justifying these results. To conclude, the simulation curves generated by the nonlinear viscoelastic–plastic creep model align well with the experimental results under two distinct confining pressures. This suggests that the nonlinear creep model, developed using FLAC3D, is capable of accurately representing the entire creep process of thin-layered carbonaceous phyllite. This validates the logic and efficiency of the secondary development strategy of the nonlinear viscoelastic–plastic creep model, as proposed in this study.

5. Engineering Verification

Through utilizing the secondary development of the GNishihara viscoelastic–plastic model in FLAC3D, the rheological properties of thin-layered tunnel can be analyzed. It also enables the prediction of deformation trends and long-term stability which can be compared with the monitoring data of surrounding rock displacement in practical construction.

5.1. Engineering Overview

Located in Ankang City, Shaanxi Province, the Xiejiapo Tunnel is a crucial connector between Langao County and Ankang City (Figure 13a). According to the engineering survey report, the lithology of the tunnel site area primarily comprises weakly weathered and medium weathered carbonaceous phyllite. In the actual construction, it is found that the carbonaceous phyllite layer is exceptionally thin over approximately 600 m of the tunnel and can be classified as thin-layered carbonaceous phyllite. The engineering geological profile is shown in Figure 13b. The support parameters for this part of the tunnel are selected based on the mechanical properties of the medium weathered carbonaceous phyllite without consideration of the substandard rock properties of the thin-layered carbonaceous phyllite. This has resulted in substantial deformation and significant tunnel intrusion during construction, as depicted in Figure 13c,d. Such issues pose a major threat to the tunnel’s safe operation. Therefore, it remains a focal point of concern.

5.2. Construction and Calculation of the Model

This paper presents a numerical simulation analysis of the IV grade surrounding rock section ZK17 + 267 ~ ZK17 + 306 in the left line of the Xiejiapo tunnel. The primary lithology is weakly weathered phyllite. Figure 14 illustrates the model grid division of the tunnel. In the figure, dark blue represents the surrounding rock, red the upper step, yellow the middle step, and green the lower step. The model consists of 23,739 nodes and 22,280 units. The X axis represents the width direction of the tunnel, the Z axis the height, and the Y axis the excavation direction. The method of simulating the surrounding rock excavation aligns with that of actual projects. The three-bench excavation method is adopted. Each step of the excavation extends 1.5 m and is divided into the upper, middle, and lower steps. Each step’s excavation length is 6 m. Figure 15 demonstrates the simulation of the specific model’s excavation process. To more accurately simulate the actual stress state of the surrounding rock, constraints are applied to the front, back, left, right, and bottom surfaces of the model. The top surface of the model is unconstrained and a uniform load is applied to mimic the pressure at the actual depth of the surrounding rock in the simulated section. The lateral pressure of the surrounding rock is simulated through uniform loading. The tunnel rheology lasts for 35 days. The secondary stress state of the tunnel is calculated first during the rheological calculation process. Then, the rheological calculations are initiated.

The constitutive model compiled above is used to assign the parameters of the surrounding rock. The cable element simulates the anchor rod of the initial support and the liner element simulates the steel arch concrete. Parameters for the initial support structure are designated in accordance with the properties of each component in the actual project. The material parameters for the cable element and the liner element are presented in

Table 2. Support parameter.

Support Structure	E/GPa	v	Cross-Sectional-Area/m2	Density/Kg·m−3	Thickness/m
Cable element	210	0.3	6.22 × 10−3	7.6	/
Liner enelemt	28	0.17	/	2.45	0.3

. During the computational process, no contact elements were added between the support element and the rock. In practical engineering, the secondary lining of the tunnel is usually constructed after the deformation of the surrounding rock is relatively stable and the strength and stiffness are large; so, only the deformation of the surrounding rock before the construction of the secondary lining is simulated.

Before simulating an excavation, monitoring points for the surrounding rock displacement are established at the vault and arch waist of the tunnel model. The displacement of the tunnel vault (point A) is primarily monitored along the z-axis direction. The displacement of the tunnel arch waist (points B and C) is primarily monitored along the x-axis direction and the displacement of the tunnel side wall (points D and E) primarily monitors its displacement along the x-axis direction. Figure 16 shows the layout of the monitoring points.

5.3. Comparative Analysis of Surrounding Rock Displacement

The rheological calculations for the tunnel reveal the distribution of the displacement and stress fields. In this study, the displacement field is used for comparison with the actual surrounding rock displacement. Additionally, the monitoring data of a typical section are used to represent the actual surrounding rock displacement.

The displacement diagrams in the Z and X directions of the model, after 35 days of creep, are shown in Figure 17 and Figure 18. The figures show that the maximum vertical displacement of surrounding rock occurs at the tunnel’s vault and inverted arch, manifesting as vault subsidence and arch bottom uplift in practical engineering. Moreover, the maximum horizontal displacement is observed at the tunnel’s arch waist and sidewall, indicating peripheral convergence in a practical engineering context.

The data recorded by the model’s monitoring points are compared with actual surrounding rock displacement data, as illustrated in Figure 19. Figure 19 shows that the numerical simulation trend aligns with that of the actual surrounding rock position. During the initial phase of tunnel excavation, the surrounding rock displacement rate is high. As time progresses, this rate gradually reduces to a minimal constant. However, the numerical simulation predicts a longer stabilization time for the surrounding rock displacement, introducing an error of about 7–10 days compared to the actual stabilization time. Furthermore, the numerical simulation’s predicted displacement of the surrounding rock closely matches the actual displacement. The final displacement predicted is larger but with a smaller error. The numerical simulation predicts a final displacement of 499.3 mm at vault point A with horizontal convergence values of 351.6 mm at spandrel BC and 397.8 mm at the side wall DE. The respective errors between these values and the actual monitored surrounding rock displacement are 3.14%, 3.75%, and 4.02%.

The comparative analysis indicates that the established GNishihara model can accurately predict the displacement of surrounding rock with minor calculation errors. Furthermore, the model’s prediction trend aligns with the actual displacement trend.

6. Discussion

The improvement to the Nishihara model proposed in this paper primarily targets the analysis of surrounding rock stability that exhibits significant rheological characteristics or pronounced softening features. It demonstrates good applicability in engineering contexts where surrounding rocks exhibit evident deformation such as in high geostress large-span tunnels, loess tunnels, expansive soil tunnels, and other soft rock tunnels. However, further discussion is still required in areas such as urban rail transit shallow buried sections or projects where strict control over ground surface deformation is essential.

7. Conclusions

Physical and mechanical properties of thin-layered carbonaceous phyllite are analyzed based on field and laboratory test data. Based on these analyses, improvements have been made to the traditional Nishihara model. Following this, the derived nonlinear viscoelastic–plastic creep model is further developed and its calculation program is validated. The following conclusions are drawn:

• Field test results indicate that the rock surrounding thin-layered carbonaceous phyllite can be divided into three damage states post-excavation. Uniaxial compression test results demonstrate that thin-layered carbonaceous phyllite softens significantly upon reaching yield strength yet maintains a certain load capacity;
• Consideration of the physical and mechanical properties of thin-layered carbonaceous phyllite has led to improvements in the traditional Nishihara model. By introducing a softening factor, the D-P yield criterion is simplified and modified;
• A detailed three-dimensional finite difference expression of the constitutive model is derived based on the finite difference theory. Utilizing the FLAC3D software’s secondary development platform, the Nishihara constitutive model is further developed and then validated via the simulation of a triaxial compression creep test. Simulation results align with the actual creep test results, thus validating the logic and efficacy of the model’s further developed calculation program;
• Engineering examples simulated by the model in this study are compared with actual monitoring results. The results reveal that the displacement error at vault point A is 3.14% with convergence errors of 3.75% at spandrel BC and 4.02% at sidewall DE. Moreover, the trend of numerical simulation results aligns with the actual situation.