1. 中国石油大学(华东)石油工程学院,山东青岛 266580；2. 山东省深地钻井过程控制工程技术研究中心, 山东青岛 266580；3. 中国石化石油工程技术研究院,北京 100101；
1. School of Petroleum Engineering in China University of Petroleum (East China), Qingdao 266580 , China；2. Shandong Ultra-deep Drilling Process Control Technology Research and Development Center, Qingdao 266580 , China；3. SINOPEC Research Institute of Petroleum Engineering, Beijing 100101 , China；
中图分类号:TE242
文献标识码:A
基于 Biot 多孔弹性力学理论,分析复杂应力条件下井底岩石的受力情况,建立考虑时间效应的渗流场与应力场全耦合模型,通过数值求解得到井底应力、孔隙压力随时间的动态变化规律,计算结果与 Warren 模型进行对比验证;结合三维莫尔库伦准则建立井底岩石强度破坏的预测模型,使用自定义损伤函数分析井底应力状况对岩石强度破坏的影响规律。 通过计算损伤函数分析不同工况、不同时刻井底岩石的强度破坏情况。 研究表明:求解所建立的流固耦合模型可以实时评估井底岩石的动态响应情况,通常在钻开井眼的瞬间三向应力达到最大值(压应力),而在数十秒后三向应力会达到最小值(部分拉应力),数小时后井底岩石应力状况趋于稳定;过平衡钻井工况下,只在井眼打开瞬间近井底面而非井底面处存在强度破坏风险;欠平衡钻井,当压差大于 1. 75 MPa 时,井底岩石在全时段出现强度破坏的趋势,且随压差增大强度破坏区域扩大。
Based on Biot̍s theory of poroelastic mechanics and the consideration of complex stress conditions of bottom hole rocks, a fully coupled model between seepage field and stress field was established, the dynamic change law of pore pressure and stress distribution at the well bottom was obtained by numerical calculation, and the results were verified by Warren̍s model. The predictive model of the rock strength failure at the bottom-hole combining with the three-dimensional Mohr-Coulomb criteria was developed. The influence of stress distribution on rock damage was discussed. By calculating the failure function, the strength condition of the rock at the bottom-hole under different working conditions and different times was analyzed. The results show that the dynamic response of the rock at the bottom-hole can be evaluated in real-time by using the established fluid-solid coupling model. Generally, the three-dimensional stress reaches the maximum value ( compressive stress)instantly as the wellbore excavated, while the stress reaches the minimum value (partially tensile stress) after tens of seconds, and the rock stress at the bottom-hole will stabilize after a few hours. Under an over-balanced drilling condition, there is a risk of strength failure at a small distance inside the rock, rather than at the bottom-hole wall only when the borehole is excavated instantaneously. When the pressure difference is greater than 1. 75 MPa under the condition of underbalanced drilling, the rock at the bottom-hole has a tendency of strength destruction during the whole period, and the strength failure area will expand with the pressure difference increasing.
钻进引起的井底岩石应力的重新分布可能会造成地层无法承受由此产生的偏应力[1-2],从而造成可能的强度破坏,因此对井底岩石应力分布的研究将有助于对井下复杂情况、岩石破坏机制等进行充分了解。 近年来相关问题的研究多集中于井壁应力分布以及井壁稳定性的影响因素分析[1,3-12],研究多采用弹塑性力学和孔隙弹性力学相关理论分析井壁失稳的力学机制,而弹塑性力学中采用位移法、应力法求解岩石的位移场与应力场[13],忽略了孔隙流体的运移,认为系统处于静态平衡,应力、应变和位移是空间坐标的函数,与时间无关,而得出的应力场分布可看做是井眼开挖一段时间以后的稳定解。 关于井底应力场分布方面的研究相对较少,Warren等[14]研究了井底岩石平均应力变化对机械钻速的影响;Ito等[15]研究了井底应力集中现象对井壁坍塌的影响; 文献[16]~[19]中考虑流固耦合作用将井底按照应力大小分为三向拉应力区、三向压应力区和拉压应力区。 以上研究忽略了井底应力的动态演化过程,时间效应是影响应力分布的重要因素,评价井底岩石强度破坏情况的指标存在不确定性。 笔者基于Biot多孔弹性力学理论,建立考虑时间效应的渗流场与应力场全耦合模型,结合三维摩尔库伦准则建立井底岩石强度破坏的预测模型,分析井底应力状况对岩石强度破坏的影响。
将井底岩石视为理想的多孔弹性介质,则岩体骨架的有效应力[20]为
其中
式中, σ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。
钻井过程中导致的井底岩石骨架变形以及开挖效应导致井眼周围岩石上覆岩层压力急剧减小都会影响孔隙压力的重新分布,同时地层中孔隙压力变化又会反过来影响岩石骨架应力重新分布[20-21]。
假设:①单相不可压缩流体;②流体渗流符合线性渗流规律;③地层岩石均质微可压缩;④弹性不稳定渗流;⑤等温渗流过程;⑥考虑重力等体积力对渗流的影响。 得到的渗流方程为
式中, km 为岩石骨架的渗透率,随骨架岩石变形而变化;μ 为流体黏度; ρ1 为流体密度; α-φKs 为岩石骨架压缩项; φK1 为流体压缩项; Ks 为岩石骨架的体积模量; K1 为流体的体积模量;φ 为孔隙度; α∂ϑ∂t 为体积应变项; qm 为源汇项。
在已知边界条件和初始条件后,就可得到井底渗流场与应力场相互影响的耦合场。
深井地层在上覆岩层压力下所产生的初始位移场是历史时期岩体沉降固结形成的,对钻进过程并没有影响。 将耦合方程的计算过程分为两步,首先求取上覆岩层压力所产生的应力场,将此应力场以第二类皮奥拉-基尔霍夫应力的形式添加到应力张量中,并以此作为第二步的初始条件计算钻进过程对现有位移场的影响[22]。
应力张量分量 σij 和渗透张量分量kij 之间的相互作用是岩石应力场和渗流场耦合的关键。 岩石的应力状态 σij 受渗流而产生的动水压力作用,从而有 σij=f 1(kij);另一方面应力张量分量 σij 变化产生的岩体变形会使渗透率改变,从而又有kij =f 2(σij), Louis [23]通过钻进实验建立了岩石平均渗透系数k与有效正应力之间的经验关系:
式中, k0 为初始渗透系数; σe 为有效应力; α'为经验系数,取决于岩石的孔隙状态。
假设应力主方向与渗流主方向一致[24],计算过程中渗透系数矩阵设为
式中,λ 为影响系数,由实验确定;σei 为有效正应力,i=1,2,3。
计算过程中考虑重力效应,所有边界载荷以及流体压力均为纵向坐标的函数。 文中分析对象仅限于井底周围小区域内,考虑孔隙中液体惯性效应所引起的压力波动。
假设:井底岩石为各向同性多孔线弹性材料, 受到均匀的水平地应力;孔隙骨架充满流体,且孔隙压力随岩石骨架变形而变化,孔隙压力变化又会对岩石骨架应力产生影响;不考虑井斜、温度、 恒定径向流动引起的压力变化;井壁具有良好的渗透性,即井筒液柱压力与地层孔隙压力能够相互传递。
将深部地层岩石的受力情况进行简化,如图1(a)所示,井周及井底岩石受到上覆岩层压力 σV、 最大及最小水平地应力 σH 和 σh、地层孔隙压力pp 和液柱压力ph 作用。 为方便分析做井筒剖面,如图1(b)所示,在模型外围设置无限元域模拟无限大地层,减小模型尺寸对计算结果的影响。
图1 深部地层岩石物理模型
Fig.1 Deep stratigraphic petro-physical model
位移边界条件:在模型底面添加法向约束。
渗流边界包括水头边界、无流动边界和透水层边界。 水头边界为井壁施加液柱压力ph;无流动边界为模型四周与底部边界,在模型四周设置远场地层孔隙压力pp;井筒井壁为透水边界层。
应力边界条件:模型顶部设置上覆岩层压力 σV;模型四周设置均匀水平地应力 σH 和 σh,并认为水平地应力由上覆岩层压力引起,即 σH=σh=ν/(1-ν)(σV-pp) +pp;考虑渗透作用所产生的渗透压力,并将其转换为渗透体力施加到对应的节点和单元之上。
为了使模型计算结果更加接近实际,模型中的计算参数使用塔里木油田某探井的岩心数据(取芯深度5300 m),岩心物性参数:弹性模量E为26308 MPa;泊松比 ν 为0.291;渗透系数k为1.1×10 -8 m/s;孔隙度 φ 为0.1;岩样密度 ρm 为2680 kg/m 3;岩石体积模量Ks 为67000 MPa;剪切模量G为9843 MPa;凝聚力S0 为25.87 MPa;摩擦角 ϕ 为25.07 °; 孔隙流体密度 ρl 为1000 kg/m 3;孔隙流体比重 γ 为1×10 4 N/m 3;流体体积模量Kl 为3270 MPa。
假设物理模型位于井深5500 m处,井眼直径为311.2 mm,钻井液密度为1100 kg/m 3。 得模型计算参数:上覆岩层压力为144.452 MPa;井底液柱压力为59.29 MPa;地层孔隙压力为53.9 MPa;均匀水平地应力(σH=σh)为91.07 MPa。
井底应力场分析的主要目的是了解井底岩石在钻进之后的强度分布情况,进而对井底可能发生的复杂情况有更清楚的认识。
使用Warren等[14] 建立的井底应力场模型对文中所建立数值模型的准确性进行验证。 Warren等通过类比温度应变与孔隙压力应变响应的相似性,建立了井底应力场计算模型,研究了井底应力特征对机械钻速的影响,并与现场试验结果进行了对比,展现了较好的相似性。 Warren模型的主要参数:井深h为3048 m;井径 Φ 为222 mm;钻井液密度 ρf 为1050 kg/m 3;上覆岩层压力 σV 为68.95 MPa;井底液柱压力ph 为32.41 MPa;孔隙压力pp 为32.41 MPa;均匀水平地应力(σH=σh) 为48.26 MPa;弹性模量E为13789.5 MPa;泊松比 ν 为0.25; 渗透率k0 为1 μm 2; 孔隙度 φ 为0.15;岩样密度 ρm 为2262 kg/m 3;岩石基体压缩系数CR 为2.697 × 10 -5 MPa -1;钻头轮廓代号为IADC5-3-7bit。
如图2 所示,文中计算得到的井底平均有效应力与Warren计算结果较为一致,说明文中所建立的模型的准确性,因此模型用于分析井底应力响应以及对岩石损伤程度的影响具有可行性。
图2 渗透性砂岩中井底附近的平均有效应力
Fig.2 Mean effective stress around a borehole in a permeable sandstone
在无孔隙、线弹性材料中,由钻进或井筒液柱压力改变所致的应力重新分布是即时的,但在多孔可渗透弹性介质中,孔隙压力的传播只能以一定速度进行,此时应力重新分布是动力学问题,具有时间效应。 多孔弹性介质的时间效应主要有两个来源:井筒钻井液侵入所导致的孔隙压力改变; 钻进过程所致地层应力的重新分布会产生体积应变所引起的孔隙压力的变化。 Detournay [1] 等和Charlez [25]详细讨论了直井、非均匀远场应力下井壁应力分布的时间效应,但是对于井底应力的研究尚未进行。
图3 为孔隙压力pp 在井底中心线上随时间的分布情况(利用孔隙压力pp 与原始地层压力p0 的比值定量衡量pp 的变化),离井底越近压力梯度越大。 在打开井眼的极短时间内(10 -2、10 -4 s),近井底处应变释放且钻井液未充分侵入井底,pp 会突然下降;随轴向距离增加,pp 变化幅度减小且其值增加。 随着时间增加( 10 2 s、1 h),应力重新分布导致地层岩石孔隙压缩且其内流体排出速度不足,pp 值上升,6 h后趋于稳定,且比原始孔隙压力要大。
Fig.3 Changes in pore pressure along the downhole centerline
计算模型采用各向同性结构、水平地应力相等且井眼轴线沿着某一主应力的方向。 使用各强度准则对岩石强度进行衡量时均需要用到有效应力,
由式(5)所给出的有效应力的概念是基于Biot常数,且方程推导时假设岩石为线弹性特性,并不能直接用于岩石损伤判别。 将Biot常数 α 设置为1, 式(5)变为广泛接受的Terzaghi有效应力[2,26-27],且能够用于岩石损伤准则,即
井底面上的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
为了方便研究井底面以下岩石的应力分布情况,定义图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 -4 s到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 -2 s),为压应力;最小值出现在约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
本文中采用三维莫尔库伦准则[28-29],其预测中间主应力与岩石强度之间有一种线性关系,并能够较好地对前人发表的多轴实验数据进行拟合。 根据该准则所设置损伤函数[30]为
其中
式中, S0 为库伦黏聚力;ϕ 为库伦摩擦角; fdmge 函数为零代表岩石开始失效; fdmge 函数小于零代表岩石已经严重损坏,其值越小代表岩石失效的可能性也就越大; fdmge 函数为正代表岩石较为稳定,其值越大代表失效的风险越小。
算例模拟过平衡钻井工况,计算 fdmge 函数得到的结果如图9 所示。 在井眼打开瞬间( t=10 -4 s), 离井底面一定距离的区域而非井底面出现强度破坏的情况,但是由于时间极短,只说明在此区域的岩石强度破坏的可能性较大,这与Detournay等[1]对井壁稳定性分析的研究结果类似。 井眼打开1 s以后, 发现井底周围岩石的 fdmge 函数值均远大于0,说明其非常稳定并没有发生强度破坏的趋势。 这也说明采用单一应力或者部分应力状态去预测岩石强度破坏存在不确定性,例如在井底不同位置处计算得到三向应力在100 s时达到最小值,且部分区域出现拉应力,但计算 fdmge 函数得到的结果显示岩石处于稳定状态。
逐步降低ph 进而模拟近平衡、欠平衡钻井工况,并对井底岩石的 fdmge 函数进行计算,结果显示(图10):近平衡钻井或当ph 小于原始孔隙压力1.75 MPa时的欠平衡钻井工况,井眼打开瞬间以及10~2500 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
(1)通过求解所建立的流固耦合模型,可以实时评估井底岩石的动态响应情况。 多孔弹性介质的时间效应可通过井眼轴线上孔隙压力随时间的分布得到体现,孔隙压力变化较大的区域限于近井底处; 井底附近岩石应力状态随时间动态变化过程明显, 通常钻开井眼的瞬间三向应力达到最大值(压应力),而在数十秒后三向应力会达到最小值(部分拉应力),数小时后井底岩石应力状况趋于稳定。
(2)使用单一应力或者部分应力状态去描述一点的岩石强度破坏情况存在较大的不确定性;过平衡钻井工况下,只在井眼打开瞬间近井底面以下区域存在强度破坏风险而非井底面;欠平衡钻井工况当压差大于1.75 MPa时,井底岩石在全时段出现强度破坏的趋势,随压差增大破坏区域也将扩大。
[20] 陈勉,金衍,张广清.石油工程岩石力学[M].北京:科学出版社,2008:33-39,55-57.
[21] 程林松.高等渗流力学[M].北京:石油工业出版社,2011:168-203.
国家重大科技专项(2016ZX05022-002)；国家自然科学基金面上项目(51674284)；中石油重大科技项目(ZD2019-183-005)；长 江学者和创新团队发展计划项目(IRT_14R58)；
呼怀刚,管志川,许玉强,等. 基于多孔弹性力学理论的深井井底应力场分析[J]. 中国石油大学学报(自然科学版),2020,44(5):52-61.
HU Huaigang, GUAN Zhichuan, XU Yuqiang ,et al. Bottom-hole stress analysis of ultra-deep wells based on theory of poroelastic mechanics[J]. Journal of China University of Petroleum(Edition of Natural Science),2020,44(5):52-61.
[20] 陈勉,金衍,张广清.石油工程岩石力学[M].北京:科学出版社,2008:33-39,55-57.
[21] 程林松.高等渗流力学[M].北京:石油工业出版社,2011:168-203.