天然气水合物开采中冲缝管套冲蚀研究

邓福成，黄斌，尹彪，李小森，易先中

(1.长江大学机械工程学院，湖北荆州，434000；2.中国科学院广州能源研究所，广东广州，510640；3.中国科学院广州天然气水合物研究中心，广东广州，510640；4.中国科学院天然气水合物重点实验室，广东广州，510640)

天然气水合物的全球储量非常巨大，是全球常规燃料总碳量的2 倍[1]，天然气水合物的资源量为1.0×1017～1.0×1018m3[2]，被视为可替代石油的重要能源。目前已有多个国家对天然气水合物的开采进行了多次试采实验，但天然气水合物储层出砂严重是其大规模商业化开采中的瓶颈问题。水合物主要分布于冻土区和海洋深水沉积物中[3]，试采中普遍使用降压法、注热法和气体置换法[4]等开采方法，以降压法最为常见，但降压开采时储层物性变化会造成大面积出砂[5-7]。开采过程中砂岩脱落的砂粒随着流体进入井筒，对井下设备造成破坏，从而制约天然气水合物安全高效开采。为减小天然气储层出砂对开采的影响，人们采用机械防砂筛管来解决油气井开发中天然气储层的出砂问题。冲缝管套作为防砂结构中的重要元件，一般用于筛管保护套或独立防砂筛管。当出现冲蚀磨损时，管套防砂性能降低，影响油井生产。因此，分析冲缝管套的冲蚀寿命成为筛管研究的重要部分。

目前，国内学者通过建立冲蚀模型，研究筛管冲蚀磨损问题，该方法已趋于成熟，并取得了丰富的研究成果。张锐等[8]建立离散颗粒筛管冲蚀模型，运用气固耦合模型模拟两相流流场分布和冲蚀特征，为气井防砂的冲蚀计算提供了研究思路；翟晓鹏等[9]用离散颗粒流数值模拟方法分析不同流速和体积分数下砂粒对金属网布速度分布和冲蚀速率影响，预测金属网布筛管冲蚀寿命；张君涛等[10]在防砂筛管模拟中引入冲蚀模型，分析不同压差下割缝筛管的冲蚀特征，得到冲蚀的主要集中区域，对于长效防砂和提高筛管使用寿命有着重要作用；陈珊珊等[11]基于CFD 的固液两相流冲蚀模拟分析，预测筛管寿命，并利用室内冲蚀试验验证数值模拟预测冲蚀寿命的可行性。但这些研究都只考虑了部分筛网结构下两相流的冲蚀模拟，未见针对水气砂三相流模拟和冲缝管套结构等冲蚀磨损问题的研究。

首先，本文作者建立冲缝管套筛缝流道模型，基于多相流模型(Mixture)和离散相模型(DPM)来模拟不同因素下砂粒流动对管套冲蚀的影响；其次，分析水合物分解后，气、液和固三相流下流体流速、颗粒体积分数、开口高度及液相体积分数对冲缝管套冲蚀速率的影响；最后，推导冲缝管套寿命预测模型，为天然气水合物储层试采及商业化开采的长效安全性提供理论支持。

1 水合物储层物性特征

国内学者根据神狐海域地质情况建立了单层水合物层和水合物-游离气层2种成藏模型[12-13]，其中第2种模型更符合现场实际，即水合物层存在于上覆层和下伏层中，且下伏层与水合物层间存在含有大量离散甲烷气体的游离气层。天然气水合物常以固态形式存在，开采时主要通过降压、注热等方法使其分解，再将分解后的天然气运送至地面。水合物分解过程如下：

降压过程中储层物质存在着相变分解与再形成，导致沉积物中包含气、液和固三相物质[14]，但水合物一般作为胶结物赋存在于沉积物中。水合物的分解导致砂岩应力变化，从而致使砂岩颗粒脱落。脱落的颗粒在甲烷气和水共同渗流的作用下发生迁移，部分未完全分解的水合物微粒也会与地层中的砂粒一起流入井筒。机械筛管是一种水合物常用的防砂手段，在生产过程中，甲烷气和水将渗流进筛管保护套，且颗粒将随流体一起进入保护套，造成筛管的冲蚀。

2 计算模型

