Transfer and dissipation of strain energy in surrounding rock of deep roadway considering strain softening and dilatancy

In this study, the transfer and dissipation of strain energy in the surrounding rock of a deep roadway were analyzed, considering the objective strain softening and dilatancy behaviors. The strain energy increment was decomposed, and its variation was analyzed based on the incremental plastic flow theory; then, a numerical simulation was conducted for verification and further analysis. The results were verified by the field monitoring data of a coal mine gateway. The results show that in the elastic stage, the volumetric elastic strain energy density Uev decreases, while the shear elastic strain energy density Ues increases. In the plastic stage, both Uev and Ues decrease. The volumetric plastic strain energy density Upv is negative, and its absolute value increases, leading to strain energy accumulation. In contrast, the shear plastic strain energy density Ups is positive and increases, leading to strain energy dissipation. Considering strain softening, the elastic strain energy decreases, the plastic strain energy increases, and the region of strain energy dissipation expands. Considering dilatancy, the plastic strain energy varies more significantly, and the effect of strain softening is amplified. The strain energy is transferred from the deep part to the shallow part of the elastic zone and then to the plastic zone. The preexisting and input strain energies in the plastic zone are transformed into considerable amounts of plastic strain energy and then dissipated. Thereafter, a significant plastic strain appears, leading to the large deformation of the surrounding rock.

According to the statistical data of strain energy storage and release in the process of rock loading and unloading, Gong et al 9,10 proposed some indices to evaluate the bursting liability. Song et al 11,12 analyzed the rock burst mechanism of the surrounding rock of roadways from the perspective of strain energy storage and dissipation. Xie et al [13][14][15] adopted the strain energy method to study the surrounding rock deformation and support design of roadways.
However, these studies mostly focused on the elastic strain energy, whereas the plastic strain energy, which is a main form of initial strain energy dissipation, 16 has rarely been studied. Moreover, the strain softening and dilatancy behaviors were either neglected or not completely considered in the above studies. Extensive laboratory tests and field investigations have shown that a strain softening model can better describe the postpeak mechanical behavior of most rock masses than perfectly elastoplastic and elasto-brittle-plastic models. [17][18][19][20][21] Dilatancy is closely related to the deformation and failure process of rock masses, [22][23][24] and the dilation angle decreases with an increase in the plastic shear strain or confining stress. [25][26][27] For a deep roadway, the surrounding rock exhibits a high initial strain energy density and large plastic zone under high in situ stress conditions. 28,29 As a result, a large amount of initial strain energy is converted into plastic strain energy and dissipated during the surrounding rock deformation process. Additionally, the consideration of strain softening and dilatancy causes a significant increase in the plastic zone and strain in the deep roadway surrounding rock, 30,31 inevitably affecting the transfer and dissipation of strain energy.
To determine the law of transfer and dissipation of the strain energy in the surrounding rock of a deep roadway, this work analyzes the variations in the elastic, plastic, volumetric, and shear strain energy densities of the surrounding rock after excavation, based on the incremental plastic flow theory and considering strain softening and dilatancy. A numerical simulation is conducted for verification and further analysis under four conditions: perfectly elastoplastic, considering dilatancy only, considering strain softening only, and considering both strain softening and dilatancy. The strain energy density components of the surrounding rock are monitored to evaluate the effects of strain softening and dilatancy. Then, considering the objective strain softening and dilatancy behaviors, the transfer and dissipation of the strain energy are analyzed to explain the large deformation phenomenon of the surrounding rock of a deep roadway. Finally, the stress monitoring, acoustic wave velocity measurement, and borehole TV imaging of a coal mine gateway are performed to verify the results.
The strain energy method (both theoretical and numerical) is widely used in roadway support design and deformation prediction. This study presents the effects of strain softening and dilatancy behaviors on the transfer and dissipation of the strain energy in the surrounding rock of a deep roadway. This can aid researchers and engineers to appropriately consider these behaviors when applying the energy method, thus improving the reliability of the roadway support design to reduce the repair rate and cost.

