基于实时温度场的井筒稳定性优化模型研究*

2022-03-22方含之孙四维张菲菲

石油机械 2022年3期

方含之王波孙四维彭涛张菲菲

(1.长江大学石油工程学院 2.青海油田钻采工艺研究院 3.中石油渤海钻探工程有限公司工程技术研究院)

0 引言

钻井液循环过程中，对地层加热可以减少井壁破裂的风险，但地层温度过高对孔隙压力的影响较大。因此如何定量地描述井下复杂温度分布情况和考虑热场影响的稳定性亟需各研究者加以重视。CHEN G.Z.等[1]开展了对比可渗透井壁和不可渗透井壁热效应的一致性分析，得到井筒温度升高井壁坍塌风险增大的结论，但并未给出定量关系。H.S.FARAHANI等[2]建立了适用于高渗透性地层中传热的3D热孔-弹性模型，同时提供了一种基于线弹性假设来耦合温度影响的地应力关系式的求解思路。曹文科等[3]结合地层强度准则与拉伸破坏准则，建立了综合地层热交换的深水井壁坍塌压力与破裂压力计算模型，分析认为钻井液循环能够抑制井壁坍塌，但过冷条件下井壁更易引发漏失。

在超深钻井的不同地层深度中，钻具钻进或划眼过程产生的机械热源和钻具与井壁接触产生的热源，会导致不同的井下温度分布情况，从而影响井筒稳定性和钻井液当量循环密度(ECD)评估计算[4]。以往的钻井风险分析多以井筒内流动为主要热源进行，对以上问题考虑较少。

本文研究了地层温度随深度和径向变化的规律，利用线性叠加原理将孔隙压力引起的应力诱导分量和地层温度变化引起的热诱导分量，组合到原位诱导孔隙弹性模型中，建立了充分考虑渗流、孔隙压力变化与地层温度变化影响的ECD计算模型。并选取了一口实例井加以验证，证明了该模型对超深钻井作业有一定指导意义。

1 热孔隙弹性模型

在定量求解适合深浅层井的各诱导应力耦合的热孔隙弹性模型时，首先应分别计算出各项应力值对于岩石变形的影响。其中各时间节点受压力和温度值，由不同循环时间步下地层温度和孔隙压力值耦合的瞬态地层传热模型给出。

1.1 液压诱导应力分量

假设沿井身轴线各层岩性满足径向对称性和均质性，则所有剪切应力和应变均为0，其中轴向没有应变，径向应力分布属于平面应变问题。因此，满足各向同性的岩石变形可用经典胡克定律描述。液压诱导应力引起的岩石变形为：

(1)

式中：σrr、σθθ、σzz分别为径向、切向、轴向诱导应力分量，MPa；αp为Biot系数；ν为泊松比，无量纲；pw为井壁处孔隙压力，MPa；rw为井径，m；r∞为远离井壁位置,m；r为径向步长，m；Δp为原位孔隙压力随时间和径向的变化量，定义地层受压为正，MPa。

式(1)充分考虑了钻井液液压影响，并且由于该式完全基于线性弹性变形，最终的热孔弹性模型可直接通过叠加原理进行组合。

1.2 热诱导应力分量

同样的，地层温度变化而引起的应力可由经典胡克定律确定，对于无外力下，取满足平面应变假设的平衡方程与位移方程式，变形有径向位移：

(2)

式中：αT为由岩石受热的线性热膨胀系数，1/K；ΔTf(r,t)为地层原始温度随时间和径向的变化量，定义地层受压为正，K；D1和D2为无量纲常数，由井壁处与远场条件下对应边界条件确定。

井壁处，径向应力主要由井壁处的温度变化(ΔTw(t)，液压诱导应力单独考虑)引起，此时有：

σrr=-αTEΔTw(t)

(3)

式中：E为岩石弹性模量，GPa。

远场条件下，径向应力为0(原始孔隙压力单独考虑)，即σrr=0。有边界条件：

(4)

代入上述边界条件，则在满足平面应变条件下有：

(5)

(6)

(7)

