A viscoelastic–plastic constitutive model incorporating creep initiation stress for surrounding rock in deep roadways and its practical application

With the increasing depth of engineering in rock masses, the issue of significant rheological deformation becomes notably prominent. To accurately characterize the deformation and failure behavior of the surrounding rock near the goaf in deep small coal pillars, we have developed a viscoelastic–plastic constitutive model that takes into account the initial stress associated with creep. This model builds upon an existing viscoelastic–plastic constitutive framework and is implemented numerically using C++ in the FLAC3D. The Oedometer test method is used to check the creep components and the calculated parameter (η/G) ensure convergence of numerical operations. We conducted both pre‐improvement and post‐improvement tests of the constitutive model in a simple numerical analysis of roadways. The enhanced model demonstrates improved accuracy in deviational stress calculations and proves suitable for modeling deep, soft roadways. Furthermore, we applied the improved model to simulate actual working conditions at Zhaozhuang Mine, where mining operations near the working face induce a fracture zone of approximately 3 m in the small coal pillar located at the goaf's edge. The new model closely aligns with the observed results, further validating its suitability for tunneling along the goaf. Tunneling and mining activities result in significant deformation at various locations, including the shoulder socket of the working face, the top angle of the coal pillar, and the bottom angle of the floor. The crushing zone at the working face extends up to 4 m. Rheological effects transform the triangular elastic core, from its original high‐pressure elastic state to a fracture and plastic zone in the deep coal. Stress in deep coal is transmitted along the triangular elastic core. Numerical simulation results closely match the excavation profile of the tunneling, providing valuable insights into the impact of mining activities. This study can provide reference for deformation of similar geological conditions and engineering situation.


| INTRODUCTION
Coal, as the fundamental and primary energy source in China, represents a valuable nonrenewable fossil energy resource. 1,2To enhance coal extraction rates and achieve safe and efficient mining practices, our country has promoted the technique of tunneling along the goaf with small coal pillars. 1,2This approach leverages the stressreduced zone at the goaf's edge to position the mining roadway within a depressurized zone. 3,4][7] As shallow coal resources gradually deplete, most coal mines are transitioning to or are already operating in the deep mining phase.In the high ground stress environment of deep mining, the plastic failure extent of coal and rock masses deepens, and the influence of rheological damage becomes pronounced.This results in significant plastic flow related to stress adjustment and viscous flow associated with time in the surrounding rock of the roadway, leading to extensive deformations. 8,9hese deformations have adverse implications for coal mine production safety and considerably hinder the application and adoption of the practice of leaving small coal pillars along the gob. 8,9o safely and efficiently implement this approach in deep mining, increase the recovery rate of deep coal resources, and elucidate the mechanisms underlying viscoplastic deterioration and significant deformations in the surrounding rock are vital prerequisites.However, existing fundamental theories of rock mechanics and support theories do not adequately guide deep mining practices. 10Consequently, there is an urgent need for indepth research into the mechanisms governing viscoplastic deterioration and substantial deformations in deep goaf excavation with small coal pillars.This research will serve as the basis for deformation prediction and support design for such roadway constructions.
Tao et al. conducted a study on the rheological characteristics of mylonite in Zhaogezhuang Coal Mine.When compared with the Burgess model, they found that the Nishihara model exhibited a better fit to the creep curve, with a higher correlation coefficient, and provided a more comprehensive description of mylonite characteristics. 11uan et al. integrated the newly proposed plastic element into the Burgess model in series, creating an improved Burgers creep model.This improved model better captures the characteristics of rock decay creep, stable creep, and accelerated creep. 12ia et al. provided a more accurate description of the creep characteristics of clayey rock and introduced a nonlinear elastic-viscoplastic constitutive model based on a new yield criterion and a modified Mohr-Coulomb model. 13iong et al. developed a composite viscoelastic rheological model suitable for hard rock by connecting a plastic element based on a strain softening model in series with a six-element viscoelastic rheological model. 14i et al. proposed a rheological constitutive model of the yield surface using the triaxial shear creep test method.This model effectively describes the viscoelasticplastic creep characteristics of soft rocks. 15alan et al. analyzed the true triaxial creep test results of mudstone and established a nonlinear empirical power function creep model that closely matches the experimental curve and field observations. 16onsidering the changing strength parameters within rock masses due to their state, Cristescu and Gioda developed a creep damage model with variable parameters.This model addresses the limitations of the inherent FLAC constitutive model. 17ellet et al. ombined damage mechanics theory with rock mechanics to establish both a damage mechanics model and a viscoelastic-viscoplastic rheological model for jointed rock masses. 18uring the study of creep characteristics in silty mudstone, granite, and mica quartz gneiss, Wang et al. noted that significant rheological deformation of rock samples only occurred when stress reached a certain level.Minimal rheological deformation was observed at lower stress levels. 19un and Pan identified a "lower rheological limit" in rock mass during rheological testing.This means that a certain minimum stress level must be exceeded for rheological deformation to occur within the rock mass. 20an et al. research confirmed the existence of an initial creep stress threshold, which linearly increases with increasing confining pressure during oil mudstone creep tests. 21heng et al. analyzed uniaxial, conventional triaxial and true triaxial creep test results, the crack initiation stress is taken as the starting point of creep and the timedependent characteristics of rock under untested true triaxial stress levels were predicted. 22ou et al. based on multi-loading creep tests of sandstone with different initial damage levels, the experimental results show that initial damage has a clear effect on the creep behavior of rock. 23hrough three-point bending tests, Wang ZY et al. derived a more accurate method of calculating toughness based on modified fracture mechanics theory, which can better explain the fracture characteristics of quasi-brittle rocks. 24an et al.'s study research explains the coupling effect of creep behavior and internal pressure in backfill stope and calculates the theoretical solution considering the time effect to verify the correctness. 25o address limitations in numerical simulation software models, many scholars frequently adopt combined element models to simulate measured stress-strain-time curves.They introduce damage mechanics and fracture mechanics to reflect the constitutive model of damage accumulation and crack propagation during the aging deformation of rock mass.
While scholars have conducted in-depth studies on the rheological properties of rock through testing and theory, improving our understanding of rheological effects, [26][27][28][29][30][31][32][33][34][35][36][37][38][39] existing results have often oversimplified the true creep behavior of rock.They often neglect the influence of the initial creep stress on the creep characteristics and do not develop models suitable for numerical calculations.This oversight may lead to the prediction of creep deformation even at lower stress levels during numerical simulations, deviating from realworld conditions.
Kang and Yi developed a viscoelastic-plastic constitutive model (BVS) based on the deformation and failure characteristics of deep soft rock roadway surrounding rock and its associated mechanical properties. 40][43] However, the existing model inadequately addresses rocks that do not exhibit creep deformation under initial creep stress, a natural occurrence that aligns with the inherent objective law. 20,21Building upon this foundation, the author enhanced the BVS model by incorporating the initial creep stress of rock to describe scenarios in which such stress does not exceed the threshold for creep deformation.This modification resulted in a new constitutive model, complete with a numerical implementation compatible with FLAC 3D , thereby producing more realistic and comprehensive simulation results.
Subsequently, leveraging rheology theory and field measurements, the model was applied to the 13,153 intake roadway of Zhaozhuang Mine in Shanxi Province.This application yielded insights into the behavior of surrounding rock during excavation, rheology, and stoping stages.Following a field survey, the suitability of the NVS constitutive model was confirmed.
Based on the phenomenon of initial creep stress found in rock creep experiments, 20,21,44,45 we takes it into account for the first time to improve the BVS model, and the engineering verification is carried out.So that the stress calculation in numerical simulation is more like a natural occurrence.

