深层页岩欠平衡钻井气液固三相瞬态流动传热模型研究

张文平，许争鸣，吕泽昊，赵 雯

（1.中石化石油工程技术研究院有限公司，北京 102206；2.中国地质大学（北京）能源学院，北京 100083；3.中国石油油气和新能源分公司，北京 100007；4.中国石油青海油田分公司采气二厂，青海茫崖 816499）

深层页岩气钻井面临高温高压、缝洞发育导致井漏事故频发、溢流和漏失同存、深部地层机械钻速低和井壁失稳等挑战，为提高深层页岩气钻进安全系数和效率，可采用欠平衡钻井技术实现井底压力精确控制和提速提效[1-8]。欠平衡钻井是一种将井筒压力保持在被钻地层压力以下，允许地层流体流入井筒并从地面循环出来的技术。使用低密度钻井液或注入惰性气体可以产生欠平衡状态。与常规钻井相比，欠平衡钻井可以减少地层损伤、提高钻速、降低漏失量、降低压差卡钻风险和快速显示产层等[9-14]。有效控制井底压力，避免大量地层流体侵入和地面环空压力过高是欠平衡钻井作业成功的关键。在欠平衡钻井时，注入液体、侵入气体和岩屑在井筒中共存，为了更好地表征流动特性和分析井下数据，需要建立钻杆内液相管流流动模型、储层侵入速率模型、环空内气液固流动模型和井筒与地层间的对流换热模型[15-17]。学者们研究了欠平衡钻井过程中的流型转变边界和不同流型下流动参数[18]、井底压力自动控制[19]、井眼清洁效率[20]、储层表征[21-22]等。然而，大多数研究忽略了岩屑对井筒多相流动的影响，直接应用气液两相流的均相模型[23]、漂移流模型[24]或双流体模型[25]模拟井筒内的多相流动特性和优化欠平衡钻井参数。欠平衡钻井技术已应用于陆上[26-27]和海上深层油藏[28-29]的开发，已经证实了流体温度沿井筒深度变化会影响流体性质（密度、黏度等），进而影响多相流动特性。因此，将环空流体与周围环境之间的对流换热嵌入多相流模型是非常重要的。

两相流模型在计算中浅层地层压力和温度时精度可满足工程需求，但随着深度和温度增加，计算精度会大幅度降低。当前模拟欠平衡钻井井筒中气液固三相流动传热特性的模型较少，且大多为等温稳态模型和不考虑岩屑的影响，影响了深部地层井底压力预测的准确性。此外，明确井筒多相流动动态演变特性，可以为现场作业者调整作业参数、实现欠平衡钻井作业性能最大化提供理论指导。因此，建立瞬态-非等温的气液固三相流动模型，可模拟欠平衡钻井作业过程中的压力、温度、相体积分数和流体性质变化规律，基于建立模型明确了关键参数对井筒中气液固三相流体动态流动和传热特性的影响规律，为控压钻井、欠平衡钻井在深层页岩油气中的高效应用提供了理论支撑。

1 气液固三相瞬态流动传热模型的建立

欠平衡钻井过程中，环空中存在钻头破岩产生的岩屑、注入的钻井液和侵入的地层气体，因此，环空内为气-液-固三相流动状态。基于以下假设建立了控制方程：1）沿井筒一维流动；2）处于同一位置的气、液、固三相具有相同的压力和温度；3）忽略不同相之间的传质；4）岩屑对气液两相流型转换准则无影响；5）岩屑在液相中分布均匀，岩屑为球形且大小均匀；6）忽略钻杆运动对多相流特性的影响。

1.1 控制方程

气相、液相和岩屑的质量守恒方程为：

式中：ρ为密度，kg/m3；α为相体积分数；v为速度，m/s；A为环空截面积，m2；下标m=g，l，s，分别表示气相、液相和岩屑，且三者的相体积分数之和为1.0(∑αm=1.0)。

气相、液相和岩屑的动量守恒方程为：

式中：p为压力，Pa；θ为井斜角，（°）；g 为重力加速度，m/s2；FWm为各相壁面剪切应力，Pa/m；FIm为各相间的作用力，Pa/m，并且各相间的作用力之和为0(∑FIm=0)。

