随着海上油气开发逐渐从浅水转向深水、超深水，海底低温高压的环境为天然气水合物（以下简称水合物）生成创造了良好条件^[1]。深水油气井打开目的层后进行避台关井时，地层中的气体会通过扩散和置换的方式侵入井筒并滑脱上升，气体到达井口后聚集形成圈闭高压，在海底附近井段内会生成水合物，严重时甚至堵塞防喷器及井筒，给避台后解封、开井作业带来严重威胁^[2]。研究关井期间井筒内含水合物相变的气泡（以下简称水合物气泡）上升规律对避台关井安全作业周期的预测至关重要。

深水油气井井筒内水合物气泡在上升过程中，水合物膜横向生长占主导地位，多采用水合物生长模型预测水合物膜厚度，但预测结果存在较大误差^[3–6]。目前，关于气泡上升动力学的研究对象主要是表面无水合物的气泡^[7–10]。气泡直径较小时，可以应用Hadamard-Rybczynski模型^[7]预测气泡的上升速度，但该模型仅适用于纯水环境条件；对于直径较大的气泡，其上升运动主要受惯性力控制，R. M. Davies等人^[8]利用势流理论建立了气泡上升速度模型，但未考虑流体的黏度和表面张力，具有一定局限性；对于中等直径的气泡，处于表面张力控制区时，上升过程中显著受到表面振荡和变形的影响，上升速度呈现分散的趋势，并建立了不同速度预测模型^[9–10]。然而，针对水合物气泡的研究相对较少，且不能完全表征全尺度气泡的上升规律^[11–12]。笔者通过实验和理论分析，建立了关井条件下考虑水合物相变的气泡上升速度综合预测模型，提出了安全作业周期的预测方法，为南海深水海域安全钻井提供了理论和技术支持。

1.   含水合物相变的气泡上升规律实验

1.1   实验装置及方法

实验装置由圆柱水槽、温度控制系统、压力控制系统、气体注入系统、数据采集系统和高速摄像系统等组成（见图1）。圆柱水槽内部直径0.30 m，高2.50 m，正面有3个观察窗，右侧有2个观测窗，便于从不同角度观察水合物气泡的上升运动。温度控制系统由冷水机、循环管线和保温层组成，冷水机通过置于水槽内的循环管线进行循环冷却，能将水冷却到3 ℃；同时，装置外层包裹有保温材料，以保持体系温度稳定。压力控制系统由手摇注水泵和注气泵组成，手摇注水泵通过压入和回吸水体来控制系统实验压力，手摇注气泵进行辅助加压。气体注入系统由气瓶和中间容器组成，气体注入之前储存在中间容器中，以便记录注入气体的压力和体积。数据采集系统主要包括温度传感器、压力传感器和数据采集终端。实验过程中所用甲烷气体纯度为99.8%，实验用水为去离子水。

图  1  实验装置示意

Figure  1.  Schematic diagram of the experimental device

下载: 全尺寸图片 幻灯片

利用冷水机组和手摇泵来控制圆柱水槽内的温度和压力，待达到水合物形成条件后，缓慢注气至喷嘴处悬浮一个甲烷气泡，然后打开高速摄像系统，拍摄气泡表面水合物膜的生长过程；当水合物膜完全覆盖气泡表面后，继续注气，使水合物气泡脱离喷嘴，向上运移，直至气泡消失在观察窗口视野之内停止拍摄，完成该直径下的甲烷水合物气泡上升实验。通过控制手摇注气泵的注气速度来控制气泡直径，该实验装置可以生成直径0.7～6.0 mm的甲烷气泡。

1.2   气泡表面水合物生长

悬浮于喷嘴处的气泡在静置一段时间后，水合物晶体开始成核，并沿着气泡表面生长形成一层水合物薄膜，直至覆盖整个气泡。在4 ℃、7 MPa下水合物气泡表面水合物膜横向生长过程如图2所示。

图  2  水合物膜沿气泡表面横向生长

Figure  2.  Lateral growth of hydrate film along the bubble surface

下载: 全尺寸图片 幻灯片

从图2可以看出，气泡表面生成水合物膜后，有水合物膜包裹的气泡表面光泽明显暗于“干净”的气泡表面，并且从开始生成水合物时，气泡底部就产生明显的褶皱现象。这是由于在注气形成气泡过程中水合物膜不仅沿气液界面生长，还随着气泡膨胀而拉伸，从而导致底部褶皱变形。同时，水合物膜前端处于同一水平线上，说明气泡表面同一高度处水合物晶体的生长速率是相同的。

