固态流化采掘海洋天然气水合物藏的多相非平衡管流特征

魏 纳 赵金洲 孙万通 周守为 张烈辉 李清平付 强, 吕 鑫 郑利军

1. “油气藏地质及开发工程”国家重点实验室·西南石油大学 2.中海油研究总院有限责任公司

0 引言

为了安全、高效开发天然气水合物（以下简称水合物）藏，笔者所在的研究团队提出了一种创新型技术思路：海洋水合物藏固态流化采掘方法[1-7]。

然而，海洋水合物藏固态流化采掘技术的研究尚处于初期阶段，对一些理论难题还需要深入研究。特别是，水合物固相颗粒在密闭管线中随海水向上输送至海面平台的过程中，管线（即井筒）中的温度不断升高、压力不断减低，水合物将会发生分解，分解产生的大量气体使井筒中的流动变为复杂多相流动。而复杂多相流动对水合物颗粒的动态分解会进一步产生影响。此过程中管流特征参数变化复杂，为多相非平衡管流，进而带来的井控、固相颗粒输送等安全风险加剧。

目前已经有很多学者开展了井筒多相流动特征研究，主要集中在流型划分与预测[8-11]、井筒压力计算与测量[12-14]方面，未有研究耦合水合物分解的井筒多相流动特征。并且，针对水合物分解的研究主要集中在分解温度和压力条件[15-19]、分解速率模型与预测[20-23]方面，未有研究水合物颗粒在多相管流上升条件下随温度不断升高、压力不断降低过程中的动态分解规律。因此，笔者通过建立井筒温度和压力场数学模型、多相上升管流中的水合物动态分解数学模型、耦合水合物动态分解的井筒多相流动数学模型，研究认识海洋天然气水合物藏固态流化采掘中的多相非平衡管流规律及相应的数值计算预测方法。

1 多相非平衡管流数学模型

1.1 井筒温度模型

基于能量守恒定律和导热问题基本方程，考虑海洋天然气水合物藏固态流化采掘工艺过程，建立了井筒中混合流体温度分布模型：

式中ρm表示井筒中混合流体密度，kg/m3；vm表示井筒中混合流体速度，m/s；cm表示井筒中混合流体比热容，J/(kg·K)；Dpi表示井筒内径，m；Tm表示井筒中混合流体温度，K；z表示井深，m；qw表示海水与井筒之间的热交换，W/m；qf表示井筒中混合流体流动摩擦产生的热量，W/m；qh表示水合物颗粒上升过程中发生分解时的相变热量，W/m；qw、qf、qh的求解参阅本文参考文献[24-26]。

1.2 井筒压力模型

基于海洋天然气水合物藏固态流化采掘工艺过程，建立了混合流体流动过程中的井筒压力模型：

式中pm表示井筒压力，MPa；pG、pF、pA分别表示重力、摩阻、流体速度变化产生的压降，MPa/m；pG、pF、pA的求解参阅本文参考文献[27-29]。

1.3 水合物相平衡模型

为了判断水合物固相颗粒在上升过程中某一温度、压力条件下是否发生分解，需要建立水合物相平衡模型。为了简化计算，笔者采用Dzyuba和Zektser[30]通过实验分析所建立的水合物相平衡模型：

式中peq表示某一温度条件下的水合物相平衡压力，MPa。

1.4 多相上升管流中的水合物动态分解模型

在上升管流中，根据水合物相平衡模型，水合物固相颗粒上升至分解临界位置前，其未发生分解，此过程中分解速率为零；水合物固相颗粒上升至分解临界位置后，其将发生分解。为了计算水合物分解情况，假设天然气水合物为甲烷水合物，以Kim模型[20]为基础，建立了多相上升管流中的水合物动态分解模型：

式中nhyd表示水合物固相颗粒中水合物物质的量，mol；thyd表示水合物分解时间，s；khyd表示水合物分解速率常数，mol/(s·m2·MPa)；Shyd表示水合物固相颗粒分解表面积，m2；feq、fm分别表示甲烷气体在Tm和peq、pm条件下的逸度，MPa。其中各参数的求解如下：

