深水气井环雾流水合物沉积和堵塞动力学模型

董 钊，刁玉乾 ，李 中，蒋东雷，张攀锋

1.中海石油（中国）有限公司湛江分公司，广东 湛江 524057

2.中海油研究总院有限责任公司，北京 朝阳 100028

3.中国石油大学（华东）石油工程学院，山东 青岛 266580

引言

深水气田在测试和开发过程中，由于高压和海底低温，泥线附近井筒和海底管线内有可能形成水合物，造成堵塞，严重时会造成重大安全事故[1-2]。目前，主要采用注入大量热力学抑制剂的方法防止水合物的形成[3-7]，同时，由于对水合物形成和堵塞动力学认识不足，对水合物防治没有得到有效运用。Dorstewitz 等[8]进行了环雾流的环路模拟实验，观察到了水合物首先在水-气界面处形成的现象。Joshi 等[9]采用环路流动实验，探究了水合物形成和堵塞机理，发现水合物颗粒在气液界面处形成后分散于液膜内，部分滞留于液膜内，部分向管壁运移沉附，沉积层逐渐增加。Camargo 等[10]在建立的水合物沉积模型中，考虑了液膜内水合物沉积的临界浓度（packing 浓度）。Wang 等[11]提出了环雾流中水合物沉积的模型，认为液膜界面生成的水合物即时沉积，没有考虑水合物颗粒的运移和临界沉积。Rao等[12]基于实验证明了在饱含水的气相中水合物形成和沉积机理，认为与霜和冰的形成机理类似。

环雾流是深水气田井筒两相流最常见的流态，包含气芯和近管壁液膜两部分，液相以液膜和液滴两种形式存在，气相携带液滴形成在管道中部流动的气芯，水合物颗粒可在气芯内和气液界面上形成。本文基于环雾流流动模式，综合考虑水合物颗粒在气芯和液膜内的形成、运移和管壁黏附特性，建立井筒水合物沉积和堵塞的动力学模型，用于预测不同工况下的井筒水合物的管壁沉附厚度及堵塞时间。

1 模型建立

1.1 液膜内水合物形成、运移和沉积动力学模型

井筒环雾流水合物沉积和堵塞动力学模型概念如图1 所示，模型考虑了水合物在液膜形成、扩散运移及在井壁上动态沉积过程和水合物在气芯形成、聚并、破碎和沉降过程。

图1 井筒环雾流水合物沉积和堵塞动力学模型概念图Fig.1 Conceptual diagram of kinetics model for hydrate deposition and blockage in wellbore annular-mist flow

Lingelem 等[13]认为，当液相充足时，水合物会在气液界面处形成，而当液相不充足时水合物会在管壁上直接形成。Englezos 等[14]提出的水合物结晶理论认为，水合物主要在气液界面处很窄的薄层内形成，薄层厚度表达式为在管壁处形成，即液膜内形成的水合物全部沉积；而当液膜厚度大于δGhyd时，水合物在气液液面处形成，随着浓度增加和颗粒生长，水合物颗粒将向管壁发生扩散运移。水合物颗粒接触到管壁后，当黏附力大于移除力和黏附力偶大于移除力偶时，水合物颗粒沉积在管壁上，否则继续滞留于液膜内。当液膜内的水合物浓度大于临界浓度时，液膜单元整体失去流动性发生沉积，即移动床层流变为固定床层[6-7]。

可将液膜划分为3 个层：水合物形成层、水合物颗粒运移层和水合物颗粒沉积层，如图1b 所示。假设水合物颗粒的成核和生长动力学行为只发生在水合物形成层内，然后在水合物运移层内运移，当水合物颗粒运移至管壁沉积层时，发生沉积动力学行为。

1.1.1 水合物颗粒形成动力学模型

水合物颗粒成核后进入快速生长期，消耗传递进液相的气体，采用Turner 模型[15]计算总水合物形成速率，见式（2）。采用Englezos 模型[14]计算单个颗粒的生长速率，见式（3），将所有颗粒的生长速率累加可得总水合物形成速率，见式（4）。