1.3   水合物气泡上升

气泡表面完全覆盖水合物膜后，通过外界扰动释放水合物气泡，观测甲烷水合物气泡上升运动轨迹（见图3（a）），并对比了表面未被水合物包裹气泡的上升运动轨迹（见图3（b））。观测结果表明：对于“干净”的气泡，当气泡直径较小时，气泡呈球形或椭球形；气泡直径增大时，由于流体动力学和表面张力的影响，气泡表面持续振荡，气泡外形开始变成扁平的椭球，然后向不规则的球冠、椭球冠转变，气泡运动轨迹呈螺旋形。对于水合物气泡，水合物膜完全覆盖气泡表面后，由于水合物膜具有抵抗变形的强度，气泡表面无表面波存在，气泡呈稳定的球形或椭球形，水合物气泡运动轨迹呈螺旋形。

图  3  气泡的上升运动轨迹

Figure  3.  The ascending trajectory of the bubble

下载: 全尺寸图片 幻灯片

2.   水合物气泡上升动力学模型

实验表明，气泡进入水合物稳定区域后表面会生成水合物，从而改变气泡的表面形态和动力学行为。水合物的生成会对气泡上升速度产生很大影响，在建立气泡上升速度模型前要对水合物生成速率和水合物厚度进行预测。通常情况下，应用水合物膜横向生长速率模型预测水合物膜的厚度。

2.1   水合物膜厚度预测

水合物膜横向生长过程中，由于气液界面没有完全被水合物覆盖，仍有部分裸露在水中，可以认为气液边界层周围有充足的气体分子供水合物生长，生长速率主要受传热过程控制。为简化计算过程，建立模型前进行如下假设：1）水合物膜在连续生长的过程中初始厚度是一个定值，即不发生增厚生长；2）水合物膜前沿向周围水相和气相中的传热为瞬变稳态过程；3）水合物膜是通过消耗前沿界面处的水和气体形成的；4）水合物膜前沿表面温度为实验压力下的三相平衡温度，并在传热过程中保持不变。水合物膜横向生长传热过程如图4所示。

图  4  水合物膜横向生长传热过程

Figure  4.  The lateral growth heat transfer process of hydrate film

下载: 全尺寸图片 幻灯片

水合物生成是放热过程，水合物膜沿气液界面生长时生成热会向周围扩散，分别包括与膜前沿接触的水相、气相及水合物膜内部传热。其中，水合物膜前沿向膜内部的传热可以忽略^[6]，因此水合物膜沿气液界面的生长速率主要受从膜前沿向水相、气相传热过程控制。根据水合物生成热与膜前沿向周围传递的热量相平衡原理，可以得到：

式中：${v _{\rm{f}}}$为水合物膜横向生长速率，m/s；${\rho _{\rm{h}}}$为水合物的密度，kg/m³；$\Delta {H_{\rm{h}}}$为水合物生成热，J/kg；${h_{\rm{w}}}$，${h_{\rm{g}}}$分别为水和气体的传热系数，W/(m²·K)；${l_{{\rm{h}} - {\rm{w}}}}$，${l_{{\rm{h}} - {\rm{g}}}}$分别为单位长度下水合物膜前沿与水、气体的接触面积（如图4黑色曲线、直线所示），m²；$\Delta T$为系统过冷度，K。

为简化计算，假设水合物膜前沿曲率半径等于水合物膜厚度，即r=$\delta$。由于水合物膜前沿处温度为实验压力下的三相平衡温度${T_{{\rm{eq}}}}$，与实验温度${T_{{\rm{ex}}}}$形成温度差，导致水合物膜前沿周围流体密度分布不均，热边界层内存在温度梯度而引起流体密度变化，从而产生局部自然对流传热。自然对流过程中水和气体的传热系数${h_{\rm{w}}}$、${h_{\rm{g}}}$可以根据特征数实验关联式求得：

式中：Nu_w/g为努塞尔特征数；Gr_w/g为格拉晓夫特征数；Pr_w/g为普朗特特征数；m和n为特征系数。

水合物膜前沿传热界面包括与水接触的水合物膜前沿界面和与气相接触的水合物膜前沿界面（见图5）。基于水合物膜前沿的几何特征，与水接触的水合物膜前沿界面的传热问题可以转化为水平圆柱在大空间水相中的自然对流传热问题（见图5（a））；与气体接触的水合物膜前沿界面的传热过程类似于水平矩形平板在大空间气相中的自然对流传热过程（见图5（b））。

图  5  水合物膜前沿传热分解示意

