*Article* **Incompatible Deformation Model of Rocks with Defects around a Thick-Walled Cylinder**

**Yingji Bao 1,2 and Binsong Jiang 1,\***


**Abstract:** Before the excavation of underground engineering, joints, fissures, and voids already exist in the rock—that is, there are defects in the rock. Due to the existence of these defects, the rock produces plastic deformation, which can lead to incompatible deformation. Therefore, the classic continuum theory cannot accurately describe the deformation of the rock. In this paper, a relationship between the strain tensor and metric tensor was studied by analyzing the three states of elastic plastic deformation, and the elasto-plastic incompatible model was built. Additionally, the stress and deformation of a thick-walled cylinder under hydrostatic pressure was investigated by using a finite element program written in the FORTRAN language. The results show that the plastic strain is associated with not only deviator stress but also the distribution of defects (represented by the incompatible parameter *R*). With the value of *R* increasing, the defects in the rock increased, but the elastic plastic stiffness matrix decreased. Thus, as more rock enters the plastic state, the deformation of the surrounding rock is enlarged.

**Keywords:** incompatibility; rock; plastic deformation; finite element

#### **1. Introduction**

The size and shape of an object will change under external forces, namely the deformation of objects. After the external force is removed, the part of the deformation that has disappeared is called elastic deformation, and the remaining deformation is plastic deformation. The important characteristic of elastic deformation is reversibility—that is, deformation occurs after the force is applied and disappears after the force is removed. This indicates that the elastic deformation is determined by the binding force between atoms. If an external force overcomes the gravitational pull between the atoms and pulls them apart, the object will break, and its strength is called breaking strength. The actual object contains defects such as dislocations and disclinations; with a small elastic deformation, the stress is capable of activating the dislocation and making it move, resulting in plastic deformation. For brittle materials, as they are more sensitive to stress concentration, when the stress is slightly higher, the concentrated stress can lead to the movement and proliferation of these defects, resulting in incompatible deformation. Therefore, the classical continuum theory is not suitable for analyzing the defected rocks, and the incompatible deformation theory can be used.

The incompatible deformation theory originates from the defect theory and has been studied by many researchers [1–3]. Kondo first introduced differential geometry into the defect theory. He established the relationship between the dislocation theory and non-Riemannian geometry [4]. Later on, Anthony pointed out the relationship between differential geometry and disclinations [5]. Then, Kroner completed the underlying theory for defects and differential geometry including dislocations and extra-matter [6]. These

**Citation:** Bao, Y.; Jiang, B. Incompatible Deformation Model of Rocks with Defects around a Thick-Walled Cylinder. *Processes* **2021**, *9*, 2215. https://doi.org/10.3390/ pr9122215

Academic Editors: Li Li and Haiping Zhu

Received: 21 October 2021 Accepted: 3 December 2021 Published: 8 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

investigations were conducted to establish a relation between the parameters of non-Euclidean geometry and the elastic strain, bend-twist, and quasi-plastic strain of the defect theory [7]. M.A. Guzev first introduced the incompatible deformation theory to geotechnical engineering and assumed that the deformed body is in a Riemannian space. The curvature tensor *R* is not zero, so it is treated as an independent thermodynamic variable. As a result, a Riemann continuous model is established by using unbalanced thermodynamics. The stress distribution of a deep circle tunnel under plane strain was analyzed [8–13]. Zhou Xiaoping and Qian Qihu considered the impact of damage on a rock, described the damage degree of surrounding rocks in a deep circular tunnel according to the damage variables defined by the energy equivalence of damage mechanics, and analyzed the relationship between the damage of a deep rock and zonal fracture [14–16]. From the perspective of differential geometry, these existing models in geotechnical engineering all assumed that the undeformed body is in a three-dimensional Euclidean space, or the deformed body is a manifold. However, there are defects such as joints and cracks in the rock before excavation; therefore, the initial configuration should also be regarded as a manifold. In this paper, three configurations are constructed to describe the elastoplastic deformation of rocks in differential geometry, and the relationship between the strain tensor and interior metric tensor is established. Accordingly, an elastoplastic incompatible model is constructed. Finally, the elastoplastic finite element program is written in the FORTRAN language, and the stress and deformation of a circular tunnel under hydrostatic pressure is also analyzed.