1.4.1 水合物分解速率常数

水合物分解速率常数存在如下关系：

式中khydc表示水合物自身分解反应速率，mol/(s·m2·MPa)；khydf表示甲烷气体在 Tm、pm下的传质速率，mol/(s·m2·MPa)。根据 Arrhenius方程，水合物自身分解反应速率计算如下：

式中khydc0表示水合物本征分解速率常数，mol/(s·m2·MPa)；Eact表示水合物分解反应活化能，J/mol；R表示通用气体常数，取8.314 J/(mol·K)。

假设分解过程中水合物固相颗粒为均质类球体，根据经验公式[31]，颗粒球体向流体中的传质关系为：

根据传质学相关理论可得：

进而得到甲烷气体在Tm、pm下的传质速率：

式中nsh表示舍伍德数；Re表示雷诺数；nsc表示施密特数；ds表示水合物固相颗粒直径，m；khydf0表示传质系数，m/s；vm表示井筒中混合流体速度，m/s；ρm表示井筒中混合流体密度，kg/m3；pm表示井筒压力，MPa；Tm表示井筒中混合流体温度，K；DAB表示甲烷气体向管流中的扩散系数，m2/s；μm表示井筒中混合流体黏度，Pa·s；Zg表示天然气压缩因子，无因次。

1.4.2 水合物固相颗粒分解表面积

假设水合物固相颗粒分解发生在界面处[32]，根据前文假设其为均质类球体。因此，分解表面积如下：

式中φ表示球形度，无量纲；Vhyd表示固相颗粒中含水合物的体积，m3；Mhyd表示水合物摩尔质量，kg/mol；ρhyd表示水合物密度，kg/m3；Ehyd表示固相颗粒中含水合物的丰度；Vs表示固相颗粒体积，m3。

1.4.3 甲烷气体的逸度

根据气体逸度定义可得[33]：

式中Ratm表示以标准大气压表示的通用气体常数，取0.082 06 atm·L/(mol·K)（1 atm=101.325 kPa）；Tx表示环境温度，K；px表示环境压力，atm；p*表示参考态压力，atm；fx表示甲烷气体在Tx、px下以标准大气压表示的逸度，atm；Vx表示甲烷气体在Tx、px下的摩尔体积，L/mol。由于R—K状态方程计算较为简单且一般具有较高的准确性。因此假设甲烷气体物态服从R—K状态方程，则：

式中a、b表示R—K常数，a提供分子间吸引力的度量，b提供分子大小的尺度。

因此，推导得到甲烷气体逸度计算公式：

1.5 耦合水合物动态分解的井筒多相流动模型

气相和液相质量守恒方程为：

井筒中气—液—固相混合动量方程为：

式中Dpi表示井筒内径，m；ρg、ρl、ρs分别表示井筒中气、液、固相密度，kg/m3；vg、vl、vs分别表示井筒中气、液、固相速度，m/s；Eg、El、Es分别表示持气率、持液率、固相含量，无量纲；mhydg、mhydl分别表示气、液相在单位长度上的质量变化速率，kg/(s·m)。

2 数值计算方法

2.1 边界和初始条件

为了简化模型，假设在水平段水合物固相颗粒均能安全运移，且运移速度等于液相速度，即

海底水平井段井筒中的混合温度等于外界海底温度，即

井口处的压力等于井口回压，即

式中zHo表示海底水平井段的井深，m；vs(zHo,t)、vl(zHo,t)分别表示海底水平段井筒中固相和液相速度，m/s；zH表示海底水平井段的垂深，m；Tm(zH,t)、Tw(zH)分别表示海底水平段井筒中混合流体温度和外界海水温度，K；pm(0,t)、pb分别表示井口压力和井口回压，MPa。

2.2 数值计算方法

为了研究海洋天然气水合物藏固态流化采掘多相非平衡管流特征，笔者采用有限差分数值计算方法。空间域为井筒，时间域为水合物固相颗粒从井底被输送至井口。假设在n时刻，井筒内任意两个节点i与i+1的参数已知，以节点i与i+1从n到n+1时刻为例说明数值计算过程（图1）。图1中α、β、γ为计算误差精度。