对于环空中向上流动的多相流，多相流温度由以下因素决定：1）沿环空向上运移的热量；2）环空流体与钻杆之间的热量交换；3）环空流体与生产套管（套管井段）或地层（裸眼段）之间的热量交换。因此，环空中混合物的能量守恒方程为[30]：

式中：u为内能，J/kg；Q为环空流体与周围环境（包括钻杆和地层）之间的热交换，W/m；Aan为环空截面积，m2；为各相单位体积的流入焓，W/m3。

环空流体与周围环境的换热计算公式为：

式中：rdp,o为钻杆外壁半径，m；hdp-an为钻杆与环空流体间的对流换热系数，W/（m2·K）；Tdp为钻杆温度，K；Tan为 环空流体温度，K；rca3,i为生产套管的内半径，m；han-ca3为环空流体与生产套管之间的对流换热系数，W/（m2·K）；Tca3为生产套管温度，K。

由于钻杆温度和生产套管温度也随时间和位置变化，因此需要额外的能量守恒方程来求解环空中混合物的能量守恒方程。钻杆内钻井液温度控制方程为：

式中：ρd为钻杆内钻井液密度，kg/m3；vd为钻杆内钻井液速度，m/s；cp,d为钻杆内钻井液比热容，J/（kg·K）；Td为钻杆内钻井液温度，K；rdp,i为钻杆内半径，m；hd-dp为钻杆内钻井液与钻杆之间的对流换热系数，W/（m2·K）。

钻杆的温度控制方程为：

式中：ρdp为钻杆密度，kg/m3；cp,dp为钻杆比热容，J/（kg·K）；λdp为钻杆导热系数，W/（m·K）。

地层的温度控制方程为：

式中：λf为地层导热系数，W/（m·K）；Tf为地层温度，K；cp,f为地层比热容，J/（kg·K）。

1.2 相间相互作用力

动量守恒方程中，FIm为各相单位体积承受的相间力，它控制着不同相之间的滑移，并高度依赖于流型。相间作用力包括虚拟质量力和拖曳力。拖曳力是由不同相之间的剪切力引起的；虚拟质量力是由不同相之间的相对加速度引起的，通常忽略虚拟质量力项，以简化计算。

对于泡状流、分散气泡流和段塞流，液相作用于气相上的拖曳力FIG通常表示为[31]：

式中：FIg为气相拖曳力，N；aIg为气相的相间面积分数；Cdg为气相的曳力系数。

假设所有流型下，岩屑在液相中均匀分布，因此，液相对岩屑施加的拖曳力FIs可表示为：

式中：FIs为岩屑拖曳力，N；ds为岩屑直径，m；Cds为岩屑的曳力系数。

由于本文液相为牛顿流体，岩屑均为球形，因此采用以下经验相关法计算岩屑曳力系数[32]：

式中：Res为岩屑雷诺数。

1.3 壁面摩擦力

对于泡状流、分散泡状流和搅拌流，将环空内气、液、固三相的混合物假设为均相流体来计算壁面摩擦力：

式中：fM为三相混合物的摩擦系数；下标 M表示环空气液固三相混合物。

fM是相对粗糙度和雷诺数的函数[33]：

式中：ReM为 环空混合物雷诺数；ε为壁面绝对粗糙度，m；A*为中间参数。

1.4 边界条件

假设钻遇的地层都是含气层，即没有液体从储层侵入井筒，因此环空底部的液体质量流速与钻杆入口的液体质量流速相等。

该模型假设侵入气体为甲烷，根据径向流动方程[34]计算气体流入速率：

式中：Qg为 气体体积流量，m3/s；pr为地层压力，Pa；pb为井底压力，Pa；K为地层渗透率，mD；h为地层厚度，m；tD为无因次时间。

井底岩屑体积流量的计算公式为：

式中：Qs为岩屑体积流量，m3/s；vpe为机械钻速，m/h。

环空井口压力设定为定值。

2 模型求解及验证

2.1 模型求解

流动方程（质量和动量守恒方程）与能量守恒方程解耦求解[35]。首先，利用上一时间步的温度剖面来更新流体性质，然后求解流动方程，得到当前时间步的相体积分数、相速度和压力；然后，利用新的流动参数更新能量守恒方程中的所有系数，求解整个井筒-地层系统的能量守恒方程，得到当前时间步的温度分布；最后，在流动方程和能量守恒方程之间进行迭代，直到满足收敛准则。

