基于多孔弹性力学理论的深井井底应力场分析

呼怀刚, 管志川, 许玉强, 韩 超, 刘永旺, 梁德阳, 路保平

(1.中国石油大学(华东)石油工程学院,山东青岛 266580; 2.山东省深地钻井过程控制工程技术研究中心,山东青岛 266580; 3.中国石化石油工程技术研究院,北京 100101)

钻进引起的井底岩石应力的重新分布可能会造成地层无法承受由此产生的偏应力[1-2],从而造成可能的强度破坏,因此对井底岩石应力分布的研究将有助于对井下复杂情况、岩石破坏机制等进行充分了解。近年来相关问题的研究多集中于井壁应力分布以及井壁稳定性的影响因素分析[1,3-12],研究多采用弹塑性力学和孔隙弹性力学相关理论分析井壁失稳的力学机制,而弹塑性力学中采用位移法、应力法求解岩石的位移场与应力场[13],忽略了孔隙流体的运移,认为系统处于静态平衡,应力、应变和位移是空间坐标的函数,与时间无关,而得出的应力场分布可看做是井眼开挖一段时间以后的稳定解。关于井底应力场分布方面的研究相对较少,Warren等[14]研究了井底岩石平均应力变化对机械钻速的影响;Ito等[15]研究了井底应力集中现象对井壁坍塌的影响;文献[16]～[19]中考虑流固耦合作用将井底按照应力大小分为三向拉应力区、三向压应力区和拉压应力区。以上研究忽略了井底应力的动态演化过程,时间效应是影响应力分布的重要因素,评价井底岩石强度破坏情况的指标存在不确定性。笔者基于Biot多孔弹性力学理论,建立考虑时间效应的渗流场与应力场全耦合模型,结合三维摩尔库伦准则建立井底岩石强度破坏的预测模型,分析井底应力状况对岩石强度破坏的影响。

1 多孔弹性介质控制方程

1.1 井底岩石应力场控制方程

将井底岩石视为理想的多孔弹性介质,则岩体骨架的有效应力[20]为

σeij=2Gεij+λϑδij-αδijpp.

(1)

其中

式中,σeij为岩石骨架的有效应力,MPa,i,j=1,2,3;G为剪切模量,MPa;E为弹性模量,MPa;ν为泊松比;εij为应变张量分量;λ为拉梅常数;ϑ为体积应变;εx、εy、εz为3个方向的正应变;u、v、w为3个方向的位移;α为Biot系数(有效应力系数);δij为Kronecher函数;pp为孔隙压力,MPa。

1.2 流体渗流场控制方程

钻井过程中导致的井底岩石骨架变形以及开挖效应导致井眼周围岩石上覆岩层压力急剧减小都会影响孔隙压力的重新分布,同时地层中孔隙压力变化又会反过来影响岩石骨架应力重新分布[20-21]。

假设:①单相不可压缩流体;②流体渗流符合线性渗流规律;③地层岩石均质微可压缩;④弹性不稳定渗流;⑤等温渗流过程;⑥考虑重力等体积力对渗流的影响。得到的渗流方程为

(2)

在已知边界条件和初始条件后,就可得到井底渗流场与应力场相互影响的耦合场。

1.3 方程求解需要处理的问题

1.3.1 初始地应力场的平衡

深井地层在上覆岩层压力下所产生的初始位移场是历史时期岩体沉降固结形成的,对钻进过程并没有影响。将耦合方程的计算过程分为两步,首先求取上覆岩层压力所产生的应力场,将此应力场以第二类皮奥拉-基尔霍夫应力的形式添加到应力张量中,并以此作为第二步的初始条件计算钻进过程对现有位移场的影响[22]。

1.3.2 交叉耦合项

应力张量分量σij和渗透张量分量kij之间的相互作用是岩石应力场和渗流场耦合的关键。岩石的应力状态σij受渗流而产生的动水压力作用,从而有σij=f1(kij);另一方面应力张量分量σij变化产生的岩体变形会使渗透率改变,从而又有kij=f2(σij),Louis[23]通过钻进实验建立了岩石平均渗透系数k与有效正应力之间的经验关系:

(3)

式中,k0为初始渗透系数;σe为有效应力;α′为经验系数,取决于岩石的孔隙状态。

假设应力主方向与渗流主方向一致[24],计算过程中渗透系数矩阵设为

(4)

式中,λ为影响系数,由实验确定;σei为有效正应力,i=1,2,3。

1.3.3 重力效应

计算过程中考虑重力效应,所有边界载荷以及流体压力均为纵向坐标的函数。文中分析对象仅限于井底周围小区域内,考虑孔隙中液体惯性效应所引起的压力波动。