虽然Turner 模型和Englezos 模型都能计算水合物形成速率，由于商业软件OLGA 中采用Turner模型计算水合物形成速率，所以本文采用Turner模型计算水合物形成速率。若，认为此时单位时间内水合物形成量为，水合物形成层的每个颗粒都发生了动力学生长行为，同时形成个临界晶核（临界晶核尺寸，m3）。若，认为此时单位时间内水合物形成量为，不是所有颗粒都发生生长行为，大尺寸的颗粒由于表面自由能低会优先生长。综上可得单位时间内，液膜内水合物形成量和每个颗粒的生长情况。

1.1.2 水合物颗粒运移动力学模型和特征

初始时，水合物颗粒尺寸和浓度较小，颗粒遵循菲克扩散定律向管壁缓慢扩散。随着水合物颗粒的生长，颗粒尺度达到湍流涡旋尺度，湍流扩散成为主要驱动作用，颗粒开始加快向管壁运移[16]，此时水合物形成层内的颗粒数量迅速减少，假设水合物颗粒离开形成层后，水合物颗粒停止生长，同时忽略水合物颗粒在液膜内的聚并和破碎行为，颗粒则以恒定尺寸向管壁运移。

水合物颗粒生长为湍流涡旋尺度le的表达式为[17]

随着水合物颗粒尺寸的生长，涡旋扩散作用逐渐增大，极大提高了颗粒的沉降通量，综合考虑对流作用和扩散作用，得到单位时间、单位面积内的沉降通量Jdep[18]，见式（6）。

当水合物颗粒尺寸小于le时，水合物的沉降通量大约10-7m3/(m2·s)，而当颗粒尺寸达到涡旋尺寸时，沉降通量达到10-3m3/(m2·s)数量级。可以假定当颗粒尺寸达到涡旋尺度后，颗粒瞬时移出水合物形成层，所以液膜内水合物颗粒的最大尺寸就是湍流涡旋尺寸le。

为了计算颗粒在液膜内的沉降通量，需要得到颗粒的径向运移速度upy、涡旋扩散系数εpt和液膜内水合物浓度分布。

通过分析颗粒在水合物形成层内的受力，得到颗粒径向运移速度upy。

根据湍流边界层理论[18]，涡旋扩散系数可表示为

由于颗粒尺寸大小不均一，而扩散运移过程会随着颗粒尺寸不同而有差异，所以对于液膜内颗粒采用多尺寸颗粒扩散理论。尺寸相同或相近的颗粒遵循相同扩散运移规律，针对不同尺寸的颗粒分别建立对流扩散方程，由于纵向无滑脱，径向速度较小，可以忽略颗粒间的碰撞和聚并等影响，所以可以近似认为不同尺寸的颗粒相互解耦。通过求解每一种尺寸颗粒的浓度分布，得到液膜内颗粒总的浓度分布。

不同尺寸颗粒扩散系数和沉降速度随时间和空间改变而改变，沉降通量Jdep是与时间和空间有关的物理量。将液膜沿径向方向划分为Ny个控制微元，（δlf—液膜厚度，m），e控制微元内某一尺寸颗粒浓度随时间变化可以表示为

1.1.3 水合物颗粒沉积动力学模型和特性

水合物颗粒运移至管壁不一定发生沉积，因为水合物颗粒受到的力（力偶）可能不足以使其黏附在管壁上。水合物颗粒运移至管壁处，受到范德华力Fvaw、液桥FaB等吸引力，同时也会受到Saffman 升力Fs、相间曳力Fd、净重力G等移除力，所以颗粒在管壁是否能够成功黏附需要分析黏附力（黏附力偶）与移除力（移除力偶）大小关系[19]，如图2 所示。

图2 颗粒在管壁黏附示意图Fig.2 Schematic diagram of particle adhesion on tube wall

（1）水合物颗粒与黏附层所有分子或原子吸引力总和为范德华力

（2）水合物颗粒具有强亲水性，颗粒与管壁或者颗粒与颗粒之间很容易形成液桥，毛细液桥力的表达式为

毛细液桥力使水合物颗粒变形，即便液桥蒸发掉，变形仍然存在，变形后的颗粒与表面接触面积增大，黏附力增大。