| Problem description
To exclude other influencing factors and considering that the lateral pressure coefficient gradually approaches 1 as the depth increases, 32 a rectangular roadway excavated in homogeneous mudstone and affected by hydrostatic in situ stress p 0 is considered as the research object, as shown in Figure 1.
Compared with the stress-strain response rate in laboratory rock tests, the surrounding rock deformation of a roadway is a slow process, and the maximum deformation typically appears at a distance of 10 m or more behind the heading face 33 as a result of the excavation spatial effect. As shown in Figure 1A, as the heading face gradually advances, F I G U R E 1 Rectangular roadway affected by hydrostatic in situ stress p 0 : (A) profile along roadway axis, and (B) profile perpendicular to roadway axis the undisturbed rock mass to be excavated in front of it gradually experiences lateral unloading and fracturing; then, the surrounding rock begins to move toward the center of the roadway section. The rock mass behind the heading face has been excavated. However, under the effect of the spatial constraint, the deformation does not instantly reach its maximum value; instead, it gradually approaches the maximum value as the heading face advances. In this process, the plastic zone of the surrounding rock is formed at the excavation boundaries and gradually extends deeper. 23,34 As shown in Figure 1B, the surrounding rock in the plastic zone overcomes the constraint within the roadway section, p c (the constraint of the rock mass to be excavated or the spatial constraint), and moves toward the center of the roadway section, indicating that the strain energy in the plastic zone is transferred to the excavation space. Moreover, the surrounding rock in the elastic zone moves toward the center of the roadway section and compresses the surrounding rock in the plastic zone, indicating that the strain energy in the elastic zone is transferred to the plastic zone. In the plastic zone, the unloading of the surrounding rock leads to a significant decrease in the elastic strain energy, which is practically completely transformed into plastic strain energy and then dissipated. 16

| Decomposition of strain energy density increment
To further explain the transfer and dissipation of strain energy and analyze the effects of strain softening and dilatancy, the strain energy density increment is decomposed based on the incremental plastic flow theory. For concision, the principal stress space is adopted, and compression and contraction are considered positive herein. Subscript i = 1, 2, 3 represents a component and applies to Einstein's summation convention. Subscript m represents the mean value of the components, α m = (α 1 + α 2 + α 3 )/3. Superscript d represents the deviatoric value of a component, d i = α i − α m . For a microunit of the roadway surrounding rock, the total elastic strain energy density U et can be written as 35 where U ev is the volumetric elastic strain energy density, U es is the shear elastic strain energy density, σ m is the mean stress (volumetric stress), d i is the deviatoric stress, and K and G are the bulk and shear moduli, respectively.
The incremental forms of U ev and U es can be derived according to the generalized Hooke's law: where ΔU ev is the volumetric elastic strain energy density increment, ΔU es is the shear elastic strain energy density increment, U evo is the volumetric elastic strain energy density of the last cycle, U eso is the shear elastic strain energy density of the last cycle, σ mo is the volumetric stress of the last cycle, d io is the deviatoric stress of the last cycle, m is the mean value of the volumetric stresses of two adjacent cycles, d i is the mean value of the deviatoric stresses of two adjacent cycles, Δ e m is the mean elastic strain increment, and Δ ed i is the deviatoric elastic strain increment.
According to the calculation methods for the total volumetric strain energy density increment ΔU tv and total shear strain energy density increment ΔU ts given by Itasca Inc, 35 the following can be derived: where Δ m and Δ d i are the mean and deviatoric strain increments, respectively.
According to the strain decomposition in the incremental plastic flow theory, the calculation methods for the volumetric plastic strain energy density increment ΔU pv and shear plastic strain energy density increment ΔU ps can be derived from Equations (2)-(5): where Δ p m and Δ pd i are the mean and deviatoric plastic strain increments, respectively.
Then, the total plastic strain energy density increment ΔU pt can be given by the sum of ΔU pv and ΔU ps as follows:

| Variations in strain energy densities and effects of strain softening and dilatancy
For a rectangular roadway, the deformation and plastic zone depth are larger at the middle (width direction) of the roof; due to symmetry, the radial stress σ r , tangential stress σ t , and axial stress σ a at these positions are always the principal stresses. Therefore, based on the stress adjustment process given by Yi et al, 30 a microunit at the middle of the roadway roof is considered for the analysis.
As shown in Figure 2, in the in situ stress stage, σ r = σ t = σ a = p 0 , σ m = p 0 , and d i = 0. The expression of the initial strain energy density U 0 can be derived from Equations (1)-(3) as After being influenced by the roadway excavation, the microunit enters the elastic stage. In this stage, σ t significantly increases, σ a slightly decreases, and σ r significantly decreases. The decrease in σ r is greater than the increase in σ t , leading to a decrease in σ m and an increase in d i . According to Equations (2) and (3), ΔU ev < 0 and ΔU es > 0, that is, U ev decreases while U es increases.
As the maximum principal stress σ 1 (σ t ) increases and the minimum principal stress σ 3 (σ r ) decreases, the microunit yields and enters the plastic stage. Only the shear failure is discussed herein because the failure mode of the surrounding rock of a deep roadway is mainly shear failure, and tensile failure only occurs near the excavation boundaries. 30,31 The M-C yield criterion f, which is widely adopted in engineering, 36,37 and the nonassociated plastic flow rule, which well describes the plastic behavior of rock masses well and is widely used in roadway deformation calculations, 38,39 are used for the analysis and given as follows: where Δ p i is the plastic principal strain increment, c is the cohesion, φ is the friction angle, ψ is the dilation angle, N α = (1 + sin(α))/(1 − sin(α)) is a transformation function, g is a plastic potential function, and λ is a constant that can be determined based on the consistency condition.
In the plastic stage, σ r , σ t , and σ a all decrease significantly, resulting in considerable decreases in σ m and d i . According to Equations (2) and (3), ΔU ev < 0 and ΔU es < 0. The plastic flow occurs with the unloading of the surrounding rock, generating a large amount of plastic strain energy. According to Equations (6) and (7), the plastic strain energy density increments are related to the stress and plastic strain increments. Δ p i can be derived from Equations (11) and (12) as follows: According to the mean and deviatoric values introduced in Section 2.1, the calculation methods for Δ p m and Δ pd i can be derived using Equation (13): When the microunit yields, λ > 0. If dilatancy is not considered, ψ = 0 and N ψ = 1. According to Equation (14), Δ p m = 0, and according to Equation (6), ΔU pv = 0. If dilatancy is considered, ψ > 0 and N ψ > 1. According to Equation (14), Δ p m < 0, and according to Equation (6), ΔU pv < 0. This indicates that the plastic volumetric expansion of the microunit occurs under the compressive volumetric stress σ m , the volumetric stress does negative work, and the elastic strain energy is accumulated. In other words, the dissipated elastic strain energy is negative via generating the volumetric plastic strain energy. According to Equation (15), Δ pd 1 > 0, Δ pd 2 ≥ 0, and Δ pd 3 < 0. According to Figure 2, σ 1 > σ m , σ 2 > σ m , and σ 3 < σ m ; then, d 1 > 0, d 2 > 0, and .
F I G U R E 2 Stress adjustment process and variations in strain energy densities of microunit at middle of roadway roof then according to Equation (7), ΔU ps > 0. This indicates that the deviatoric stress does positive work, and the elastic strain energy is dissipated. If dilatancy is considered, according to Equations (14) and (15), the absolute values of Δ pd i and Δ p m increase with an increase in ψ, as shown in Figure 3. Consequently, ΔU ps and the absolute value of ΔU pv increase. In other words, the plastic strain energy varies more drastically.
If strain softening is considered, the cohesion c decreases with an increase in the plastic shear strain in the plastic stage, and the yield surface enclosed by f in the principal stress space continually shrinks, 31 resulting in the decreases in σ i , σ m , and d i . Therefore, according to Equations (2) and (3), the absolute values of ΔU es and ΔU ev increase. Simply put, the elastic strain energy decreases. However, the variations in ΔU pv and ΔU ps cannot be determined according to Equations (6) and (7). Moreover, considering that strain softening leads to an expansion of the plastic zone of the surrounding rock, and the region of strain energy dissipation expands accordingly.
If both the strain softening and dilatancy are considered, the dilatancy will accelerate the strain softening process 31 and amplify its effect. Therefore, the elastic strain energy further decreases, and the range of strain energy dissipation in the surrounding rock further expands.
In summary, during the deformation process after roadway excavation, the strain energy in the elastic zone of the surrounding rock is transferred into the plastic zone, and the strain energy in the plastic zone is transferred into the excavation space. For a microunit of the surrounding rock, in the elastic stage, U ev decreases while U es increases. In the plastic stage, both U ev and U es decrease, and U pv is negative, leading to an elastic strain energy accumulation; U ps is positive, leading to an elastic strain energy dissipation. If strain softening is considered, the elastic strain energy decreases, and the region of strain energy dissipation expands. If dilatancy is considered, the plastic strain energy in the surrounding rock varies more drastically, and the effect of strain softening is amplified.
However, the above theoretical analysis is only qualitative, and cannot determine the variations in the total elastic strain energy U et and total plastic strain energy U pt , as well as the effect of strain softening on the variations in the plastic strain energy densities. These problems are solved, and the results of the theoretical analysis are verified in the following numerical simulation.