2 模型建立

2.1 几何模型

假设:井底岩石为各向同性多孔线弹性材料,受到均匀的水平地应力;孔隙骨架充满流体,且孔隙压力随岩石骨架变形而变化,孔隙压力变化又会对岩石骨架应力产生影响;不考虑井斜、温度、恒定径向流动引起的压力变化;井壁具有良好的渗透性,即井筒液柱压力与地层孔隙压力能够相互传递。

将深部地层岩石的受力情况进行简化,如图1(a)所示,井周及井底岩石受到上覆岩层压力σV、最大及最小水平地应力σH和σh、地层孔隙压力pp和液柱压力ph作用。为方便分析做井筒剖面,如图1(b)所示,在模型外围设置无限元域模拟无限大地层,减小模型尺寸对计算结果的影响。

图1 深部地层岩石物理模型Fig.1 Deep stratigraphic petro-physical model

2.2 边界条件

位移边界条件:在模型底面添加法向约束。

渗流边界包括水头边界、无流动边界和透水层边界。水头边界为井壁施加液柱压力ph;无流动边界为模型四周与底部边界,在模型四周设置远场地层孔隙压力pp;井筒井壁为透水边界层。

应力边界条件:模型顶部设置上覆岩层压力σV;模型四周设置均匀水平地应力σH和σh,并认为水平地应力由上覆岩层压力引起,即σH=σh=ν/(1-ν)(σV-pp)+pp;考虑渗透作用所产生的渗透压力,并将其转换为渗透体力施加到对应的节点和单元之上。

2.3 岩心物性参数及模型计算参数

为了使模型计算结果更加接近实际,模型中的计算参数使用塔里木油田某探井的岩心数据(取芯深度5 300 m),岩心物性参数:弹性模量E为26 308 MPa;泊松比ν为0.291;渗透系数k为1.1×10-8m/s;孔隙度φ为0.1;岩样密度ρm为2 680 kg/m3;岩石体积模量Ks为67 000 MPa;剪切模量G为9 843 MPa;凝聚力S0为25.87 MPa;摩擦角φ为25.07 °;孔隙流体密度ρl为1 000 kg/m3;孔隙流体比重γw为1×104N/m3;流体体积模量Kl为3 270 MPa。

假设物理模型位于井深5 500 m处,井眼直径为311.2 mm,钻井液密度为1 100 kg/m3。得模型计算参数:上覆岩层压力为144.452 MPa;井底液柱压力为59.29 MPa;地层孔隙压力为53.9 MPa;均匀水平地应力(σH=σh)为91.07 MPa。

3 深井井底应力场分析

井底应力场分析的主要目的是了解井底岩石在钻进之后的强度分布情况,进而对井底可能发生的复杂情况有更清楚的认识。

3.1 模型验证

使用Warren等[14]建立的井底应力场模型对文中所建立数值模型的准确性进行验证。Warren等通过类比温度应变与孔隙压力应变响应的相似性,建立了井底应力场计算模型,研究了井底应力特征对机械钻速的影响,并与现场试验结果进行了对比,展现了较好的相似性。Warren模型的主要参数:井深h为3 048 m;井径Φ为222 mm;钻井液密度ρf为1 050 kg/m3;上覆岩层压力σV为68.95 MPa;井底液柱压力ph为32.41 MPa;孔隙压力pp为32.41 MPa;均匀水平地应力(σH=σh)为48.26 MPa;弹性模量E为13 789.5 MPa;泊松比ν为0.25;渗透率k0为1 μm2;孔隙度φ为0.15;岩样密度ρm为2 262 kg/m3;岩石基体压缩系数CR为2.697×10-5MPa-1;钻头轮廓代号为IADC5-3-7bit。

如图2所示,文中计算得到的井底平均有效应力与Warren计算结果较为一致,说明文中所建立的模型的准确性,因此模型用于分析井底应力响应以及对岩石损伤程度的影响具有可行性。

图2 渗透性砂岩中井底附近的平均有效应力Fig.2 Mean effective stress around a borehole in a permeable sandstone

3.2 多孔弹性介质的时间效应

在无孔隙、线弹性材料中,由钻进或井筒液柱压力改变所致的应力重新分布是即时的,但在多孔可渗透弹性介质中,孔隙压力的传播只能以一定速度进行,此时应力重新分布是动力学问题,具有时间效应。多孔弹性介质的时间效应主要有两个来源:井筒钻井液侵入所导致的孔隙压力改变;钻进过程所致地层应力的重新分布会产生体积应变所引起的孔隙压力的变化。Detournay[1]等和Charlez[25]详细讨论了直井、非均匀远场应力下井壁应力分布的时间效应,但是对于井底应力的研究尚未进行。