（3）黏性流体流经静止的颗粒时，由于相间摩擦作用，对颗粒施加曳力和附加力偶影响，使颗粒具有移除的潜能，O′Neill[20]推导出了黏性流体流经与固定平面接触的固定球体时产生的相间曳力和附加力偶

（4）Saffman 升力

（5）净重力

结合颗粒在管壁附近的受力分析，竖直管中颗粒在壁面处的移除机制包括径向拉升、轴向运移和滚动移除3 种。

（1）径向拉升FS≥(Fvaw+FaB)

当颗粒所受的Saffman 升力大于黏附力时，颗粒沿径向远离壁面。

（2）轴向运移Fd≥G+fs(Fvaw+FaB-FS)

当颗粒所受的曳力大于重力和与壁面之间的静摩擦力之和时，颗粒发生轴向运移（fs—静摩擦阻力系数，无因次，fs=0.2）。

（3）滚动移除Md+Fdl1+Fsl2≥Gl1+Fvawl2

当颗粒受到的移除力偶大于黏附力偶时，颗粒发生滚动移除。

根据3 种移除机制得到颗粒在壁面的临界移除粒径，即大于此粒径的颗粒不发生黏附，而是滞留在液膜内，从而导致液膜内水合物颗粒浓度持续增大。

图3 是产气量20×104m3/d，产液量5 m3/d 时不同粒径的颗粒受到的力和力矩对比图，可以发现水合物颗粒被移除的主要原因是滚动移除，临界移除粒径为30µm。

图3 不同粒径颗粒受力和力矩图Fig.3 Diagram of force and torque of different particle sizes

1.2 气芯内水合物形成、聚并、破碎和沉降动力学模型

气芯中液滴被气相包裹，未形成水合物前，液滴尺寸满足对数正太分布，假设液滴一旦进入水合物风险区域便瞬时转变为水合物颗粒，则初始尺寸为dp的水合物颗粒的数量密度函数为

水合物颗粒尺寸分布的对数期望

当液滴进入水合物风险区域后即刻形成水合物颗粒，所以初始控制单元内的水合物体积为：VChyd=（Qld—气芯内液滴体积流量，m3/s；uc—气芯速度，m/s）。

控制单元内初始水合物颗粒的总数量：NChyd=，所以，尺寸为dp的水合物颗粒数量为Ndp=fn(dp)NChyd。

水合物颗粒在气芯中运移时会发生多种动力学行为，如聚并、破碎、形成和沉降等，如图1c 所示。每种动力学行为均服从统计学规律，控制单元内的各个动力学事件发生的概率都有一个定量描述模型（核模型），水合物颗粒在风险区域内流动时，尺寸为dp的颗粒的数量满足群平衡方程，见式（20），研究重点是明确颗粒尺寸分布函数随时空的演变规律。本文研究的水合物颗粒形态均为球形，颗粒聚并后形成的聚集体为球形，颗粒破碎后形成的子颗粒也为球形。

由于气芯速度远大于扩散速度，所以扩散项可以忽略。对流项表征井筒不同深度处气芯流速差异的影响。聚并项包含两部分，尺寸的两种颗粒聚并形成dp尺寸颗粒，尺寸dp的颗粒与其他尺寸的颗粒聚并消耗dp尺寸的颗粒数量。破碎项包括dp尺寸的颗粒破碎消耗其数量和（大于dp）尺寸颗粒破碎形成dp颗粒。沉降项是水合物颗粒沉降入液膜，形成项是液膜雾化的小液滴转化为水合物颗粒。求解群平衡模型的关键是得到描述水合物颗粒动力学行为的核模型，包括形成核Ge、聚并核，破碎核和沉降核DEP。

1.2.1 水合物颗粒形成

气液接触面呈波浪状，较高的气芯流速对液膜有较强的剪切力，使液团从液膜处剥离，在气芯中进一步被分散，从而使液膜发生雾化。由于液膜雾化是气芯中液滴的来源，所以可根据液膜雾化率得到气芯中水合物颗粒的形成核模型，同时液膜雾化形成的液滴尺寸应该与气芯中液滴尺寸满足同样的对数正太分布规律，则雾化形成的液滴的数量密度函数可用初始液滴的数量密度函数，且认为在水合物风险区域内液滴瞬时转化为水合物颗粒。

水合物形成核模型为

