A Unified Solution for Surrounding Rock of Roadway Considering Seepage, Dilatancy, Strain-Softening and Intermediate Principal Stress

Abstract

The analytic solution for surrounding rock of roadway is of significance for stability analysis and roadway support. However, analytical solution for surrounding rock of roadway which took influences of water seepage, strain softening, dilatancy and intermediate principal stress all into account did not receive much reporting. To promote research in this aspect, a mechanical model simultaneously considering water seepage, strain softening, rock dilatancy, and intermediate principal stress was established based on porous elastoplastic mechanics, and then unified analytical solution for surrounding rock of roadway was obtained. Based on an example, influences of water seepage, strain softening, rock dilatancy, residual cohesion and intermediate principal stress on surrounding rock of roadway were thoroughly investigated using single factor analysis. The obtained results are as follows: radii of plastic zones and surface displacement of roadway would increase exponentially with water pressure increasing and their magnitudes are greater than corresponding values without water seepage considered; with softening modulus increasing, peak circumferential stress location would slightly shift to deeper surrounding rock, while broken zone radius and surface displacement of roadway would increase in a decay velocity; rock dilatancy has little effect on peak circumferential stress and plastic softening zone radius, while broken zone radius and surface displacement of roadway increase linearly with dilatancy coefficient α1 increasing indicating their magnitudes are overestimated if associated flow rule is adopted; with weighted coefficient increasing, stress components in plastic zones at the same distance from roadway center would increase, while radii of two plastic zones and surface displacement of roadway are reduced, i.e., self-bearing capacity of rock is enhanced considering intermediate principal stress effect compared to Mohr–Coulomb criterion; with residual cohesion increasing, peak circumferential stress remains unchanged, while stress components in plastic zones at the same distance from roadway center would increase and radii of two plastic zones decrease significantly. The above results implicated that water seepage effect should be carefully considered for roadway stability under groundwater environment; strain-softening and flow rule of rock should be reasonably analyzed and chosen to accurately predict surface displacement and broken zone radius of roadway; rock bolt length should be increased with softening modulus increasing, while it can be decreased with intermediate principal stress effect considered; grouting measure is an effective measure to improve roadway stability. In short, the research provides a theoretical basis and some practical engineering implication for roadway support.

1. Introduction

Coal is of significance for the development of China’ s society and economy. With exhaustion of shallow coal resources, the coal exploitation has gradually entered deeper rock mass with mechanical characteristics of “three highs and one disturbance” [1], which affects the stability of underground roadway. Determination of deformation and plastic zone of roadway is of significance for stability analysis and roadway support. Analytic solution for surrounding rock of roadway, which provides direct understanding of plastic zones and deformation, is a hot issue studied by scholars [2,3,4,5,6]. There are many factors to be considered in deriving analytic expressions for surrounding rock of roadway. Firstly, groundwater always exists in underground rock mass of mining engineering. When the roadway is excavated, the original stress state would be destroyed, thus simultaneously generating a seepage field because of groundwater flow. The water seepage would then influence stress distribution of surrounding rock [7,8,9]. Hence, an accurate theoretical solution for surrounding rock should take the influence of groundwater into account. Besides, when circumferential stress of surrounding rock reaches its compressive strength during stress redistribution process, its stress would be reduced to its residual strength gradually; namely, surrounding rock of roadway would present post-peak strain-softening behavior [3,4,5,6]. Thus, elastoplastic analysis of surrounding rock should consider strain-softening characteristics. Furthermore, related studies have demonstrated that there is little volume change of rock in the pre-peak process, while rock dilatancy exists at the post-peak process, which would affect the deformation characteristics of surrounding rock [3,5,10,11]. Therefore, rock dilatancy should be also considered. Last but not least, surrounding rock of roadway is generally in three-dimensional stress state, while related studies showed that intermediate principal stress would affect surrounding rock strength to an extent [10,11,12]. Hence, to accurately describe failure or deformation characteristics for surrounding rock, influences of water seepage, strain softening, dilatancy and intermediate principal stress should be all considered [11,12,13,14,15,16].

Recently, regarding solutions for surrounding rock of roadway, numerous achievements have been obtained. Zhang [9] and Gao [15] obtained elastoplastic solution which considered effects of dilatancy, strain softening and seepage, while influence of second principal stress is ignored. Based on unified strength theory, Liu [17] obtained analytical expressions considering dilation characteristics and intermediate principal stress, but effect of strain softening and seepage were not considered. Considering effects of strain softening, dilatancy and intermediate principal stress, scholars [11,18,19,20,21] obtained analytic solutions based on different strength criteria, such as generalized spatially mobilized plane criterion, generalized three-dimensional Hoek–Brown criterion, modified Lade criterion and Drucker–Prager criterion, which indicated that intermediate principal stress had influenced on stress distributions, but the seepage effect was not considered. In comparison of four different strength criteria, Pan [12] and Fan [16] analyzed influences of seepage, rock dilatancy and intermediate principal stress, which indicated that the rock dilatancy was related to plastic potential function, but they did not distinguish the strain-softening zone and broken zone. For bolted roadway, Gu [10] had analyzed roadway stability in consideration of intermediate principal stress and dilatancy, but seepage effect was not considered; Zhou [22] proposed a new analytical solution with seepage, dilatancy, intermediate principal stress considered, but an elastic–brittle model is used, which is not suitable for describing strain-softening characteristics of surrounding rock without rock bolt or other support measures. Furthermore, scholars [23,24,25,26] obtained numerical solutions for surrounding rock of roadway considering influences of confining pressure and seepage on rock parameters, in which extra relatively complex calculation procedures were needed and not all factors were considered. In conclusion, both analytic and numerical solutions for circular roadway have been extensive studied recently, and the previous literature has analyzed influences of factors including water seepage, strain softening, dilatancy and intermediate principal stress on roadway stability, respectively. However, to our best knowledge, there is little reporting regarding analytical solution for surrounding rock of roadway, which takes influences of seepage, strain softening, dilatancy and intermediate principal stress all into account, and how these factors influence stress, plastic zone radius or displacement was not fully understood. The goal of the present study is to promote and improve research in this aspect. Therefore, based on existing studies, a mechanical model considering water seepage, strain softening, rock dilatancy, and intermediate principal stress was established according to porous elastoplastic mechanics in Section 3, and then unified analytical solutions for surrounding rock of roadway were obtained, which can represent analytical solutions for an ordered series of solutions. Based on an example, influences of each factor on stress distribution, radii of plastic zones and displacement were further studied using single factor analysis. The results provide a theoretical basis and some practical engineering implication for roadway support.