| Numerical model and schemes
A deep roadway problem can be simplified as a plane strain problem, 40,41 and the excavation spatial effect ( Figure 1) is equivalent to the slow reduction in the constraint within the roadway section p c . 35 Therefore, the thickness of the model is set to 1 m, and p c is set to gradually decrease from p 0 to 0, with a 5% reduction in each calculation before p c is reduced to 0.2p 0 (the slower the variation in p c , the higher the calculation precision in the plastic stage of the surrounding rock 35 ; however, the calculation time increases. At this stage, all or practically all of the surrounding rock is in the elastic stage, so a fast variation rate of p c is selected to improve the calculation efficiency); a 2% reduction in each calculation is set after p c is reduced to 0.2p 0 (at this stage, a large area of the surrounding rock enters the plastic stage, so the variation in p c is decelerated for higher precision). As shown in Figure 4, the roadway section is set to 3 m × 3 m, the model boundaries are located 30 m away from the excavation boundaries (10 times the dimension of the roadway to reduce the error caused by the artificial boundaries), and all the model boundaries are fixed in the normal directions (if the artificial boundaries are sufficiently far, the results obtained using this boundary condition and a stress boundary condition are almost the same). Almost uniform meshing is adopted to minimize the impact of strain localization problems 42 ; there are 16 129 zones and 32 768 grids in the model. The hydrostatic in situ stress condition is the same as that implemented in the above theoretical analysis, and the excavation depth is set to 1500 m (p 0 = 37.5 MPa). The surrounding rock of the roadway is set to be homogeneous mudstone, and the constitutive model established in our previous study 31 is adopted. This constitutive model can accurately consider the strain softening and dilatancy behaviors of rock masses based on the equations of cohesion c and dilation angle ψ given by Pourhosseini and Shabanimashcool 26 : where c 0 is the peak cohesion, ψ 0 is the peak dilation angle, γ p is the plastic shear strain, σ c is the uniaxial compressive strength (which can be determined by c 0 and φ), and a, b, A, and B are the fitting parameters depending on the rock type (a = 0.340, b = 0.125, A = 6.794, and B = 15.308 for mudstone 26 ).
The remaining parameters are assigned the common values of mudstone. 31 Young's modulus is set to 6.0 GPa, Poisson's ratio is set to 0.2, tensile strength is set to 1.0 MPa, peak cohesion c 0 is set to 6.37 MPa, and friction angle φ is set to 25° (σ c = 20 MPa). The numerical simulation is conducted under four conditions: perfectly elastoplastic (scheme A), considering only dilatancy (scheme B), considering only strain softening (scheme C), and considering both strain softening and dilatancy (scheme D). The material parameters are shown in Table 1.