水合物储层生产过程中的冲蚀过程是气、液和固三相流作用下的冲蚀问题。根据储层中各物质特性，计算中将气体和液体视为连续相，颗粒为离散相。因此，水合物储层中流体携带颗粒冲蚀计算主要分为流体控制方程、颗粒运动方程及冲蚀模型[15]。

2.1 气、液流体控制方程

控制方程包括连续性方程、动量方程和能量方程，流体分为气体与液体，双流体的连续性方程为[16]：

式中：α为体积分数，%；ρ为密度，kg/m3；t为时间，s；∇为拉普拉斯算子；u为速度矢量，m/s；下标g和l分别表示气相和液相；αl+αg= 1。

动量方程为[17]：

式中：pl和pg分别为液相和气相压力，Pa；g为重力加速度，通常取9.8 m/s2；和分别为液相和气相的应力张量，Pa；M为相间动量交换系数，kg/(m·s)2；μl，μs和ξ分别为液体、气体和流体的剪切黏度，Pa·s；I为单位向量。

能量方程为[18]：

式中：σl和σg分别为液相和气相剪应力张量；el和eg分别为液相和气相的比热力学能，J/kg；K为导热系数，W/(m·K)；T为热力学温度，K；Q为内部热源。

在水合物储层开采中，水合物分解产生的气体与液体混合，从而形成水气两相流。模拟中采用混合流模型，忽略水气间的作用力。

2.2 固相颗粒运动方程

颗粒随流体冲蚀冲缝管套过程中，流体(气体和液体)视为连续相，颗粒为离散相。根据牛顿第二定律，颗粒的受力方程为[19]

式中：G为颗粒所受重力；FD，FΔp，FF，FV和FL分别为颗粒所受的拖曳力、压力梯度力、升力、附加质量力以及浮力。在颗粒运动中还存在许多力，如Besset 力、Magnus 力和Saffman 力等。Besset 力是瞬时的流动阻力，与流体的稳定性有关，当颗粒加速度较小时，可忽略该力。Magnus和Saffman力等都只在颗粒高速旋转的特定情况下才有意义，故受力分析中可忽略。

颗粒在流体中所受重力G为

式中：d为颗粒粒径，m；ρs为颗粒密度，kg/m3；Vs为颗粒体积，m3。

压力梯度力FΔp为

式中：p为压力，Pa。

流体内部压强差形成的浮力FL为

颗粒表面流体流速差导致的升力FF为

式中：CL为升力系数；v为流体作用在颗粒上的速度，m/s；As为颗粒平均投影面积，m2。

流体相对于颗粒运动，流体产生的曳力(气动阻力)FD为

式中：CD为阻力系数；νf为颗粒速度，m/s。

颗粒加速所需要的高于颗粒本身的附加作用力FV。

颗粒在冲蚀过程中不断撞击壁面，颗粒对壁面的作用力是颗粒与壁面的碰撞力和摩擦力，颗粒间碰撞力忽略不计。根据GRANT 等[20]的研究，得到颗粒冲击-反弹系数为：

式中：eN为切向反弹系数；eT为垂向反弹系数；θ为颗粒入射角度，rad。

2.3 颗粒冲蚀模型

在多相流模型计算中，欧拉-欧拉法将不同的相视为连续介质；而欧拉-拉格朗日法则将流体视为连续相，颗粒视为离散相。本文将颗粒视为离散相，运用动能、质量和能量传递作为桥梁来计算模型。在气、液和固三相模型计算中，需要用欧拉-拉格朗日法进行计算，离散相冲蚀模型计算公式[21]为

式中：Rerosion为壁面磨损速率，kg/(m2·s)；N为碰撞颗粒数；C(dp)为颗粒直径函数，C(dp)=1.8×10-9；f(ω)为冲蚀角度函数，如式(19)所示；ω为颗粒轨迹与壁面的冲击角度，rad；νsl为颗粒的相对速度，m/s；b(νsl)为颗粒相对速度函数，取2.6；Aface为冲蚀面积，m2。

3 冲缝管套冲蚀模型

3.1 冲缝管套筛缝模型

