# Comparative Study of Wellhole Surrounding Rock under Nonuniform Ground Stress

Authors: Zhongyu Jiang; Guoqing Zhou
Journal: Advances in Civil Engineering
 (2019)

---

## Abstract

The stress analysis of the wellhole surrounding rock and the regular failure of the wellhole has always been a concern for the well builders. Firstly, the Hamilton canonical equations are obtained by using the Hamiltonian variational principle in the sector domain, and the zero eigensolution and nonzero eigensolutions of the homogeneous equation are solved. According to the Hamiltonian operator matrix with the orthogonal eigenfunction system, the special solution form of the nonhomogeneous boundary condition equation is obtained. Then, according to the principle of the same coefficient being equal, the relationship equation between the direction eigenvalue and the angle coefficient is obtained, from which the specific expression of the special solution of the equation can be determined. Furthermore, the analytical solution of the wellhole surrounding rock problem under nonuniform ground stress is obtained by using the linear elastic accumulative principle. Finally, a concrete example is given to compare the finite element method and the symplectic algorithm. The results are consistent, which ensures the accuracy and the reliability of the symplectic algorithm. The relationship between the circumferential stress distribution around the hole and the lateral pressure coefficient is further analyzed.


---

## Body

## 1. Introduction

In the mine construction, the stress analysis of the wellhole surrounding rock and the regular failure of the wellhole are always been a concern for the well builders [1]. This has not only important theoretical significance but also practical application value. In the past, the axisymmetric plane strain analysis of a circular wellhole subjected to two-dimensional isobaric original ground stress field in the distance is familiar to all [2]. However, the nonaxisymmetric problem under the two-dimensional unequal pressure condition is more common in engineering. Because of the complexity of boundary conditions, it is more difficult to solve these kinds of problems.At present, although there are many accurate and easy-to-use finite element calculation programs, these methods make it easier to analyze the stress and deformation of the wellhole surrounding rock [3]. However, the pursuit of basic theoretical solutions is still the goal of engineers and technicians unremitting efforts. This is because the theoretical solution can provide a valuable means for design engineers to evaluate the results of numerical analysis. Moreover, the theoretical solution can reveal the relationship between variables more intuitively. Verruijt mapped the half-plane problem of a circular tunnel to a circular ring using complex variables and gave a detailed analytical solution to the case of a uniform radial displacement of the circular tunnel boundary [4]. Exadaktylos and Stavropoulou gave the closed-form plane strain solution of stresses and displacements around a semicircular or “D” shaped cross-sectional tunnel periphery based on complex potential function and conformal mapping representation [5]. Zhu used the theory of complex variable function to analyze the surrounding rock stress of the underground excavation tunnel. By conformal transformation, the original boundary problem of the excavated section is mapped to a new problem with a unit circle as its boundary. Then, two analytic functions are used to express the surrounding rock stress of the excavated tunnel [6]. However, it is difficult to determine the analytic function in the derivation process, and it needs to be recalculated for different external boundary conditions, lacking generality.Zhong [7] and Yao et al. [8] developed symplectic elasticity theory in the Hamilton system. It shows great advantages in revealing the physical essence of elastic solutions and calculating local mechanical properties. Hamilton state space [6–9] breaks through the limitations of the traditional separating variable method [5] and provides a broad prospect for solving existing problems of elasticity. In state space, according to the expansion theorem of eigenfunction, the arbitrary integrable state vector can be represented by orthogonal eigenvectors. The determination of eigensolution is very important for the explicit solution of a problem.Leung applied the Hamiltonian method to study the two-dimensional elastic plane problems. It is a method of separable variables for partial differential equations using displacements and their conjugate stresses as unknowns. The analytical solution of stress distribution under arbitrary boundary and without considering the assumption of the beam is obtained [9]. Tseng and Tarn based on the Hamiltonian state space method, by means of variable separation and symplectic eigenfunction expansion, discussed the theoretical solution of stress field around a circular hole in an elastic plate with a hole subjected to unidirectional tension [10]. Tseng and Tarn used the state space method to analyze the deformation and stress distribution of circular elastic plates under axisymmetric loads. Based on this, the effect of thickness on the bending of circular plates was evaluated [11].