| Variations in strain energy densities at measurement point
For a rectangular roadway, the deformation and the plastic zone depth at the middle (width direction) of the roof are generally larger, 31 so the strain energy densities vary more drastically at these positions. In addition, the deep part of the surrounding rock is not significantly disturbed, whereas tensile failure may occur at the excavation boundary. Therefore, the measurement point is selected at 1.25 m away from the excavation boundary and at the middle of the roadway roof, as shown in Figure 4. The energy calculation keywords and their combinations in the FISH language of FLAC 3D35 are adopted to monitor the strain energy densities, and the results are shown in Figure 5.
As shown in Figure 5, in the elastic stage, U ev tends to decrease with a decrease in p c (except that U ev slightly increases when the microunit nearly yields), and U es increases. In the plastic stage, U ev and U es decrease, U pv is negative, and its absolute value increases; U ps is positive and increases. This is consistent with the conclusions of the theoretical analysis above. In addition, Figure 5 also shows the variations in U et and U pt . In the elastic stage, U et first decreases and then increases; in the plastic stage, U et decreases while U pt increases.
Schemes C and D, which consider strain softening, enter the plastic stage earlier than schemes A and B, which not consider strain softening. In the plastic stage, based on the comparison between schemes A and B, considering dilatancy has almost no effect on U ev , U es , and U et , but causes U es , U et , and the absolute value of U pv to all increase, which is consistent with the conclusion in Section 2.3 that considering dilatancy leads to a more drastic variation in the plastic strain energy. The comparison between schemes A and C indicates that considering strain softening reduces U ev , U es , and U et , which is consistent with the conclusion in Section 2.3 that the elastic strain energy decreases when strain softening is considered. In addition, considering strain softening leads to increases in U ps , U pt , and the absolute value of U pv , which could not be determined in Section 2.3. By comparing schemes A and D, on the basis of considering strain softening, considering dilatancy causes U ev , U es , and U et to decrease, and U ps and the absolute value of U pv to increase, which is consistent with the In each scheme, ΔU et + ΔU pt > 0, that is, the increment in the total plastic strain energy density is larger than the reduction in the total elastic strain energy density. This indicates that the strain energy input occurs at the measurement point during the roadway deformation process, as discussed below.

| Transfer and dissipation of strain energy in surrounding rock
The objective strain softening and dilatancy behaviors of the surrounding rock are considered in scheme D, which is selected for analysis. The reduced elastic strain energy in the surrounding rock is almost completely transformed into the plastic strain energy and dissipated. 16 Therefore, it is assumed herein that the reduced elastic strain energy in the surrounding rock is completely transformed into plastic strain energy (ie, ΔU et + ΔU pt = 0, ΔU et + ΔU pt > 0, and ΔU et + ΔU pt < 0, representing the conservation, input, and output of the strain energy, respectively). Figures 6-11 show ΔU et + ΔU pt , ΔU et , ΔU pt , plastic zone, stress tensor, and cohesion of the surrounding rock, respectively. The black contours in Figures 6-11 represent ΔU et + ΔU pt = 0, which combined with the color of a region in Figure 6 can be used to determine whether the region is a strain energy output region (ΔU et + ΔU pt < 0) or an input region (ΔU et + ΔU pt > 0).
According to Figures 6-11, when p c = 0.5p 0 , all of the surrounding rock is in the elastic zone, and ΔU pt = 0. The strain energy is only transferred within the surrounding rock and not dissipated. The value of ΔU et + ΔU pt is only related to the stress adjustment. The stress decreases at the middle of each excavation boundary and increases at each corner of the roadway; consequently, a strain energy output region forms at the middle of each excavation boundary, with an output strain energy density U out of up to 0.04 MJ/m 3 (considering the maximum value in the region, the same below). A strain energy input region also forms at each corner of the roadway, with an input strain energy density U in of up to 0.10 MJ/m 3 . When p c = 0.2p 0 , a plastic zone appears near each excavation boundary of the roadway. However, both the range and deformation of the plastic zone are small, resulting in ΔU pt being considerably small, and the strain energy is almost not dissipated. The value of ΔU et + ΔU pt still depends on the stress adjustment. The stress shifts to the depth of the surrounding rock, resulting in U out of the strain energy output region increasing to 0.06 MJ/m 3 . The strain energy input region expands, and U in increases to 0.22 MJ/m 3 .
When p c = 0.04p 0 , the range and deformation of the plastic zone significantly increase, resulting in the strain energy being significantly dissipated. The tangential contraction and radial expansion at the middle of a roadway excavation boundary are mainly plastic deformation, 31 leading to an increase in the deviatoric plastic strain d i . Therefore, the total plastic strain energy density increases, and its increment exceeds the reduction in the total elastic strain energy density. Then, the strain energy output region at the middle of each roadway excavation boundary is transformed into a strain energy input region. The stress continues to shift to the depth of the surrounding rock, resulting in the continuous expansion of the strain energy input region, and U in continuously increases to 0.43 MJ/m 3 . The strain localization of the surrounding rock results in the discontinuity of strain softening in the plastic zone, and areas with high strength appear. These areas with low total plastic strain energy density belong to the strain energy output region, and U out reaches 0.15 MJ/m 3 .
When p c = 0, the range and deformation of the plastic zone further increase, and the strain energy continues to be considerably dissipated. After the stress shifts to the depth of the surrounding rock, the ultimate equilibrium is reached, the strain energy input region expands to the final state, and U in finally increases to 0.81 MJ/m 3 . The discontinuity of strain softening is more significant, resulting in both the number and range of the areas with high strength increasing, and U out finally increases to 0.16 MJ/m 3 .
As p c decreases to 0, the strain energy is transferred from the deep to the shallow part of the elastic zone, and then to the plastic zone. The existing and input strain energies in the plastic zone are transformed into plastic strain energy and dissipated, and the plastic strain energy is generated with a significant plastic strain according to Equations (6) and (7). The surrounding rock of a deep roadway has a higher initial strain energy density and larger plastic zone. After the roadway excavation, a large amount of plastic strain energy is generated, and a significant plastic strain appears, leading to a large deformation of the surrounding rock.

