西南石油大学学报(自然科学版)  2020, Vol. 42 Issue (5): 135-144
考虑岩石弹塑性变形的井周应力场计算模型研究    [PDF全文]
李小刚1, 李昀昀2, 胡秋萍3, 张翔3, 易良平1    
1. 油气藏地质及开发工程国家重点实验室·西南石油大学, 四川 成都 610500;
2. 中海油能源发展股份有限公司工程技术分公司非常规院, 天津 滨海新区 300452;
3. 中联煤层气有限责任公司晋城分公司, 山西 晋城 048000
摘要: 针对煤岩、页岩等岩石具有弹塑性这一特性,考虑岩石非线性硬化与软化变形,并耦合工作液渗滤作用,沿井筒径向方向将地层分为塑性软化、塑性硬化和弹性区,根据弹性区与塑性区交界处应力连续的特点,建立了耦合工作液渗滤和岩石非线性硬化与软化变形的井壁应力场模型。分析表明,塑性区大小与井眼流体压力、泊松比和地层初始孔隙流体压力呈正相关关系,与岩石强度呈负相关关系。与弹性状态相比,岩石塑性变形降低了井眼应力集中效应,岩石塑性变形对井周径向应力影响较小,但会使井周周向应力增加。
关键词: 煤岩    页岩    弹塑性    非线性本构    井周应力    
A Study on Calculation Model of Stress Field Around Wells Considering Rock Elastic-plastic Deformation
LI Xiaogang1, LI Yunyun2, HU Qiuping3, ZHANG Xiang3, YI Liangping1    
1. State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation, Southwest Petroleum University, Chengdu, Sichuan 610500, China;
2. CNOOC EnerTech-Drilling & Production Co., Binhai New Area, Tianjin 300452, China;
3. Jincheng Branch, China United Coalbed Methane Co. Ltd., Jincheng, Shanxi 048000, China
Abstract: Aiming at the elastic-plastic characteristics of coal, shale and other rocks, a calculation method for stress field of elastic-plastic formation is established, which divide the formation into the plastic softening, the plastic hardening and the elastic zone along the radial direction. In this new model the nonlinear hardening and softening deformation of rock and the percolation effect of the working fluid are considered. The analysis shows that the plastic zone range is positively correlated with the borehole fluid pressure, Poisson's ratio and the initial pore fluid pressure, and is inversely related to the strength of the rock. Compared with the elastic state, the plastic deformation of rock decreases the stress concentration effect, and the plastic deformation of the rock has a little influence on the radial stress of the formation while it can increase the circumferential stress of the formation.
Keywords: coal    shale    elastic-plastic    nonlinear constitutive equation    stress field around wellbore    
引 言

在油气田开发过程中,井周应力场的准确预测非常重要。因为这是钻井液安全密度窗口预测和压裂破裂压力预测的基础[1-2]。传统的破裂压力和坍塌压力计算模型主要是建立在弹性力学基础上的[3-5],但煤岩、页岩等岩石实际上具有弹塑性[6-7]。目前,关于弹塑性地层井周或巷道周围应力场的研究得到了国内外石油工程领域和岩土工程领域学者的广泛关注。Aadnoy等根据弹塑性理论,建立了弹塑性地层井周应力场模型,并根据最大拉应力准则,建立了弹塑性地层裂缝起裂模型[8-9],但其将地层考虑为理想弹塑性,且没有考虑孔隙流体压力的作用。Wang等研究了非渗透性地层和渗透性地层井眼屈服和裂缝起裂[10-11],使用摩尔库伦应变软化准则估算塑性区域扩展,针对渗透性地层耦合稳态孔隙压力分布,建立了井眼裂缝起裂模型,但采用的岩石本构模型过于理想,没有考虑破坏前的非线性变形特征。潘岳等考虑岩石的非线性硬化与软化变形,采用全量岩体本构方程研究了圆孔隧道围岩应力场[12-13],但忽略了孔隙流体压力作用。郭建春等采用塑性全量理论,考虑岩石破坏前的非线性变形特征,建立了弹塑性地层井周应力场和起裂分析模型[14],定量研究了弹塑性地层起裂模式和起裂压力,但并没有考虑孔隙流体压力的作用。Shao等采用有限元数值方法研究岩石弹塑性变形对井周应力场和井壁稳定性的影响[15-16],但是有限元数值模拟方法计算复杂,编程能力要求高,不便于推广应用。

本文在前人研究的基础之上,全面考虑岩石非线性变形特性,耦合工作液渗滤诱导孔隙压力,建立了一种考虑岩石塑性变形的井周应力场计算模型。

1 基本假设及岩石本构 1.1 基本假设及井筒受力分析

假设水平方向上两个主应力差异不大,为方便方程推导,定义远场水平地应力$\sigma _{{\rm h}}$=($\sigma_{{ x}}+\sigma_{{ y}}$)/2,井周岩石同时受到远场水平地应力和井筒内流体压力作用,如图 1所示。由于地层岩石为多孔介质,因此井周岩石骨架还受到孔隙流体压力的作用。由于应力集中效应,井眼附近区域应力值比远井地带高,因此,井眼附近最容易发生屈服[17]。根据岩石是否屈服,将井周分为弹性区和塑性区,其中,塑性区可以进一步分为塑性硬化区和塑性软化区[13]

图1 井眼受力分析图 Fig. 1 Schematic of stress around the borehole
1.2 岩石本构 1.2.1 岩石应力—应变关系

在岩石破坏之前,岩石变形分为线弹性阶段和非线性塑性阶段,岩石破坏前的两个阶段的变形分别表示为[13]

$ \sigma =E\varepsilon , {\kern 25pt} 0\leqslant \varepsilon \leqslant {{\varepsilon }_{\rm s}} $ (1)
$ \sigma ={{E}_{0}}\left( \varepsilon {\rm e}^{ -\frac{\varepsilon }{{{\varepsilon }_{0}}} }-C\dfrac{{{\varepsilon }_{0}}}{\varepsilon } \right), {\kern 5pt} \varepsilon > {{\varepsilon }_{\rm s}} $ (2)

由式(2)对$\varepsilon$求导,可得

