引 言

在油气田开发过程中，井周应力场的准确预测非常重要。因为这是钻井液安全密度窗口预测和压裂破裂压力预测的基础^[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.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所示。

由应力平衡条件，有

$ \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可知，井眼流体压力越大，塑性硬化区半径和塑性软化区半径均越大，即井眼内外压力差越大，塑性区越大。

3.1.2 岩石强度对塑性区范围影响分析

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

3.1.3 泊松比对塑性区范围影响分析

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

3.1.4 地层初始孔隙流体压力对塑性区范围影响

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

3.2 井周应力场影响因素研究 3.2.1 井眼流体压力对井周应力影响分析

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

从图 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~图 12可知，岩体强度为25~MPa时，井周未发生塑性变形，当发生塑性变形时，岩体强度越小，塑性区周向应力和径向应力越大，应力强度越小。

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

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

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

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

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

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

4 结论

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

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

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

符号说明