| Field monitoring scheme
To exclude other influencing factors, the above analyses were conducted under ideal conditions. Although such ideal conditions do not exist in reality and the strain energy cannot be directly monitored, a roadway with similar geological conditions can be selected to monitor the tangential stress σ t , and observe the acoustic wave velocity and damage degree of the surrounding rock to verify the obtained results. The excavation depth of the 050502 rail gateway of the Qingyun coal mine in Shanxi Province, China, is 802 m; the roof is a mudstone, the floor is a sandy mudstone, and the sidewall is a coal seam. The cross section of the gateway is rectangular with a width of 4.0 m and height of 3.1 m. As shown in Figure 12, the GYW25 type surrounding rock stress meter was used to monitor σ t ; the borehole diameter was 42 mm, and the initial pressure was set to 4.9-5.1 MPa. The CLC1000-Z type ultrasonic detector was used to measure the acoustic wave velocity; the single-borehole test method and water coupling were adopted; the borehole diameter was 38 mm, and readings were recorded at intervals of 200 mm from the borehole bottom outward. The CXK12 type borehole TV was used to observe the damage degree of the surrounding rock; the borehole diameter was 28 mm. Because the heading face was gradually advancing and each test work required some time, the positions of the boreholes in Figure 12 were relative to the heading face and did not coincide.