在天然气水合物的开采防砂过程中，冲缝管套用于保护筛管或者独立防砂，其与地层流体直接接触，易受到冲蚀磨损。冲缝管套结构如图1所示，根据其流道结构简化出单位流道模型结构，建立冲缝管套筛缝的单位流体域模型。入口面为矩形截面，长为22 mm，宽为3 mm；壁面高度为11 mm，出口高度根据冲缝管套的设计条件设置为不同尺寸。图2所示为单位筛缝流道模型及网格模型，图2中蓝色面网格(网格模型上平面)为流场入口面，绿色面网格为筛缝壁面，黄色面网格为流场出口面(开口面)。

图1 冲缝管套示意图Fig.1 Schematic diagram of punched screen

图2 单位筛缝流道模型及网格模型图Fig.2 Flow channel model and grid model of unit screen

3.2 模型条件设置

离散相采用DPM 模型，对壁面冲蚀进行气、液和固三相三维湍流稳态计算。计算中湍流模型采用RNGk-ε模型，运用DPM 模型对颗粒运动进行追踪，固体球形颗粒由入口面注入。流体为气相和液相混合，压力和速度耦合用Coupled 算法。流体模型为多相流模型，液相与气相的入口速度相同。当液相体积分数较小时，液相均匀分散在气相中，忽略水气间的相间作用力，将流道区域内流体拟定为气液混合流体(等效为在气液混合均匀的单相流体)，两者共享同一压力场和速度场，且忽略相对滑移速度。

气相为甲烷，考虑低温、高压对气体的影响，密度为1.225 kg/m3；液相为水，密度为998.2 kg/m3；固相为砂粒颗粒，密度为1 500 kg/m3，沉积物中大部分颗粒是粉砂，其颗粒粒径为4～63 μm[22]，取颗粒平均粒径为40 μm，粒径正态分布。由于分析的流道模型是管套内实际流道模型的基本单元，并沿轴线轴向分布，且分析模型尺寸较小，故不考虑重力影响。入口处边界条件设置速度入口、流体和颗粒的速度及方向。流动中不涉及热力学，即不考虑流动过程中水合物二次生成等现象，流动过程视为等温绝热。管壁条件采用无滑移壁面，出口边界设置outflow。

2017年我国在南海神狐海域第一轮试采中累计产气量超过3.0×104m3，日平均产量为5.0×103m3，最高日产量为3.5×104m3；第二次试采创造了产气总量为8.614×105m3，日平均产量为2.87×104m3的世界纪录[23]。假定该井的水平段长度为100 m，按日平均产量2.87×104m3可计算得到流体流速为0.125 m/s，计算时初始流体流速取0.100 m/s。根据申志聪等[24]建立的模型，水合物储层下方存在游离气层。若水平井处于水合物层中，开采初始阶段压差使得水合物分解和游离气层逸散出的混合气体流入管道。此时，气体流量将远大于初始流量。考虑到各因素对管道破坏的影响，合理选取不同因素的取值，如表1所示。

表1 模拟因素水平表Table 1 Simulation factor level table

在冲蚀计算中，由于流体流速及组成比例不同，颗粒数量根据体积来确定更适宜，颗粒体积分数(含砂体积分数)按在入口面处所占面积比进行计算，因流体流速发生变化，单位时间内进入筛缝流道的颗粒总体积随之改变。计算得到不同颗粒体积分数和不同流体流速下的颗粒总质量流量：

式中：Qm为颗粒总质量流量，kg/s；φs为颗粒体积分数，%；vs为颗粒在入口处速度，m/s；Aintel为入口面积，m2。

4 结果分析与讨论

4.1 流场特征

分析计算在入口流速为0.1 m/s、颗粒体积分数为5%、液相体积分数为0.25%和开口高度为0.30 mm 条件下颗粒对筛缝结构的冲蚀磨损情况。图3所示为上述条件下冲蚀速率分布云图。从图3可见，颗粒进入流道对底部壁面和斜端壁面造成冲蚀磨损，且底部壁面靠近开口处的冲蚀磨损速率最大。

图3 模型冲蚀速率分布云图Fig.3 Cloud map of model erosion rate distribution

冲蚀过程中颗粒会在底部与壁面发生碰撞[25]后经出口面流出，造成筛缝开口处的冲蚀磨损。图4所示为颗粒轨迹图。从图4可看出，部分颗粒撞击壁面后，颗粒随流体向出口运动。分析颗粒运动轨迹，发现颗粒在筛缝壁面和底部中线区域的停留时间长、速度小、轨迹复杂。