| COMPOSITION OF THE NVS MODEL AFTER ENHANCEMENTS TO THE BVS MODEL
The BVS model effectively captures transient strain softening, plastic dilatancy characteristics, and various rheological effects, including three stages of creep, differentiation between stable and unstable creep, rheological damage, and dilatancy. 40To address situations where rocks exhibit no rheological deformation under low stress, 20,21,44,45 the model incorporates a mechanism to discriminate the initial creep stress of the rock, as illustrated in Figure 1.
The model is composed of the enhanced Nishihara body, improved viscoplastic and strain softening components from the BVS model.The NVS model is briefly described below.

| Enhanced Nishihara body considering the theoretical limitations of creep termination curve
When a certain stress level is applied to a rock, an instantaneous elastic strain is generated immediately.Chen argued that if the stress remains below a certain level, the rock sample will not undergo creep deformation. 44Experiments conducted by Xia et al. revealed that the initial creep stress for red sandstone is 12% of the uniaxial compressive strength, indicating that low-stress levels do not induce creep deformation. 45n the BVS model, the acceleration creep is initiated by the damage resulting from the accumulation of deceleration and constant velocity creep in the three stages of creep.The deceleration creep is described by the Kelvin body, while the constant velocity creep is described by the limited Newtonian body (referred to as the limited Newtonian body, constrained by the creep termination curve).Essentially, the initial creep stress serves as the switch to activate the limited Newtonian body and the Kelvin body.Consequently, when the stress level is lower than the initial creep stress, no creep deformation occurs.To achieve this, a series and parallel connection of certain basic elements is employed.Specifically, a creep initial stress switch (lacking plastic properties with a switch value of 12% of the uniaxial compressive strength) is connected in parallel with the improved Newton body and Kelvin body.In addition, the Hooke body, improved Bingham body, and Kelvin body are connected in series to create an enhanced Nishihara body, as depicted in Figure 2.