可以看出，当地层温度升高时，井壁附近地层环向应力增大，对地层加热会产生较小的拉应力，从而降低井壁破裂风险；同样的，对井壁循环降温或减小受压应力，可降低井壁坍塌风险。

1.3 瞬态地层传热模型

在求解地层温度及孔隙压力耦合的温度场时，可以综合井筒传热模型和孔隙压力剖面计算模型，进而计算出地层温度随深度和径向的变化规律。

1.3.1 深井井筒传热

对于给定的井眼轨迹与钻具组合，地层性质和钻井液性质，井筒传感计算可以使用A.R.HASAN等[5-6]提出的分析模型。基于迎风差分格式与能量守恒，推导出瞬态井筒传热模型。各段控制体传热过程示意图如图1所示。

图1 控制体传热过程示意图

基本传热模型建立过程中有如下假设：

(1)钻柱与环空内钻井液均属不可压缩流体；

(2)钻井液的径向温度梯度可忽略；

(3)地温梯度恒定；

(4)无辐射传热；

(5)地层性质与温度无关；

(6)地层孔隙流体无热交换。

图1中：rpi、rpo和rw分别表示控制体内径、外径和环空直径，m；qp和qa分别表示流入/出钻具内和流入/出环空热量，J；qpa表示钻具内热源经钻具壁去环空热量，J；qaf表示环空热源经井壁去地层热量，J。

整理有各段热平衡表达式，具体包含环空中耗散热源、进入环空热源、环空向地层的导热以及钻柱内热源向环空的传热[7]，于是有：

(8)

其中：

(9)

式中：Uap为钻柱内与环空传热机制下等效热阻；ρm为流入钻井液密度，kg/m3；m为流入单元控制体流量，kg/s；ha为环空换热系数，W/(m2·K)；Aa为环空过流面积，m2；cp,m为钻井液比热容，J/(kg·K)；Ta、Tp、Tf分别为环空、钻杆内和地层各点温度，K；Δz为控制体轴向步长，m。

热平衡表达式为：

qa(z+Δz,t)+qap+qs=qa(z,t)+qaf+qEa

(10)

式中：qa、qs、qap和qaf分别为环空段、钻柱内产生的热流量，钻柱到环空、环空到地层传递的热流量，W/m2；qEa为内能变化率。

同样的，钻杆内热平衡表达式为：

(11)

其中：

(12)

(13)

式中：Ap,i为钻柱内过流面积，m2。

式(11)～式(13)对应初始条件为初始温度剖面与地层温度相同。

井口处边界条件为钻井液温度已知，即：

Tp(0,t)=Tin

(14)

井底处边界条件为钻杆与环空流体换热为0，即有：

Tp(L,t)=Ta(L,t)

(15)

式中：Tin为井筒入口温度，K；Ta(L,t)为环空井底任一时刻温度，K。

1.3.2 深井地层传热

地层温度曲线可由扩散方程得出，包括地层孔隙流体的对流传热：

(16)

(17)

式中：αf为地层扩散系数；ρf为地层岩性密度，kg/m3；Cp,f为地层岩性比热容，J/(kg·K)；ρf1为孔隙流体密度，kg/m3；Cp,f1为孔隙流体比热容，J/(kg·K)；vr为孔隙流体流速，m/s；p为孔隙压力，MPa；φ为孔隙度；κf为地层岩石渗透率，mD；μf1为孔隙流体黏度，cP。

为了简化计算，利用牛顿流体假设孔隙内流体流动性质，且地层中仅沿径向流动，最终有地层温度剖面表达式[8]：

(18)

对应有初始条件：地层温度剖面即为初始地热温度。于是有：

Tf(r,z,0)=Tsf+GtDz

(19)

式中：Tsf为地表温度，K；Gt为地温梯度，K/m；Dz为垂深，m。

此时，边界条件如下。

(1)远场条件下，地层温度分布与地热温度相同：

Tf(r→∞,z,t)=Ts+Gtz

(20)

(2)井口与井底地层温度恒定不变；

(3)各边界节点，如钻杆内、外壁及井壁等热流密度可由傅里叶定理直接写出热平衡。