图4 颗粒轨迹图Fig.4 particle trajectory

4.2 冲蚀速率的影响因素

4.2.1 流体速度

由式(20)可知：当流速改变时，进入流道内的颗粒质量流量增大，单位时间内撞击壁面的颗粒数量增加、颗粒撞击频率升高，冲蚀速率增大。根据表1中数据，分析计算得到冲缝管套冲蚀速率。图5所示为不同颗粒体积分数下冲蚀速率随流速变化折线图，对数据进行拟合，冲蚀速率y和流速x的拟合关系为

式中：ɑ和b均为幂函数参数。图5中各组拟合线与相应数据的重合程度较好。筛缝上冲蚀速率随流速增加呈现指数增长趋势[26]。当流速增加时，颗粒获得更大动能，对壁面的冲击磨损程度增大，冲蚀速率增大。

图5 不同颗粒体积分数下冲蚀速率随流速变化折线图Fig.5 Line chart of erosion rate change with velocity at different sand volume fractions

4.2.2 颗粒体积分数

图6所示为不同流速下冲蚀速率随颗粒体积分数变化。从图6可得：冲蚀速率随颗粒体积分数增加，呈现线性增长趋势。由式(20)可知，随颗粒体积分数增大，颗粒质量流量增加。颗粒数量增多，颗粒与壁面碰撞频率上升，造成更严重的冲蚀磨损。SUN 等[13]研究发现：当颗粒质量流量达到一定值后，冲蚀速率随质量流量增加出现增长缓慢的趋势。颗粒增多导致颗粒间相互碰撞更加剧烈，动能减小，冲蚀速率降低。在本次分析中，未考虑颗粒间碰撞的因素。

4.2.3 液相体积分数

图7所示为不同液相体积分数下冲蚀速率随流速变化折线图。由图7可知：冲蚀速率随着液相体积分数增加而增大，呈现出线性增长的趋势。由于流道开口面积远小于入口面积，相同流量流体通过开口面的流速将远大于入口面的速度。流体中液相黏性大，影响颗粒运动，增加颗粒动能；而气相黏性小，颗粒速度的增加不明显[27]。因此，液相体积分数增加时，出口处颗粒动能增加，冲蚀速率增大。在计算中，未考虑重力的影响，故气液两相的密度对冲蚀速率影响不大，而两者黏度的差异造成气相更易从流道中逸出，残余的液相对颗粒运动产生影响。

图7 不同液相分数下冲蚀速率随流速变化折线图Fig.7 Line chart of erosion rate change with flow velocity at different liquid phase volume fractions

在流速为0.2 m/s 时，冲蚀速率随液相体积分数增加出现先减小后增大的现象。由于模型计算中2种流体为均匀混合，但两者黏度的差异导致混合流体在流经筛缝底部时，液相更容易附着于底部。气体逸出筛缝后，在底部冲蚀区域液体随着体积分数增加更易形成液膜。液相体积分数增加对冲蚀速率有着增大和抑制2个方面作用：

1)由于液相对颗粒动能变化起主导作用，随着液相体积分数增加，冲蚀速率增大；

2)当液相体积分数增加时，由于液相更易附着于壁面，部分冲蚀区域将形成液膜，减小颗粒动能，从而减弱颗粒对管套壁面的冲蚀作用。

4.2.4 出口面高度

图8所示为不同开口高度下冲蚀速率随流速变化折线图。由图8可见：当流速一定时，开口高度减小，壁面冲蚀速率增大。开口高度减小为0.10 mm时，随流体流速增加，冲蚀速率的增长速率增加得更为明显。当开口高度变小时，出口面积减小使得出口处流速增加、动能增大，碰撞壁面的颗粒数量增多。颗粒动能和碰撞数量增加造成冲蚀区域的冲蚀程度加剧。

图8 不同开口高度下冲蚀速率随流速变化折线图Fig.8 Line diagram of erosion rate change with flow velocity at different opening heights

开口高度分别为0.10 mm和0.15 mm时，堵塞现象导致冲蚀速率增长幅度明显。分析计算中颗粒粒径为0.02～0.06 mm，当粒径大于开口的1/3时，易形成砂桥，发生堵塞现象[28]。开口高度再次减小或颗粒数量增加时，该现象发生的概率将增大。在现场应用中，若部分筛管发生堵塞，由于生产井产量一定，将增加未堵塞筛管的过流速度，增大冲蚀速率。