## 2. Description of the Problem

Consider the wellhole with radiusR0 in the x and y direction at infinity is being affected by nonuniform ground stress σx and σy (Figure 1). The lateral pressure coefficient is given by the relation m=σy/σx, where assume σx≥σy, then 0 ≤ m ≤ 1.Figure 1
Force diagram of the surrounding rock wall of nonuniform pressure.Basic assumptions are as follows:(1)
The entire stress state is in the elastic zone, and the wellhole and surrounding rock wall are isotropic linear elastic materials(2)
The stress boundary conditions at infinity remain constant along the axial direction of the wellhole, and this problem is reduced to the plane strain problem in the infinite domainFor the square boundary conditions of this problem, it is advisable to use Cartesian coordinates. But here the stress near the wellhole is mainly discussed, so it is more appropriate to use polar coordinates. Firstly, the square boundary is transformed into a circular boundary. The center of the wellhole is taken as the center of the circle, a lengthR∞ is taken which is much larger than radius R0 of the wellbore as the radius, and great circle (shown in the dotted line in Figure 1) is drawn. Considering that the radius of the great circle is much larger than the radius of the wellhole, the stress of point A at the great circle is the same as that without hole. According to the stress state formula of a point, the stress components in the polar coordinates of that point can be obtained. Thus, the original problem is transformed into a ring with an inner diameter of R0 and outer diameter of R∞. The boundary conditions are as follows.Whenρ=R0,(1)σρ=0,τρφ=0.Whenρ=R∞,(2)σρ=1+m2σx+1−m2σxcos2φ,τρφ=1−m2σxsin2φ.


## 3. Hamiltonian Variational Principle of the Polar Coordinate System

Define new variables:Sρ=ρσρ, Sφ=ρσφ, and Sρφ=ρτρφ, in a typical sector region (Figure 2) R1≤ρ≤R2, α≤φ≤β. Then, perform variable substitution ξ=lnρ, that is, ρ=eξ, and ξ1=lnR1 and ξ2=lnR2. In the variational principle, φ is simulated as time coordinates, ξ is horizontal, and Sρ lateral force factors are eliminated. Then, the mixed energy variational principle of the Hamiltonian system under polar coordinates is obtained:(3)δ∫αβ∫ξ1ξ2Sρφ∂uρ∂φ+Sφ∂uφ∂φ+Sφuρ+v∂uρ∂ξ−Sρφuφ−∂uφ∂ξ+12E∂uρ∂ξ2−12E1−v2Sφ2+21+vSρφ2dξdφ=0.Figure 2
Polar coordinates sector region.IntroducingqT=uρuφ for the displacement vector, pT=SρφSφ for the dual vector of q, and ⋅=∂/∂φ for the derivative of φ, equation (3) can be rewritten as follows:(4)δ∫αβ∫ξ1ξ2pTq˙−Hq,pdξdφ=0,where the Hamiltonian density function H is(5)Hq,p=Sρφuφ−∂uφ∂ξ−Sφuρ+v∂uρ∂ξ−12E∂uρ∂ξ2−12E1−v2Sφ2+21+vSρφ2.The variational principle is developed to obtain the Hamiltonian regular equations:(6)q˙p˙=ADB−A∗qp,where(7)A=01−∂∂ξ−1−v∂∂ξ0,B=−E∂2∂ξ2000,D=21+vE001−v2E,A∗=0−1+v∂∂ξ1+∂∂ξ0.The mutual adjoint operator matrix withA∗ is A, and the total-state vector υ and the matrix H are introduced; Hamilton's regular equations (6) become(8)υ˙=Hυ,where(9)υ=qp=uρuφSρφSφT,H=ADB−A∗.It is not difficult to prove [8] that the matrix H is a Hamilton operator matrix of symplectic geometric space.The homogeneous boundary conditions for the sidesξ=ξ1 and ξ=ξ2 are(10)E∂uρ∂ξ+vSφ=0,Sρφ=0.Separate the variable to solve the homogeneous equation, and letυφ,ξ=eμφψξ; substituting it into the Hamiltonian dual equation, obtain the following eigenequation:(11)Hψξ=μψξ,where μ is the eigenvalue and ψξ is the eigenfunction vector.