1.3.3 钻井附加热源

由于流体黏度效应，钻柱与地层/套管之间的接触会与运动的钻井液产生接触阻力摩擦。但是，与其他情况下的机械摩擦相比，黏性摩擦的重要性通常较小。因此，本文忽略了由黏性阻力产生的热量。对应的，有钻头位置以及喷嘴节流效应产生机械热源：

(21)

式中：J为焦耳系数；β为钻头效率，%；FWOB为钻压，kN；vROP为机械钻速，m/s；N为转速，r/min；Mbit为钻头扭矩，kN·m。

钻进或划眼过程中，由于钻杆与井壁(或套管)接触产生热源，考虑到目前国内外学者并未深入研究钻杆接触对于井筒温度的影响，本文基于R.MITCHELL等[9]提出的3D摩阻扭矩模型，来估算其接触阻力产生的热源qdrag。

钻进工况：

(22)

划眼工况：

(23)

式中：μf为摩擦因数，无量纲；rp为钻杆外径，m；Tc为单位接触力，kN/m；s1、s2为钻杆移动间距，m。

1.3.4 孔隙压力剖面

低渗透地层中，容易在井壁上形成泥饼，可能不会发生水力扩散[10]。因此，本文所建立的孔隙压力剖面模型忽略了泥饼的存在，考虑了可能发生的从井筒到地层的水力扩散。

一般情况下，假设井筒周围地层为饱和单相孔隙流体，孔隙流体只存在径向流动，并认为孔隙度φ不变，对于柱坐标系，有地层孔隙中流体质量平衡，即：

(24)

Ct=Cf1+Cfg

(25)

(26)

(27)

式中：Ct为总的等温压缩系数，无量纲；Cf1为孔隙流体的等温压缩系数；Cfg为地层岩石的等温压缩系数。

由钻井液在环空的水力扩散和孔隙流体变形耦合引起的瞬时孔隙压力为：

(28)

由于式(28)最后二阶偏导项值较小，大多数计算中或可忽略不记，或可表征流体微团变形导致的压力影响。

假定地层为半无限长圆柱，式(28)对应的初始条件与边界条件为：

(1)近井壁处，孔隙压力等于井眼内的钻井液静水压；

(2)远场条件下，孔隙压力等于初始孔隙压力值。

结合式(5)、式(6)、式(7)、式(10)、式(17)、式(28)及附加热源，即可得到随着深度和径向位移变化的瞬态地层传热模型。该模型充分耦合了钻井液在井筒中流动的状态和地层孔隙流体水力扩散的影响。

2 井壁稳定性模型

2.1 围岩应力分布

通过叠加原理分别将基于线性弹性方法描述的孔隙压力引起的应力诱导分量和热诱导分量组合到原位诱导孔隙弹性模型中，建立起热应力与孔隙压力影响的围岩应力场的热孔隙弹性模型，并可推导出基于热应力影响的当量循环密度模型。再利用破坏准则来确定变形极限并评估钻井过程中的井筒稳定性。

(29)

(30)

(31)

各节点温度值与压力值由上述瞬态地层温度分布模型计算，采用迎风差分格式来离散，并可以通过改变径向步长来减少远场条件下计算节点数量，以加快模型运算速率。

2.2 井壁失效准则

在校核饱和流体的多孔岩石层稳定性时，由于砂岩地层易发生脆性劈裂，考虑到库仑-摩尔(M-C)剪切破坏准则，其定义失效基准面为六边形，忽略了中间主应力的影响，并未考虑井壁渗透作用，因此给出的预测可能过度保守[11-14]。最终基于Drucker-Prager(D-P)模型对Mises模型的六角形基准面光顺修改，能克服Mises准则不适用拉压强度差异较大的岩土材料。

通过对受外力岩石控制体的三个主要应力分量校核[15]，又由于大多数强度标准用主应力表示，则有最终式：

(32)

(33)

(34)

(35)

(36)

不同失效准则下的钻井液安全密度窗口预测值对比如图2所示。

图2 不同准则的ECD对比示意图