#### **2. The Incompatible Plastic Deformation Theory of Rocks**

For rocks with elastoplastic deformation, we need to analyze three configurations. As shown in Figure 1 below, in the initial configuration Ω0, the rock is undeformed, and *X*, the micro element at any point, is denoted as *dX*. The deformed rock is the current configuration Ω, which has undergone much plastic deformation. The micro element at any point, namely *x* of the deformed rock, is denoted as *dx*. In addition, the intermediate configuration Ω\* is also introduced. The intermediate configuration is assumed to be obtained from the initial configuration through pure plastic deformation or from the current configuration through the process of elastic unloading, where the micro element at any point *ζ* is denoted as *dζ*.

**Figure 1.** Three configurations of plasticity deformation.

A point *X* on Ω<sup>0</sup> is mapped to a point *ζ* on the intermediate configuration Ω\* under the mapping *χp*:

$$
\chi^p: X \to \mathbb{Q} \tag{1}
$$

Due to the presence of these defects, the initial configuration is a material of manifold. In local coordinates, the neighborhood of point *X* has coordinates (*X*1, ... , *X*3); the

neighborhood of point *ζ* has coordinates (*ζ*1, ... , *ζn*); *χ<sup>p</sup>* is a continuous mapping; and the transformation matrix between *X* and *ζ* is *∂ζ*/*∂X*.

In differential geometry, *χ<sup>p</sup>* is an immersion mapping of manifold *M* to manifold *L*, and the intermediate configuration is assumed to be a Riemannian manifold with a metric structure of a symmetric nonsingular covariant tensor field of the second order. In an analysis of local properties, it can be considered that *M* is embedded in *L*, inducing the push from the tangent space:

$$\chi\_{\ast}^{p}: \eta^{I}(X)\exists\_{I} \to \eta^{I}(X)A\_{I}^{a}\partial\_{a} = \overline{\eta}^{a}\partial\_{a} \tag{2}$$

A metric is defined in a Riemannian manifold as

$$ds\_\*^2 = \mathcal{g}\_{a\beta} d\zeta^a d\zeta^\beta \,\tag{3}$$

The metric tensor field *gαβ* on *L* is dragged back to the metric field *GIJ* on *M* by

$$ds\_\*^2 = \mathcal{g}\_{a\beta} A^a\_I A^\beta\_{I\dot{J}} dX^I dX^J = \mathcal{G}\_{IJ} dX^I dX^J \tag{4}$$

The metric tensor *g* of manifold *L* is determined by the distribution of rock defects and the degree of plastic deformation, while the metric tensor *G* of manifold *M* is determined by the distribution of these rock defects.

The element *dX* at a point *X* in the initial configuration can be related to the element *dζ* at a point *ζ* in the intermediate configuration through nonholonomic transformation:

$$d\mathbb{Z} = A dX\tag{5}$$

*A* is the nonholonomic coefficient. The components satisfy

$$d\mathbb{S}^{a} = A\_{I}^{a}dX^{I}, \mathbf{e}\_{a} = A\_{I}^{a}\mathbf{e}\_{I} \tag{6}$$

Point *ζ* on Ω\* is mapped to point *x* on Ω under mapping *χ<sup>e</sup>* :

$$\chi^{\mathfrak{x}} : \mathscr{J} \to \mathbf{x} \tag{7}$$

There is a tangent mapping at *ζ*:

$$
\psi\_\* : T\_{\tilde{\zeta}}(\Omega\_\*) \to T\_x(\Omega) \tag{8}
$$

The local coordinate of point *ζ* is (*ζ*1, ... , *ζn*), the corresponding coordinate of point *x* is (*x*1, ... , *x*3), *χ<sup>e</sup>* is a continuous mapping, and the transformation matrix between *ζ* and *x* is *∂x*/*∂ζ*. The metric tensor of the current configuration *N* is determined by the metric tensor *G* of the intermediate configuration and the elastic deformation *ε<sup>e</sup>* .

Because the elastic deformation is reversible, we have

$$d\mathfrak{x} = Fd\mathfrak{y} \tag{9}$$

*F* is the deformation gradient, and its components satisfy