## 4. General Solution of Homogeneous Boundary Conditions

First, the eigenequation (11) is expanded:(12)−μuρ+uφ−duφdξ+21+vESρφ+0=0,−uρ−vduρdξ−μuφ+0+1−v2ESφ=0,−Ed2uρdξ2+0−μSρφ+Sφ−vdSφdξ=0,0+0−Sρφ−dSρφdξ−μSφ=0.This is the simultaneous ordinary differential equations forξ. The solution should first find the direction eigenvalue λ of ξ. The characteristic equation is(13)det−μ1−λ21+vE0−1−vλ−μ01−v2E−Eλ20−μ1−vλ001−λ−μ=0.Expand its determinant, and solve(14)λ1,2=±1+μi,λ3,4=±1−μi.The general expression of equation (12) is(15)uρ=A1eλ1ξ+A2eλ2ξ+A3eλ3ξ+A4eλ4ξ,uφ=B1eλ1ξ+B2eλ2ξ+B3eλ3ξ+B4eλ4ξ,Sρφ=C1eλ1ξ+C2eλ2ξ+C3eλ3ξ+C4eλ4ξ,Sφ=D1eλ1ξ+D2eλ2ξ+D3eλ3ξ+D4eλ4ξ.According to the difference in eigenvalues ofμ, the expression of the general solution has changed. When the eigenvalues μ=0 and ±i, the direction eigenvalue λ has multiple roots, which means that the eigenvalue μ has Jordan-type eigensolutions and has its special physical meaning. The eigensolution does not decay with the increase or decrease of the ordinate φ and is the basic eigensolution.Whenμ=0, the eigensolution and its Jordan-type eigensolution are expressed as(16)ψ00=0eξ00,ψ01=c1eξ+c2e−ξ+1−v2ξeξ00E1−vc1eξ+E1+vc2e−ξ+E2eξξ+2−v1−v,where(17)c1=−12−1−v2R22lnR2−R12lnR1R22−R12,c2=−1+v2R22R12R22−R12lnR2R1.Whenμ=±i, the eigensolution and its Jordan-type eigensolution are expressed as(18)ψi0=1i00T,ψ−i0=1−i00T,ψi1=uρ1iuφ1Sρφ1Sφ1iT,ψ−i1=−uρ1iuφ1Sρφ1−Sφ1iT,where(19)uρ1=121−vξ+α1−3ve2ξ+β1+ve−2ξuφ1=−121+v+1−vξ+α5+ve2ξ+β1+ve−2ξSρφ1=E12+2αe2ξ−2βe−2ξSφ1=E12+6αe2ξ+2βe−2ξ,α=−14R12+R22,β=−αR12R22.The original problem corresponding to the eigensolution and its Jordan eigensolution is solved as follows:(20)υμn=∑m=0n1m!φmψμn−m,n=0,1.Since the eigenvalues of the current problem are complex numbers, the corresponding eigenvectors are also complex forms. The emergence of complex numbers has brought great trouble to operation. The original problem was of the real type, so the complex eigensolution can be converted into a real type. According to the properties of symplectic orthogonal conjugation of the Hamiltonian matrix [8], the real and imaginary parts of the complex eigenvectors correspond to the symplectic conjugate eigenvalue solution of the original problem, respectively. It can be seen that υi0 and υ−i0 are a pair of complex conjugate eigensolutions with each other. The real part and imaginary part are written as follows:(21)υiR0=cosφ−sinφ00T,υiI0=sinφcosφ00T,υi1 and υ−i1 are another pair of complex conjugate eigensolutions. The real part and imaginary part are written as follows:(22)υiR1=φcosφ−uρ1sinφ−φsinφ+uφ1cosφSρφ1cosφ−Sφ1sinφ,υiI1=φsinφ+uρ1cosφφcosφ+uφ1sinφSρφ1sinφSφ1cosφ.These real solutions are all solutions of the original problem.The general solution of homogeneous boundary conditions for the original problem equation (8) can be expressed as(23)υh=a1υ00+a2υ01+a3Reυi0+a4Imυi0+a5Reυi1+a6Imυi1,where the coefficients a1∼a6 are determined by the boundary conditions at both ends.