| Improved Bingham body
The Bingham body consists of a Hooke body and an improved viscoplastic body (friction sheet without plastic properties) connected in parallel.When σ σ < ct , the element undergoes purely elastic changes, the constitutive equation is equivalent to Hooke's body, which is not detailed here.When σ σ > ct , the improved viscoplastic body exhibits mechanical response and undergoes creep deformation.Its one-dimensional (1D) constitutive equation is as follows: ( )

| Enhanced Nishihara body
The Nishihara body is composed of the Hooke body, enhanced Kelvin body, and enhanced viscoplastic body in series.Unlike the traditional Nishihara body, when σ σ < ct , the Kelvin body does not exhibit a response, the enhanced Nishihara body undergoes only elastic changes.When σ σ > ct , as only a friction plate switch is introduced, the stress in the model must merely exceed the resistance of the friction plate.Therefore, there is no need for a detailed derivation of the creep equation, 46 but it is presented directly here.The creep equation is as follows:

| Operation of creep termination curve
In the BVS model, to simplify and expedite the calculation of the constitutive equation, a linear creep termination curve, which does not consider the creep starting stress, is employed.Therefore, when considering the creep starting stress, the slope of the creep termination curve in this model needs modification based on the BVS model.Consequently, the maximum creep value of the rock mass under a certain stress level also requires adjustment.In addition, improvements were made to the Newtonian body in the Bingham body.
The creep deformation process of rock mass under various viscoelastoplastic stress states primarily manifests as differences in the viscoplastic stage.Accurate expression of the change in viscoplastic creep rate in different time intervals can yield the fundamental equation for the creep deformation process. 47,48With this understanding, determining the boundary points for different deformation stages on the entire creep process curve of rock becomes possible using the total stress-strain curve of rock.Corresponding constitutive relationships for these deformation stages can be established.Simultaneously, after accounting for the creep starting stress at a specific stress level and based on the strain space, 49 the quantitative description of the rock's creep deformation process can be formulated.Figure 3 illustrates this, with σ c representing uniaxial compressive strength; σ ct denoting the initial creep stress; σ ∞ indicating the long-term strength of rock; σ f is the residual stress; k representing the ratio of σ ∞ to σ c ; E c being the slope of the creep termination trajectory line; and E f representing the post-peak deformation modulus.For ease of calculation, the creep's initial stress is defined as where the stress exceeds σ ct and remains constant until reaching σ ∞ , at which point the creep will intersect the creep termination curve, and creep will cease.The endpoint of creep corresponding to such a stress level nearly lies on the same straight line.If the stress level remains constant at σ c1 , then the rock produces the maximum creep strain value ε e at the current stress level, along with the instantaneous elastic strain ε c .To enhance the accuracy of calculations and facilitate subsequent development, an adjustment was made based on the traditional "three-line model" to account for the initial creep stress and include the trajectory of creep termination.
To obtain the expression of the upper limit of creep ε c , the slope of the creep termination curve is first determined. First: Merge Equation ( 4) with Equation ( 5), Then where the proportion coefficient k is 0.12, then: According to the basic assumption of elastoplastic mechanics, it is assumed that the rock mass is a continuous and isotropic medium.The rheological deformation of the rock mass does not cause volume deformation; therefore, the volume stress of the viscoplastic part, including the Bingham body and Kelvin body, is neglected, focusing solely on the deviatoric strain response, which is the fundamental source of creep deformation.This can be expressed as: The changes in the creep deviatoric stress formula due to variations in the slope of the creep termination curve and the derivation of the 3D central difference scheme for the constitutive equation remain consistent with the BVS model and are detailed in the literature. 40

| CONSTITUTIVE MODEL VERIFICATION AND COMPARISON
During the development of this model, the "compilation separately" approach was employed. 50In this approach, functions, variable declarations, macro definitions, and structure definitions are written into the.hheader file, while the specific logic implementation of functions, variable initialization, constitutive model algorithms, stress-strain feedback, and elastic-plastic corrections are written into the .cppsourcefile.This C++- developed constitutive model can be called by FLAC 3D .
F I G U R E 3 Slope of creep termination curve and maximum creep strain value.