4.3 寿命预测分析

4.3.1 预测模型计算

由图3可知，最大冲蚀速率出现在冲缝管套筛缝底部边缘区域。分析计算入口流速为0.10 m/s、颗粒体积分数为5%、液相体积分数为1%和开口高度0.30 mm条件下筛管结构的冲蚀速率，得出筛缝底部的冲蚀速率云图如图9所示。从图9可知，底部冲蚀区域面积占比约为2/3。

图9 筛缝底部冲蚀云图Fig 9 Cloud map of bottom of screen seam

拟用最大冲蚀速率ER表示底部区域的平均冲蚀速率E(kg·m-2·s-1)，但由于冲蚀速率会随时间而变化，故需引入修正系数k。从图9可看出：部分冲蚀区域冲蚀速率与最大冲蚀速率相差较大，冲蚀区域的冲蚀速率的最小值与最大值的比值为0.1，综合考虑后取k=0.1。

冲蚀区域的平均冲蚀速率E为

颗粒冲蚀筛管将造成筛缝开口高度增大，导致筛管过滤精度发生改变，防砂失效。冲蚀深度按壁面材料的质量损失进行计算。该冲蚀模型的冲蚀主要发生在开口底部一侧，计算冲缝管套的质量损失MR为

式中：S为筛缝处冲蚀的面积，m2。

根据体积公式，将质量损失转化为体积损失后，可计算冲蚀深度为

式中：L为冲蚀深度，mm；VR为体积损失，m3；ρb为管套材料的密度，kg/m3。

冲缝管套预测使用寿命T为

式中：D为筛缝允许冲蚀深度，mm。

管套常用材料304 不锈钢密度为7 930 kg/m3。根据GILLESPIE 等[27]进行的冲蚀磨损试验，缝宽增加0.05 mm时，出砂量会严重增多。考虑到冲蚀发生在底部，故D=0.05 mm。

4.3.2 寿命预测结果

根据式(25)，当允许冲蚀深度为0.05 mm 时，预测使用寿命与冲蚀速率呈反比例关系，分析计算后可得到管套寿命预测数据。在天然气产出的初始阶段，流速大于初始条件值。按最大日产量3.5×104m3计算，流速取0.2 m/s。将此条件下的数据绘制为冲缝管套预测寿命折线图如图10所示。

图10 不同因素下冲缝管套预测寿命折线图Fig.10 Line diagram of predicted life of punched screen for different factors

从图10可以得出：冲缝管套的预测使用寿命随流体速度增加而减小，在速度增加的初始段，预测使用寿命下降趋势明显，该情况与实际符合。预测使用寿命随着开口高度增加而增大，在开口高度为0.15～0.30 mm 间减少趋势为线性关系。开口高度为0.15 mm和0.10 mm时有所偏差的原因是堵塞现象，导致冲蚀速率骤增、预测寿命减小幅度大。预测使用寿命随流体中液体体积分数和颗粒体积分数上升而减小。在液相体积分数为0.5%时，由于底部区域残余的液相较少，无法形成液膜，导致冲蚀速率大、预测寿命小。冲蚀速率与颗粒体积分数之间的关系为线性关系，由式(25)和图10(c)可知，预测使用寿命与颗粒体积分数之间的关系为反比关系。

5 结论

1)利用CFD 方法分析了冲缝管套筛缝的流场分布情况，管套筛缝受到的冲蚀磨损主要集中在底部边缘区域。

2)在水合物储层开采中，流体流速、液相体积分数增加和开口高度减小都是通过直接或间接方式增大颗粒动能，进而增大冲蚀速率，而颗粒体积分数增加通过增加颗粒碰撞壁面的数量来增大冲蚀速率。

3)流体中液相体积分数增加对冲蚀有着增大和抑制2 种作用。当液相体积分数增大时，一方面，会增加颗粒动能，增大冲蚀速率；另一方面，在冲蚀区域形成液膜，削弱颗粒动能、减少摩擦磨损。

4)冲缝管套预测使用寿命在流速增大的初始阶段下降较为明显。当液相体积分数增大时，形成的液膜将起到一定的减小冲蚀、延长管套使用寿命的作用。