2. Unified Strength Theory

As mentioned above, the Mohr–Coulomb and Hoek–Brown criteria had widely been applied in rock mechanics. However, both two criteria have only considered the influence of minimum and maximum principal stresses, whereas intermediate principal stress has not been considered. Existing studies [10,11,12] revealed that intermediate principal stress would influence rock failure and deformation in a certain extent, which demonstrated intermediate principal stress should not be ignored. Hence, to obtain a reasonable analytical solution for surrounding rock of roadway, the unified strength theory proposed by Yu et al. [27] is adopted:

$$\{\begin{array}{c}\frac{2c\cos \phi }{\sin \phi -1}=-{\sigma }_1+\frac{\left({\sigma }_3+{\sigma }_{2}b\right)\left(\sin \phi +1\right)}{\left(b+1\right)\left(1-\sin \phi \right)} \left({\sigma }_2\le \frac{{\sigma }_3-{\sigma }_1}{2}\sin \phi +\frac{{\sigma }_3+{\sigma }_1}{2}\right)\\ \frac{2c\cos \phi }{\sin \phi -1}=-\frac{{\sigma }_1+{\sigma }_{2}b}{b+1}-\frac{{\sigma }_3\left(1+\sin \phi \right)}{\sin \phi -1} \left({\sigma }_2>\frac{{\sigma }_3-{\sigma }_1}{2}\sin \phi +\frac{{\sigma }_3+{\sigma }_1}{2}\right)\end{array},$$

(1)

where φ and c represent internal friction angle and cohesion; σ₃, σ₂ and σ₁ represent minor, intermediate and major principal stresses, respectively; b represents weighted coefficient accounting for influencing of σ₂ which ranges from zero to one. In the present paper, positive value is for compressive stress.

3. Mechanical Model for Surrounding Rock of Roadway

3.1. Distribution of Pore Water Pressure

A mechanical model, which is used to calculate deformation and stress calculation for surrounding rock of roadway, is established based on following assumptions: (1) the roadway is approximately considered as an equivalent circular cavern with infinite length, so the surrounding rock is in a plane strain state; (2) surrounding rock is assumed as elastoplastic porous media with homogeneous, isotropic characteristics; (3) the weight of surrounding rock, which is subjected to hydrostatic stress and pore water pressure, is ignored; (4) the permeability of surrounding rock is a constant in all directions. As illustrated in Figure 1, the roadway radius is r₀, and it is imposed by a hydrostatic stress of σ₀ before excavation. If the support pressure p_i is not great enough, two plastic zones, namely broken zone with radius being R_b and plastic softening zone with radius being R_p, would be commonly formed around the roadway, while the surrounding rock with radius greater than R_p is still under elastic state.

Before excavation of roadway, an initial pore water pressure p₀ is existed throughout surrounding rock (Figure 1). After roadway excavation, internal water pressure is decreased to zero. As is well known, the subsurface water will flow due to pore water pressure difference between inner and outer domain. According to above assumptions, the water seepage would only occur in radial direction, which results in the axisymmetric redistribution of pore water pressure. To obtain analytical solution of water pressure distribution, a boundary on which water pressure is unaltered should be determined. Related studies [2] showed that when a point is located on boundary with radius of R₀ ≥ 30 r₀, the alternation of pore water pressure p_w and stress can be negligible compared to infinite condition. Thus, the mechanical model is established with a maximum radial dimension of R₀ = 30r₀ (shown in Figure 1). In this study, neglecting the buoyancy part of seepage body force and compressibility of water, the governing equation for steady water seepage can be obtained [11]:

$$r\frac{{\partial }^2{p}_w\left(r\right)}{\partial r^2}+\frac{\partial p_w\left(r\right)}{\partial r^2}=0,$$

(2)

where r denotes radial distance from roadway center.

Based on seepage boundary conditions of p_w(r)|_r=r0 = 0 and p_w(r)|_r=R0 = p₀, the distribution of water pressure around roadway surrounding rock is given:

$$p_w\left(r\right)=s\ln \frac{r_0}{r},$$

(3)

where s = p₀/ln(r₀/R₀).

3.2. Elastoplastic Solution of Mechanical Model

3.2.1. Basic Equations under Polar Coordinate System

The stresses in circumferential and radial directions under polar coordinate system are denoted by σ_θ and σ_r, respectively. Under the assumption of (1), the strain in axial direction of roadway is zero, and then stress σ_z in axial direction can be expressed by σ_z = a(σ_r + σ_θ)/2, in which a is a constant. If the surrounding rock is in elastic zone, then a = 2v in which v represents Poisson’s ratio with value interval of [−1, 0.5]; while if the surrounding rock is in the two plastic zones, then a = 1 [13]. Hence, the relations among minor, intermediate and major principal stresses under polar coordinate system are σ₃ = σ_r ≤ σ₂ = σ_z ≤ σ₁ = σ_θ, and then intermediate principal stresses in surrounding rock satisfy the relation of σ₂ ≥ (σ₁ + σ₃)/2 − (σ₁ − σ₃)sin φ/2. Then unified strength theory of Equation (1) can be rewritten as an alternative form:

$${\sigma }_{\theta }-m{\sigma }_r-n=0,$$

(4)

where $m=\left[0.5b\left(1-\sin \phi \right)-\left(\sin \phi +1\right)\left(b+1\right)\right]/\left[\left(0.5b+1\right)\left(\sin \phi -1\right)\right]$, and $n=-2c\cos \phi \left(1+b\right)/\left[\left(1+0.5b\right)\left(\sin \phi -1\right)\right]$. When weighted coefficient b is zero, unified strength theory will be simplified as Mohr–Coulomb criterion.

According to elastoplastic mechanics and principle of effective stress [15,28], the equilibrium equation for surrounding rock of roadway under seepage effect is

$$\frac{{\sigma }_r-{\sigma }_{\theta }}{r}+\frac{d{\sigma }_r}{dr}+\eta \frac{d{p}_w}{dr}=0,$$

(5)

and the displacement-strain equations of surrounding rock is

$${\epsilon }_{\theta }=\frac{u}{r},{\epsilon }_r=\frac{du}{dr},$$

(6)

where η represents effective pore pressure coefficient; ε_θ and ε_r denote circumferential strain and radial strain, respectively; and u denotes displacement in radial direction.

3.2.2. Analytical Solutions of Displacement and Stress in Elastic Zone

The constitutive relationship of surrounding rock in elastic zone under plane strain condition is given by:

$$\{\begin{array}{c}{\sigma }_r^e=\frac{Ev}{\left(1-2v\right)\left(1+v\right)}\left({\epsilon }_{\theta }^e+\frac{1-v}{v}{\epsilon }_r^e\right)\\ {\sigma }_{\theta }^e=\frac{Ev}{\left(1+v\right)\left(1-2v\right)}\left(\frac{1-v}{v}{\epsilon }_{\theta }^e+{\epsilon }_r^e\right)\end{array},$$

(7)

where ${\sigma }_{\theta }^e$ and ${\sigma }_r^e$ denote circumferential and radial stresses of elastic zone; E is elastic modulus; and ${\epsilon }_{\theta }^e$ and ${\epsilon }_r^e$ are circumferential strain and radial strain in elastic zone, respectively.

Substituting Equations (6) and (7) into (5), the displacement and stress components in elastic zone are obtained:

$$u^e=\frac{C_0}{2}r\ln r+C_{1}r+\frac{C_2}{r},$$

(8)

$$\{\begin{array}{c}{\sigma }_r^e=\frac{C_{1}E}{\left(-2v+1\right)\left(v+1\right)}-\frac{C_2}{r^2}\frac{E}{v+1}+\frac{E}{\left(v+1\right)\left(1-2v\right)}\frac{C_0}{2}\left[\left(-v+1\right)+\ln r\right]\\ {\sigma }_{\theta }^e=\frac{E{C}_1}{\left(1-2v\right)\left(v+1\right)}+\frac{E}{r^2}\frac{C_2}{1+v}+\frac{E}{\left(1-2v\right)\left(1+v\right)}\frac{C_0}{2}\left(\ln r+v\right)\end{array},$$

(9)

where $C_0=\eta s\left(2v-1\right)\left(v+1\right)/\left(\left(v-1\right)E\right)$, and $C_1$, $C_2$ are constants given below in Equation (11).

At the plastic–elastic interface, r = R_p, the surrounding rock is on the verge of yielding, which satisfies yield condition of Equation (4); and the boundary condition of radial stress at r = R₀ is σ₀ + p₀, namely

$$\{\begin{array}{c}{{\sigma }_{\theta }^e|}_{r=R_p}=m{{\sigma }_r^e|}_{r=R_p}+n\\ {{\sigma }_r^e|}_{r=R_0}={\sigma }_0+p_0\end{array},$$

(10)

Substituting Equation (9) into Equation (10), the constants of $C_1$ and $C_2$ are given by

$$\{\begin{array}{c}{C}_1=\frac{\left(1+v\right)\left(1-2v\right)}{E}\left({\sigma }_0+p_0\right)+\frac{1-2v}{R_0^2}\times C_2-\frac{C_0}{2}\left[\left(1-v\right)+\ln R_0\right]\\ C_2=\frac{\left[n+\left({\sigma }_0+p_0\right)\left(m-1\right)+\frac{E{C}_0}{2\left(1-2v\right)\left(1+v\right)}\left[\left(m-1\right)\left(\ln \frac{R_p}{R_0}+1-2v\right)\right]\right]\left(v+1\right)}{E\left(\frac{1+m}{R_p^2}+\frac{1-m}{R_0^2}\right)}\end{array},$$

(11)

As the strains induced by initial stress σ₀ have been accomplished before excavation of roadway, so the actual strain components are given by [29,30]:

$$\{\begin{array}{c}{\epsilon }_r^e=\frac{ 1-v^2}{E}\left[\left({\sigma }_r^e-{\sigma }_0\right)-\frac{v}{1-v}\left({\sigma }_{\theta }^e-{\sigma }_0\right)\right]\\ {\epsilon }_{\theta }^e=\frac{1-v^2}{E}\left[\left({\sigma }_{\theta }^e-{\sigma }_0\right)-\frac{v}{1-v}\left({\sigma }_r^e-{\sigma }_0\right)\right]\end{array},$$

(12)

For elastic zone, the actual radial displacement u^e* is solved by substituting Equation (12) into left-hand side equation of Equation (6):

$$u^{e\ast }=r{C}_1+\frac{C_2}{r}+\frac{C_0}{2}r\ln r+\frac{{\sigma }_{0}r\left(2v^2+v-1\right)}{E},$$

(13)

3.2.3. Analytical Solutions of Displacement and Stress in Plastic Softening Zone

Related studies [18] demonstrated that surrounding rock has characteristics of strain softening after the peak stress is reached. Figure 2 shows that compressive strength of rock would gradually decrease to the residual strength after peak stress reaches. The residual strength represents the constant stress when the strain softening of rock would no longer occur after failure. The experimental strain–stress curve of rock [25,29] could be simplified and subdivided into a three-line segment as illustrated in Figure 2.

Generally, elastic strains occurring in plastic softening zone would be much smaller in contrast to plastic strains.

Hence, elastic strains in softening zone are assumed to the same as corresponding strains at the plastic–elastic interface [6]. For plastic softening zone, the radial and circumferential strains are expressed by the following [12,18]:

$$\{\begin{array}{c}{\epsilon }_{\theta }^p={\left({\epsilon }_{\theta }^e\right)}_{r=R_p}+\mathsf{\Delta }{\epsilon }_{\theta }^p\\ {\epsilon }_r^p={\left({\epsilon }_r^e\right)}_{r=R_p}+\mathsf{\Delta }{\epsilon }_r^p\end{array},$$

(14)

where ${\epsilon }_r^p$ and ${\epsilon }_{\theta }^p$ represent radial and circumferential strains; $\mathsf{\Delta }{\epsilon }_{\theta }^p$ and $\mathsf{\Delta }{\epsilon }_r^p$ represent increments of circumferential and radial plastic strains, respectively.

According to plastic potential theory [18], the increment of plastic strain is given by

$$d{\epsilon }_{ij}^p=d\lambda \frac{\partial \mathrm{g}}{\partial {\sigma }_{ij}},$$

(15)

where $d{\epsilon }_{ij}^p$ is increment of plastic strain; $d\lambda$ is the plastic factor increment; g is the plastic potential function; and ${\epsilon }_{ij}$ is the stress tensor. Assuming the plastic potential function g has the same form of unified strength theory of Equation (4):

$$\mathrm{g}={\sigma }_{\theta }-m^{\prime }{\sigma }_r-n^{\prime },$$

(16)

where $m^{\prime }$ is $\left[0.5b\left(1-\sin \psi \right)-\left(\sin \psi +1\right)\left(b+1\right)\right]/\left[\left(0.5b+1\right)\left(\sin \psi -1\right)\right]$, $n^{\prime }$ is $-2c\cos \psi \left(1+b\right)/\left(\left(1+0.5b\right)\left(\sin \psi -1\right)\right)$, and $\psi$ is dilatancy angle.

Substituting (16) into (15), the plastic strain increment could be expressed by

$$\{\begin{array}{c}d{\epsilon }_{\theta }^p=d\lambda \frac{\partial \mathrm{g}}{\partial {\sigma }_{\theta }}=d\lambda \\ d{\epsilon }_r^p=d\lambda \frac{\partial \mathrm{g}}{\partial {\sigma }_r}=-m^{\prime }d\lambda \end{array},$$

(17)

Existing studies [3,5] have demonstrated that there is little volume change of rock in the pre-peak process, while there is a volumetric dilatancy of rock at the post-peak process. Therefore, rock dilatancy occurs both at broken zone and plastic softening zone. Relationships between strains and dilatancy coefficients are illustrated in Figure 3.

Based on flow rule of volumetric dilatancy in plastic softening zone [12,18]:

$$\mathsf{\Delta }{\epsilon }_r^p+{\alpha }_1\mathsf{\Delta }{\epsilon }_{\theta }^p=0,$$

(18)

where α₁ (α₁ ≥ 1) denotes dilatancy coefficient, which is slope between increments of radial and circumferential plastic strains shown in Figure 3. When α₁ = 1, it represents a condition of surrounding rock without volumetric dilatancy.

Combined with Equations (17) and (18), the dilatancy coefficient is obtained by plastic potential function:

$${\alpha }_1=-\frac{\mathsf{\Delta }{\epsilon }_r^p}{\mathsf{\Delta }{\epsilon }_{\theta }^p}=-\frac{d{\epsilon }_r^p}{d{\epsilon }_{\theta }^p}=m^{\prime },$$

(19)

Similarly, for broken zone, the circumferential and radial strains are expressed as [12,18]:

$$\{\begin{array}{c}{\epsilon }_{\theta }^b=\mathsf{\Delta }{\epsilon }_{\theta }^b+{\left({\epsilon }_{\theta }^p\right)}_{r=R_b}\\ {\epsilon }_r^b=\mathsf{\Delta }{\epsilon }_r^b+{\left({\epsilon }_r^p\right)}_{r=R_b}\end{array},$$

(20)

where ${\epsilon }_{\theta }^b$ and ${\epsilon }_r^b$ represent the total circumferential and radial strains; $\mathsf{\Delta }{\epsilon }_{\theta }^b$ and $\mathsf{\Delta }{\epsilon }_r^b$ represent increments of radial and circumferential plastic strains, respectively.

Based on flow rule of volumetric dilatancy in the broken zone [12,18],

$$\mathsf{\Delta }{\epsilon }_r^b+{\alpha }_2\mathsf{\Delta }{\epsilon }_{\theta }^b=0,$$

(21)

where α₂ (α₂ ≥ 1) denotes dilatancy coefficient, which is slope between increments of radial and circumferential plastic strains in broken zone as illustrated in Figure 3, and it can be calculated similarly to Equation (19). When α₂ = 1, it means that no volumetric dilatancy occurs in the broken zone.

For the plastic softening zone, substituting Equation (14) into Equation (6), we obtain

$$\{\begin{array}{c}\frac{u^{p\ast }}{r}={\epsilon }_{\theta }^p={\left({\epsilon }_{\theta }^e\right)}_{r=R_p}+\mathsf{\Delta }{\epsilon }_{\theta }^p\\ \frac{d{u}^{p\ast }}{dr}={\epsilon }_r^p={\left({\epsilon }_r^e\right)}_{r=R_p}+\mathsf{\Delta }{\epsilon }_r^p\end{array},$$

(22)

Multiplying the first equation in Equation (22) by ${\alpha }_1$, which is then added to the second Equation in Equation (22), the following expression are obtained:

$${\alpha }_1\frac{u^{p\ast }}{r}+\frac{d{u}^{p\ast }}{dr}={\alpha }_1{\left({\epsilon }_{\theta }^e\right)}_{r=R_p}+{\left({\epsilon }_r^e\right)}_{r=R_p}+\mathsf{\Delta }{\epsilon }_r^p+{\alpha }_1\mathsf{\Delta }{\epsilon }_{\theta }^p,$$

(23)

According to Equation (18), Equation (23) can be simplified as

$${\alpha }_1\frac{u^{p\ast }}{r}+\frac{d{u}^{p\ast }}{dr}={\alpha }_1{\left({\epsilon }_{\theta }^e\right)}_{r=R_p}+{\left({\epsilon }_r^e\right)}_{r=R_p},$$

(24)

where ${\left({\epsilon }_{\theta }^e\right)}_{r=R_p}$ and ${\left({\epsilon }_r^e\right)}_{r=R_p}$ are circumferential and radial strains at the interface of plastic softening zone and elastic zone, which can be obtained by solving Equation (12) with r being R_p.

Combining with the continuity condition for displacement u^p* = u^e* at r = R_p (see Equation (13)), we obtain the displacement u^p* and circumferential strain ${\epsilon }_{\theta }^p$ through solving ordinary differential Equations (24) and (6):

$$u^{p\ast }=\left(\frac{1}{2}{C}_0\ln R_p+\frac{C_0}{2\left({\alpha }_1+1\right)}+C_1+\frac{{\sigma }_0\left(2v^2+v-1\right)}{E}+\frac{C_2\left({\alpha }_1-1\right)}{\left({\alpha }_1+1\right)R_p^2}\right)r+\frac{R_p^{{\alpha }_1+1}}{\left({\alpha }_1+1\right)r^{{\alpha }_1}}\left(\frac{2C_2}{R_p^2}-\frac{C_0}{2}\right),$$

(25)

$${\epsilon }_{\theta }^p=\frac{1}{2}{C}_0\ln R_p+\frac{C_0}{2\left({\alpha }_1+1\right)}+C_1+\frac{{\sigma }_0\left(2v^2+v-1\right)}{E}+\frac{C_2\left({\alpha }_1-1\right)}{\left({\alpha }_1+1\right)R_p^2}+\frac{R_p^{{\alpha }_1+1}}{\left(1+{\alpha }_1\right)r^{1+{\alpha }_1}}\left(\frac{2C_2}{R_p^2}-\frac{C_0}{2}\right),$$

(26)

The rock within plastic softening zone is yielding, which satisfies the following relation:

$${\sigma }_{\theta }^p-m{\sigma }_{\theta }^r-n^p=0,$$

(27)

where ${\sigma }_{\theta }^p$ and ${\sigma }_r^p$ denote circumferential stress and radial stress, respectively; n^p denotes compressive strength of surrounding rock as illustrated in Figure 4.

In Figure 4, S denotes the compressive strength parameter of the last term on the left side in Equation (4), M = tan θ represents the softening modulus of rock. For rock in strain softening region, its compressive strength parameter of cohesion would be significantly influenced, while internal friction angle is minimally affected [25], which means that the internal friction angle can be assumed as a constant. Therefore, the evolution model for compressive strength parameter S, i.e., the evolution of cohesion, is subdivided into three-line segment (Figure 4). The greater of plastic strain increment, the smaller is n^p. The compressive strength of surrounding rock, n^p, is given by

$$n^p=-M\mathsf{\Delta }{\epsilon }_{\theta }^p+n=-M\left[{\epsilon }_{\theta }^p-{\left({\epsilon }_{\theta }^e\right)}_{r=R_p}\right]+n,$$

(28)

Substituting Equations (12) and (26) into (28) gives

$$n^p=n-\frac{M}{1+{\alpha }_1}\left(\frac{2C_2}{R_p^2}-\frac{C_0}{2}\right)\left({\left(\frac{R_p}{r}\right)}^{{\alpha }_1+1}-1\right),$$

(29)

Furthermore, substituting Equation (22) into (20) leads to

$${\sigma }_{\theta }^p=m{\sigma }_r^p+n-\frac{M}{1+{\alpha }_1}\left(\frac{2C_2}{R_p^2}-\frac{C_0}{2}\right)\left({\left(\frac{R_p}{r}\right)}^{{\alpha }_1+1}-1\right),$$

(30)

In combination with continuity condition of stress ${\sigma }_r^e$ = ${\sigma }_r^p$ at r = R_p, and substituting Equation (30) into Equation (5), the stress components in plastic softening zone are expressed as

$${\sigma }_r^p=\frac{n+\eta s}{1-m}+\frac{M}{1+{\alpha }_1}\left(-\frac{C_0}{2}+\frac{2C_2}{R_p^2}\right)\left(\frac{1}{1-m}+\frac{R_p^{{\alpha }_1+1}}{\left(m+{\alpha }_1\right)r^{{\alpha }_1+1}}\right)+E_1{\left(\frac{r}{R_p}\right)}^{m-1},$$

(31)

$${\sigma }_{\theta }^p=m\left(\frac{n+\eta s}{1-m}+\frac{M}{1+{\alpha }_1}\left(-\frac{C_0}{2}+\frac{2C_2}{R_p^2}\right)\left(\frac{1}{1-m}+\frac{R_p^{{\alpha }_1+1}}{\left(m+{\alpha }_1\right)r^{{\alpha }_1+1}}\right)+E_1{\left(\frac{r}{R_p}\right)}^{m-1}\right)+n-\frac{M}{1+{\alpha }_1}\left(\frac{2C_2}{R_p^2}-\frac{C_0}{2}\right)\left({\left(\frac{R_p}{r}\right)}^{{\alpha }_1+1}-1\right),$$

(32)

$$E_1=\frac{C_{0}E}{2\left(1-2v\right)\left(1+v\right)}\left(1-v+\ln R_p\right)+\frac{C_{1}E}{\left(1-2v\right)\left(1+v\right)}-\frac{E}{1+v}\frac{C_2}{R_p^2}-\frac{n+\eta s}{1-m}-\frac{M}{\left(1-m\right)\left(m+{\alpha }_1\right)}\left(\frac{2C_2}{R_p^2}-\frac{C_0}{2}\right),$$

(33)

3.2.4. Analytical Solutions for Displacement and Stress in Broken Zone

For surrounding rock of broken zone, it satisfies unified strength theory:

$${\sigma }_{\theta }^b-m{\sigma }_r^b-n^{\ast }=0,$$

(34)

where ${\sigma }_{\theta }^b$ and ${\sigma }_r^b$ denote circumferential and radial stresses, respectively; n* denotes the residual strength.

Combining with the stress condition ${\sigma }_r^b$|_r=r0 = p_i, and substituting Equation (30) into Equation (5), stress components in broken zone are then given as:

$${\sigma }_r^b={\left(\frac{r}{r_0}\right)}^{m-1}\left(p_i-\frac{n^{\ast }+s\eta }{1-m}\right)+\frac{\eta s+n^{\ast }}{1-m},$$

(35)

$${\sigma }_{\theta }^b=\left[{\left(\frac{r}{r_0}\right)}^{m-1}\left(p_i-\frac{\eta s+n^{\ast }}{1-m}\right)+\frac{\eta s+n^{\ast }}{1-m}\right]m+n^{\ast },$$

(36)

Similarly, combining with Equations (6), (20) and (21) and the displacement continuity condition of $u^{p\ast }=u^{b\ast }$ at $r=R_b$, then $u^{b\ast }$ is obtained:

$$u^{b\ast }=Ar+\left(\frac{{\alpha }_2-{\alpha }_1}{{\alpha }_2+1}\right)Br+\left(\frac{{\alpha }_1+1}{{\alpha }_2+1}\right)B\frac{R_b^{{\alpha }_2+1}}{r^{{\alpha }_2}},$$

(37)

where A and B are constants given by

$$A=\frac{1}{2}{C}_0\ln R_p+\frac{C_0}{2\left({\alpha }_1+1\right)}+C_1+\frac{{\sigma }_0\left(2v^2+v-1\right)}{E}+\frac{C_2\left({\alpha }_1-1\right)}{\left({\alpha }_1+1\right)R_p^2},$$

(38)

$$B=\frac{1}{\left({\alpha }_1+1\right)}{\left(\frac{R_p}{R_b}\right)}^{{\alpha }_1+1}\left(\frac{2C_2}{R_p^2}-\frac{C_0}{2}\right),$$

(39)

3.2.5. Radii of Plastic Zones

According to the interface condition n^p = n* at r = R_b and Equation (29), we obtain

$$n^{\ast }=n-\frac{M}{1+{\alpha }_1}\left(\frac{2C_2}{R_p^2}-\frac{C_0}{2}\right)\left({\left(\frac{R_p}{R_b}\right)}^{{\alpha }_1+1}-1\right),$$

(40)

Then relationship between softening zone radius and broken zone radius can be given by

$$R_b=R_p{\left(1+\frac{\left(1+{\alpha }_1\right)\left(n-n^{\ast }\right)}{M\left(\frac{2C_2}{R_p^2}-\frac{C_0}{2}\right)}\right)}^{\frac{-1}{1+{\alpha }_1}},$$

(41)

Utilizing the stress continuity condition ${\sigma }_r^p$ = ${\sigma }_r^b$ at r = R_b, the radius of softening zone is be given by

$${\left(\frac{R_p}{r_0}\right)}^{m-1}=\left[\frac{-{\left[\left(1+{\alpha }_1\right)\left(n-n^{\ast }\right)+M\left(\frac{2C_2}{R_p^2}-\frac{C_0}{2}\right)\right]}^{\frac{{\alpha }_1+m}{{\alpha }_1+1}}}{\left(m-1\right)\left({\alpha }_1+m\right){\left[M\left(\frac{2C_2}{R_p^2}-\frac{C_0}{2}\right)\right]}^{\frac{m-1}{1+{\alpha }_1}}}+E_1\right]/\left(p_i+\frac{\eta s+n^{\ast }}{m-1}\right),$$

(42)

3.3. Validation of Unified Analytical Solution

The proposed analytical solution for surrounding rock of roadway has taken influences of water seepage, strain softening, dilatancy and intermediate principal stress all into account. Therefore, it can be deemed as a unified solution for an ordered series of solutions. Obviously, the series of solutions obtained should involve analytic solution for elastic and perfectly plastic surrounding rock. The simplification process of unified solution is as follows: $\eta$ = 0 and p₀ = 0, the water seepage influence is not regarded; $n^{\ast }=n$, the strain softening of surrounding rock is ignored; ${\alpha }_1={\alpha }_2=1$, the rock dilatancy is not considered; b = 0, the unified strength theory is simplified as Mohr–Coulomb criterion. For the analytic solution of plastic zone radius (i.e., strain-softening zone radius or broken zone radius in present study, as the two radii are the same under such condition), substituting the above parameters into Equation (42), and with combination of Equations (33), (11), and (4), the plastic zone radius is given:

$$R_p=r_0{\left[\frac{\left(1-sin\phi \right)\left(c·cot\phi +{\sigma }_0\right)}{c·cot\phi +p_i}\right]}^{\frac{1-\sin \phi }{2\sin \phi }},$$

(43)

Let $r=r_0$, substituting the above parameters into Equation (37), and with combination of Equations (38) and (39), then surface displacement of roadway can be obtained:

$$\begin{gathered} u^{b\ast }=\frac{\left(1+v\right)sin\phi }{r_{0}E}\left(c·cot\phi +{\sigma }_0\right)R_{p }^2 \\ =\frac{r_0\left(1+v\right)sin\phi }{E}\left(c·cot\phi +{\sigma }_0\right){\left[\frac{\left(1-sin\phi \right)\left(c·cot\phi +{\sigma }_0\right)}{c·cot\phi +p_i}\right]}^{\frac{1-\sin \phi }{\sin \phi }}, \end{gathered}$$

(44)

As can be seen, both Equations (43) and (44) are consistent with the classical Kastner’s solution in literature [30], which validates the correctness of the unified solution for surrounding rock of roadway. Thus, the unified solution obtained can indeed represent an ordered series of solutions.

4. Example Analysis

To further investigate influences of water seepage, strain softening, rock dilatancy and intermediate principal stress on radii of plastic zone, displacement and stress in surrounding rock of roadway, an engineering example is thoroughly studied using the single factor analysis. Namely, one of the influencing factors, such as pore water pressure, softening modulus, dilatancy coefficient and weighted coefficient are changed each time. The radius of the circular roadway r₀ is 2 m with uniform original stress σ₀ being 15 MPa.

Table 1. Calculation parameters for mechanical model.

Parameters	Values	Parameters	Values
Poisson’s ratio v	0.25	Dilatancy coefficient in plastic softening zone α1	2
Elastic modulus E/MPa	2000	Dilatancy coefficient in broken zone α2	1.5
Residual Cohesion c*/MPa	1.0	Softening modulus M/MPa	2000
Original Cohesion c/MPa	3.0	Internal friction angle φ/°	30
Weighted coefficient b	0.5	Initial pore water pressure p0/MPa	2

lists the calculation parameters for mechanical model of surrounding rock of roadway.

4.1. Influence of Water Seepage

Figure 5a illustrates the distributions of stresses with water seepage considered (η = 1) and without water seepage considered (η = 0). As can be seen, variation laws of stresses under both conditions are coincident with each other; namely, the circumferential stresses present tendency of increasing firstly and decreasing next, radial stresses present a gradually increasing trend, and peak circumferential stresses occur at plastic–elastic interface; however, water seepage plays a significant role in distributions of two stress components; the peak circumferential stress with water seepage considered is greater than that of without water seepage considered, and location of peak stress would shift to the righthand side; with initial pore water pressure increasing, the peak circumferential stresses and its corresponding distances from roadway center would both increase. The results obtained agree well with existing study in literatures [22,25]. When p₀ is increased from 1 MPa to 3 MPa, peak circumferential stress increases by 18% from 28.7 MPa to 33.8 MPa. Figure 5b,c shows variations of radii of plastic zones and of roadway surface displacements under different water pressure. As can be seen, radii of broken and plastic softening zones and surface displacement with water seepage considered are greater than that of without water seepage considered. The reason is that effective stress in surrounding rock would increase with seepage effect considered [2,6]; moreover, radii of two plastic zones and surface displacement of roadway all present an increase trend as water pressure increases with water seepage considered and without water seepage considered although their tendencies are different. The former increases exponentially with water pressure, while the latter increases linearly with water pressure. For instance, with p₀ increasing from 1 MPa to 3 MPa, plastic softening zone radius increases by 22% from 3.2 m to 3.9 m, surface displacement of roadway increases by 80% from 0.05 m to 0.09 m. In short, the obtained results indicate that water seepage plays a significant role in radii of plastic zones and surface displacement of roadway, and its influence should be carefully considered for the stability of roadway under groundwater occurrence environment.

4.2. Influence of Strain Softening

Three softening modulus M, i.e., 1000 MPa, 2000 MPa and 3000 MPa, are chosen to investigate effect of strain softening on stress distribution. As illustrated in Figure 6a, with increase of softening modulus, the magnitude of peak circumferential stress is kept unchanged, which is consistent with analysis in literature [30] indicating that stresses at interface of elastic zone and plastic softening zone would not be influenced by strain softening in plastic zones. Besides, location of peak circumferential stress would slightly shift to the righthand side when softening modulus increases from 1000 MPa to 3000 MPa; and the greater of softening modulus, the smaller of two stress components in plastic zones at the same distance from roadway center. The reason is that the softening modulus denotes the compressive strength reduction per unit of plastic strain increment, the greater of softening modulus, the greater of strength reduction per unit of plastic strain increment.

Furthermore, softening modulus M has a greater influence on R_b but less influence on radius of strain-softening zone as illustrated in Figure 6b. With the increase in M, R_b would gradually expand with expansion velocity decreasing. For instance, when softening modulus M is increased from 1000 MPa to 2000 MPa, radius of broken zone increases by 21% from 2.4 m to 2.9 m; while softening modulus M is increased from 2000 MPa to 3000 MPa, radius of broken zone increases by 7% from 2.9 m to 3.1 m. Ratio of R_p and R_b decreases in a decay form with softening modulus increasing, which further indicates that softening modulus has a greater influence on R_b than R_p. The variation curve between softening modulus and surface displacement of roadway. As seen from Figure 6c, when softening modulus M is increased from 1000 MPa to 3000 MPa, surface displacement of roadway only increases from 51.8 mm to 62.3 mm with its slope between softening modulus, and surface displacement of roadway is decreasing gradually, which showed that softening modulus M would affect surface displacement of roadway to an extent. In short, the obtained results indicate that strain-softening characteristics of surrounding rock should be analyzed to accurately predict broken zone radius and surface displacement of roadway; for bolting support roadway, the length of rock bolt is usually determined by the radius of plastic zones [1], so the length of rock bolt should be increased with softening modulus increasing.

4.3. Influence of Dilatancy

The dilatancy characteristics of rock are described by dilatancy coefficient α₁ and α₂ in this paper (see Figure 3). As described in theoretical analysis of Section 3, dilatancy coefficient α₂ has no influence on radii of broken and softening zones. Noting that when rock dilatancy coefficient is 3.4, it is corresponding to associated flow rule (i.e., internal friction angle is equal to dilatancy angle); while rock dilatancy coefficient is set to 1, it is corresponding to the condition without dilatancy considered (i.e., dilatancy angle is set to zero).

Figure 7a illustrates the distributions of stresses around roadway. As can be seen, the rock dilatancy has no effect on magnitude of peak circumferential stress, which is consistent with the results in the literature [22,30]. Figure 7b shows that the dilatancy coefficient has little influence on strain-softening zone radius, but the greater of dilatancy coefficient α₁, the greater of radius of broken zone, which agrees well with the related literature [18]. The broken zone radius increases linearly with dilatancy coefficient α₁ increasing. Dilatancy coefficients play an important role in surface displacements of roadway (Figure 7c). As can be seen, surface displacements of roadway increase linearly with dilatancy coefficient α₁. When dilatancy coefficient α₁ is increased from 1 to 3.4, surface displacement of roadway increases by 32.4% from 50.0 mm to 66.2 mm under condition of α₂ = 1; while surface displacement of roadway increases by 37.8% from 52.6 mm to 72.5 mm under condition of α₂ = 1.6, the increment of surface displacement of roadway is about 7.52 mm with dilatancy coefficient increasing by one. Relevant studies showed that the general value of dilatancy coefficient is 1.5 [22]. Hence, the radius of broken zone and surface displacement of roadway is underestimated with dilatancy not considered, while it is overestimated using the associated flow rule of rock.

4.4. Influence of Intermediate Principal Stress

As mentioned in Section 2, intermediate principal stress effect is manifested by weighted coefficient b. Figure 8a illustrates the distribution of stresses under different weighted coefficient.

As can be seen, the greater the weighted coefficient, the larger the two stress components in plastic zones at the same distance from roadway center and corresponding peak circumferential stress; however, with a smaller weighted coefficient, the circumferential stress is increased while radial stress is decreased, which is consistent with results in the literature [17]. Moreover, the location of peak circumferential stress would move to the lefthand side with weighted coefficient increasing. The farther away from roadway center, the less influence of weighted coefficient on stresses. It can be explained that the compressive strength of rock would be increased with effect of intermediate principal stress considered as shown in unified strength theory of Equation (4). Figure 8b shows the relation between weighted coefficient and radii of plastic zones. Both radii of broken and plastic softening zones are reduced exponentially with weighted coefficient increasing. When weighted coefficient is increased from 0 to 1, peak circumferential stress increases by 6% from 29.5 MPa to 31.2 MPa, radius of plastic softening zone decreased by 28% from 4.4 m to 3.15 m, surface displacement of roadway decreased by 48% from 92 mm to 48 mm. Results demonstrate that limit strength of surrounding rock has been enhanced with influence of intermediate principal stress considered. Noting when b is set to zero, unified strength criterion is equivalent to Mohr–Coulomb criterion, which is more conservative in predicting radii of plastic zones and surface displacement of roadway. Therefore, in engineering practice, self-bearing capacity enhancement of surrounding resulting from intermediate principal stress must be considered. Due to length of rock bolt is mainly determined by plastic zone radius, so rock bolt length can be decreased with effect of intermediate principal stress considered. If grouting measure is adopted, the grouting range can also be decreased. The results obtained are helpful for reasonable measure selection of surrounding rock, which can save support cost.

4.5. Influence of Residual Cohesion

Figure 9a presents the distribution of stresses under different residual cohesion c*. As can be seen, the greater the residual cohesion, the larger the two stress components in plastic zones at the same distance from roadway center, and the smaller the distance of plastic–elastic interface from roadway center. However, the peak circumferential stress remains unchanged, which is in good agreement with the literature [30].

Variations of radii of plastic zones under different residual cohesions are illustrated in Figure 9b. The broken zone and softening zone decreased exponentially with the increase of residual cohesion. For instance, the radius of broken zone decreases by 20% from 3.3 m to 2.65 m with residual cohesion changes from 0.6 MPa to 2.0 MPa. Hence, the grouting measure, which will improve residual cohesion, can be adopted in roadway to ensure its stability.

5. Conclusions

In this paper, a mechanical model which can simultaneously consider water seepage, strain softening, rock dilatancy, and intermediate principal stress was established, and then unified analytical solutions for pore water pressure, stress, displacement, radii of plastic zones were obtained for surrounding rock of roadway. Based on an example, influences of water seepage, strain softening, rock dilatancy, residual cohesion and intermediate principal stress on surrounding rock of roadway were thoroughly investigated using single factor analysis. The main conclusions are as follows:

(1)

With water pressure increasing, peak circumferential stress, radii of two plastic zones and surface displacement of roadway with water seepage considered would increase, and their magnitudes are all greater than their corresponding conditions with water seepage not considered. Particularly, radii of plastic zones and surface displacement of roadway increases exponentially with water pressure increasing, water seepage influence should be carefully considered to ensure roadway stability under groundwater environment.

(2)

With softening modulus increasing, the magnitude of peak circumferential stress is kept unchanged, and location of peak circumferential stress would slightly shift to the deeper surrounding rock. Softening modulus has a greater influence on broken zone radius than that of strain-softening zone radius, and it also affects surface displacement of roadway to an extent. To accurately predict broken zone radius and surface displacement of roadway, strain-softening characteristics of surrounding rock should be analyzed; and length of rock bolt should be increased with softening modulus increasing for bolting support roadway.

(3)

Rock dilatancy coefficient has little effect on magnitude of peak circumferential stress and plastic softening zone radius, while radius of broken zone and surface displacement of roadway increase linearly with dilatancy coefficient α1 increasing. For actual engineering, surface displacement of roadway would be underestimated with rock dilatancy not considered, while they are overestimated if associated flow rule is adopted, indicating that flow rule of rock should be reasonably chosen for calculating roadway surface displacement and broken zone radius of roadway.

(4)

With weighted coefficient increasing, stress components in plastic zones at the same distance from roadway center would increase; while the distance of peak circumferential stress location from roadway center, radii of broken and plastic softening zones, and surface displacement of roadway are reduced. Self-bearing capacity enhancement of surrounding rock resulting from effect of intermediate principal stress should be considered. The traditional Mohr–Coulomb criterion is more conservative. Rock bolt length and grouting range can be decreased after considering effect of inter-mediate principal stress, which is helpful for reasonable measure selection of surrounding rock.

(5)

With residual cohesion increasing, peak circumferential stress remains unchanged, stress components in plastic zones at the same distance from roadway center would increase, and the distance of plastic–elastic interface from roadway center decreases, which implicates that grouting measure can be adopted to improve roadway stability effectively.