图3为孔隙压力pp在井底中心线上随时间的分布情况(利用孔隙压力pp与原始地层压力p0的比值定量衡量pp的变化),离井底越近压力梯度越大。在打开井眼的极短时间内(10-2、10-4s),近井底处应变释放且钻井液未充分侵入井底,pp会突然下降;随轴向距离增加,pp变化幅度减小且其值增加。随着时间增加(102s、1 h),应力重新分布导致地层岩石孔隙压缩且其内流体排出速度不足,pp值上升,6 h后趋于稳定,且比原始孔隙压力要大。

图3 井底中心线上孔隙压力随时间变化Fig.3 Changes in pore pressure along the downhole centerline

3.3 井底应力场

计算模型采用各向同性结构、水平地应力相等且井眼轴线沿着某一主应力的方向。使用各强度准则对岩石强度进行衡量时均需要用到有效应力,

σeij=σij-αpp.

(5)

由式(5)所给出的有效应力的概念是基于Biot常数,且方程推导时假设岩石为线弹性特性,并不能直接用于岩石损伤判别。将Biot常数α设置为1,式(5)变为广泛接受的Terzaghi有效应力[2,26-27],且能够用于岩石损伤准则,即

σeij=σij-pp.

(6)

3.3.1 井底面岩石应力分布

井底面上的pp与ph保持一致;岩石强度由有效应力控制,从而衡量岩石内一点的应力状态需要同时用到pp、σr、σθ和σz。如图4所示,井底面上有效径向σer、周向σeθ及轴向应力σez随离井眼中心线距离的增加,其值逐渐上升至极大值,且三者极大值出现的位置会逐渐向井壁移动。σer和σeθ随时间变化其值波动较大(图4(a)、(b)),σez随时间波动较小;各有效应力值在100 s时达到最大值,在较长时间(24 h)后逐渐变为最小值。井底面上靠近井眼轴线的位置,σer和σeθ出现负值(表示拉应力),井眼钻开后随着时间的增长,拉应力会越加明显;σez随时间的变化规律比较一致,都保持在约0 MPa,拉应力并不明显。

图4 井底有效应力分布Fig.4 Distribution of effective stress on bottom hole

3.3.2 井底面以下岩石应力分布

为了方便研究井底面以下岩石的应力分布情况,定义图5所示4个位置:以井底面作为对象分别向下偏移50、100、156和250 mm,形成了位置1、2、3、4。对以上位置处的σer、σeθ和σez分布进行计算,结果如图6～8所示。

图5 井底面各分析位置示意图Fig.5 Schematic diagram of each analysis position under bottom-hole

(1)有效径向应力σer。从分析位置1到位置4,σer随轴向距离增加而波动减小,并由图6(a)、(b)的全部压应力变为图6(c)、(d)部分拉应力;分析位置1因距井底较近且受井底形状的影响较大,σer(压应力)会出现波动(图6(a));其余分析位置处,σer随径向距离的增加,其值逐步减小。从时间效应角度看:井眼钻开后10-4s到24 h,σer都会先减小后增加,且在井眼打开后100 s时达到最小值。

(2)有效周向应力σeθ。从分析位置1到位置4,随距离井底轴向距离的增加,σeθ逐渐减小(图7(a)～(d)),并且由全部压应力变为位置4处部分拉应力,这与σer的变化规律类似。位置1处随径向距离的增加,各时刻的σeθ逐渐增大;位置2、3和4处,σeθ对于径向位置的变化不敏感。从时间上来看:σeθ的最大值出现在井眼打开的短时间内(10-4、10-2s),为压应力;最小值出现在约100 s,且在位置4处出现了拉应力;同一分析位置处,随着时间延长,σeθ会先减小而后增大。

(3)有效轴向应力σez。从空间上看:从分析位置1到4处(图8(a)～(d)),σez距井底面轴向距离的增加使其值越来越大;同一分析位置处,随距井眼中心线径向距离的增加,σez增大;分析位置1、2、3处,在靠近井眼中心线附近会出现σez为负的拉应力状态区域。各个分析位置处σez的最小值都出现在井眼钻开后约100 s,且此时近井眼轴线出现拉应力状态;σez在井眼钻开瞬间的极大值、24 h后的稳定值都会随分析位置加深而增大。