$$d\boldsymbol{x}^{i} = F\_{\boldsymbol{a}}^{i} d\zeta^{\boldsymbol{a}}, \boldsymbol{e}\_{i} = F\_{i}^{\boldsymbol{a}} \boldsymbol{e}\_{\boldsymbol{a}} \tag{10}$$

The deformation of rocks can be considered as occurrences in a three-dimensional Euclidean space in a small range. The total deformation from the initial configuration to the current configuration is

$$
\boldsymbol{\chi} = \boldsymbol{\chi}^c \cdot \boldsymbol{\chi}^p \tag{11}
$$

Since the rock is in a three-dimensional Euclidean space before and after the deformation, the total deformation is compatible, so *χ* is bijective, and according to the composite mapping theorem, if *χ<sup>p</sup>* is injective, then *χ<sup>e</sup>* is surjective.

The initial configuration is in a three-dimensional Euclidean space, and the length square of the element at point *X* is

$$ds\_0^2 = \delta\_{I\bar{I}} dX\_I X\_{\bar{I}} \tag{12}$$

The intermediate configuration is a Riemannian manifold, whose element length squared is

$$ds\_\*^2 = \mathcal{g}\_{a\beta} d\zeta^a d\zeta^\beta = \mathcal{g}\_{a\beta} A^a\_I A^\beta\_{I} dX^I dX^I \tag{13}$$

The current configuration is also in a three-dimensional Euclidean space, and its element length squared is

$$ds^2 = \delta\_{i\bar{j}} dx^i dx^{\bar{j}} = \delta\_{i\bar{j}} F^i\_I F^{\bar{j}}\_{\bar{I}} dX^I dX^{\bar{I}} = F^K\_I F^K\_{\bar{I}} dX^I dX^{\bar{I}} \tag{14}$$

Therefore, the change in the squared length of the element from the intermediate configuration to the current configuration is

$$ds^2 - ds\_\*^2 = (F\_I^K F\_I^K - g\_{IJ})dX^I dX^J = 2\varepsilon\_{IJ}^c dX^I dX^J \tag{15}$$

where *ε<sup>e</sup> I J* is the elastic strain tensor. If we substitute the displacement equation into Equation (15), we will have

$$\begin{array}{l} \mathfrak{e}\_{I\!I}^{\varepsilon} = \frac{1}{2} (\mathbf{F}\_{I}^{K} \mathbf{F}\_{I}^{K} - \mathbf{g}\_{II}) = \frac{1}{2} [ (\delta\_{k\!I} + \frac{\partial u\_{k}}{\partial \mathbf{x}\_{I}})(\delta\_{k\!I} + \frac{\partial u\_{k}}{\partial \mathbf{x}\_{I}}) - \mathbf{g}\_{II}] \\ = \frac{1}{2} [\delta\_{I\!\!I} + \frac{\partial u\_{I}}{\partial \mathbf{x}\_{I}} + \frac{\partial u\_{I}}{\partial \mathbf{x}\_{I}} + \frac{\partial u\_{k}}{\partial \mathbf{x}\_{I}} \frac{\partial u\_{k}}{\partial \mathbf{x}\_{I}} - \mathbf{g}\_{II}] \\ = \frac{1}{2} (\frac{\partial u\_{l}}{\partial \mathbf{x}\_{l}} + \frac{\partial u\_{l}}{\partial \mathbf{x}\_{l}} + \frac{\partial u\_{k}}{\partial \mathbf{x}\_{l}} \frac{\partial u\_{k}}{\partial \mathbf{x}\_{l}}) + \frac{1}{2} (\delta\_{II} - \mathbf{g}\_{II}) \end{array} \tag{16}$$

where *εIJ* is the strain tensor in classic elasticity strain. When the deformation of the rock is small, the quadratic term in the above equation can be ignored, i.e.,

$$
\varepsilon\_{I\!I}^{\varepsilon} = \varepsilon\_{II} + \frac{1}{2}(\delta\_{II} - \mathbb{g}\_{II}) \tag{17}
$$

The last two terms of Equation (17) are incompatible deformation terms—that is, the existence of plastic deformation makes the elastic deformation itself uncoordinated. Generally, the local coordinate basis {*EI*} of the initial configuration is different from the local coordinate basis {*ei*} of the current configuration. When the incremental deformation is small, the same coordinate basis vector {*EI*}={*ei*} can be selected. Then, the above equation can be written as