Figure  5.  Schematic diagram of heat transfer decomposition of hydrate film front

下载: 全尺寸图片 幻灯片

M. W. Fishenden等人^[13]进行了矩形平板的自然对流传热过程实验研究，对实验数据拟合得到特征数n=0.2，大部分学者普遍采用该值进行计算。W. M. Lewandowski等人^[14]在研究超薄恒温水平板的自然对流传热过程中发现，当$L \ll {\rm{1}}$时，特征数m取值0.982，因此，对于微米级的水合物膜，特征数m可近似为1.0。同理，对于水平圆柱自然对流传热，特征长度$L = 2\delta$，m=0.50，n=0.25^[15]。把这些特征系数代入式（4），可以得到气相和水相传热系数；然后，将式（4）代入式（1），可以得到水合物膜横向生长速率与过冷度、水合物膜厚度、体积膨胀系数的关系：

式中：${\lambda _{\rm{g}}}$，${\lambda _{\rm{w}}}$分别为气体和水的导热系数，W/(m·K)；${a_{{\rm{v}},{\rm{g}}}}$，${a_{{\rm{v}},{\rm{w}}}}$分别为气体和水的体积膨胀系数，T^–1；c_{p, g}，c_{p, w}分别气体和水的比热容，J/(kg·K)；ρ_g，ρ_w分别为气体和水的密度，kg/m³；$g$为重力加速度，m/s²；μ_g，μ_w分别为气体和水的黏度，mPa·s。

式（5）中等式右边第一项代表水合物膜前沿向周围水相的局部自然对流传热对水合物膜横向生长速率的影响；第二项表示水合物膜前沿向周围气相的局部自然对流传热对水合物膜横向生长速率的影响。研究表明，水合物膜的厚度与过冷度成反比^[16]，则水合物膜厚度用过冷度表示为：

式中：k为反比例系数。

将式（8）代入式（5），可以得到水合物膜横向生长速率与过冷度的关系：

2.2   考虑水合物相变的气泡上升速度

气泡上升速度是预测避台关井期间钻井安全作业周期的关键参数之一。在深水钻井作业中，气泡上升速度不仅会受到钻井液物性参数（密度、黏度和表面张力等）的影响，同时与井筒温度和压力密切相关。地层气体从裸眼段附近岩石孔隙进入井筒后，以气泡的形式向上迁移，并迅速从受浮力作用显著的不稳定状态过渡到重力、浮力和流动阻力相平衡的稳定状态。气泡达到平衡状态后，可以根据受力平衡求解气泡上升速度：

式中：v_b为气泡上升速度，m/s；C_D为气泡的拖曳力系数；D_eq为气泡的当量直径，m；ρ_l为流体的密度，kg/m³；ρ_b为气泡的密度，kg/m³；E为气泡变形率。

气泡在水合物稳定区域外迁移时，气泡直径会随着周围温压环境而变化，根据气体状态方程，可以得到不同井深处的气泡直径：

式中：D_eq,H为井深H处气泡的气泡当量直径，m；Z_H为井深H处的气体压缩因子；p_H为井深H处的井筒压力，Pa；T_H为在井深H处的井筒内温度，K；D_eq,0为井底处气泡的气泡当量直径，m；Z₀为井底处气体压缩因子；p₀为井底压力，Pa；T₀为井底温度，K。

地层气体通过井壁周围的岩石孔隙侵入井筒，气泡在裸眼段的初始直径可以利用Wilkinson模型^[17]计算得到：

式中：σ_l为钻井液的表面张力，N/m；μ_l为钻井液黏度，mPa·s。

除气泡直径外，气泡上升速度还与气泡密度、变形率和拖曳力系数等3个参数密切相关。

在水合物稳定区域外，由于气体密度远小于水的密度，可以忽略气泡重力的影响，气泡运动主要受浮力和拖曳力控制；当气泡进入水合物稳定区域后，由于气泡表面覆盖了一层水合物薄膜，且水合物的密度略低于水的密度，不可忽略水合物气泡的重力，因此，气泡在钻井液中受到浮力、重力和拖曳力的共同作用，此时，气泡的密度可以用水合物气泡的当量密度表示：

式中：ρ_eq为水合物气泡的当量密度，kg/m³；ρ_g为气体密度，kg/m³；δ为水合物膜厚度，m。

3.   模型重要参数的确定及应用

水合物气泡的密度与水合物膜厚度相关，可以利用水合物膜横向生长速率模型确定；气泡变形率和拖曳力系数可通过实验数据回归的经验公式得到。