图1 固态流化采掘多相非平衡管流数值计算流程图

2.3 数值计算验证

为了验证海洋天然气水合物藏固态流化采掘多相非平衡管流数值计算的准确性，采用所建立的多相非平衡管流数值计算方法，得到了井筒温度、压力、持气率、固相含量、气—液—固相速度曲线如图2所示。将数值计算结果与周守为等[5,34,35]所做的海洋天然气水合物固态流化采掘大型物理模拟实验数据进行对比，文献中实验基础数据如下：海水深度为1 315 m，井深为1 500 m，实验管段井筒直径为0.076 m，模拟海面温度为308 K，固相颗粒中含水合物的丰度为70%，模拟施工排量为0.008 m3/s，液相密度为1 030 kg/m3，固相输送量为0.54 m3/h。

图2 数值计算值与文献中的实验值对比图

由图2对比可知，所得到的数值计算值与文献中的实验值误差较小，变化趋势一致。为了更进一步定量分析数值计算值与文献中的实验值的误差大小，采用如下误差定义式：

式中σ表示误差，无量纲；xs表示数值计算值；xr表示文献中的实验值。

通过计算得到平均井筒温度误差为0.022 ℃，平均井筒压力误差为0.023 MPa，平均持气率误差为0.010%，平均固相含量误差为0.004%，平均气相速度误差为0.005 m/s，平均液相速度误差为0.003 m/s，平均固相速度误差为0.004 m/s。可以看出，笔者所建立的数值模型与计算方法具有较高的准确性，为下文开展多相非平衡管流特征的分析奠定了基础。

3 固态流化采掘多相非平衡管流特征

海洋天然气水合物藏固态流化采掘如图3所示，井深为1 400 m，垂深为1 000 m，垂直井段深0～830 m，斜井段井深830～1 100 m，水平段井深1 100～1 400 m；井筒外径为0.508 m，井筒内径为0.476 m，管输液相密度为1 030 kg/m3。南中国海天然气水合物藏特征：海水深度为1 000 m；海面温度为293 K；固相颗粒中含水合物的丰度为70%。

图3 海洋水合物藏固态流化采掘多相非平衡管流示意图

为了得到水合物固相颗粒在密闭管线中随海水向上输送至海面平台过程的多相非平衡管流特征规律，采用所建立的数学模型及数值计算方法，下面将开展不同液相排量、固相输送量、井口回压条件下的固态流化采掘多相非平衡管流特征分析。

3.1 不同液相排量下的管流特征

数值计算参数如下：液相排量分别为0.03 m3/s、0.05 m3/s、0.07 m3/s，固相输送量为1.5 m3/h（对应日产气量约为4 133 m3），井口回压为0.1 MPa。通过数值计算得到了海洋天然气水合物藏固态流化采掘中，不同液相排量下的井筒温度、井筒压力、多相上升管流中的水合物动态分解规律、耦合水合物动态分解的井筒多相流动特征，如图4所示。

从图4-a、b中可以看出，水合物固相颗粒向上输送至海面平台的过程中，在海底水平段1 400 m至1 100 m，井筒温度等于外界海水温度。在斜井段和垂直段1 100 m至0 m，井筒温度受外界海水温度升高的影响而升高；且随着液相排量的增大，井筒中的混合流体与海水换热时间缩短，海水向井筒中的传热量减小，因此井筒温度降低，进而通过水合物相平衡模型计算得到的相平衡压力降低（图4-b）。随着液相排量的增大，相同固相输送量下，井筒中固相含量（Es）降低（图4-i），造成井筒中混合流体密度降低，进而有使井筒压力降低的趋势；然而井筒中各相速度的增大（图4-d～f）会使摩阻压降增大，又有使井筒压力升高的趋势，因此，井筒压力变化不明显（图4-b）。