| Component verification
To ensure the correctness of the Kelvin body and Maxwell body components and prevent them from causing calculation errors in the constitutive model, verification is necessary.The SingleZone Oedometer test method is used to check the creep components in the constitutive model by applying single-axis lateral confinement and a dead load, 50 as shown in Figure 4A.Since FLAC 3D 's builtin Burgers and Burgers-Mohr models contain Kelvin body and Maxwell body components, the viscosity of the Kelvin body in the NVS model is set to infinity to disable its influence during the Maxwell body test.For Kelvin's verification, the viscosity of the Maxwell body in the NVS model is set to infinity.A comparison with FLAC 3D 's built-in Burgers and Burgers-Mohr models reveals consistent stress-strain responses of the creep element, confirming that the viscoelastic part of the constitutive model is correctly programmed, as illustrated in Figure 4B,C.According to literature data, the stress in this model is lower than 12% of the uniaxial compressive strength (below 4.5 MPa), while the NVS model exhibits only elastic strain response, as shown in Figure 4D.

| Constitutive model comparison
As the NVS model is an improvement upon the BVS model, it is essential to compare and analyze the numerical performance of both models in FLAC 3D .
F I G U R E 4 Loading method and verification results.

| Model size and mesh division
The numerical model, depicted in Figure 5, was created using FLAC 3D .The roadway had dimensions of 5.5 m in width and 4.5 m in height, with the model boundary situated 35 m away from the excavation boundary.Initially, without considering support, this problem can be simplified as a plane strain problem, and the model can be set to a unit thickness.However, when considering bolt support, it does not fall under the category of a plane strain problem.Therefore, the bolt (cable) row distance should be taken into account when determining the model thickness.In this simulation, an on-site support scheme was employed, with a bolt (cable) row distance of 0.95 × 1.2 m, anchoring section length of 1 m, a yield load of 200 kN, and a preload of 100 kN.For the side anchors, the row distance was 2.4 m, the anchoring section was 2 m in length, the yield load was 300 kN, and the preload was 200 kN.Analyzing the middle section of the model thickness, the model thickness was set at 4.8 m for comparative analysis between schemes, and the thickness of the model without support was also 4.8 m.The mesh was refined within 20 m from the roadway's excavation boundary, while the remainder was relatively sparse, resulting in a total of 57,596 units being divided.

| Initial conditions and boundary conditions
The tunnel is buried at a depth of 700 m, with a vertical ground stress σ z = 17.5 MPa.To facilitate a more comprehensive comparison between schemes, a larger side pressure coefficient is chosen, leading to more severe deformation and damage of the roadway surrounding rock.Therefore, the side pressure coefficient is set to 1.2, and the horizontal ground stress is σ x = σ y = 21 MPa.The model's base is fixed, and the normal displacement is restricted on the sides, while the overlying rock weight is accounted for as a stress boundary on the top.

| Material parameters
Due to the complexity of rock mass properties such as joints, cracks and faults, to eliminate the influence of other factors, only the response process of bias stress simulation is analyzed.In subsequent engineering verification, the roof and floor of roadway are mainly sandy mudstone.
So this simulation is based on the analysis of homogeneous argillaceous sandstone roadway.Considering that the strength and elastic modulus of the on-site rock mass are significantly lower than those of the rock sample, the parameters of the rock sample should be adjusted for the calculation.The surrounding rock is considered as common soft rock with the following properties: uniaxial compressive strength of 20 MPa, elastic modulus of 6 GPa, and Poisson's ratio of 0.2.The reduction methods for other parameters in the constitutive model are detailed in the literature, and the parameters in Table 1 are derived accordingly.

| Simulation analysis
After tunnel excavation, the surrounding rock does not exhibit sudden and severe deformation but rather undergoes gradual rheological deformation as the mining face advances.The average tunneling speed at Zhaozhuang Mine is 7 m/day.In this numerical simulation of tunneling, due to the short duration of this process, rheology is ignored, and only instantaneous elastic-plastic calculations are performed.Subsequently, the rheological switch is activated, and the simulation continues for 50 days.Since deviational stress primarily influences the deformation and failure of the surrounding rock in the roadway and controls the creep process, a comparative analysis of the deviational stress in the surrounding rock for both BVS and NVS models, with and without support, is conducted.
It can be observed from Figure 6 that in the absence of roadway support, the deviator stress values distributed at the roof and floor of the roadway in the NVS model are 21.7 MPa, which is smaller than those in the BVS model.
The distribution of deviator stress is more gradual, indicating improved calculation precision after accounting for the initial creep stress in the numerical calculation process.When the roadway is supported, the distribution area of partial stress around the roadway in the NVS model is significantly reduced compared to the BVS model.It shows higher sensitivity to supporting components, leading to an improved stress state of the roadway surrounding rock.The partial stress is mainly concentrated at the corner of the roof, floor, and near the shoulder socket of the roof.The high partial stress area in the middle of the roof and floor is also significantly reduced, aligning better with the actual on-site conditions.