式（1）和式（2）可以表示为以下形式：

式中：W、F(W)、Q(W)为密度、速度、相体积分数和压力的函数，如表1 所示。

表1 式（15）中的变量计算公式Table 1 Variables in Eq.(15)

用有限差分法对式（15）进行离散化[36]，得到如下方程：

采用有限差分法，求解整个井筒-地层系统的能量守恒方程。

2.2 模型验证

为了验证所提出的气-液-固三相流动模型的准确性，将计算结果与实验数据[37]进行了比较。采用空气和水作为实验工质，以陶瓷球为固体颗粒。空气的密度由理想气体的状态方程决定。实验过程中，圆管直径为0.026 m，圆管长度为6.74 m，固相密度为2 540 kg/m3，岩屑直径为6.1 mm，液相密度为1 000 kg/m3，液相黏度为1.0 mPa·s。计算结果与实测值的对比如图1 所示。

图1 液相速度实测值与模型预测值间的对比Fig.1 Comparison between measured and predicted liquid velocities

由图1 可知，平均相对误差为3.66%，最大误差为16.80%，计算得到的液相速度与实测值吻合较好。此外，还将预测的井筒压力与墨西哥Agave 301 井的现场数据进行了比较。Agave 301 井的压力计在井深2 259 m 处，地表温度为20.85 ℃，氮气注入速率为10 m3/min，环空回压为0.069 MPa，钻井液密度为949 kg/m3，钻井液黏度为10.0 mPa·s，机械钻速为3.0 m/h，地温梯度为1.745 ℃/100 m，钻井液注入速率为0.45 m3/min，钻杆和环空压力的预测与现场测量的平均相对误差为1.76%和1.29%（见图2）。

图2 Agave 301 井实测压力与预测压力比较Fig.2 Comparison between measured and predicted pressure of Well Agave 301

3 井筒多相流动特性影响因素分析

研究了岩屑及井筒与地层间对流换热对井筒流动传热特性的影响，利用控制变量法，分析了敏感性参数（欠平衡压差、机械钻速、地温梯度）的影响规律。模拟井的井深为7 000 m，钻杆内径为70.2 mm、钻杆外径为88.9 mm，钻头直径为152.4 mm，液相质量流量为9.0 kg/s，地层渗透率为60 mD，地层孔隙度为50%，地层厚度为50 m，地表温度为25 ℃，钻井液注入温度为10 ℃，液体黏度为10.0 mPa·s，节流阀压力为1.0 MPa，岩屑密度为2 650 kg/m3，液相注入压力为11.5 MPa。

3.1 岩屑

以前的欠平衡钻井模型，通常忽略岩屑的影响，而岩屑影响井筒中的压力分布，进而影响井筒多相流动特性。基于本文所建立的模型，比较了模型中考虑和不考虑岩屑存在时环空多相流动特性差异。欠平衡压差和地温梯度分别设定为1.0 MPa 和2.0 ℃/100m。计算中考虑岩屑影响时，岩屑直径和机械钻速分别设为2.5 mm 和15 m/h。考虑和不考虑岩屑存在时，利用所建模型计算出的井底压力曲线、岩屑体积分数/岩屑重力压降、气体侵入速率/气体侵入总质量、气相体积分数剖面如图3 所示。从图3（a）可以看出，不考虑岩屑时，井底压力被低估了约2 MPa，60 000 s 时，有岩屑存在时为57.043 MPa，无岩屑存在时为55.041 MPa，2 种工况下的井底压力差为2.002 MPa，相对误差为3.64%。分析认为其原因是：首先，岩屑的存在增加了环空中混合物的重力压降，如60 000 s时井筒底部岩屑重力压降为1.219 MPa（见图3（b））。其次，由于井底压力升高，环空底部与储层之间的压差减小，从储层侵入井筒气体的速率减小，考虑和不考虑岩屑存在时，60 000 s时储层侵入井筒气体的速率分别为1.018 kg/s 和1.137 kg/s（见图3（c））；从储层侵入井筒的气体越少，井底压力降低程度越小。此外，井筒压力升高会使气体密度升高，从而进一步降低环空中的气相体积分数；考虑和不考虑岩屑存在时，60 000 s 沿井筒的最大气体分数差为0.031 6（见图3（d））。