图4-b中，任一确定的液相排量条件下的井筒压力和相平衡压力曲线交点，即为水合物固相颗粒随海水向上输送中发生分解的临界位置。之后，随着水合物固相颗粒在井筒中不断上升，其物质的量不断减小（图4-c），至0时的位置即为水合物完全分解的位置。对比不同液相排量的影响，可以看出，随着液相排量的增大，水合物分解临界位置上移，完全分解位置上移。

从图4-d～f中可以看出，在海底水平段1 400 m至1 100 m，根据前文中假设，水平段水合物固相颗粒速度等于液相速度；井筒中不含气相，气相速度为0。在斜井段1 100 m至830 m，由于井筒中流道截面积不变（井筒内径不变），液相速度不变；固相速度逐渐降低。在直井段830 m至0 m，结合图4-b和4-c中的水合物分解临界位置，水合物未发生分解时，气相速度依然为0，液相速度和固相速度保持稳定；当水合物固相颗粒随海水向上输送至分解临界位置时，固相颗粒中的水合物分解为气体和水，因此，井筒中的液—固两相流动变为复杂的气—液—固多相流动，随着多相流不断上升，井深不断减小，气体不断膨胀体积变大，气相速度升高，进而其携带作用使液相速度升高，然而由于井筒中混合流体密度的降低，引起携固能力下降，因此，固相速度随井深减小而降低。同时，对比不同液相排量的影响，可以看出，随着液相排量的增大，井筒中的气、液、固相速度均增大。

图4 不同液相排量下的井筒温度压力、水合物动态分解与多相流动特征图

图4-g～i中，结合图4-d～f，可以看出，在海底水平段1 400 m至1 100 m，持气率为0，持液率和固相含量保持稳定。在斜井段1 100 m至830 m，固相速度逐渐降低，引起固相含量逐渐升高，对应井筒中持液率降低。在直井段830 m至0 m，结合图4-c，水合物未发生分解时，持气率为0，持液率和固相含量保持稳定；当水合物固相颗粒随海水向上输送至分解临界位置时，固相颗粒中的水合物分解为气体和水，因此，水合物分解过程中，持气率升高、持液率升高、固相含量降低；当水合物固相颗粒随海水向上输送至完全分解的位置时，持液率不再升高、固相含量不再降低，且随着井深的减小，气体膨胀体积变大，持气率升高，造成井筒中混合流体密度降低，进而引起携固能力下降，因此，固相含量升高，对应持液率则降低。同时，对比不同液相排量的影响，可以看出，随着液相排量的增大，井筒中的持气率减小、持液率增大、固相含量减小；当液相排量最小为0.03 m3/s时，在井口附近（井深接近0），固相速度突然升高（图4-f）、固相含量突然降低（图4-i），说明其他条件一定、排量较小的情况下，井口处由于持气率过高，造成多相流流型发生了转变，进而造成混合流体携固能力的转变，同时，流型转变对井控安全造成了较高的挑战。

综合图4可以看出，海洋天然气水合物藏固态流化采掘现场施工中，为了防止井控安全问题、提高水合物固相颗粒输送能力，应适当提高液相排量。

3.2 不同固相输送量下的管流特征

数值计算参数如下：固相输送量分别为1.5 m3/h、3.0 m3/h、4.5 m3/h（对应产气量约为4 133 m3/d、8 266 m3/d、12 399 m3/d），液相排量为0.03 m3/s，井口回压为0.1 MPa。通过数值计算得到了海洋天然气水合物藏固态流化采掘中，不同固相输送量下的井筒温度、井筒压力、多相上升管流中的水合物动态分解规律、耦合水合物动态分解的井筒多相流动特征，如图5所示。

图5-a、b中，随着固相输送量的增大，海水向井筒中的传热量变化不大，因此井筒温度变化不明显，进而计算得到的相平衡压力变化不大。图5-b中，更高的固相输送量会使单位时间内有更多的固相进入井筒，有使井筒中混合流体密度升高的趋势。在水合物未发生分解的井段（井筒的下部），井筒中为液—固两相流动，因此井筒中混合流体密度升高，使重力压降增大，进而使井筒压力在下部井段升高；而在水合物发生分解后的井段（井筒的上部），气相从水合物中分解出，井筒中变为气—液—固多相流动，而气相有使井筒中混合流体密度降低的趋势，综合固相含量的影响，因此井筒压力在上部井段随着固相输送量的增大变化不明显。同时， 水合物固相颗粒随海水向上输送中发生分解的临界位置随着固相输送量的增大变化不明显，完全分解位置同样变化不明显（图5-c）。