| Engineering background
The research focuses on the intake roadway of the 1315 working face in Zhaozhuang Coal Industry.It analyses the deformation and failure behavior of surrounding rock in a goaf excavation roadway with deep cutting and unloading small coal pillars, and verifies the accuracy of the NVS model.
The layout of the working face and roadway for this project is depicted in Figure 7.The 1315 working face is situated in the 1 pan area and extracts coal from the 3# coal seam.The working face's ground elevation ranges from +1013.4 to +1154.3 m, while the coal seam floor elevation is between +341 and +412 m.The cover hill thickness above the working face is approximately 700 m, with an average coal seam thickness of 4.5 m.During mining at the 1312 working face, in alignment with production requirements, top pressure relief was implemented.The narrow coal pillar has a width of 6 m, and the top and bottom of the coal seam primarily consist of weak mudstone and sandy mudstone.
The 13,153 intake roadway is positioned within the coal seam, with a cross-section width × height of 5.5 × 4.5 m and a net fault area of 24.75 m 2 .The supporting parameters are illustrated in Figure 8. Upon field observation, it was noted that despite the adoption of the support scheme, the surrounding rock of the roadway still exhibited partial fragmentation and severe deformation.

| Numerical simulation
Based on the geological data of the 13,153 intake roadway in Zhaozhuang Mine, a numerical model, as shown in Figure 9, was created.The model dimensions are 382.5 m (X) × 200 m (Y) × 72.5 m (Z), where X = 1312 half the length of the working face + width of the narrow coal pillar + width of the 13,153 intake roadway lane +1315 half the length of the working face.The value of Y is determined according to the mine's monitoring requirements, and Z is based on geological data.To expedite the numerical simulation, a portion of the model was refined and equipped with the NVS constitutive model based on the research focus, while the rest was assigned the Mohr-Coulomb model.Assuming the tunnel is buried at a depth of 700 m and the average bulk density of the rock mass is 25 kN/m 3 , the vertical ground stress σ z = 17.5 MPa.A higher lateral pressure coefficient was chosen to induce more significant roadway deformation and damage, so a lateral pressure coefficient of 1.2 was set, resulting in an initial horizontal stress of 21 MPa.The model base was fixed, lateral displacement was constrained, and stress boundaries were employed to simulate the overlying rock load.To account for gravity-induced stress gradients, a gradient value of 0.025 MPa was applied using the "gradient" keyword.Due to the use of the INITIAL command for element stress initialization, there was a mismatch between strain and stress in the Maxwell and Kelvin bodies of the NVS model.To rectify this, a Fish function was programmed to predefine k_exx, k_eyy, k_ezz (principal strains in x, y, and z directions of the Kelvin body), k_exy, k_exz, and k_eyz (shear strain components of the Kelvin body) to match their corresponding strain values.The NVS model takes into account the initial creep stress, necessitating modification of the formula for calculating deviational stress in the Maxwell body to align with its corresponding strain value.The parameters required for the NVS constitutive model and the parameter reduction algorithm were obtained from the literature, 40 and the rock mechanics parameters for each rock stratum are presented in Table 2. Accounting for the goaf subsidence following working face excavation, the stope's pressuring step was set as the backfilling step to fill the goaf.
Considering the dilatancy of caving rock mass in goaf, most of them are low-cement granular materials, therefore, the null model is replaced by the Mohr-Coulomb model with low rock mechanical strength parameters, and all stresses are removed from the backfilled zones.Ini keyword is used to set the zone stress, gridpoint velocity and displacement to 0 to complete the stress simulation of the backfilled materials.A cable structural unit was employed to simulate tunnel anchor (cable) support, as depicted in and mining activities, in addition to calculating the instantaneous elastoplastic equilibrium, time-dependent rheological calculations were performed.The excavation and calculation steps of the model are as follows: 1312 mining face → conducting rheological calculations consistent with on-site conditions for half a year → driving 13,153 intake roadway → 1315 mining face.Since FLAC 3D employs the finite difference method as its core algorithm, the time-related stress change in creep calculations should not significantly exceed the stress change based on strain conversion.Otherwise, the excessive unbalanced force generated by the difference between the two can affect the model and lead to the failure of numerical convergence.According to the FLAC 3D manual, the iteration time increment is explicitly calculated.In the creep calculation of the NVS model, the maximum time step cannot exceed η N 0 /G H or η K 0 /G K .As shown in Table 2, the minimum value of the two is 0.7 h, so Δt = 0.1 h is used to ensure the model's calculation convergence.

