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

: Before the excavation of underground engineering, joints, ﬁssures, 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 ﬁnite 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.


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][2][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 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][9][10][11][12][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][15][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.

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ζ. These 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][9][10][11][12][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][15][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.

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ζ.  A point X on Ω 0 is mapped to a point ζ on the intermediate configuration Ω * under the mapping χ p : 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 Processes 2021, 9, 2215 3 of 12 neighborhood of point ζ has coordinates (ζ 1 , . . . , ζ n ); χ p is a continuous mapping; and the transformation matrix between X and ζ is ∂ζ/∂X.
In differential geometry, χ p 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: A metric is defined in a Riemannian manifold as The metric tensor field g αβ on L is dragged back to the metric field G IJ on M by 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: A is the nonholonomic coefficient. The components satisfy Point ζ on Ω * is mapped to point x on Ω under mapping χ e : There is a tangent mapping at ζ: The local coordinate of point ζ is (ζ 1 , . . . , ζ n ), the corresponding coordinate of point x is (x 1 , . . . , x 3 ), χ e 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 ε e .
Because the elastic deformation is reversible, we have F is the deformation gradient, and its components satisfy 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 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 χ p is injective, then χ e is surjective. The initial configuration is in a three-dimensional Euclidean space, and the length square of the element at point X is The intermediate configuration is a Riemannian manifold, whose element length squared is The current configuration is also in a three-dimensional Euclidean space, and its element length squared is Therefore, the change in the squared length of the element from the intermediate configuration to the current configuration is where ε e I J is the elastic strain tensor. If we substitute the displacement equation into Equation (15), we will have 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., 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 {E I } of the initial configuration is different from the local coordinate basis {e i } of the current configuration. When the incremental deformation is small, the same coordinate basis vector {E I } = {e i } can be selected. Then, the above equation can be written as

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: 1.
The internal metric tensor g ij ; 5.
The plastic potential w; 9.
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 By substituting Equation (18) into Equation (19), we will obtain The relationship between the elastic stress increment and strain increment can be determined by Hooke's law as where C ijkl 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 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 with the definitions of The vector a is termed the flow vector. Equation (22) can be immediately rewritten as where D is the usual matrix of elastic constants. Setting dg kl /dε kl = 2 R, according to Equation (18), we can obtain where E e 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 T D = a T D and eliminating a T dσ by the use of Equation (23), we can obtain Substituting Equation (28) into Equation (26) we can obtain the complete incompatible elastic-plastic incremental Stress-Strain relation as where 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 σ 1 , σ 2 , σ 3 are given as the roots of the cubic equation [18] Noting the trigonometric identity and setting t = rsinθ and substituting into Equation (31), we can obtain Comparing Equation (32) and Equation (33), we can obtain 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 = rsinθ 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 with σ 1 > σ 2 > σ 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 Substituting Equation (38) in Equation (37) and using Equation (35), we can write a as where and The C 1 , C 2 , and C 3 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) and

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 Processes 2021, 9, 2215 8 of 12 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. 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.  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. 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. 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. 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  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. 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    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. 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. 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.

Conclusions
(1) The existence of defects is the cause of incompatible deformation, and plastic deformation is the result of incompatible deformation as well. Although the deformed object has undergone much plastic deformation, it still exists in a threedimensional Euclidean space, so the rock before and after the deformation is in a Euler space. The incompatible plastic deformation leads to the incompatible elastic deformation, while the total deformation can be coordinated.
(2) When the deformation of the rock is small, the deformed and undeformed body can use the same coordinate base vector, and assuming the rock is an ideal plastic material, we can obtain the incremental plastic stress-strain relationship of the rock. The plastic strain increment is not only associated with the deviatoric stress state but also with the defect distribution. This distribution can be expressed as an incompatible parameter R. The larger the value of R is, the more defects there are in the rock and the smaller the elastic-plastic stiffness matrix will be. (3) The excavation of a deep underground tunnel can be regarded as a problem of twodimensional plane strain. An elastic-plastic finite element program is written in the FORTRAN language to analyze circular tunnels under hydrostatic pressure. As the existence of defects decreases the stiffness of the surrounding rock, the larger the value of parameter R is, the faster the rock enters the plastic state and the larger the deformation of the surrounding rock will be.

Conclusions
(1) The existence of defects is the cause of incompatible deformation, and plastic deformation is the result of incompatible deformation as well. Although the deformed object has undergone much plastic deformation, it still exists in a three-dimensional Euclidean space, so the rock before and after the deformation is in a Euler space. The incompatible plastic deformation leads to the incompatible elastic deformation, while the total deformation can be coordinated. (2) When the deformation of the rock is small, the deformed and undeformed body can use the same coordinate base vector, and assuming the rock is an ideal plastic material, we can obtain the incremental plastic Stress-Strain relationship of the rock. The plastic strain increment is not only associated with the deviatoric stress state but also with the defect distribution. This distribution can be expressed as an incompatible parameter R. The larger the value of R is, the more defects there are in the rock and the smaller the elastic-plastic stiffness matrix will be. (3) The excavation of a deep underground tunnel can be regarded as a problem of twodimensional plane strain. An elastic-plastic finite element program is written in the FORTRAN language to analyze circular tunnels under hydrostatic pressure. As the existence of defects decreases the stiffness of the surrounding rock, the larger the value of parameter R is, the faster the rock enters the plastic state and the larger the deformation of the surrounding rock will be. Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

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