3.1   水合物膜生长速率和初始厚度

在不同过冷度（2～10 K）下水合物膜沿甲烷气泡表面的横向生长速率为0.05～0.85 mm/s。目前，经典传热模型包括Uchida的稳态热传导模型^[18]、Mori的对流传热模型^[19]、Freer的传热和本征动力学耦合模型^[20]及Mochizuki的瞬态热传导模型^[21]。笔者结合实验数据，对自然对流传热生长模型与各经典模型水合物膜横向生长速度进行了对比（见图6）。

图  6  自然对流传热生长模型与经典模型模拟对比

Figure  6.  Comparison on the simulation results of natural convection heat transfer growth model and conventional model

下载: 全尺寸图片 幻灯片

从图6可以看出，Freer模型和Mochizuki模型的模拟结果相对偏小；Mori模型的模拟结果相对偏大；Uchida模型虽然能较好反映水合物生长速率随过冷度的变化规律，但在较低或较高过冷度下预测结果偏差较大；水合物自然对流传热生长模型与实验数据吻合较好。这说明与自然对流传热相比，强制对流理论高估了膜前沿向周围流体的传热速率，而热传导理论却低估了热扩散速率，模拟结果验证了水合物自然对流传热生长模型的正确性。

由于水合物膜厚度与过冷度成反比，水合物膜厚度可以根据反比例系数k来估算。根据实验数据回归得到参数k=3.10835×10^–5 m·K，计算得到过冷度0.5～3.0 K时甲烷气泡表面水合物膜厚度为8.5～60.0 μm。Li Shengli等人^[22]利用光学显微镜测量得到相同过冷度下甲烷水合物膜的初始厚度为8.0～85.0 μm。水合物膜厚度的预测结果与实验结果之间的误差分析结果表明，预测的水合物膜厚度相对误差在±20%以内，进一步验证了模型的准确性（见图7）。

图  7  水合物膜厚度预测的误差分析

Figure  7.  Error analysis on the predicted thickness of hydrate film

下载: 全尺寸图片 幻灯片

3.2   水合物气泡变形率和拖曳力系数

水合物气泡的上升速度随气泡直径的变化情况如图8所示（Ⅰ为黏性力控制阶段，Ⅱ为惯性力控制阶段）。

图  8  水合物气泡上升速度与气泡当量直径的相关性

Figure  8.  Correlation between the ascending velocity and equivalent diameter of hydrate bubble

下载: 全尺寸图片 幻灯片

从图8可以看出，水合物气泡上升速度的变化规律与N. K. Bigalke等人^[11-12]的研究结果基本一致。同时，与“干净”气泡上升经历的控制阶段^[9]（包括黏性力控制阶段、表面张力控制阶段和惯性力控制阶段）对比，可以发现：1）在黏性力控制阶段，水合物气泡的上升速度与“干净”气泡的上升速度均随气泡当量直径呈线性增大，这是因为气泡直径较小时这2种情况的气泡呈球形，表面波的振荡影响较小；2）在表面张力控制阶段，水合物气泡仍然随气泡直径呈线性趋势增大，且小于相同直径下“干净”气泡的上升速度，然而“干净”气泡的上升速度随气泡直径增加而先增大，达到极值点后开始逐渐下降，并且数据点呈现分散的趋势，这是因为“干净”气泡在该阶段受表面张力影响显著，气泡呈现不规则的球冠或椭球冠形，表面流动性较强，受到的阻力分布不均匀，而水合物气泡由于水合物的形成压制了表面波振荡的影响，气泡表面更加稳定；3）在惯性力控制阶段，水合物气泡和“干净”气泡的上升速度均随气泡直径的增加缓慢增大，分散程度一致，这是因为在该阶段气泡受黏性力和表面张力影响较小，主要受惯性力控制。

由于水合物气泡表面生成了水合物膜，表面张力控制阶段消失。因此，水合物气泡随气泡直径变化时存在2个阶段：首先会经历黏性力控制阶段（Ⅰ），然后直接进入惯性力控制阶段（Ⅱ）。这2个阶段之间存在临界点，该实验条件下的分界点为气泡直径1.68 mm。

为了研究气泡变形率与钻井液密度、黏度和表面张力的关系，将气泡变形率与莫顿数M_o关联。实验中气泡变形率与莫顿数M_o之间的关系如图9所示。

图  9  水合物气泡变形率与莫顿数的相关性

Figure  9.  Correlation between the bubble deformation rate and Morton number of hydrate bubble

下载: 全尺寸图片 幻灯片