| Influence of adjacent working face (NVS model compared to BVS model)
After the completion of mining in the 1312 working face and subsequent goaf backfill, the goaf's edge enters the stress reduction zone.The coal wall facing the goaf transitions from its pre-mining three-way stress state to a two-way stress state due to stress redistribution caused by the secondary stress field.This concentration of stress in the coal wall leads to local areas bearing stress that exceeds the strength of the coal body, causing it to gradually deteriorate and enter a state of plastic flow.In this context, k c is defined as the strength deterioration ratio, which is the cohesion c divided by the peak cohesion c 0 of the surrounding rock.As shown in Figure 10A, within 3 m to the left of the 1312 working face after backfilling, the coal pillar's strength deterioration ratio is less than 0.3, and in its middle region, the deterioration ratio is 0.1.This indicates that the edge strength deterioration ratio of the coal pillar, acting as a protective barrier during goaf excavation, is about 0.2.Under the influence of high stress, a fracture zone forms within 3 m of the goaf side, which aligns closely with the observations made by miners on the wall of 13,121 lane, as shown in Figure 10B.However, in Figure 10C, the strength deterioration ratio of the coal body and goaf after backfilling within 1 m of the hole depth is about 0.4, which somewhat differs from actual observations.Therefore, the NVS model proves more suitable for simulating the surrounding rock of deep small coal pillar excavation along the goaf.This further validates the assertion made in Section 3.2.4 that the "NVS model is suitable for deep soft rock roadway."The simulation results of the NVS model will be the focus of the following discussion.Figure 10 illustrates that the edge coal body undergoes unloading after fracturing, and the stress is transferred to the coal pillar.The deterioration ratio of the deeper coal pillar is approximately 0.85, and it exhibits an elastic core in a triangular shape.
Based on the actual operational conditions at Zhaozhuang Mine, the coal body's strength deterioration occurred  approximately 6 months after the cessation of mining in 1312 and before the excavation of the ventilation channel in 13,153, as illustrated in Figure 11.A comparison between the two figures reveals that the future coal pillar's strength deterioration ratio, which remains below 0.3, stabilizes over time, and the range of its fracture zone does not significantly expand.As the plastic zone remains stable and not enlarged, the deformation rate of the coal wall at the edge decreases, gradually stabilizing and approaching rheological deformation.Meanwhile, the elastic core area on the left contracts significantly, indicating that rheology has a more pronounced impact on the coal pillar as time progresses.
During this period, the strength of coal and rock mass at the limit equilibrium in the deterioration area gradually deteriorates, leading to rheological damage.This causes the elastic-plastic interface to move deeper with time, intensifying the fracture of coal and rock mass.The direction of deterioration shifts towards the roof at an angle of about 105°.

| Impact of tunneling
As depicted in Figure 12, following the design specifications of lane 13,151 in Zhaozhuang Mine, the narrow coal pillar roadway is driven under peak residual abutment pressure.The roof rock of the roadway is a weak layer with low strength, comprising black mudstone and black sandy mudstone.Approximately 2 m away from the tunneling face, the original balance is disrupted by goaf excavation, leading to stress concentration in the surrounding rock of the roadway.During the formation of a new plastic zone, the narrow coal pillar undergoes destruction and unloading, causing it to strongly shift toward the roadway.The deterioration at the shoulder socket of the roadway working face is substantial, with severe fracturing occurring within a 1 m length and 2.8 m height range.Displacement in this area measures approximately 300 mm, while displacement near the bottom plate is around 180 mm.The upper top corner of the coal pillar side experiences significant degradation, extending to about 2.8 m in height, with displacement approaching 380 mm.
The lateral stress concentration and long-term rheological effects resulting from the adjacent working face make the coal pillar side more severely damaged than the working face side.The pressure on both sides of the roadway is higher, leading to stress concentration in the middle of the roof's surrounding rock, resulting in roof subsidence in the roadway area, with displacement of around 300 mm.The roadway floor comprises black mudstone and sandy mudstone.Due to the absence of support measures for the bottom angle and floor of the roadway, increased roof and side pressure in the roadway leads to stress concentration in the floor, causing plastic deformation and shear failure.This results in floor heave of approximately 200 mm.
At a distance of 4 m from the tunneling face, the area with a deterioration strength ratio of 0.1 on the working face wall of the roadway gradually expands, moving towards the bottom plate.The bottom angle experiences transferred stress, and the deterioration area connects with the coal pillar wall.The bottom angle displacement is approximately 220 mm, about 20 mm higher than the floor heave.The triangular elastic core area on the left side of the roadway significantly reduces, extending the broken area of the coal wall to 3.3 m and deepening the damage range.The deterioration strength ratio near the bottom side of the coal pillar wall shifts from 0.3 to 0.15, indicating further deepening of the breakage.The plastic area of the roof remains relatively stable, with no significant increase in displacement.
As the tunnel advances to 6 m, the area with a deterioration strength ratio of 0.1 on the working face wall slowly expands, and the elastic core range further reduces.This suggests that the compressive shear stress resulting from raw rock stress in different directions exceeds the rock strength in the deep coal face wall.The elastic area under high pressure transitions into the fracture zone and plastic zone, with the coal pillar wall's crushing area shifting toward the top angle.The displacement of the roadway's surrounding rock remains relatively unchanged.
When the roadway is 15 m away from the tunneling face, the strength deterioration ratio of the surrounding rock shows no significant change.Over time, the deformation rate of the surrounding rock gradually stabilizes, indicating that the area has maintained stability.
Field observations during tunnel excavation corroborate the accuracy of the numerical simulation results, which closely match actual working conditions.