从图2可以看出：基于M-C失效准则的ECD预测较为保守，若对超深钻井进行预测，该准则可能完全失效；Lade准则做数值积分耗时过长且不易收敛，不建议采用。修正的D-P准则分别通过对受外力岩石控制体的各向主要应力分量校核，并充分考虑控制体内摩擦角和岩石黏聚力影响，结合本文热孔隙弹性模型，也能很好地考虑热应力与地层孔隙压力以及渗流的耦合影响，所以基于该准则的ECD预测值更为准确，且经济效益更好。

3 实例分析

3.1 实例模型建立

基于井筒传热模型和热孔隙弹性井眼稳定性模型，对一口实际井钻井过程中钻柱、环空和地层各节点温度值进行模拟计算，并对井筒稳定性进行了分析。具体井身结构、层位信息和计算域分网方法如图3所示。钻具组合：ø215.9 mm牙轮钻头×0.37 m+ø172.0 mm螺杆×8.40 m+ø209.0 mm扶正器×1.51 m+浮阀×0.50 m+ø165.0 mm无磁钻铤×8.86 m+近钻头测量短节×3.07 m+ø127.0 mm加重钻杆×15根+震击器×3.55 m+ø127.0 mm加重钻杆×12根+ø127.0 mm钻杆(若干)。首先输入井眼轨迹与钻具组合等参数，并结合其他钻井操作参数，如井口处传热系数、钻井液性能参数以及各段地层岩性参数，模型自主建立由井口至井底的自适应域网格，并判别各个计算域最适宜的方法，迭代给出有效应力分布，并建立基于D-P失效准则的井眼稳定性模型。测试参数为：循环流量12.6 L/s，转速10.5 rad/s，钻头扭矩1 360 N·m，钻压22 200 N，地表温度289.0 K。

图3 实例井井身结构、层位信息及计算域分网

该盆地源岩成熟度和古地温梯度拟合计算表明：晚古生代(石炭系和二叠纪)和中生代早期古地温梯度较高，地温梯度取0.032 ℃/m。白垩纪后古地温梯度逐步降低，现以成为典型“冷盆”，地温梯度取0.020 ℃/m。

3.2 温度场与稳定性分析

井筒内不同循环时间下各井深位置热应力与地层温度分布情况如图4～图6所示。

图4 循环5 h对地层温度的影响

由图4可以看出，钻井液在钻杆内逐渐加热，进入环空后，由于地层传递热源的叠加及喷嘴的节流效应，在井底往上某一处达到最大温度值，而非井底位置。由图5可知，随着循环时间延长，近井底位置冷却效果更好。由图6可知，循环效果与循环时间并非呈线性关系，且近井壁位置温度越低，热诱导应力值越大，因此需考虑适当的循环时间对于地层热应力的影响，而非越久越好。综合分析认为，若停钻期间需要井眼清洁操作，则该口实例井建议循环时间为5 h左右。

图5 井底处切向热应力分布

图6 近井壁位置各值与循环时间的关系

一般认为，井壁失稳主要受钻井操作参数(如井径变化、ROP、泵送流量及岩石各项异性参数)的影响[16-17]。其中地层受力分布不均匀以及渗流作用对井壁压力影响较大[18]。由于本文旨在考虑循环时间对井壁失稳的影响，所以分析孔隙流体流动状态变化(Biot系数)及钻井液流动状态(流性指数)对井壁压力的影响，结果如图7所示。

图7 钻井液不同流性指数对安全密度窗口影响曲线

如图7a所示，Biot系数作为孔隙压力对于失效准则中有效应力的直接评判依据，其值越低，说明岩石可压缩性越高，钻井液安全密度窗口范围增大，并且由于钻井液对地层的冷却作用，地层坍塌风险较破裂风险趋势更大，也表明与地层中其他流体渗流作用(如钻井液渗流作用，见图7b)相比，孔隙流体流动状态对于地层稳定性影响最大。因此钻井循环过程中，钻遇不同地层时，务必准确获取当前进尺下岩石特性参数，这对于安全密度窗口范围评价有重要参考价值。