液膜雾化率[22]为

1.2.2 水合物颗粒聚并

聚并事件是涉及两个颗粒的动力学事件，两个颗粒发生碰撞的概率用碰撞核模型βij表示，碰撞使颗粒可能发生聚并也可能不发生聚并，发生聚并的概率用聚并效率aij表示。所以颗粒i和j的聚并核模型为{aijβij}（βij— 碰撞核，m3/s）。

根据碰撞的驱动力不同，颗粒碰撞机理主要分为两种：表征布朗扩散作用的布朗碰撞核和表征湍流作用的湍流碰撞核。根据文献[23]中整理的模型可得βij。

由于颗粒的布朗作用和气流的湍动作用使控制单元内的颗粒发生碰撞，发生碰撞的两个颗粒能否聚并成一个整体取决于受到的黏附力和剪切力大小关系。当液滴刚进入水合物风险区域时，液滴外围迅速形成水合物壳，具有颗粒性质，此时水合物颗粒之间的黏附力为液桥力FaB[24]。黏附力使颗粒聚并为整体，而气流对颗粒的流动剪切使颗粒碰撞后分离，颗粒所受的剪切力为

通过比较发生碰撞的两个颗粒所受到的黏附力和流动剪切力大小得到聚并效率

1.2.3 水合物颗粒破碎

破碎核模型包括尺寸为dpi的颗粒i发生破碎的概率为Bi、破碎产生的子颗粒数目bi和子颗粒群中尺寸为dpj水合物颗粒数量比例γij。颗粒发生破碎的机理主要有两种，剪切破碎和碰撞破碎。根据Serra 等[25]研究，颗粒i的破碎概率可表示为

假设颗粒破碎后子颗粒尺寸呈等分二元分布，即

1.2.4 水合物颗粒沉降

假设颗粒沉降过程是随机的，颗粒尺寸对沉降事件无影响，则沉降颗粒的尺度分布与气芯中颗粒的尺度分布相同。沉降率[22]为

气芯中液滴转化为水合物颗粒后，水合物颗粒的沉降核模型为DEP=2.5πrcDepΔz/ρl。水合物颗粒沉降入液膜后，采用液膜内水合物颗粒运移和沉降规律分析。

2 模型求解及与实验结果对比验证

2.1 模型求解方法和计算流程

首先，利用商业软件（如Pipesim）及相关成熟模型计算环雾流井筒温压场、液膜和气芯的基础参数，与水合物相态曲线耦合，得到水合物形成的过冷度沿泥线井口以下井深变化规律。然后，根据1.1 所述的模型和公式计算单位时间内液膜内水合物形成量、每个颗粒生长、运移和沉积参数；根据1.2 所述的模型和公式计算单位时间内气芯内每个水合物颗粒形成、聚并、破碎和沉降参数。依次得到当前时刻（步骤n）液膜内水合物浓度、气芯内水合物颗粒尺寸分布、水合物管壁沉积量（水合物层厚度），然后开始下一时刻（步骤n+1）的计算。具体计算流程如图4 所示。

图4 环雾流水合物沉积和堵塞动力学模型计算流程图Fig.4 Flow chart for calculation of kinetic model of annular-mist flow hydrate deposition and blockage

由于液膜和气芯内水合物颗粒数量庞大，为简化计算，采用虚拟颗粒理念，即将尺寸相同或相近的多个实际颗粒整合为一个或一组虚拟颗粒，减少实际颗粒的计算数量。

2.2 模型预测结果的验证

Lorenzo[26]利用环路流动装置进行了环雾流下水合物沉积堵塞实验，实验参数见表1。

表1 环路流动装置进行的环雾流水合物沉积堵塞实验参数表[20]Tab.1 Experimental parameter table of annular-mist flow hydrate deposition and blockage of flow loop device

水合物在管壁上沉积，造成管壁缩径，增大流动压降。通过对比管路压降随时间变化的实验结果，对本文所建立的模型进行验证。图5 为模型模拟结果与实验结果的对比，二者的相对误差小于10%。

图5 模型模拟结果与环路水合物沉积堵塞实验结果对比图Fig.5 Comparison between the simulation results of the model and the experimental results of the hydrate deposits in the loop