图3 考虑与不考虑岩屑时井底压力、岩屑体积分数/重力压降、气体流入速率/质量和气体分数剖面Fig.3 BHP,cuttings fraction/hydrostatic pressure,gas influx rate/mass,and gas fraction profile with and without considering cuttings

3.2 井筒-地层传热

井筒与地层间的传热影响井筒温度场，进而影响井筒中流体性质，最终影响井筒多相流动特性。为此，比较了考虑与不考虑井筒-地层间对流换热时的环空多相流动特性。欠平衡压差为1.0 MPa，岩屑直径为2.5 mm，机械钻速为15 m/h，地温梯度为2.0 ℃/100 m。在不考虑环空流体与周围环境之间换热的情况下，利用所建模型计算考虑与不考虑井筒-地层对流换热时的井底压力、温度、气体密度和岩屑体积分数，结果见图4。计算时将环空流体温度设为常数，并采用了2 种环空流体温度分布方式：1）环空流体温度等于钻井液注入温度（条件1）；2）随着井深的增加，环空流体温度从钻井液注入温度开始，随着地温梯度增大而升高。

图4 考虑与不考虑换热效应时的井底压力、环空流体温度、气体密度和岩屑体积分数分布曲线Fig.4 BHP,annulus fluid temperature,gas density,and cuttings fraction distribution curve with and without considering heat transfer effects

从图4（a）可以看出，不考虑井筒-地层对流换热时，井底压力被高估了。在100 000 s 时，考虑对流换热、条件1 和条件2 等3 种条件下的井底压力分别为56.996，61.932 和57.694 MPa。不考虑换热效应时井底压力的相对误差分别为8.6 6% 和1.21%。图4（b）为不同工况下环空流体温度分布。考虑传热效应的情况下，环空流体温度随时间推移而升高。考虑换热效应、条件1 和条件2 等3 种条件下，在100 000 s 时，考虑换热效应的环空流体温度均大于条件1，而其先大于后小于条件2；与考虑换热效应的计算结果相比，条件1 和条件2 的最大温差分别为101.12 和43.67 ℃。环空流体温度不同也会使气体密度等沿井筒流体性质的不同，从而使井筒压力也不同（见图4（c））；在不考虑换热效应的情况下，井筒岩屑体积分数被高估（见图4（d））。

3.3 欠平衡压差

欠平衡压差影响气相进入井筒中的速率，进而影响井筒多相流动特性。岩屑直径为2.5 mm，机械钻速为15 m/h，地温梯度为2.0 ℃/100m，欠平衡压差分别为0.5，1.0 和1.5 MPa，利用所建模型分别计算得到井筒压力随井深变化剖面和气体侵入速率随时间的变化，结果见图5。从图5 可以看出，随着欠平衡压差增大，井筒压力升高。在60 000 s 时，欠平衡压差为0.5，1.0 和1.5 MPa 时的井底压力分别为60.790，57.044 和54.136 MPa，这主要是因为气体侵入速率随着欠平衡压差增大而升高（见图5（b））。

图5 不同欠平衡压差下井筒压力分布及气体侵入速率随时间的变化曲线Fig.5 Variation curve of wellbore pressure distribution and gas influx rate with time under different underbalanced pressure difference

不同欠平衡压差条件下，利用所建模型分别计算得到9 520 和60 000 s时气相体积分数随深度的变化关系，结果见图6。从图6 可以看出，由于气体流入速率较高，沿井筒的气体分数随欠平衡压差增大而增大；60 000 s 时，欠平衡压差为0.5 MPa 与1.0 MPa和1.0 MPa 与1.5 MPa 的最大气相体积分数差分别为0.108 0 和0.076 6。从图6（a）可以看出，由于气体速度较大，随着欠平衡压差增大，气相到达地面的时间会更早。

图6 不同欠平衡压差下的气相体积分数分布曲线Fig.6 Gas fraction distribution curve under different underbalanced pressure differences

3.4 机械钻速