综上所述,在不同分析位置处会同时存在σer、σeθ为正而σez为负(分析位置1、2、3)或σer、σeθ为负而σez为正(分析位置4)的情况,说明井底面以下岩石应力状态复杂,同时此现象又与井眼开挖后应力演化过程相关。以上计算结果与井底形状密切相关。

图6 不同分析位置处有效径向应力分布Fig.6 Effective radial stress distribution on different analysis positions

图7 不同分析位置处有效周向应力分布Fig.7 Effective hoop stress distribution on different analysis positions

图8 不同分析位置处有效轴向应力分布Fig.8 Effective axial stress distribution on different analysis positions

3.4 井底岩石的强度破坏

3.4.1 破坏准则

本文中采用三维莫尔库伦准则[28-29],其预测中间主应力与岩石强度之间有一种线性关系,并能够较好地对前人发表的多轴实验数据进行拟合。根据该准则所设置损伤函数[30]为

(7)

其中

式中,S0为库伦黏聚力;φ为库伦摩擦角;fdmge函数为零代表岩石开始失效;fdmge函数小于零代表岩石已经严重损坏,其值越小代表岩石失效的可能性也就越大;fdmge函数为正代表岩石较为稳定,其值越大代表失效的风险越小。

3.4.2 井底岩石强度破坏

算例模拟过平衡钻井工况,计算fdmge函数得到的结果如图9所示。在井眼打开瞬间(t=10-4s),离井底面一定距离的区域而非井底面出现强度破坏的情况,但是由于时间极短,只说明在此区域的岩石强度破坏的可能性较大,这与Detournay等[1]对井壁稳定性分析的研究结果类似。井眼打开1 s以后,发现井底周围岩石的fdmge函数值均远大于0,说明其非常稳定并没有发生强度破坏的趋势。这也说明采用单一应力或者部分应力状态去预测岩石强度破坏存在不确定性,例如在井底不同位置处计算得到三向应力在100 s时达到最小值,且部分区域出现拉应力,但计算fdmge函数得到的结果显示岩石处于稳定状态。

逐步降低ph进而模拟近平衡、欠平衡钻井工况,并对井底岩石的fdmge函数进行计算,结果显示(图10):近平衡钻井或当ph小于原始孔隙压力1.75 MPa时的欠平衡钻井工况,井眼打开瞬间以及10 ～2 500 s时间段井底面靠近井壁处出现了岩石强度破坏的情况,其余时间段井底面岩石较为稳定。这一结果与传统理论认为井底与井壁交界面处岩石强度较大而难于破碎有所不同,有研究[28,31]认为井眼外围岩石可钻性降低多是由于井眼形状所导致的岩屑清除需要更多的机械能量;同时钻头切削井眼外围岩层的工作量较大、钻头外围旋转角速度较大也进一步对钻头质量提出了更高的要求。

当压差大于1.75 MPa时的欠平衡钻井工况(图11),井底岩石强度在全时段内均表现出不稳定甚至强度破坏的趋势;随压差继续增加,强度破坏区域随之增大,且破坏区域从井底逐渐上升至井壁处;近平衡与欠平衡条件下的岩石强度破坏区域更加接近于井底面,而非在井眼以下的岩石内部。

图9 过平衡钻井工况下井底fdmge函数值Fig.9 Values of the fdmge calculated with overbalanced drilling condition

图10 欠平衡(压差1.75 MPa)钻井工况下井底fdmge函数值Fig.10 Values of fdmge calculated with underbalanced drilling condition

图11 欠平衡(压差10 MPa)钻井工况下井底fdmge函数值Fig.11 Values of the fdmge calculated with underbalanced drilling condition

4 结 论

(1)通过求解所建立的流固耦合模型,可以实时评估井底岩石的动态响应情况。多孔弹性介质的时间效应可通过井眼轴线上孔隙压力随时间的分布得到体现,孔隙压力变化较大的区域限于近井底处;井底附近岩石应力状态随时间动态变化过程明显,通常钻开井眼的瞬间三向应力达到最大值(压应力),而在数十秒后三向应力会达到最小值(部分拉应力),数小时后井底岩石应力状况趋于稳定。

(2)使用单一应力或者部分应力状态去描述一点的岩石强度破坏情况存在较大的不确定性;过平衡钻井工况下,只在井眼打开瞬间近井底面以下区域存在强度破坏风险而非井底面;欠平衡钻井工况当压差大于1.75 MPa时,井底岩石在全时段出现强度破坏的趋势,随压差增大破坏区域也将扩大。

(3)文中的计算结果与岩石物性、所受载荷、流体性质等密切相关;强度破坏准则的正确性需要建立在区域岩石物性测量和大量强度破坏试验校正的基础之上。