| Analysis of results
The monitoring data of σ t , measurement results of acoustic wave velocity, and borehole TV images are shown in Figures 13-15, respectively. Figure 13 shows that the stress meter reading at depth of 1 m remains the initial pressure 5.1 MPa, indicating that at a position 0.5 m away from the heading face, the depth of 1 m of the roof has entered the plastic stage, and σ t ≤ 5.1 MPa. After the installation, the stress meter readings at depths of 2-5 m rapidly rise as the boreholes deform. With the advance of the heading face, the stress meter readings at depths of 2 and 3 m successively increase, then rapidly decrease, and finally tend to stabilize. The stress meter readings at depths of 4 and 5 m continuously rise and finally stabilize. This verifies the result presented in Section 3.3 that the stress continues to shift to the depth of the surrounding rock. The plastic zone depth at the measuring station continuously increases and finally stabilizes between 3 and 4 m. After the surrounding rock enters the plastic stage, σ t significantly decreases, indicating that a large amount of elastic strain energy is transformed into plastic strain energy and dissipated, which is consistent with the conclusions in Sections 2.3 and 3.3.  The higher the damage degree of a rock mass, the lower its acoustic velocity. 43 As shown in Figure 14, borehole nos. 1-3 for acoustic wave velocity measurement have higher acoustic wave velocities (1.91-2.53 km/s) and tend to be stable at depths exceeding 1.8, 3.0, and 4.2 m, respectively. Therefore, the plastic zone depths at 0.5, 5, and 15 m from the heading face are 1.8, 3.0, and 4.2 m, respectively, indicating that the plastic zone depth of the surrounding rock increases with the advance of the heading face. In the plastic zone, the acoustic wave velocity is low (0.53 km/s at the lowest), indicating that the surrounding rock is severely damaged and has significant plastic deformation. Thus, a large amount of strain energy is dissipated, which is consistent with the conclusions in Sections 2.3 and 3.3 that the strain energy is dissipated in the plastic zone. Moreover, the variation in acoustic wave velocity in the plastic zone is nonmonotonic, indicating that the strain softening is discontinuous, which is consistent with the conclusion in Section 3.3 "that the strain localization of the surrounding rock results in the discontinuity of strain softening in the plastic zone." As shown in Figure 15, borehole nos. 1-3 for TV imaging have many open or slip fractures at depths less than 1.57, 2.71, and 4.03 m, respectively. Therefore, the plastic zone depths at 0.5, 5, and 15 m from the heading face are 1.57, 2.71, and 4.03 m, respectively (these boreholes do not coincide with those used for acoustic wave velocity measurement, so the plastic zone depths differ). This indicates that the plastic zone depth increases with the advance of the heading face. The opening or slipping of these fractures in the plastic zone will dissipate a considerable amount of strain energy, 44  is consistent with the conclusions in Sections 2.3 and 3.3 that the strain energy is dissipated in the plastic zone. Strain softening is closely associated with the microcrack growth in the surrounding rock, as well as the upscaling fracturing. [45][46][47][48] These fractures are not uniformly distributed in the plastic zone, indicating that the strain softening of the surrounding rock is discontinuous, verifying the conclusion in Section 3.3.

| CONCLUSIONS
In this work, based on the incremental plastic flow theory and considering strain softening and dilatancy, theoretical analysis and numerical simulation were adopted to analyze the transfer and dissipation of the strain energy in the surrounding rock of a deep roadway; the results were verified through the field monitoring data. The main conclusions are summarized as follows.
For a microunit of the surrounding rock, in the elastic stage, the volumetric elastic strain energy density U ev decreases, whereas the shear elastic strain energy density U es increases, and the total elastic strain energy density U et first decreases and then increases. In the plastic stage, U ev , U es , and U et all decrease; the volumetric plastic strain energy density U pv is negative, and its absolute value increases, leading to the accumulation of the elastic strain energy; the shear plastic strain energy density U ps is positive and increases, leading to the dissipation of the elastic strain energy; moreover, the total plastic strain energy density U pt increases. When strain softening is considered, the elastic strain energy decreases, the plastic strain energy increases, and the region of strain energy dissipation expands. When dilatancy is considered, the plastic strain energy in the surrounding rock varies more drastically, and the effect of strain softening is amplified.
As the heading face gradually advances, the stress also gradually shifts to the depth of the surrounding rock, and the plastic zone and strain energy input region gradually expand. The shallow part of the elastic zone and most of the plastic zone belong to the strain energy input region. The deep part of the elastic zone and areas with a lower strain softening level in the plastic zone belong to the strain energy output region. The strain energy in the surrounding rock of a deep roadway is transferred from the deep to the shallow part of the elastic zone, and then to the plastic zone. The existing and input strain energies in the plastic zone are transformed into a large amount of plastic strain energy and dissipated. A significant plastic strain then appears, leading to a large deformation of the surrounding rock.