$ \dfrac{{\rm d}\sigma }{{\rm d}\varepsilon }={{E}_{0}}\left[ \left( 1-\dfrac{\varepsilon }{{{\varepsilon }_{0}}} \right){\rm e}^ { -\frac{\varepsilon }{{{\varepsilon }_{0}}} }+C\dfrac{{{\varepsilon }_{0}}}{{{\varepsilon }^{2}}} \right] $ (3)

式(1)和式(2)光滑连接的条件为

$ \dfrac{{\rm d}\sigma }{{\rm d}\varepsilon }=\dfrac{\sigma }{{{\varepsilon }_{\rm s}}}=E $ (4)

在式(2)和式(3)中,令$\varepsilon=\varepsilon_{{\rm s}}$,代入式(4),可得

$ C=\dfrac{\varepsilon _{\rm s}^{3}}{2\varepsilon _{0}^{2}}{\rm e}^{ -\frac{{{\varepsilon }_{\rm s}}}{{{\varepsilon }_{0}}} } $ (5)

由式(3)~(5)可得岩体的弹性模量为[13]

$ E={{E}_{0}}\left( 1-\dfrac{{{\varepsilon }_{\rm s}}}{{{\varepsilon }_{0}}} \right){\rm e}^{ -\frac{{{\varepsilon }_{\rm s}}}{{{\varepsilon }_{0}}} }=\dfrac{{{\sigma }_{\rm c}}{\rm }}{{{\varepsilon }_{\rm c}}}\left( 1-\dfrac{{{\varepsilon }_{\rm s}}}{{{\varepsilon }_{0}}} \right){\rm e}^{-\frac{{{\varepsilon }_{\rm s}}}{{{\varepsilon }_{0}}} } $ (6)

当岩石变形达到极限应变时有[13]

$ \dfrac{{\rm d}\sigma }{{\rm d}\varepsilon }\left| _{\varepsilon ={{\varepsilon }_{\rm c}}} \right.={{E}_{0}}\left[ \left( 1-\dfrac{{{\varepsilon }_{\rm c}}}{{{\varepsilon }_{0}}} \right){\rm e}^{-\frac{{{\varepsilon }_{\rm c}}}{{{\varepsilon }_{0}}}}+C\dfrac{{{\varepsilon }_{0}}}{\varepsilon _{\rm c}^{2}} \right]=0 $ (7)

计算中暂取$\varepsilon_{{\rm s}}$/$\varepsilon_{{\rm 0}}$=1/3,由式(5),可得

$ C={{\varepsilon }_{0}}/\left( 54\sqrt[3]{\rm e} \right) $ (8)

将式(8)代入式(7),可得

$ \begin{array}{l} {{{\varepsilon }_{\rm c}}}/{{{\varepsilon }_{0}}} =1.03487 \end{array} $ (9)
$ \begin{array}{l} C=0.01280{{\varepsilon }_{\rm c}} \end{array} $ (10)
$ \begin{array}{l} \dfrac{{{\varepsilon }_{\rm c}}-{{\varepsilon }_{0}}}{{{\varepsilon }_{0}}}=0.03487 \end{array} $ (11)

经过上述推导,考虑岩石弹性变形和非线性塑性变形的岩石本构为[13]