通过曲线拟合，得到了水合物气泡变形率的计算公式：

从图9可以看出，水合物气泡的变形率随莫顿数的增加而减小，这说明了随着钻井液黏度的增大，表面张力的减小，气泡受到的拖曳力越大，水合物气泡更容易变形。该公式与实验数据吻合较好，水合物气泡变形率预测的相对误差在±15%以内。

间接影响水合物气泡上升速度的另一个重要因素是拖曳力系数。由于水合物生成的影响，目前拖曳力系数经验关系式都不适于描述水合物生成区域内气泡的上升运动。为了计算水合物气泡的拖曳力系数，将其与雷诺数进行关联。实验中气泡拖曳力系数与雷诺数之间的关系如图10所示。

图  10  水合物气泡拖曳力系数与雷诺数的相关性

Figure  10.  Correlation between the drag coefficient and Reynolds number of hydrate bubble

下载: 全尺寸图片 幻灯片

通过曲线拟合，得到了水合物气泡的拖曳力系数的计算公式：

同样地，水合物气泡拖曳力系数也经历了2个阶段：黏性力控制阶段和表面张力控制阶段。雷诺数较小（50～150）条件下，拖曳力系数随过冷度增大而减小；雷诺数较大（150～1 000）时，拖曳力系数具有相反的变化趋势；雷诺数等于150时，水合物气泡拖曳力系数最小。

4.   实例分析

南海某深水油气井基本参数：井深4 800.00 m，气侵段深度4 650.00 m，海水深度1 500.00 m，海水表面温度20.0 ℃，储层温度61.5 ℃，储层压力72.8 MPa，地层孔隙度0.13，地层渗透率10 mD，地温梯度2.5 ℃/100m，钻井液密度1.38 kg/L，钻井液黏度40 mPa·s，压缩系数0.000 2 MPa^–1。应用气体上升速度模型及相关参数的子模型，结合南海某深水油气井的实际钻井参数，计算得到分别考虑和不考虑水合物相变2种情况下的侵入气体上升至海底井口装置的时间，即避台关井的安全作业周期，结果如图11所示。

图  11  侵入气体上升井深与时间的关系

Figure  11.  Relationship between the gas ascending depth with time in wellbore

下载: 全尺寸图片 幻灯片

从图11可以看出，随着井深的减小，曲线斜率逐渐增大，气体上升速度不断增加，但进入水合物生成区域（1 500.00～2 213.00 m井段）后，气泡上升速度突然减小，然后再逐渐增大。这进一步说明了水合物的生成显著降低了气泡上升速度，从而延长了气体上升至井口装置的时间；相应地，避台关井期间钻井安全作业周期从166 h增加至221 h。因此，考虑水合物相变的影响，可以准确预测钻井安全作业周期。

避台时间小于安全作业周期时，侵入井筒内气体还未上升至海底井口，在开井作业时相对安全；避台时间略微超过安全作业周期时，侵入气体已上升至井口，在顶部形成一定高压，在开井作业时应注意排气减压；避台时间远大于安全作业周期时，井口附近的圈闭气体由于水合物的形成、生长与聚集，极有可能造成防喷器及井筒堵塞，进一步增加了解封、开井作业时的井控风险。

上述研究表明，水合物生成虽然延缓了气体的上升，延长了安全作业周期，但是在气体到达井口后水合物生成增大了井控风险。因此，针对不同的避台时间，应采取不同的井控方法，例如当预计的避台时间略大于安全作业周期时，在撤离平台前应向井底注入一段高黏液体，降低气体上升速度，延长钻井安全作业周期；当预计的避台时间远大于安全作业周期时，在撤离平台前应向井口附近注入一定量的水合物抑制剂，阻止水合物的生成，在开井作业时立即关闭防喷器，进行节流循环减压，确保井控安全。

5.   结　论

1）基于自然对流传热理论，提出了水合物生成速率和初始厚度预测模型，与经典传热模型模拟结果对比发现，该模型与实验数据吻合较好，水合物膜厚度的预测误差在±20%以内。

2）水合物气泡上升过程中先进入黏性力控制阶段，然后进入惯性力控制阶段；水合物气泡变形率随莫顿数增加而减小；拖曳力系数随雷诺数的增大先逐渐减小后开始增大。基于这些特点，在关井前可向井筒内注入一段高黏液体，以延缓气体上升。

3）在解封、开井作业时，应根据避台时间判断海底井口附近是否存在圈闭高压气体和水合物生成堵塞等问题，并采取合理的井控措施，以延长安全作业周期和预防水合物的生成。