| Effect of rheology
In a deep, high-stress environment, the surrounding rock of a roadway typically does not experience plastic failure shortly after excavation.Instead, it often exhibits strong rheological characteristics over time.As indicated in Figure 12, after moving 15 m away from the excavation face, the surrounding rock failure caused by stress redistribution from excavation gradually stabilizes.Generally, coal and rock masses exhibit rheological properties, leading to slow deformation of the surrounding rock over time.However, the rate of deformation is significantly lower than that observed during the early stages of excavation.
Figure 13 illustrates that due to stress adjustments in the coal pillar wall of the roadway, the coal body's strength deteriorates.Consequently, the area with a deterioration strength ratio of approximately 0.1 continues to expand over time, gradually increasing in scope.Building upon the damage to black mudstone, the damaged area progressively moves downward and deepens, exacerbating the damage to sandy mudstone.The displacement cloud map reveals that the displacement of the fractured coal pillar wall increased from 380 to 480 mm within 50 days, with rheological effects accounting for 20% of this displacement.

| Effects of mining
As illustrated in Figure 14, mining activities at working face 1315 induce secondary disturbances in the stress distribution of the roadway's surrounding rock.Stress redistributes to establish a new equilibrium state, further deteriorating the roadway's surrounding rock.The elastic-plastic limit equilibrium zone continues to deepen as a result.
When the 1315 working face advances to within 30 m of the roadway section, a sharp increase in surrounding rock stress occurs due to the combination of advancing abutment pressure and residual abutment pressure in the goaf area adjacent to the working face.The primary degradation area of the coal pillar on the working face remains concentrated in the middle of the wall near the roof, with a larger broken area near the shoulder, approximately 4 m wide.The strength deterioration ratio of the surrounding rock ranges between 0.1 and 0.3.
Stress originating from the working face side propagates through the deep triangular-shaped elastic core area at an angle of 105°to the top of the elastic core.It then follows a direction of approximately 45°around the triangle, causing the plastic zone on the working face side to expand once more.The surrounding rock deformation increases by about 900 mm.On the coal pillar side, which is less affected, the larger displacement area is primarily concentrated around the top corner, resulting in approximately 750 mm of displacement.The bottom corner area near the bottom plate exhibits greater displacement, approximately 350 mm.The overall roof subsidence reaches 600 mm.
When the working face advances to within 10 m of the roadway section, the superposition of advancing support pressure and residual abutment pressure reaches its maximum.Consequently, the deterioration zone of the working face wall expands again, extending to approximately 6.5 m and primarily concentrated near the bottom corner of the floor.
As mining progresses and reaches a point 2 m from the mining face, stress redistribution becomes relatively stable.Deformation of the roadway's surrounding rock becomes less pronounced than the initial pressure from the front abutment stress.Stress on the working face side transfers to the triangular elastic core, causing the deterioration area to expand.The displacement at the shoulder socket on the working face side measures approximately 1200 mm, while the displacement of the top and bottom, as well as the coal pillar side, remains relatively stable.
In summary, the roof and floor of the roadway consist mainly of black mudstone and sandy mudstone, lithologically weak with low strength and poor self-stability.Under the influence of long-term structural stress and gravity stress, the surrounding rock is relatively loose.Prolonged exposure of the surrounding rock due to delayed grouting exacerbates its deterioration and weakens its load-bearing capacity.This results in severe deformation in the shoulder sockets of the surrounding rock and significant floor heave in the bottom corner of the floor, ultimately compromising the overall self- | 473 stability of the roadway.The roof is squeezed into the roadway space near the narrow coal pillar side.