机械钻速影响岩屑进入井筒的速率，因此分析机械钻速对井筒多相流动特性影响规律。在欠平衡压差1.0 MPa、岩屑直径2.5 mm、地温梯度2.0 ℃/100 m，机械钻速分别为10，15 和20 m/h 的条件下，利用所建模型计算井底压力随时间的变化，结果见图7。

图7 不同机械钻速下井底压力随时间的变化曲线Fig.7 Changes curve of BHP with time under different ROPs

从图7 可以看出，井筒压力随着机械钻速增大而升高。60 000 s 时，机械钻速为10，15 和20 m/h 时的井底压力分别为56.320，57.044 和57.811 MPa。这是因为机械钻速越大，环空混合物的静液柱压力也越大。

绘制不同机械钻速下环空流体温度分布曲线，结果如图8 所示。由图8 可以看出，空流体温度随着机械钻速增大而升高。随着机械钻速增大，环空中的岩屑体积分数也随之增大。由于岩屑比热容小于液相和气相比热容，环空流体的比热容随着岩屑体积分数增大而减小。因此，由于混合物的比热容较小，环空流体温度随着机械钻速增大而升高。

图8 60 000 s 时不同机械钻速下的环空流体温度分布曲线Fig.8 Annulus fluid temperature distribution curve under different ROP at 60 000 s

3.5 地温梯度

地温梯度是影响井筒温度的重要因素。在欠平衡压差1.0 MPa，岩屑直径2.5 mm，机械钻速15 m/h，地温梯度分别为1.0，2.0 和3.0 ℃/100m 的条件下，利用所建模型计算60 000 s 时的环空流体温度分布，结果见图9。从图9 可以看出，环空流体温度随着地温梯度增大而升高。这是因为随着地温梯度增大，环空流体与周围环境的温差增大，可以从周围环境吸收更多的热量。此外，环空流体的最高温度并不出现在井底。这种情况下，环空流体最高温度出现在井深5 795 m处，是井底以上井深的六分之一。这种现象的原因可以从环空流体与周围环境的热交换来解释：首先，钻杆中的钻井液吸收环空流体的热量，其温度在井底达到最高；然后，钻井液随岩屑和侵入气体沿环空向上流动，由于环空中混合物的温度低于地层温度，混合物继续从地层中吸收热量；在环空下部，从储层吸收的热量大于钻杆流体损失的热量，因此环空流体温度升高；当环空流体的净换热为零时，环空流体温度达到最高。

图9 60 000 s 时不同地温梯度下的环空流体温度分布曲线Fig.9 Annulus fluid temperature distribution curve under different geothermal gradients at 60 000 s

绘制不同地温梯度下井底压力随时间变化的曲线和气相密度随深度的变化曲线，如图10 所示。从图10 可以发现，地温梯度增大会导致沿井筒气体的密度降低，从而导致井底压力降低。

图10 不同地温梯度下的井底压力和气体密度分布曲线Fig.10 BHP and gas density distribution curve profiles under different geothermal gradients

综合上述研究结果可以看出，考虑岩屑影响时，井筒中的气体侵入速率和气相体积分数都会相应降低；不考虑井筒-地层对流换热时，井筒压力会被高估。随着机械钻速增大，井筒压力、井筒温度和岩屑体积分数都会相应增大；且随着机械钻速增大，欠平衡钻井过程中气相和岩屑到达地面的时间会更早。

4 结论

1）针对深层页岩欠平衡钻井过程，建立了井筒内气-液-固三相瞬态流动传热模型；与欠平衡钻井作业的现场数据进行了对比验证，验证了模型的可靠性和准确性。

2）欠平衡钻井过程中不考虑岩屑影响时，井底压力被低估。这是因为岩屑会产生额外的岩屑重力压降，使井底压力升高，从而导致侵入气体速率减小，并且进一步降低环空中的气相体积分数。

3）欠平衡钻井过程中不考虑井筒-地层对流换热影响时，井底压力被高估，因此需要准确预测井筒温度场及其对井筒流体性质和压力分布的影响。

4）本文所建立的气-液-固三相流动模型拓展了传统的气-液两相流动模型，为页岩欠平衡钻井过程中井筒压力准确预测和高效调控提供了理论基础。