## 5. Special Solution for Nonhomogeneous Boundary Conditions

Above all, the general solution of the eigenequation under the homogeneous boundary conditions is solved, but the special solution under the nonhomogeneous boundary condition is much more complicated. According to the Hamiltonian matrix property [12, 13], block operator A has orthogonal eigenfunction systems in Hilbert space X×X and the eigenvalues and eigenfunction systems of A can be expressed as follows:(24)Uk=coskφsinkφ,K=0,±1,±2,⋯.Block operator−A∗ has orthogonal eigenfunction systems in Hilbert space X×X, and the eigenvalues and eigenfunction systems of −A∗ can be expressed as follows:(25)U^k=sinkφcoskφ,K=0,±1,±2,⋯.The boundary conditions of the original problem are relatively complex, so it is difficult to solve it directly. Since the rock deformation is small elastic deformation, the linear elastic accumulative theory is applicable. Therefore, the complex boundary conditions can be decomposed into a plurality of simple stress situations. The problem can be solved by accumulating each simple stress situation result [14]. Any one complicated boundary function can be decomposed into a simple regular trigonometric function system by Fourier decomposition. Therefore, the two sides only subjected to the boundary load of normal cosine angle should be taken into consideration firstly. The specific form is as follows.Whenξ=ξ1,(26)E∂uρ∂ξ+vSφ=0,Sρφ=0.Whenξ=ξ2,(27)E∂uρ∂ξ+vSφ=P0coskφ,Sρφ=0,where P0 is the normal load coefficient and k is the angle coefficient, taking an integer. Especially when k=0, the boundary load is a constant term.According to the general solution of the homogeneous equation (15) and the eigenfunction equations (24) and (25) of the block operators A and −A∗, the special solution of the equation satisfying the nonhomogeneous boundary conditions equation (26) can be written as(28)u˜ρ=P0A1eλ1ξ+A2eλ2ξ+A3eλ3ξ+A4eλ4ξcoskφ,u˜φ=P0B1eλ1ξ+B2eλ2ξ+B3eλ3ξ+B4eλ4ξsinkφ,S˜ρφ=P0C1eλ1ξ+C2eλ2ξ+C3eλ3ξ+C4eλ4ξsinkφ,S˜φ=P0D1eλ1ξ+D2eλ2ξ+D3eλ3ξ+D4eλ4ξcoskφ.These constants in the above formula are not completely independent, and they should also satisfy equation (8) so that the relationship between these constants can be obtained:(29)Ai=1−vλi1+λi−k2Ekλi2Ci,Bi=λi3+λi2+vk2−1λi+k2−1Ek2λi2Ci,Bi=1−k2+1−vλi+2+vλi2Eλi−1λi2Ci,Di=1+λikCi.According to the principle that the same coefficient must be equal, a quartic equation about the direction eigenvalueλ and the angle coefficient k is obtained.(30)λ4−2k2+1λ2+k2−12=0.Thus, the direction eigenvalueλ can be calculated according to the angle coefficient k, and the specific expression of the special solution of the nonhomogeneous equation can be determined. The values of the integral constants Ai, Bi, Ci, and Di are then determined according to the boundary conditions equation (26).Similarly, the boundary loads in the form of tangential direction sine multiple times angle on both sides are considered. The specific form isWhenξ=ξ1,(31)E∂uρ∂ξ+vSφ=0,Sρφ=0.Whenξ=ξ2,(32)E∂uρ∂ξ+vSφ=0,Sρφ=T0sinkφ,where T0 is the tangential load coefficient and k is the angle coefficient, taking an integer.By the same method, a special solution of the equation that satisfies the boundary condition equation (31) and (32) can be calculated. Finally, according to the linear elastic accumulative theory, the theoretical solution which satisfies the boundary condition equations (1) and (2) of the original problem can finally be obtained.