$ \left\{ \begin{array}{l} \sigma =E\varepsilon , {\kern 88pt} 0\leqslant \varepsilon \leqslant {{\varepsilon }_{\rm s}} \\[5pt] \sigma ={{E}_{0}}{{\varepsilon }_{\rm c}}\left( \dfrac{\varepsilon }{{{\varepsilon }_{\rm c}}}{\rm e} { -g\dfrac{\varepsilon }{{{\varepsilon }_{\rm c}}} }-h\dfrac{{{\varepsilon }_{\rm c}}}{\varepsilon } \right), {\kern 5pt} \varepsilon > {{\varepsilon }_{\rm s}} \\[5pt] g=\dfrac{{{\varepsilon }_{\rm c}}}{{{\varepsilon }_{0}}} \\[8pt] h=C\dfrac{{{\varepsilon }_{0}}}{\varepsilon _{\rm c}^{2}} \end{array}\right. $ (12)
1.2.2 塑性区全量理论本构关系

由于地层厚度大,因此,井周应力场计算可作平面应变问题处理,则$\varepsilon_{{ z}}$=0。

设岩石体积应变$\theta=3\varepsilon_{{ m}}=0$,有

$ \left\{ \begin{array}{l} {{\varepsilon }_{z}}=0 \\ 3{{\varepsilon }_{m}}={{\varepsilon }_{\theta }}+{{\varepsilon }_{r}}=0 \end{array}\right. $ (13)

在平面应变和体积应变为0的条件下,轴向应力$\sigma_{{ z}}$$\sigma_{\theta}$$\sigma_{r}$有如下关系[17]

$ \begin{array}{l} {{\sigma }_{z}}=\dfrac{1}{2}\left( {{\sigma }_{\theta }}+{{\sigma }_{r}} \right) \end{array} $ (14)

由此,可得井周岩石应力

$ \begin{array}{l} {{\sigma }_{i}}\!=\!\dfrac{1}{\sqrt{2}}\sqrt{{{\left( {{\sigma }_{\theta }}\!-\!{{\sigma }_{z}} \right)}^{2}}\!+\!{{\left( {{\sigma }_{z}}\!-\!{{\sigma }_{r}} \right)}^{2}}\!+\!{{\left( {{\sigma }_{r}}\!-\!{{\sigma }_{\theta }} \right)}^{2}}}\!=\!\\{\kern 40pt}\dfrac{\sqrt{3}}{2}\left( {{\sigma }_{r}}-{{\sigma }_{\theta }} \right) \end{array} $ (15)

由于井筒在水平方向上所受地应力相等,因此井筒受力问题为轴对称问题,其几何方程为

$ \left \{ \begin{array}{l} {{\varepsilon }_{\theta }}=\dfrac{u}{r}\\[6pt] {{\varepsilon }_{r}}=\dfrac{{\rm d}u}{{\rm d}r} \end{array} \right . $ (16)

将式(16)代入式(13)第二式并积分可得

$ \left \{ \begin{array}{l} u\left( r \right)=\dfrac{A}{r} \\[6pt] {{\varepsilon }_{\theta }}\left( r \right)=\dfrac{A}{{{r}^{2}}}\\[6pt] {{\varepsilon }_{r}}\left( r \right)=-\dfrac{A}{{{r}^{2}}} \end{array} \right . $ (17)

结合式(17)和式(13)可得井周岩石应变

$ {{\varepsilon }_{i}}\, =\, \dfrac{\sqrt{2}}{3}\sqrt{{{\left( {{\varepsilon }_{\theta }}\, -\, {{\varepsilon }_{z}} \right)}^{2}}\, +\, {{\left( {{\varepsilon }_{z}}\, -\, {{\varepsilon }_{r}} \right)}^{2}}\, +\, {{\left( {{\varepsilon }_{r}}\, -\, {{\varepsilon }_{\theta }} \right)}^{2}}}\, =\, \dfrac{2}{\sqrt{3}}\dfrac{A}{{{r}^{2}}} $ (18)

$r$等于软化区半径$R$$\varepsilon_{i}$=$\varepsilon_{\rm c}$,可得积分常数$A$

$ \begin{array}{l} A=\dfrac{\sqrt{3}}{2}{{R}^{2}}{{\varepsilon }_{\rm c}} \end{array} $ (19)

将式(19)分别代入式(17)和式(18),可得

$ \left \{ \begin{array}{l} u\left( r \right)=\dfrac{\sqrt{3}}{2}{{\varepsilon }_{\rm c}}\dfrac{{{R}^{2}}}{r}\\[6pt] {{\varepsilon }_{\theta }}=\dfrac{\sqrt{3}}{2}{{\varepsilon }_{\rm c}}\dfrac{{{R}^{2}}}{{{r}^{2}}}\\[6pt] {{\varepsilon }_{r}}=-\dfrac{\sqrt{3}}{2}{{\varepsilon }_{\rm c}}\dfrac{{{R}^{2}}}{{{r}^{2}}} \end{array} \right . $ (20)
$ \begin{array}{l} {{\varepsilon }_{i}}={{\varepsilon }_{\rm c}}\dfrac{{{R}^{2}}}{{{r}^{2}}}=\dfrac{2}{\sqrt{3}}\dfrac{u\left( r \right)}{r} \end{array} $ (21)

由式(13)知,$\varepsilon_{ z}$:$\varepsilon_{r}$:$\varepsilon_{\theta}$= 0:1:—1,围岩应变分量保持固定比例,属于简单加载。从而得到围岩岩体应力$\sigma_{i}$与应变$\varepsilon_{i}$之间的物理关系式[13]

$ \left \{ \begin{array}{l} {{\sigma }_{i}}=E{{\varepsilon }_{i}}, {\kern 85pt}0\leqslant {{\varepsilon }_{i}}\leqslant {{\varepsilon }_{\rm s}} \\ {{\sigma }_{i}}={{E}_{0}}{{\varepsilon }_{\rm c}}\left( \dfrac{{{\varepsilon }_{i}}}{{{\varepsilon }_{\rm c}}}{\rm e}^{ -g\frac{{{\varepsilon }_{i}}}{{{\varepsilon }_{\rm c}}}} -h\dfrac{{{\varepsilon }_{\rm c}}}{{{\varepsilon }_{i}}} \right), {\kern 7pt} {{\varepsilon }_{i}}> {{\varepsilon }_{\rm s}} \end{array} \right . $ (22)
2 井周应力场模型建立 2.1 弹性区应力场 2.1.1 不考虑渗滤作用的井周应力场

岩石达到塑性屈服阶段时,井周存在塑性区,塑性区外为弹性区,此时弹性区的应力分布为[18]

$ \left\{ \begin{array}{l} \sigma _{r}^{\rm e}={{\sigma }_{\rm h}}\left( 1-\dfrac{R_{\rm s}^{2}}{{{r}^{2}}} \right)+\dfrac{R_{\rm s}^{2}}{{{r}^{2}}}{{\sigma }_{{{R}_{\rm s}}}} \\[8pt] \sigma _{\theta }^{\rm e}={{\sigma }_{\rm h}}\left( 1+\dfrac{R_{\rm s}^{2}}{{{r}^{2}}} \right)-\dfrac{R_{\rm s}^{2}}{{{r}^{2}}}{{\sigma }_{{{R}_{\rm s}}}} \\ \end{array} \right. $ (23)
2.1.2 流体渗滤诱导应力场

单相不可压缩稳定渗流极坐标微分方程为[19]

$ \begin{array}{l} \dfrac{1}{r}\dfrac{\rm d}{{\rm d}r}\left( r\dfrac{{\rm d}p}{{\rm d}r} \right)=0 \end{array} $ (24)

边界条件

$ \begin{array}{l} \left\{ \begin{array}{l} p\left| _{r={{r}_{\rm w}}} \right.={{p}_{\rm w}}\\ p\left| _{r={{r}_{\rm e}}} \right.={{p}_{\rm s}} \end{array} \right . \end{array} $ (25)

对式(24)进行积分,并将边界条件式(25)代入,可得孔隙流体压力和压降沿半径分布为

$ \begin{array}{l} p={{p}_{\rm s}}-\left( {{p}_{\rm s}}-{{p}_{\rm w}} \right)\dfrac{\ln {{{r}_{\rm e}}}-\ln{r}}{\ln {{{r}_{\rm e}}}-\ln{{{r}_{\rm w}}}}, {\kern 8pt} {{r}_{\rm w}}\leqslant r\leqslant {{r}_{\rm e}} \end{array} $ (26)
$ \begin{array}{l} \dfrac{{\rm d}p}{{\rm d}r}=\dfrac{ {{p}_{\rm s}}-{{p}_{\rm w}} }{\ln {{{r}_{\rm e}}}-\ln{{{r}_{\rm w}}}}\dfrac{1}{r}, {\kern 20pt} {{r}_{\rm w}}\leqslant r\leqslant {{r}_{\rm e}} \end{array} $ (27)

渗流造成的孔弹性应力为

$ \left\{ \begin{array}{l} {{\sigma }_{r}}\!=\!F\left( \dfrac{{{r}^{2}}\!-\!r_{\rm w}^{2}}{r_{\rm e}^{2}\!-\!r_{\rm w}^{2}} {\int}_{{{r}_{\rm w}}}^{{{r}_{\rm e}}}{\Delta pr{\rm d}r\!-\! {\int}_{{{r}_{\rm w}}}^{r}{\!\!\Delta pr{\rm d}r}} \right)\\[8pt] {{\sigma }_{\theta }}\!=\!F\left( \dfrac{{{r}^{2}}\!+\!r_{\rm w}^{2}}{r_{\rm e}^{2}\!-\!r_{\rm w}^{2}} {\int}_{{{r}_{\rm w}}}^{{{r}_{\rm e}}}{\Delta pr{\rm d}r\!-\! {\int}_{{{r}_{\rm w}}}^{r}{\!\!\Delta pr{\rm d}r\!-\!\Delta p{{r}^{2}}}} \right)\\[8pt] F = \dfrac{1\!-\!2\upsilon }{\left( 1\!-\!\upsilon \right){{r}^{2}}} \end{array} \right . $ (28)

由式(26),有

$ \begin{array}{l} \Delta p\!=\!p\!-\!{{p}_{\rm s}}\!=\!\left( {{p}_{\rm w}}\!-\!{{p}_{\rm s}} \right)\dfrac{\ln {{{r}_{\rm e}}}\!-\!\ln{r}}{\ln {{{r}_{\rm e}}}\!-\!\ln{{{r}_{\rm w}}}}, {\kern 5pt}{{r}_{\rm w}}\!\leqslant \!r\!\leqslant\!{{r}_{\rm e}} \end{array} $ (29)

将式(29)代入式(28),积分,%由于$r_{\rm w} \leqslant r \leqslant r_{\rm e}$,可得

$ \left\{ \begin{array}{l} {{\sigma }_{r}}=\dfrac{\left( 1-2\upsilon \right)\left( {{p}_{\rm w}}-{{p}_{\rm s}} \right)}{2\left( 1-\upsilon \right)}\left[ 1-{{\left( \dfrac{{{r}_{\rm w}}}{r} \right)}^{2}} \right]\\[12pt] {{\sigma }_{\theta }}=\dfrac{\left( 1-2\upsilon \right)\left( {{p}_{\rm w}}-{{p}_{\rm s}} \right)}{2\left( 1-\upsilon \right)}\left[ 1+{{\left( \dfrac{{{r}_{\rm w}}}{r} \right)}^{2}} \right] \end{array} \right . $ (30)

因此,考虑渗滤作用耦合孔隙压力条件下的弹性区应力场为

$ \left\{ \begin{array}{l} \sigma _{r}^{\rm e}={{\sigma }_{\rm h}}\left( 1-\dfrac{R_{\rm s}^{2}}{{{r}^{2}}} \right)+\dfrac{R_{\rm s}^{2}}{{{r}^{2}}}{{\sigma }_{{{R}_{\rm s}}}}+\\{\kern 40pt}\dfrac{\left( 1-2\upsilon \right)\left( {{p}_{\rm w}}-{{p}_{\rm s}} \right)}{2\left( 1-\upsilon \right)}\left[ 1-{{\left( \dfrac{{{r}_{\rm w}}}{r} \right)}^{2}} \right] \\[12pt] \sigma _{\theta }^{\rm e}={{\sigma }_{\rm h}}\left( 1+\dfrac{R_{\rm s}^{2}}{{{r}^{2}}} \right)-\dfrac{R_{\rm s}^{2}}{{{r}^{2}}}{{\sigma }_{{{R}_{\rm s}}}}+\\{\kern 40pt}\dfrac{\left( 1-2\upsilon \right)\left( {{p}_{\rm w}}-{{p}_{\rm s}} \right)}{2\left( 1-\upsilon \right)}\left[ 1+{{\left( \dfrac{{{r}_{\rm w}}}{r} \right)}^{2}} \right] \end{array} \right . $ (31)
2.2 塑性区应力场 2.2.1 应力平衡方程

假设井周岩石为岩石力学性质各向同性的连续介质,考虑其渗透性。围岩中取一微分体,其应力状态如图 2所示。

图2 微分体应力状态 Fig. 2 Stress state of differential body

由应力平衡条件,有

$ \left( {{\sigma }_{r}}+\dfrac{\partial {{\sigma }_{r}}}{\partial r}{\rm d}r \right)\left( r+{\rm d}r \right){\rm d}\theta -{{\sigma }_{r}}r{\rm d}\theta +\left( {{\sigma }_{r\theta }}+\dfrac{\partial {{\sigma }_{r\theta }}}{\partial \theta }{\rm d}\theta \right){\rm d}r\cos \dfrac{{\rm d}\theta }{2}- \\{\kern 40pt}{{\sigma }_{r\theta }}{\rm d}r\cos \dfrac{{\rm d}\theta }{2}-{{\sigma }_{\theta }}{\rm d}r\sin \dfrac{{\rm d}\theta }{2}-\left( {{\sigma }_{\theta }}+\dfrac{\partial {{\sigma }_{\theta }}}{\partial \theta }{\rm d}\theta \right){\rm d}r\sin \dfrac{{\rm d}\theta }{2}+\dfrac{{\rm d}{{p}_{r}}}{{\rm d}r}r{\rm d}r{\rm d}\theta =0 $ (32)
$ \begin{array}{l} \left( {{\sigma }_{\theta }}+\dfrac{\partial {{\sigma }_{\theta }}}{\partial \theta }{\rm d}\theta \right){\rm d}r\cos \dfrac{{\rm d}\theta }{2}-{{\sigma }_{\theta }}{\rm d}r\cos \dfrac{{\rm d}\theta }{2}+\left( {{\sigma }_{r\theta }}+\dfrac{\partial {{\sigma }_{r\theta }}}{\partial \theta }{\rm d}\theta \right){\rm d}r\sin \dfrac{{\rm d}\theta }{2} + \\{\kern 40pt}{{\sigma }_{r\theta }}{\rm d}r\sin \dfrac{{\rm d}\theta }{2}+\left( {{\sigma }_{r\theta }}+\dfrac{\partial {{\sigma }_{r\theta }}}{\partial r}{\rm d}\theta \right)\left( r+{\rm d}r \right){\rm d}\theta -{{\sigma }_{r\theta }}r{\rm d}\theta +\dfrac{{\rm d}{{p}_{\theta }}}{{\rm d}\theta }r{\rm d}r{\rm d}\theta =0 \end{array} $ (33)

由于${\rm d}\theta$非常小,因此

$ \begin{array}{l} \left\{ \begin{array}{l} \sin \dfrac{{\rm d}\theta }{2}\approx \dfrac{{\rm d}\theta }{2} \\[8pt] \cos \dfrac{{\rm d}\theta }{2}\approx 1 \end{array} \right . \end{array} $ (34)

经整理并略去高阶无穷小项,则式(32)和式(33)可简化为

$ \begin{array}{l} \left\{ \begin{array}{l} \dfrac{\partial {{\sigma }_{r}}}{\partial r}+\dfrac{1}{r}\dfrac{\partial {{\sigma }_{r\theta }}}{\partial \theta }+\dfrac{{{\sigma }_{r}}-{{\sigma }_{\theta }}}{r}+\dfrac{{\rm d}{{p}_{r}}}{{\rm d}r}=0 \\[8pt] \dfrac{1}{r}\dfrac{\partial {{\sigma }_{\theta }}}{\partial \theta }+\dfrac{\partial {{\sigma }_{r\theta }}}{\partial r}+\dfrac{2{{\sigma }_{r\theta }}}{r}+\dfrac{{\rm d}{{p}_{\theta }}}{{\rm d}\theta }=0 \end{array} \right . \end{array} $ (35)

由于应力呈轴对称分布,应力分量与$\theta$无关。式(35)可简化为轴对称条件下考虑流体渗流耦合的应力平衡方程

$ \begin{array}{l} \dfrac{{\rm d}{{\sigma }_{r}}}{{\rm d}r}+\dfrac{{\rm d}p}{{\rm d}r}+\dfrac{{{\sigma }_{r}}-{{\sigma }_{\theta }}}{r}=0 \end{array} $ (36)
2.2.2 塑性区应力分布规律

令式(12)第二项中$\varepsilon=\varepsilon_{\rm s}$$\varepsilon=\varepsilon_{\rm c}$,可得

$ \begin{array}{l} \left\{ \begin{array}{l} {{\sigma }_{\rm s}}={{E}_{0}}{{\varepsilon }_{\rm c}}\left( \dfrac{{{\varepsilon }_{\rm s}}}{{{\varepsilon }_{\rm c}}}{\rm e}^ {-g\frac{{{\varepsilon }_{\rm s}}}{{{\varepsilon }_{\rm c}}} } -h\dfrac{{{\varepsilon }_{\rm c}}}{{{\varepsilon }_{\rm s}}} \right) \\ {{\sigma }_{\rm c}}={{E}_{0}}{{\varepsilon }_{\rm c}}\left( {{\rm e}^{-g}}-h \right) \end{array} \right . \end{array} $ (37)

利用式(21)和式(37),可将式(22)中的$\sigma_{i}$写为

$ \begin{array}{l} {{\sigma }_{i}}=\dfrac{{{\sigma }_{\rm c}}}{ {{\rm e}^{-g}}-h }\left( \dfrac{{{R}^{2}}}{{{r}^{2}}}{\rm e}^ { -g\frac{{{R}^{2}}}{{{r}^{2}}} }-h\dfrac{{{r}^{2}}}{{{R}^{2}}} \right) \end{array} $ (38)

结合式(27)、式(15)、式(36)和式(38),可得

$ \begin{array}{l} \dfrac{{\rm d}{{\sigma }_{r}}}{{\rm d}r}=-\dfrac{{\rm d}p}{{\rm d}r}-\dfrac{{{\sigma }_{r}}-{{\sigma }_{\theta }}}{r} =-\dfrac{ {{p}_{\rm s}}-{{p}_{\rm w}} }{\ln {{{r}_{\rm e}}}-\ln{{{r}_{\rm w}}}}\dfrac{1}{r}- \dfrac{2}{\sqrt{3}}\dfrac{{{\sigma }_{\rm c}}}{\left( {{\rm e}^{-g}}-h \right)}\left( \dfrac{{{R}^{2}}}{{{r}^{3}}}{\rm e}^ { -g\frac{{{R}^{2}}}{{{r}^{2}}} }-h\dfrac{r}{{{R}^{2}}} \right) \end{array} $ (39)

对式(39)进行积分,且在井壁处$\sigma_{{\rm r}}=p_{{\rm w}}$,可得

$ \begin{array}{l} \sigma _{r}^{p}={{p}_{\rm w}}+\dfrac{ {{p}_{\rm s}}-{{p}_{\rm w}} }{\ln {{{r}_{\rm e}}}-\ln{{{r}_{\rm w}}}}\ln \dfrac{{{r}_{\rm w}}}{r} -\dfrac{1}{\sqrt{3}}\dfrac{{{\sigma }_{\rm c}}}{\left( {{\rm e}^{-g}}-h \right)}\left( \dfrac{1}{g}{\rm e}^{ -g\frac{{{R}^{2}}}{{{r}^{2}}} }-\dfrac{1}{g}{\rm e}^ { -g\frac{{{R}^{2}}}{r_{\rm w}^{2}} }+h\dfrac{r_{\rm w}^{2}-{{r}^{2}}}{{{R}^{2}}} \right) \end{array} $ (40)

将式(40)代入式(15),可得

$ \begin{array}{l} \sigma _{\theta }^{p}={{p}_{\rm w}}+\dfrac{ {{p}_{\rm s}}-{{p}_{\rm w}} }{\ln {{{r}_{\rm e}}}-\ln{{{r}_{\rm w}}}}\ln \dfrac{{{r}_{\rm w}}}{r} -\dfrac{1}{\sqrt{3}}\dfrac{{{\sigma }_{\rm c}}}{\left( {{\rm e}^{-g}}-h \right)}\left[\dfrac{1}{g}\left( {\rm e}^ { -g\frac{{{R}^{2}}}{{{r}^{2}}} }-{\rm e}^ { -g\frac{{{R}^{2}}}{r_{\rm w}^{2}} } \right)+h\dfrac{r_{\rm w}^{2}-3{{r}^{2}}}{{{R}^{2}}}+\dfrac{2{{R}^{2}}}{{{r}^{2}}}{\rm e}^ { -g\frac{{{R}^{2}}}{{{r}^{2}}} } \right] \\ \end{array} $ (41)

塑性硬化区和软化区应力可同时由式(40)、式(41)来表示。在$r=R_{{\rm s}}$处,塑性硬化区和弹性区应力连续,即

$ \begin{array}{l} \left\{ \begin{array}{l} \sigma _{\theta }^{\rm e}+\sigma _{r}^{\rm e}=\sigma _{\theta }^{p}+\sigma _{r}^{p} \\ \sigma _{\theta }^{\rm e}=\sigma _{\theta }^{p}\\ \sigma _{r}^{\rm e}=\sigma _{r}^{p} \end{array} \right. \end{array} $ (42)

结合式(31)、式(40)、式(41)和式(42),可得

$ \begin{array}{l} 2{{\sigma }_{\rm h}}\!+\!\dfrac{\left( 1\!-\!2\upsilon \right)\left( {{p}_{\rm w}}\!-\!{{p}_{\rm s}}\right) }{ 1\!-\!\upsilon }\!=\!2\left( {{p}_{\rm w}}\!+\!\dfrac{ {{p}_{\rm s}}\!-\!{{p}_{\rm w}} }{\ln \dfrac{{{r}_{\rm e}}}{{{r}_{\rm w}}}}\ln \dfrac{{{r}_{\rm w}}}{{{R}_{\rm s}}} \right)\!-\\[8pt]{\kern 40pt} \left[ \dfrac{1}{g}\left( {\rm e}^{ -g\frac{{{R}^{2}}}{R_{\rm s}^{2}} }\!-\!{\rm e}^ { -g\frac{{{R}^{2}}}{r_{\rm w}^{2}} } \right)\!+\!h\dfrac{r_{\rm w}^{2}\!-\!2R_{\rm s}^{2}}{{{R}^{2}}}\!+\!\dfrac{{{R}^{2}}}{R_{\rm s}^{2}}{\rm e}^ { -g\frac{{{R}^{2}}}{R_{\rm s}^{2}} } \right]\cdot \\[8pt]{\kern 40pt} \dfrac{2}{\sqrt{3}}\dfrac{{{\sigma }_{\rm c}}}{\left( {{\rm e}^{-g}}-h \right)} \end{array} $ (43)

由式(21),可知

$ \begin{array}{l} \dfrac{{{\varepsilon }_{\rm s}}}{{{\varepsilon }_{\rm c}}}=\dfrac{{{R}^{2}}}{R_{\rm s}^{2}} \end{array} $ (44)

将式(44)代入式(43),可得

$ \begin{array}{l} 2{{\sigma }_{\rm h}}\!+\!\dfrac{\left( 1\!-\!2\upsilon \right)\left( {{p}_{\rm w}}\!-\!{{p}_{\rm s}}\right) }{ 1\!-\!\upsilon }\!=\!2\left( {{p}_{\rm w}}\!+\!\dfrac{ {{p}_{\rm s}}\!-\!{{p}_{\rm w}} }{\ln \dfrac{{{r}_{\rm e}}}{{{r}_{\rm w}}}}\ln \dfrac{{{r}_{\rm w}}}{{{R}_{\rm s}}} \right)\!-\\[8pt]{\kern 40pt} \left[ \dfrac{{\rm e}^ { -g\frac{{{\varepsilon }_{\rm s}}}{{{\varepsilon }_{\rm c}}} }}{g}\!-\!\dfrac{{\rm e}^ { -g\frac{{{\varepsilon }_{\rm s}}}{{{\varepsilon }_{\rm c}}}\frac{R_{\rm s}^{2}}{r_{\rm w}^{2}} }}{g} \!+\!h\left( \dfrac{{{\varepsilon }_{\rm c}}}{{{\varepsilon }_{\rm s}}}\dfrac{r_{\rm w}^{2}}{R_{\rm s}^{2}}\!-\!\dfrac{2{{\varepsilon }_{\rm c}}}{{{\varepsilon }_{\rm s}}} \right)\!+\!\dfrac{{{\varepsilon }_{\rm s}}}{{{\varepsilon }_{\rm c}}}{\rm e}^ { -g\frac{{{\varepsilon }_{\rm s}}}{{{\varepsilon }_{\rm c}}} } \right] \cdot \\[8pt]{\kern 40pt} \dfrac{2}{\sqrt{3}}\dfrac{{{\sigma }_{\rm c}}}{\left( {{\rm e}^{-g}}-h \right)} \end{array} $ (45)

由式(45)可求得塑性硬化区半径$R_{{\rm s}}$,再结合式(44)可求得塑性软化区半径$R$,再将$R_{{\rm s}}$$R$代入式(31),式(40)和式(41)则可求出弹性区、塑性硬化区和塑性软化区的应力场。

3 井周应力场分布规律研究 3.1 塑性区范围影响因素研究 3.1.1 井眼流体压力对塑性区范围影响分析

表 1中的参数代入本文所建立模型,得到塑性区范围随井眼流体压力的变化关系,如图 3所示。由于水平方向上远场地应力相等,井周任意角度上的变形和应力状态是相同的,因此,后文图中可表示任意井周角上岩石应力应变状态。从图 3可知,井眼流体压力越大,塑性硬化区半径和塑性软化区半径均越大,即井眼内外压力差越大,塑性区越大。

表1 井眼流体压力对塑性区范围影响分析计算参数 Tab. 1 Calculation parameters for the influence of the wellbore fluid pressure on the range of plastic zone
图3 井眼流体压力对塑性区范围影响 Fig. 3 The effects of wellbore fluid pressure on the plastic zone
3.1.2 岩石强度对塑性区范围影响分析

表 2中的参数代入本文建立的模型,得到塑性区范围随岩石强度的变化关系,如图 4所示。从图 4可知,相同条件下,岩石强度越大,塑性硬化区半径和塑性软化区半径都越小,即塑性区越小。由此可知,岩石强度越小,越容易发生塑性变形。

表2 岩体强度对塑性区范围影响分析计算参数 Tab. 2 Calculation parameters for the influence of rock strength on the range of plastic zone
图4 岩石强度对塑性区范围影响 Fig. 4 The effects of rock strength on the plastic zone
3.1.3 泊松比对塑性区范围影响分析

表 3中的参数代入本文所建立的模型,得到塑性区范围随泊松比的变化关系,如图 5所示。从图 5可知,相同条件下,泊松比越大,塑性硬化区半径和塑性软化区半径都越大,即塑性区越大,岩石越容易发生塑性变形。

表3 泊松比对塑性区范围影响分析计算参数 Tab. 3 Calculation parameters for the influence of Poisson ratio on the range of plastic zone
图5 泊松比对塑性区范围影响 Fig. 5 The effects of Poisson ratio on the plastic zone
3.1.4 地层初始孔隙流体压力对塑性区范围影响

表 4中的参数代入本文所建立的模型,得到塑性区范围随地层初始孔隙流体压力的变化关系,如图 6所示。从图 6可知,相同条件下,地层初始孔隙流体压力越大,塑性硬化区半径和塑性软化区半径都越大,即塑性区越大,岩石越容易发生塑性变形。

表4 地层初始孔隙流体压力对塑性区范围影响计算参数 Tab. 4 Calculation parameters for the influence of initial formation pore pressure on the range of plastic zone
图6 地层初始孔隙流体压力对塑性区范围影响 Fig. 6 The effects of initial formation pore pressure on the range of plastic zone
3.2 井周应力场影响因素研究 3.2.1 井眼流体压力对井周应力影响分析

表 1中的数据代入本文建立的模型,计算井眼流体压力分别为35,40和45~MPa时井周应力分布,如图 7~图 9所示。

图7 井眼流体压力对井周周向应力分布影响 Fig. 7 The effects of wellbore fluid pressure on the circumferential stress
图8 井眼流体压力对井周径向应力分布影响 Fig. 8 The effects of wellbore fluid pressure on the radial stress
图9 井眼流体压力对井周应力强度分布影响 Fig. 9 The effects of wellbore fluid pressure on the stress intensity

图 7可知,井眼压力为35~MPa时井周还未发生塑性变形,随着$r/r_{{\rm w}}$增加,周向应力增加,周向受压程度增加,即弹性变形条件下,周向应力最小值在井壁处。井眼压力为40和45~MPa时,井周发生了塑性变形,岩石塑性变形使得井周周向应力发生了变化,随着$r/r_{{\rm w}}$增加,周向应力先减小后增加,且最小值在塑性软化区外边界处。

图 8可知,随着$r/r_{{\rm w}}$增加,径向应力减小。且3条曲线均未出现明显拐点,即岩石塑性变形对径向应力分布影响不明显。

图 9可知,井眼压力为35和40~MPa时,随着$r/r_{{\rm w}}$增加,岩石应力强度减小。井眼压力为45~MPa时,随着$r/r_{{\rm w}}$增加,岩石应力强度先增加后减小。和岩石弹性变形相比,岩石弹塑性变形减小了井壁岩石应力强度。

3.2.2 岩体强度对井周应力影响分析

表 2中的数据代入本文建立的模型,计算岩体强度分别为15,20和25~MPa时井周应力分布,如图 10~图 12所示。

图10 岩体强度对井周周向应力分布影响 Fig. 10 The effects of rock strength on the circumferential stress
图11 岩体强度对井周径向应力分布影响 Fig. 11 The effects of rock strength on the radial stress
图12 岩体强度对井周应力强度分布影响 Fig. 12 The effects of rock strength on the stress intensity

图 10~图 12可知,岩体强度为25~MPa时,井周未发生塑性变形,当发生塑性变形时,岩体强度越小,塑性区周向应力和径向应力越大,应力强度越小。

3.2.3 泊松比对井周应力影响分析

表 3中的数据代入本文建立的模型,计算泊松比分别为0.1,0.2和0.3时井周应力分布,如图 13~图 15所示。

图13 泊松比对井周周向应力分布影响 Fig. 13 The effects of Poisson ratio on the circumferential stress
图14 泊松比对井周径向应力分布影响 Fig. 14 The effects of Poisson ratio on the radial stress
图15 泊松比对井周应力强度分布影响 Fig. 15 The effects of Poisson ratio on the stress intensity

图 13~图 15可知,泊松比对井周应力分布影响很大,这是由于考虑了孔隙流体压力作用,从而导致泊松比对弹性区应力场造成影响,进而对塑性区应力场造成影响,但在井壁处,泊松比对应力场没有影响。泊松比对井周应力场的影响具体表现为:泊松比越大,井周周向应力和径向应力越小,井周应力强度越大。

3.2.4 地层初始孔隙流体压力对井周应力影响

表 4中的数据代入本文建立的模型,计算地层初始孔隙流体压力分别为10,15和20~MPa时井周应力分布,如图 16~图 18所示。

图16 地层初始孔隙流体压力对井周周向应力分布影响 Fig. 16 The effects of initial formation pore pressure on the circumferential stress
图17 地层初始孔隙流体压力对井周径向应力分布影响 Fig. 17 The effects of initial formation pore pressure on the radial stress
图18 地层初始孔隙流体压力对井周应力强度分布影响 Fig. 18 The effects of initial formation pore pressure on the stress intensity

图 16~图 18可知,地层初始孔隙流体压力对井周应力场分布影响很大。具体表现为:地层初始孔隙流体压力越大,井周周向应力和径向应力都越小,即受压程度减小。这是由于孔隙压力越大,岩石骨架所承受的压应力越小。而初始孔隙压力越大,离井壁5倍井眼半径距离内井周应力强度越大,超过5倍井眼半径距离过后,初始孔隙压力对岩石应力强度影响不大。

4 结论

(1) 根据弹塑性地层力学特征,考虑工作液渗滤作用和岩石非线性硬化与软化变形,建立了耦合工作液渗滤和岩石非线性硬化与软化变形的井周应力场模型。

(2) 计算分析表明,井周塑性区大小是受岩石力学性质以及岩石所受外力所影响的,具体而言,塑性区大小与井眼流体压力、泊松比和地层初始孔隙流体压力呈正相关关系,与岩石强度呈反相关关系。因此现场应用时应根据实际地层参数估计塑性区大小。

(3) 岩石塑性变形对井周径向应力影响较小,对周向应力和岩石应力强度影响较大,具体表现为,与弹性状态相比,塑性变形后周向应力增大(受压程度增加),岩石发生塑性变形会降低岩石的应力强度,即塑性变形降低了井眼应力集中效应。

符号说明

参考文献
[1]
MOHAMMAD E Z. Mechanical and physico-chemical aspects of wellbore stability during drilling operations[J]. Journal of Petroleum Science and Engineering, 2012, 82-83: 120-124. doi: 10.1016/j.petrol.2012.01.006
[2]
杨兆中, 刘云锐, 张平, 等. 煤层气直井地层破裂压力计算模型[J]. 石油学报, 2018, 39(5): 578-586.
YANG Zhaozhong, LIU Yunrui, ZHANG Ping, et al. A model for calculating formation breakdown pressure in CBM vertical wells[J]. Acta Petrolei Sinica, 2018, 39(5): 578-586. doi: 10.7623/syxb201805009
[3]
马天寿, 陈平. 页岩层理对水平井井壁稳定的影响[J]. 西南石油大学学报(自然科学版), 2014, 36(5): 97-104.
MA Tianshou, CHEN Ping. Influence of shale bedding plane on wellbore stability for horizontal wells[J]. Journal of Southwest Petroleum University (Science & Technology Edition), 2014, 36(5): 97-104. doi: 10.11885/j.issn.1674-5086.2013.06.30.03
[4]
曾凡辉, 唐波涛, 王涛, 等. 考虑渗滤效应的压裂裸眼井破裂压力预测模型[J]. 天然气地球科学, 2019, 30(4): 549-556.
ZENG Fanhui, TANG Botao, WANG Tao, et al. Prediction model of fracture initiation pressure of open hole well considering penetration effect[J]. Natural Gas Geoscience, 2019, 30(4): 549-556. doi: 10.11764/j.issn.1672-1926.2019.01.010
[5]
HUANG J, GRIFFITHS D V, WONG S W. Initiation pressure, location and orientation of hydraulic fracture[J]. International Journal of Rock Mechanics and Mining Sciences, 2012, 49(2): 59-67. doi: 10.1016/j.ijrmms.2011.11.014
[6]
李玉梅, 苏中, 张涛, 等. 欠平衡钻井煤层井壁稳定有限元数值计算研究[J]. 系统仿真学报, 2018, 30(11): 226-232.
LI Yumei, SU Zhong, ZHANG Tao, et al. Numerical simulation of borehole stability in underbalanced drilling coal seam[J]. Journal of System Simulation, 2018, 30(11): 226-232. doi: 10.16182/j.issn1004731x.joss.201811025
[7]
张鹏伟, 孙峰, 叶贵根, 等. 基于塑性损伤理论的软泥页岩地层井壁稳定性分析[J]. 力学与实践, 2018, 40(3): 273-280.
ZHANG Pengwei, SUN Feng, YE Guigen, et al. Analysis of the wellbore stability based on plastic damage theory in soft shale formation[J]. Mechanics in Engineering, 2018, 40(3): 273-280. doi: 10.6052/1000-0879-17-369
[8]
AADNOY B S. Geomechanical analysis for deep-water drilling[C]. SPE 39339-MS, 1998. doi: 10.2118/39339-MS
[9]
AADNOY B S, BELAYNEH M. Elasto-plastic fracturing model for wellbore stability using non-penetrating fluids[J]. Journal of Petroleum Science and Engineering, 2004, 45(3): 179-192. doi: 10.1016/j.petrol.2004.07.006
[10]
WANG Y, DUSSEAULT M B. Borehole yield and hydraulic fracture initiation in poorly consolidated rock strata-part Ⅰ. Impermeable media[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1991, 28: 235-246. doi: 10.1016/0148-9062(91)90591-9
[11]
WANG Y, DUSSEAULT M B. Borehole yield and hydraulic fracture initiation in poorly consolidated rock strata-part Ⅱ. Permeable media[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts,, 1991, 28: 247-260. doi: 10.1016/0148-9062(91)90592-A
[12]
潘岳, 王志强, 吴敏应. 非线性硬化与非线性软化的巷、隧道围岩塑性分析[J]. 岩土力学, 2006(7): 10381042.
PAN Yue, WANG Zhiqiang, WU Minying. Plastic analysis of surrounding rock of tunnel based on strain nonlinear hardening and nonlinear softening[J]. Rock and Soil Mechanics, 2006(7): 1038-1042. doi: 10.3969/j.issn.1000-7598.2006.07.004
[13]
潘岳, 王志强, 王在泉. 非线性硬化与软化的巷道围岩应力分布与工况研究[J]. 岩石力学与工程学报, 2006(7): 1343-1351.
PAN Yue, WANG Zhiqiang, WANG Zaiquan. Study on stress distribution of surrounding rock and work condition of tunnel based on strain nonlinear hardening and softening[J]. Chinese Journal of Rock Mechanics and Engineering, 2006(7): 1343-1351. doi: 10.3321/j.issn:1000-6915.2006.07.008
[14]
郭建春, 何颂根, 邓燕. 弹塑性地层水力压裂起裂模式及起裂压力研究[J]. 岩土力学, 2015, 36(9): 2494-2550.
GUO Jianchun, HE Songgen, DENG Yan. Study of hydraulic fracturing initiation mode and initiation pressure of elastoplastic formation[J]. Rock and Soil Mechanics, 2015, 36(9): 2494-2550. doi: 10.16285/j.rsm.2015.09.008
[15]
SHAO J F, HOMAND S. Influences of heat convection and plastic deformation on hydraulic fracture initiation and reservoir compaction[J]. SPE/ISRM Rock Mechanics in Petroleum Engineering, 1998. doi: 10.2118/47387-MS
[16]
XIA H W, MOORE I D. Estimation of maximum mud pressure in purely cohesive material during directional drilling[J]. Geomechanics and Geoengineering:An International Journal, 2006, 1(1): 3-11. doi: 10.1080/17486020600604024
[17]
何颂根.具有弹塑性特征的地层水力压裂起裂模型研究[D].成都: 西南石油大学, 2014.
HE Songgen. Study of hydraulic fracturing initiation mode of elastoplastic formation[D]. Chengdu: Southwest Petroleum University, 2014. http://cdmd.cnki.com.cn/Article/CDMD-10615-1014409749.htm
[18]
徐芝纶. 弹性力学(上)[M]. 北京: 高等教育出版社, 2006.
XU Zhilun. Elasticity (Part 1)[M]. Beijing: Higher Education Press, 2006.
[19]
李晓平. 地下油气渗流力学[M]. 北京: 石油工业出版社, 2008.
LI Xiaoping. Percolation mechanics of underground oil and gas[M]. Beijing: Petroleum Industry Press, 2008.