$$
\varepsilon\_{ij}^{\varepsilon} = \varepsilon\_{ij} + \frac{1}{2}\delta\_{ij} - \frac{1}{2}g\_{ij} \tag{18}
$$

#### **3. Stress-Strain Relationship of the Incompatible Plastic Deformation**

In the present section, we consider a Stress-Strain relationship of the incompatible plastic deformation. A number of variables and symbols are listed as follows:


After the rock enters the state of plastic deformation, the strain at any point is composed of elastic strain and plastic strain. When the external load has a small increment, the total strain also has a small increment. The total strain increment consists of an elastic strain increment and a plastic strain increment, namely

$$d\varepsilon\_{i\dot{j}} = d\varepsilon\_{i\dot{j}}^{\varepsilon} + d\varepsilon\_{i\dot{j}}^{p} \tag{19}$$

By substituting Equation (18) into Equation (19), we will obtain

$$d\varepsilon\_{ij}^{\varepsilon} = d\varepsilon\_{ij} + \frac{1}{2}d\delta\_{ij} - \frac{1}{2}dg\_{ij} - d\varepsilon\_{ij}^{p} = d\varepsilon\_{ij} - \frac{1}{2}dg\_{ij} - d\varepsilon\_{ij}^{p} \tag{20}$$

The relationship between the elastic stress increment and strain increment can be determined by Hooke's law as

$$d\sigma\_{i\bar{j}} = \mathbb{C}\_{i\bar{j}k\bar{l}} d\varepsilon\_{k\bar{l}}^{\varepsilon} \tag{21}$$

where *Cijkl* is the elastic stiffness tensor. Since the plastic strain can be obtained through the flow law, the Stress-Strain relationship for ideal elastoplastic materials can be expressed as

$$d\sigma\_{i\bar{j}} = \mathbb{C}\_{i\bar{j}kl} \left( d\varepsilon\_{i\bar{j}} - d\lambda \frac{\partial w}{\partial \sigma\_{kl}} - \frac{1}{2} d\mathfrak{g}\_{kl} \right) \tag{22}$$

where *dλ* is an undetermined nonnegative scalar termed the plastic multiplier. During the plastic deformation, the stress point stays on the yield surface, and this supplementary condition is called the consistency condition, which is expressed by the following formula as

$$dF = \frac{\partial F}{\partial \sigma} d\sigma + \frac{\partial F}{\partial \kappa} d\kappa = a^T d\sigma - A d\lambda = 0 \tag{23}$$

with the definitions of

$$a^T = \frac{\partial F}{\partial \sigma} = \left[ \frac{\partial F}{\partial \sigma\_x}, \frac{\partial F}{\partial \sigma\_y}, \frac{\partial F}{\partial \sigma\_z}, \frac{\partial F}{\partial \tau\_{yz}}, \frac{\partial F}{\partial \tau\_{zx}}, \frac{\partial F}{\partial \tau\_{xy}} \right] \tag{24}$$

$$A = -\frac{1}{d\lambda} \frac{\partial F}{\partial \kappa} d\kappa \tag{25}$$

The vector *a* is termed the flow vector. Equation (22) can be immediately rewritten as

$$(1 - R)d\varepsilon = \left[D\right]^{-1} d\sigma + d\lambda \frac{\partial F}{\partial \sigma} \tag{26}$$