## 6. Examples

In order to verify the correctness of the results obtained by the symplectic elastic method, the finite element software was used for the comparative calculation. The relevant data used in the calculation are as follows: the wellhole radiusR0=5m, calculating boundary radius at a distance R∞=50m, elastic modulus of surrounding rock E=30GPa, Poisson’s ratio v=0.2, nonuniform ground stress conditions: σx=−50MPa and σy=−30MPa, and the lateral pressure coefficient m=0.6. According to the boundary condition equations (1) and (2) of the original problem, the boundary condition can be divided into three parts: (1) the normal direction is affected by the uniform load; (2) the normal direction is affected by the 2 times cosine load; and (3) the tangential direction is affected by the 2 times sinusoidal load.Firstly, the boundary is subjected to normal uniform load, that is, the normal load coefficientP0=−40MPa and the angle coefficient k=0, and substituting them into the formula equation (30), the direction eigenvalue λ1,2,3,4=±1 can be obtained. Then, the special solution equation (28) is substituted into the boundary condition to obtain the undetermined coefficients, and a special solution of corresponding boundary conditions is obtained. In the same method, special solutions of corresponding other boundary conditions can be obtained. Finally, by accumulating the special solutions under various boundary conditions, we can get the displacement and stress cloud diagram of the surrounding rock of the wellhole under nonuniform ground stress.By comparing and analyzing the result cloud diagrams in Figures3 to 7, it can be seen that the distribution of each result cloud diagrams obtained by the two methods is consistent, and there is little difference in the value of extreme points. These show the correctness and reliability of the symplectic method. The difference between the two methods is that the symplectic method is an analytical method, while the finite element method is a numerical method.Figure 3
Cloud diagram of the radial displacementuρ.
(a)
(b)Figure 4
Cloud diagram of the circumferential displacementsuφ.
(a)
(b)Figure 5
Cloud diagram of the radial stressesσρ.
(a)
(b)Figure 6
Cloud diagram of the circumferential stressesσφ.
(a)
(b)Figure 7
Cloud diagram of the shear stressτρφ.
(a)
(b)Comparing the stress extremum in Figures5 and 6, it can be seen that the circumferential stress extremum around the hole is much larger than the radial stress extremum, which is consistent with the engineering practice. In Figure 6, the circumferential stresses are all compressive stresses, and the distribution is extremely inhomogeneous at the edge of the hole. The value of circumferential stress in the horizontal position is relatively small, and it reaches the maximum value in the vertical direction. Tensile stress may occur in the horizontal direction under certain conditions. Considering the brittleness of geotechnical materials, tensile failure should be avoided as far as possible.The influence of different lateral pressure coefficientsm on the circumferential stress value around the hole is detailed in Table 1. From the data, it can be seen that the circumferential stress around the hole changes periodically along the azimuth angle. Especially when the azimuth angle is equal to 0°, 90°, 180°, 270°, and 360°, the circumferential stress reaches the extreme value. The maximum occurs at 0°, 180°, and 360°, and the minimum occurs at 90° and 270°. The data show that when the azimuth angle is 0°, 180°, and 360°, the maximum circumferential stress around the hole gradually changes from compressive stress to tensile stress with the decrease of the lateral pressure coefficient. When the lateral pressure coefficient m is between 0.3 and 0.4, the circumferential stress will appear on the inflection point; that is, the circumferential stress is zero.Table 1
Circumferential stressσφ around the hole corresponding to different m values (unit: MPa).
Angle
m
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0°
−101.01
−85.76
−70.5
−55.25
−40
−24.74
−9.49
5.76
21.02
36.27
51.53
30°
−101.01
−90.86
−80.71
−70.55
−60.4
−50.25
−40.1
−29.95
−19.79
−9.64
0.51
60°
−101.01
−101.06
−101.11
−101.16
−101.21
−101.27
−101.32
−101.37
−101.42
−101.47
−101.52
90°
−101.01
−106.16
−111.32
−116.47
−121.62
−126.77
−131.93
−137.08
−142.23
−147.38
−152.54
120°
−101.01
−101.06
−101.11
−101.16
−101.21
−101.27
−101.32
−101.37
−101.42
−101.47
−101.52
150°
−101.01
−90.86
−80.71
−70.55
−60.4
−50.25
−40.1
−29.95
−19.79
−9.64
0.51
180°
−101.01
−85.76
−70.5
−55.25
−40
−24.74
−9.49
5.76
21.02
36.27
51.53
210°
−101.01
−90.86
−80.71
−70.55
−60.4
−50.25
−40.1
−29.95
−19.79
−9.64
0.51
240°
−101.01
−101.06
−101.11
−101.16
−101.21
−101.27
−101.32
−101.37
−101.42
−101.47
−101.52
270°
−101.01
−106.16
−111.32
−116.47
−121.62
−126.77
−131.93
−137.08
−142.23
−147.38
−152.54
300°
−101.01
−101.06
−101.11
−101.16
−101.21
−101.27
−101.32
−101.37
−101.42
−101.47
−101.52
330°
−101.01
−90.86
−80.71
−70.55
−60.4
−50.25
−40.1
−29.95
−19.79
−9.64
0.51
360°
−101.01
−85.76
−70.5
−55.25
−40
−24.74
−9.49
5.76
21.02
36.27
51.53Figure8 depicts the influence of different lateral pressure coefficients m on the circumferential stress around the hole. When the value of m changes from 1 to 0, the stress distribution around the hole changes from circular to dumbbell-shaped, then to “8” shape, and finally to double “8” shape (where the vertical “8” shape is much larger than the horizontal “8” shape). The circumferential stress around the hole in the horizontal direction gradually changes from compressive stress to tensile stress with the decrease of the m value. Especially when the value of m is near 0.3, the stress critical point appears. When the value of m is less than 0.3, the tensile stress will appear in the horizontal direction of the hole edge. This is consistent with the data in Table 1. The circumferential stress around the hole in the vertical direction increases slowly with the decrease of the m value. Especially when m = 0, the maximum tensile stress appears in the horizontal direction and the maximum compressive stress appears in the vertical direction. According to the properties of wellhole materials, the circumferential stress of the sidewall should be avoided as far as possible in engineering practice, which is particularly important.Figure 8
Circumferential stress curves around the hole corresponding to differentm values.