图5 不同固相输送量下的井筒温度压力、水合物动态分解与多相流动特征图

图5-d～i中，固相输送量越大，固相含量越高（图5-i）；当水合物固相颗粒随海水向上输送至分解临界位置后，水合物分解产生的气体量越多，持气率越高（图5-g），对应持液率越低（图5-h）；随着井深不断减小，固相输送量越大，气体膨胀后的体积越大，气相速度越高（图5-d），进而其携带作用使液相速度更高（图5-e）。同时，高固相输送量条件下，井口处持气率更高，造成多相流流型转变更剧烈，进而引起携固能力转变越大，因此，井口处固相速度升高越多（图5-f）、固相含量降低越多（图5-i），同时，对井控安全造成了更大的挑战。

综合图5可以看出，海洋天然气水合物藏固态流化采掘现场施工中，提高固相输送量可以增大单位时间内输送至海面平台的水合物量，进而提高天然气产量，然而所带来的井控安全问题会更加严重。因此，应该在井控安全的前提下才可提高固相输送量。

3.3 不同井口回压下的管流特征

数值计算参数如下：井口回压分别为0.1 MPa、1.0 MPa、2.0 MPa，液相排量为0.03 m3/s，固相输送量为1.5 m3/h（对应产气量约为4 133 m3/d）。通过数值计算得到了海洋天然气水合物藏固态流化采掘中，不同液相排量下的井筒温度、井筒压力、多相上升管流中的水合物动态分解规律、耦合水合物动态分解的井筒多相流动特征，如图6所示。

图6-a～c中，随着井口回压的增大，海水向井筒中的传热量变化不大，因此井筒温度变化不明显，进而相平衡压力变化不大；而井筒压力明显升高；水合物发生分解的临界位置和完全分解位置均上移。

图6 不同井口回压下的井筒温度压力、水合物动态分解与多相流动特征图

图6-d～i中，当水合物固相颗粒随海水向上输送至分解临界位置之前，井筒中为液—固两相流动，随着井口回压的增大，各相速度和各相含量均不变；至分解临界位置之后，水合物分解产生气体，井口回压越高，气体膨胀越小，气相速度和持气率越低（图6-d、g），液相速度越低（图6-e），持液率越高（图6-h），且井筒中混合流体密度越高，携固能力越高，固相速度越大（图6-f），固相含量越低（图6-i）。同时，较高的井口回压条件下，在井口处固相速度和固相含量（图6-f、i）变化越小，说明越有助于保障井控安全。

综合图6可以看出，提高井口回压可以显著提高井筒压力，降低井口处气体膨胀带来的井控安全风险。因此，海洋天然气水合物藏固态流化采掘现场施工中，应该适当施加井口回压。

4 结论

1）增大液相排量，可显著降低井筒温度和水合物相平衡压力，使水合物分解的临界位置上移；可显著提高气、液、固相速度，使持气率和固相含量降低。

2）随固相输送量的增大，井筒温度和水合物相平衡压力不变，井筒压力在下部井段升高、在上部井段变化不明显，水合物分解的临界位置不变，固相含量显著增大，井口处持气率显著增大、持液率显著减小、气液固相速度均增大。

3）增大井口回压，可显著提高井筒压力，使水合物分解的临界位置上移；可显著降低水合物分解后的持气率，且气、液相速度降低，固相速度增大，固相含量降低，持液率增大。

4）海洋天然气水合物藏固态流化采掘现场施工中，适当提高固相输送量可以提高天然气产量，然而井控风险等问题将会加剧；为了防止井控风险等井筒流动安全问题，并提高水合物固相颗粒输送能力，应同时适当提高液相排量、施加井口回压。