# Numerical Simulation Analysis of Wellbore Stability in Weak Interlayer of Igneous Rock
## 3 Numerical Simulation Model for Wellbore Stability
### Determination of the Extent of the Plastic Zone  
The plastic deformation region of the wellbore envelope is calculated using a two-dimensional planar model that is subject to asymmetric loading, as shown in Figure 5.  
When the borehole is opened, the stresses are redistributed, leading to stress concentrations near the wellbore. If the concentrated stress is too high, causing the rock near the wall to exceed its elastic limit and enter a plastic state, it can be assumed that the stresses \(U_{1}(x,y)\) in the plastic strain range follow hyperbolic functions. Meanwhile, the elastic stresses \(U_{2}(x,y)\) can be described using the biharmonic equation [24]:  
\[\frac{\partial^{4}U_{2}}{\partial x^{4}}+2\frac{\partial^{4}U_{2}}{\partial x ^{2}\partial y^{2}}+\frac{\partial^{4}U_{2}}{\partial y^{4}}=0 \tag{4}\]  
The boundary conditions for stresses in the vicinity of the borehole are as follows:  
\[\sigma_{n}=f_{1}(s),\ \tau_{nt}=f_{2}(s) \tag{5}\]  
After the borehole is opened, although it will result in stress concentration, its influence range is limited. It can be considered that the rock formation infinitely far from the borehole remains unaffected and in an elastic state; the boundary conditions at infinity can be expressed as follows:  
\[ \begin{split}\sigma_{x}^{(\infty)}&=P_{1}(x,y)\\ \sigma_{y}^{(\infty)}&=P_{2}(x,y)\\ \tau_{xy}^{(\infty)}&=P_{3}(x,y)\end{split}  \tag{6}\]  
Suppose \(U_{1}(x,y)\) satisfy the functional relationship as:  
\[F_{1}\left(x,y,\frac{\partial^{2}U_{1}}{\partial x^{2}},\frac{\partial^{2}U_{ 1}}{\partial y^{2}},\dots,\frac{\partial U_{1}}{\partial x},\frac{\partial U_ {1}}{\partial y}\right)=0 \tag{7}\]  
Figure 5: Plastic deformation diagram of circular hole under stress. Based on the characteristics of the hyperbolic function, the equation expressing the region where the plastic change occurs be defined as follows:  
\[\left(\frac{\partial^{2}U_{1}}{\partial x^{2}}-\frac{\partial^{2}U_{1}}{\partial y ^{2}}\right)^{2}+4\left(\frac{\partial^{2}U_{1}}{\partial x\partial y}\right)^{ 2}=4k^{2} \tag{8}\]  
where \(k\) is related to the rock yield limit, \(k=\sigma_{r}/2\). On the borehole, \(\sigma_{r}=p\), \(\tau_{r\theta}=0\); and at infinity, \(\sigma_{x}^{(\infty)}=A\), \(\sigma_{y}^{(\infty)}=B\). The solution to Equation (8) is expressed as:  
\[U_{1}(x,y)=kr^{2}\ln\frac{r}{R}-\frac{k-p}{2}r^{2} \tag{9}\]  
Combined with the stress conditions on the borehole:  
\[\left\{ \begin{array}{l}\sigma_{x}^{(1)}+\sigma_{y}^{(1)}=2k+2p+2k\ln\frac{ \Xi\Xi}{R^{2}}\\ \sigma_{x}^{(1)}+\sigma_{y}^{(1)}+2i\tau_{xy}^{(1)}=2k\frac{\Xi}{\Xi}\end{array} \right. \tag{10}\]  
If the boundary of the region of plastic change is denoted as \(L\), it is assumed that the region of elastic deformation is related to functions \(\phi_{2}^{*}(z)\) and \(\psi_{2}^{*}(z)\) as follows:  
\[\frac{\partial^{2}U_{2}}{\partial y^{2}}+\frac{\partial^{2}U_{2}}{\partial x ^{2}}=4\text{Re}[\phi_{2}^{*}(z)] \tag{11}\]  
\[\frac{\partial^{2}U_{2}}{\partial x^{2}}-\frac{\partial^{2}U_{2}}{\partial y ^{2}}-zi\frac{\partial^{2}U_{2}}{\partial x\partial y}=2\left[\Xi\phi_{2}^{* *}(z)+\psi_{2}^{*}(z)\right] \tag{12}\]  
At the boundary, the stress function should satisfy the requirement for both plastic and elastic deformation to occur:  
\[\left\{ \begin{array}{l}4\text{Re}\phi_{2}^{*}(z)=2k+2p+2k\ln\frac{z\Xi}{R^ {2}}\\ 4\text{Re}\phi_{2}^{*}(z)=A+B\\ 2[\Xi\phi_{2}^{*}(z)+\psi_{2}^{*}(z)]=2k\frac{\Xi}{\Xi}\\ 2[\Xi\phi_{2}^{*}(z)+\psi_{2}^{*}(z)]=B-A\end{array} \right. \tag{13}\]  
Introduce functions to solve the above problems:  
\[z=\omega(\zeta)=c\zeta+\frac{c_{1}}{\zeta}+\frac{c_{2}}{\zeta^{2}}+\cdots+ \frac{c_{n}}{\zeta^{n}} \tag{14}\]  
Assuming a hypothetical surface \(\zeta\) with a unit circle \(\gamma\) on it, we can mathematically project the region outside the plastic boundary onto the exterior of the hypothetical circle, resulting in the following:  
\[\phi_{2}^{*}[\omega(\zeta)]=\phi_{2}(\zeta) \tag{15}\]  
\[\psi_{  
The function \(\phi_{2}^{*}(z)\) is obtained in the form expressed as:  
\[\left\{