## 7. Conclusion

The symplectic method realizes the conversion from Euclidean space to Hilbert space, breaking through the limitations of the traditional variable separation method. In this paper, the eigenvalues and eigenfunctions of Hamiltonian block operators are written according to the symplectic elastic mechanics method in the polar coordinate system. According to the characteristics of Hamiltonian block operators, the problem of nonhomogeneous boundary conditions is solved. And the relationship equation (30) between the direction eigenvalue λ and the angle factor k is found, which provides a theoretical foundation for further calculation of arbitrary boundary conditions.The symplectic elastic mechanics method is used to accurately analyze the wellhole surrounding rock under nonuniform ground stress. By comparing the results of the symplectic algorithm and finite element method, we can see that the distribute regulation of the two methods is completely consistent, and the numerical values are almost the same. The results show that the circumferential stress distribution of the wellhole rock is closely related to the lateral pressure coefficientm. With the change of the lateral pressure coefficient m from 1 to 0, the circumferential stress of the wellhole rock changes differently in horizontal and vertical directions. In the horizontal position, the circumferential compressive stress gradually decreases with the decrease of the m value until it becomes tensile stress and finally reaches the maximum of tensile stress. In the vertical direction, the circumferential compressive stress increases slowly with the decrease of the m value and finally reaches the maximum of compressive stress. The research results not only provide theoretical guidance for wellhole design but also provide technical support for preventing wellhole fracture.



---

*Source: 7424123-2019-06-25.xml*