where *D* is the usual matrix of elastic constants. Setting *dgkl*/*dεkl* = 2 *R*, according to Equation (18), we can obtain

$$R = \frac{d\mathbf{g}\_{kl}}{d\varepsilon\_{kl}} = 1 - \frac{d\varepsilon\_{kl}^{\varepsilon}}{d\varepsilon\_{kl}} = 1 - \frac{E}{E\_c} \tag{27}$$

where *Ee* is the elastic modulus of an intact rock, and *E* is the elastic modulus of a rock with defects. The parameter *R* reflects the influence of microscopic defects on macroscopic deformation of the rock.

Premultiplying both sides of Equation (26) by using *d<sup>T</sup> <sup>D</sup>* = *<sup>a</sup>TD* and eliminating *<sup>a</sup>Td<sup>σ</sup>* by the use of Equation (23), we can obtain

$$d\lambda = \frac{(1 - R)}{A + a^T D a} a^T d\_D d\varepsilon \tag{28}$$

Substituting Equation (28) into Equation (26) we can obtain the complete incompatible elastic-plastic incremental Stress-Strain relation as

$$d\sigma = D\_{cp} d\varepsilon \tag{29}$$

where

$$D\_{\varepsilon p} = (1 - R)(D - \frac{d\_D d\_D^T}{A + d^T a});\\ d\_D = Da \tag{30}$$

For numerical computation, it is convenient to write the yield function in terms of stress invariants. The advantage is that it permits the computer coding of yield function and the flow rule in a general form [17]. The principal deviatoric stresses *σ* <sup>1</sup>, *σ* <sup>2</sup>, *σ* <sup>3</sup> are given as the roots of the cubic equation [18]

$$t^3 - l\_2't - l\_3' = 0 \tag{31}$$

Noting the trigonometric identity

$$
\sin^3\theta - \frac{3}{4}\sin\theta + \frac{1}{4}\sin3\theta = 0\tag{32}
$$

and setting *t* = *r*sin*θ* and substituting into Equation (31), we can obtain

$$
\sin^3 \theta - \frac{l\_2'}{r^2} \sin \theta - \frac{l\_3'}{r^3} = 0 \tag{33}
$$

Comparing Equation (32) and Equation (33), we can obtain

$$r = \frac{2}{\sqrt{3}} (l\_2')^{\frac{1}{2}} \tag{34}$$

$$\sin 3\theta = -\frac{I\_3'}{r^3} = -\frac{3\sqrt{3}}{2} \frac{I\_3'}{\left(I\_2'\right)^{3/2}}\tag{35}$$

By noting the cyclic nature of sin (*θ* + 2 *nπ*), we have immediately made the three possible values of sin*θ* which define the three principal stresses. The deviatoric principal stresses are given by *t* = *r*sin*θ* upon substitution of the three values of sin*θ* in turn. Substituting for *r* from Equation (34) and adding the mean hydrostatic stress component give the total principal stresses as

$$\left\{ \begin{array}{c} \sigma\_1\\ \sigma\_2\\ \sigma\_3 \end{array} \right\} = \frac{2}{\sqrt{3}} (l\_2')^{\frac{1}{2}} \left\{ \begin{array}{c} \sin(\theta + \frac{2\pi}{3})\\ \sin\theta\\ \sin(\theta + \frac{4}{3}\pi) \end{array} \right\} + \frac{l\_1}{3} \left\{ \begin{array}{c} 1\\ 1\\ 1 \end{array} \right\} \tag{36}$$

with *σ*<sup>1</sup> > *σ*<sup>2</sup> > *σ*3. For geomaterials, they typically have a frictional strength and different strengths in tension and in compression, so we can choose the Mohr–Coulomb yield criterion and write it in terms of *J*1, *J*2, and *θ* as

$$\frac{1}{3}l\_1\sin\phi + \left(l\_2'\right)^{\frac{1}{2}}(\cos\theta - \frac{1}{\sqrt{3}}\sin\theta\sin\phi) = c\cos\phi\tag{37}$$

In order to calculate the *Dep* matrix in (30), we need to express vector *a* in a form suitable for numerical computation as

$$a^T = \frac{\partial F}{\partial \sigma} = \frac{\partial F}{\partial f\_2} \frac{\partial f\_2}{\partial \sigma} + \frac{\partial F}{\partial \left(f\_2'\right)^{1/2}} \frac{\partial \left(f\_2'\right)^{1/2}}{\partial \sigma} + \frac{\partial F}{\partial \theta} \frac{\partial \theta}{\partial \sigma} \tag{38}$$

Differentiating Equation (35), we can obtain

$$\frac{\partial\theta}{\partial\sigma} = \frac{-\sqrt{3}}{2\cos 3\theta} \left[ \frac{1}{\left(l\_2'\right)^{3/2}} \frac{\partial l\_3}{\partial \sigma} - \frac{\partial l\_3}{\left(l\_2'\right)^2} \frac{\left(l\_2'\right)^{1/2}}{\partial \sigma} \right] \tag{39}$$

Substituting Equation (38) in Equation (37) and using Equation (35), we can write *a* as

$$a = \mathbb{C}\_1 a\_1 + \mathbb{C}\_2 a\_2 + \mathbb{C}\_3 a\_3 \tag{40}$$

where

$$a\_1^T = \frac{\partial I\_1}{\partial \sigma} = \{1, 1, 1, 0, 0, 0\} \tag{41}$$

$$a\_2^T = \frac{\partial \left(l\_2'\right)^{1/2}}{\partial \sigma} = \frac{1}{2\left(l\_2'\right)^{1/2}} \left\{ s\_{x,\prime} s\_{y,\prime} s\_{z,\prime} 2\tau\_{yz,\prime} 2\tau\_{zx,\prime} 2\tau\_{xy} \right\} \tag{42}$$

$$a\_3^T = \frac{\partial f\_3}{\partial \sigma} = \left\{ \begin{array}{c} (s\_y s\_z - \tau\_{yz}^2 + \frac{l\_2'}{3}), (s\_x s\_z - \tau\_{xz}^2 + \frac{l\_2'}{3}), (s\_x s\_y - \tau\_{xy}^2 + \frac{l\_2'}{3}), \\ 2(\tau\_{zx} \tau\_{xy} - s\_x \tau\_{yz}), 2(\tau\_{yz} \tau\_{xy} - s\_y \tau\_{zx}), 2(\tau\_{yz} \tau\_{zx} - s\_z \tau\_{xy}) \end{array} \right\} \tag{43}$$

and

$$C\_1 = \frac{\partial F}{\partial f\_1} \tag{44}$$

$$C\_2 = \frac{\partial F}{\partial \left(l\_2'\right)^{1/2}} - \frac{\tan 3\theta}{\left(l\_2'\right)^{1/2}} \frac{\partial F}{\partial \theta} \tag{45}$$

$$\mathcal{C}\_3 = \frac{-\sqrt{3}}{2\cos 3\theta} \frac{1}{\left(l\_2'\right)^{3/2}} \frac{\partial F}{\partial \theta} \tag{46}$$

The *C*1, *C*2, and *C*<sup>3</sup> are the constants related to the yield surface. The tunnel excavated in underground engineering can be regarded as a problem of two-dimensional plane strain. According to the Mohr-Coulomb yield criterion, the equation from Equation (40) to Equation (46) can be written as

$$a\_1^T = \{1, 1, 0, 1\} \tag{47}$$

$$a\_2^T = \frac{1}{2\left(l\_2'\right)^{1/2}} \{s\_{x'}s\_{y'}2\tau\_{xy'}s\_z\} \tag{48}$$

$$a\_3^T = \frac{\partial l\_3}{\partial \sigma} = \left\{ (s\_y s\_z + \frac{l\_2'}{3}), (s\_x s\_z + \frac{l\_2'}{3}), 2s\_z \tau\_{xy'} (s\_x s\_y + \frac{l\_2'}{3}), \right\} \tag{49}$$

and

$$C\_1 = \frac{1}{3}\sin\phi\tag{50}$$

$$C\_2 = \cos\theta \left[ (1 + \tan\theta \tan 3\theta) + \sin\phi (\tan 3\theta - \tan \theta) / \sqrt{3} \right] \tag{51}$$

$$\mathcal{C}\_3 = \frac{\sqrt{3}\sin\theta + \cos\theta\sin\phi}{2f\_2'\cos3\theta} \tag{52}$$

#### **4. Numerical Analysis**

The FORTRAN language was used to write the elastic-plastic finite element program for the problem of two-dimensional plane strain. The parameters are listed as follows. It is assumed that the rock satisfies the associated flow rule, which is the Mohr–Coulomb yield criterion: the material constants include a Young's modulus *E* = 21 GPa, a Poisson's ratio *v* = 0.3, an internal friction angle *φ* = 30◦, and a cohesion *c* = 5.85 MPa. The problem chosen is a thick-walled cylinder subjected to the external pressure *P*, with plane strain conditions assumed in the axial direction. Axisymmetry normally only requires analyzing a segment of the cylinder. For example, a 90◦ segment (see Figure 2) is used here to take advantage of the symmetry conditions along the *x* and *y* axes.

**Figure 2.** Different mesh divisions for the problem.

Before analyzing the problem, a sensitivity analysis of mesh is necessary. A number of different mesh divisions are shown in Figure 2, where an eight-node quadrangle element was chosen.

Referring to Figure 3, it can be seen that the results of meshes are quite reasonable since the results of radial stress and tangential stress are coincident with the results in the theoretical curve.

**Figure 3.** The distribution of radial stress and tangential stress in the thick cylinder; (**a**) Radial stress distribution under *P* = 14 MPa; (**b**) Tangential stress distribution under *P* = 14 MPa.

Thus, mesh(c) will be used because of the higher accuracy and relatively small amount of computation involved. After having decided on the element type and mesh division, it is now possible to prepare the data before putting them into the program. Such data are extracted directly from the drawing of the structure (see Figure 4) in which all information concerns nodes and element numbers, and load and boundary conditions are available.

**Figure 4.** Details of the example problem.

When the surrounding rock is intact, or with no defects, the distribution of stress is shown in Figure 5 below. When the pressure level of the surrounding rock is at 10 MPa, the surrounding rock is in an elastic state, and the distribution of hoop stress and radial stress is completely consistent with the elastic solution. With the increase in the surrounding rock stress, when the pressure of the surrounding rock is at 12 MPa, the plastic state appears 0.5 m closer to the inside of the tunnel, and the hoop stress increases first and decreases with the radius. The rock to the left of the peak point is in a plastic state, while the rock to the right of the peak point is in an elastic state. As the pressure of the surrounding rock continues to increase (14 MPa, 16 MPa, 18 MPa), the plastic zone of the surrounding rock also expands, and the radius of the plastic zone can be 3.35 m, 4.25 m, and 5.25 m. When the pressure of the surrounding rock reaches 19.4 MPa, all rocks in the region come to a plastic state.

When the surrounding rock is a defect mass, the stiffness of the surrounding rock decreases due to the existence of the defect, per se, and the changes are shown in Figure 6 below. When *R* = 0, the rock has no defects, and the surrounding rock reaches a plastic state at 19.4 MPa. When *R* = 0.2, the surrounding rock reaches the plastic state at 19.2 MPa, indicating that the existence of defects reduces the stiffness of the surrounding rock, but the defects have little influence on the stress of the surrounding rock. When *R* = 0.4, the surrounding rock reaches a plastic state at 17.4 MPa, and the influence of the defects of the surrounding rock stress is further increased. When *R* = 0.5, the defects have a significant effect on the surrounding rock, and the surrounding rock reaches plasticity at 12.8 MPa.

Due to the existence of defects, the stiffness of the surrounding rock decreases, and the deformation increases. The variation rule of deformation with incompatible parameters when the stress of the surrounding rock is at 19 MPa is assessed, as shown in Figure 7 below. When *R* increases from 0.1 to 0.3, the displacement of the surrounding rock increases linearly; when *R* reaches 0.4, the displacement increases more greatly than what it was before.

**Figure 5.** The distribution of radial stress and tangential stress varies with the confining pressure; (**a**) *P* = 10 MPa; (**b**) *P* = 12 MPa; (**c**) *P* = 14 MPa; (**d**) *P* = 16 MPa; (**e**) *P* = 18 MPa; (**f**) *P* = 19.4 MPa.

**Figure 6.** Variation in the confining pressure with incompatible parameter *R* when the rock reaches plasticity.

**Figure 7.** Variation in the deformation of the surrounding rock with incompatible parameters; (**a**) displacement of the surrounding rock; (**b**) displacement of the cylinder surface.

#### **5. Conclusions**


value of parameter *R* is, the faster the rock enters the plastic state and the larger the deformation of the surrounding rock will be.

**Author Contributions:** Conceptualization, Y.B. and B.J.; methodology, B.J.; software, Y.B.; validation, Y.B. and B.J.; formal analysis, Y.B.; investigation, Y.B.; resources, Y.B.; data curation, Y.B.; writing original draft preparation, Y.B.; writing—review and editing, Y.B.; visualization, Y.B.; supervision, B.J.; project administration, B.J.; funding acquisition, B.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was financially supported by General Project of National Natural Science Foundation of China (Grant No. 51174196).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The study did not report any data.

**Acknowledgments:** This work was financed by the General Project of National Natural Science Foundation of China (Grant No. 51174196).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