| Analysis of monitoring data
Due to severe deformation near the top and bottom angles of the 13,153 intake roadway, a cross-layout method was employed to position measurement points at the top and bottom of the roadway, as well as the shoulder sockets on both sides, to monitor surface displacements.The results are depicted in Figures 15  and 16.
Within the initial 8 days following tunnel excavation, the deformation rate of the roadway's roof and sides exceeded 20 mm/day, with the maximum deformation  rate observed on the roof and coal pillar side reaching 70 mm/day.As time progressed, the deformation rate gradually decreased and stabilized.The displacement increment attributable to continuous rheological effects was approximately 130 mm.
Upon completion of mining activities in this working face, the displacements measured on the two sides and the top plate were 700, 900, and 500 mm, respectively.This indicates that the displacement caused by rheology accounts for roughly 19%, 15%, and 26% of the total displacement observed on the two sides and the top plate.In contrast, the bottom plate exhibited minimal rheological deformation.The deformation of the roadway's surrounding rock suggests that stoping disturbance is the primary influence on rheology, with other factors playing a secondary role.
The comparison between the numerical simulation results and the measured data is shown in Table 3.It can be seen that the numerical simulation results closely align with the observed deformation of the surrounding rock during the excavation of the 13,153 intake roadway.Given that the working face of 1315 has not yet been mined, the simulated working face's mining results can serve as a reference.

| CONCLUSION
Using the BVS constitutive model as a foundation, we derived the creep equation for each element, taking into account creep starting stress, and subsequently constructed the NVS model.We quantitatively described the maximum creep strain of rock at a specific stress level by calculating the slope of the creep ending curve using different deformation stages' boundary points on the total stress-strain curve.
In a further development phase, we conducted a comprehensive comparison and verification of the internal constitutive model of FLAC 3D with the modified model to ensure accuracy when incorporating creep starting stress.The results validated the correctness of our programming.A comparative analysis between the BVS and NVS models revealed that, in scenarios without roadway support, the NVS model provided more accurate calculations of deviational stress.Conversely, when roadway support was present, the simulation results closely matched real-world conditions, enhancing the NVS model's applicability to soft rock roadways.
We applied the NVS model to lane 13,153 and achieved numerical simulation results that closely resembled field observations.Mining activities near the working face led to stress redistribution, resulting in fracture zones within small coal pillars, while rheological processes caused deep shifts and gradual degradation of the elastic-plastic limit equilibrium zone in the deep rock mass.This phenomenon primarily occurred along a 105°direction from the roof, causing severe deformation at the shoulder socket of the working face, the top angle of the coal pillar, and the bottom angle of the floor.Rheology transformed the coal body, particularly the triangular elastic core, from its initial state as an elastic zone under high pressure to a fracture zone and plastic zone after deformation and failure.Consequently, the surrounding rock of the roadway became weak and fractured, with prominent coal softening and reduced load-bearing capacity.
Mining operations significantly increased displacement, particularly in the two sides of the roadway, reaching up to 1600 mm.The shoulder socket of the working face experienced pronounced damage.Stress in the surrounding rock after unloading transferred deep into the elastic core, with the extent of plastic deformation primarily reliant on the scope of the plastic zone and fracture zone within the roadway.The expansion of the plastic zone exhibited a noticeable time-dependent effect, with deformation moderating once this zone ceased to expand.Rheology contributed to 15%, 19%, and 26% of the displacement in the working face, pillar, and roof of the surrounding rock, respectively, while having minimal impact on the floor.After conducting underground field surveys, we confirmed that the numerical simulation results aligned with the changes in surrounding rock behavior in the excavation section.Given that mining activities had yet to commence in the working face of 1315 at the time, our simulated deformation patterns still held significance for understanding the potential impact of mining on the roadway working face.
In the process of numerical calculation, the velocity and displacement are calculated by the equation of motion, and the new stress is obtained from the constitutive equation by the known strain rate, which is the core of the secondary development.Then, the KANG ET AL.| 465

Figure 9 .
Based on on-site construction operations, the average mining speed of the 1312 working face was 7.8 m/day, the average mining speed of the 1315 working face was 4.3 m/day, and the average tunneling speed of the 13,153 intake roadway was 7 m/day.After simulating excavation F I G U R E 7 Layout of the working face and roadway.F I G U R E 6 Deviational stress cloud diagram.

F
I G U R E 8 Roadway section support diagram.F I G U R E 9 Numerical model diagram.

F
I G U R E 10 Strength deterioration ratio (c/c 0 ) of the surrounding rock and peep result.

F
I G U R E 11 Strength deterioration ratio (c/c 0 ) of the surrounding rock after half-year stoppage of 1312 working face.F I G U R E 12 Displacement and deterioration of roadway surrounding rock during tunneling.

F
I G U R E 13 Displacement and deterioration of roadways surrounding rock during rheology.KANG ET AL.

F
I G U R E 14 Displacement and deterioration of roadway surrounding rock during mining.F I G U R E 15 Changes in surface displacements of the surrounding rock.

F
I G U R E 16 Deformation curve of roadway during stoping stage.
Parameters for model comparison.
T A B L E 1 Parameters in the numerical model.
T A B L E 2 T A B L E 3 Comparison between numerical simulation and measured data of roadway deformation.properties of the model elements are verified by numerical examples, which has certain reference value for the development of research ideas for users. mechanical