3 模型应用案例分析

利用本文建立的水合物沉积和堵塞动力学模型，可以模拟深水气井开发过程中低产量生产时水合物在井筒管壁的沉积堵塞情况，由此可以明确水合物沉积堵塞风险程度，若是可以保证在井筒未被水合物完全堵塞时，产量恢复高产，随后使沉积的水合物热融解，则整个过程就可以少注甚至不注热力学抑制剂。

目标气田为南海LS17–2 深水气田，水深为1 336 m，泥线温度为3°C，气藏埋深为2 073 m、气藏压力为36 MPa。目标气田气体组分C1=93.33%，C2=4.33%，C3=0.98%，其他组分1.36%。

图6 是产气量5×104m3/d、水气比0.010%（产水量5 m3/d）时，不同时间水合物层厚度沿井筒分布曲线。由于液膜内部分水合物会随液膜流动，当水合物颗粒浓度大于临界沉积浓度时，该控制单元内液膜整段沉积；同时气芯内水合物运移和管壁黏附是一个随机的过程，导致水合物层厚度沿井筒分布的曲线不是光滑的。水合物层最大生长速率为8.28 mm/h，最大沉积厚度发生在1 524 m 处（海底井口下190 m），以此计算大约需要5.82 h 会造成井筒完全堵塞。

图6 不同时刻水合物层厚度随深度变化曲线Fig.6 Curve of hydrate layer thickness with depth at different times

表2 和图7 为水气比一定（0.010%）时，不同产气量下水合物沉积堵塞风险。

表2 不同产气量下水合物沉积堵塞风险Tab.2 Hydrate blockage risk under different gas production rates

图7 不同产气量下40 min 时水合物层厚度随深度变化曲线Fig.7 The variation curve of hydrate layer thickness with depth under different gas production at 40 min

随着产气量的增大，堵塞井筒的时间越长，发生完全堵塞的位置越靠近海底井口。因为产气量越大，水合物风险范围减小，在水合物稳定区内过冷度小，生成的水合物量少，同时部分水合物容易被带出井筒，所以水合物沉积的少，堵塞井筒的风险较小。当产气量为12×104m3/d 时，水合物在井筒内的沉积厚度为0，堵塞风险消失。

产气量一定（Qg=5×104m3·d-1）时，水气比对水合物层沉积厚度影响的模拟结果如表3 和图8所示。

表3 不同水气比下水合物堵塞风险Tab.3 Hydrate blockage risk under different WGR

图8 不同水气比下40 min 时水合物层厚度随深度变化曲线Fig.8 The variation curve of hydrate layer thickness with depth under different WGR at 40 min

水气比增加，水合物层增长速率有所增加，但由于产水量越大（携带热量越大），井筒内水合物形成的范围越小。模拟结果表明，当产水量较小（水气比为0.001%）或者较大（水气比为0.100%）时，堵塞位置靠近泥线井口。水气比为0.001%时，水合物层生长速率较小，大于12.30 h 可堵塞井筒，水气比为0.010%和0.100%时，堵塞时间大约为5.85 h。

4 结论

（1）基于深水气田井筒环雾流模式，在液膜内考虑了水合物颗粒的生长、运移和管壁黏附现象，同时在气芯内考虑了水合物颗粒的形成、聚并、破碎、沉降行为，建立了环雾流水合物动态沉积和堵塞的动力学模型，可用于模拟计算不同工况下的井筒内水合物层增长速率、水合物沉积厚度分布和堵塞时间，模型模拟结果与环路实验结果的相对误差小于10%，验证了模型的可靠性。

（2）典型深水气田（水深1 336 m、泥线温度3°C、气藏埋深2 073 m、压力36 MPa）案例分析表明，在产气量5×104m3/d、水气比0.010% 工况下，水合物沉积层生长速率高达8.28 mm/h，泥线井口下190 m 处水合物层最厚，井筒完全堵塞所需时间为5.82 h。水气比降至0.001%时，堵塞时间增加到12.30 h。产气量增加，海底井口下水合物形成的风险区域范围减小，水合物沉积厚度降低，堵塞位置向海底井口处靠近。当产气量为12×104m3/d 时，水合物沉积厚度降为0，堵塞风险消除。