临界CO2传输中水合物对弯管冲蚀规律分析

袁显宝， 刘 超， 谭 伟， 张永红*， 张彬航， 程子阳

(1.三峡大学机械与动力学院， 宜昌 443002； 2.湖北省水电机械设备设计与维护重点实验室，宜昌 443002)

随着温室效应、气候变化等环保问题的日益突出，碳排放的有效控制是当今社会在发展中需要考虑的重要问题。利用碳捕获、利用与封存(carbon capture，utilization and storage,CCUS)技术将CO2资源化的技术受到各国重视[1-5]。成熟的CCUS技术可将含有杂质的CO2气体经过分离、净化和压缩处理为超临界CO2，并通过管道运输。当输送流质中混杂有游离水时，在温压产生波动条件下易形成CO2固态水合物，其物理性质与冰类似。在输送过程中， CO2水合物会与管壁发生碰撞将对管壁造成冲蚀破坏，影响管道传输安全性；冲蚀坑为水合物提供生长附着点，水合物沿管壁堆积将影响输送效率，甚至可能堵塞整个管道[6-7]。弯管是维系管道输送系统结构完整的必要结构，而据大量研究表明，含有固体颗粒的两相流对于弯管的侵蚀危害程度相对于直管段更加严重[8-11]。

目前，根据应用方向的差异，对管道壁面冲蚀研究主要从混杂有固体颗粒物的水、油、气三种材料进行开展，且相关研究成果较多。但对于超临界CO2输送管道内水合物对管道冲蚀的研究较少。王博等[12]研究了气固混合流集输管道的弯管段受颗粒冲蚀影响的因素；刘丹[13]利用实验方法研究了超临界CO2流动状态对碳钢表面的腐蚀效应的影响；李乐等[14]利用实验方法研究了在流量变化对高压环境下CO2水合物生成过程的影响；巨熔冰[15]利用FLUENT软件研究了超临界CO2流体产生的应力和湍流动能对管壁的影响；刘洪亮等[16]使用FLUENT软件求解了压裂工况下T形三通冲蚀影响因素。超临界CO2输送管道中水合物对管道的冲蚀效应未见更详细的研究报道。因此，根据超临界CO2输送管线的输送工艺和固液两相物质的属性，采用数值模拟方法，选取输送管路中的弯管段的冲蚀情况进行计算，对于输送管道的防护和安全性评估有重要参考价值，有利于CCUS技术的发展。现基于COMSOL多物理场耦合软件研究超临界CO2输送弯管的冲蚀效应，求解不同流速、颗粒直径和不同弯角对弯管冲蚀效应的影响。

1 理论模型及控制方程

1.1 颗粒离散相模型

使用COMSOL Multiphysics软件的流体流动颗粒追踪模块与湍流模块耦合求解流场中微粒的运动轨迹。通过对每个分量建立相应的常微分方程，计算每个时间步内，粒子受连续相的作用力与粒子间相互作用力的合力，求解粒子运动状态。离散相颗粒的受力与惯性力保持平衡[17]，流体对粒子的曳力FD占主导作用，粒子的受力方程可表示为

(1)

(2)

(3)

(4)

式中:Vp为离散相粒子速度，m/s；dp为粒子直径，m；ρp为粒子密度，kg/m3;Rep为粒子雷诺数；CD为阻力系数；FD为流体拖曳力；Vf、ρf分别为流体速度、密度大小，m/s1、kg/m3；μf为流体的动力黏度，Pa·s；a1、a2、a3为经验常数，与粒子雷诺数和颗粒形状相关。

除受到曳力外，粒子在相载流体中受到浮力FB、压力梯度FP、附加质量力FA、旋转力FR，各力的表达式为

(5)

(6)

(7)

(8)

式中:∇p为粒子容积范围内压强梯度，Pa；ωp为粒子的旋转角速度，rad/s。

引入Stokes数St来描述粒子运动对流场变化的响应特点，其定义为粒子相应时间与流体传播时间的比率[18]。其表达式为

(9)

式(9)中:D为外部容器的管径，m。

1.2 湍流模型及壁反射函数

参考Menter[19]对于k-ω湍流模型应用于强流线曲率流动准确性的研究工作，采取标准k-ω湍流模型来求解弯管冲蚀过程中的湍流流场。其控制方程为

(10)

(11)

式中:k为湍动动能；ω为湍流动能比耗散率；u为流体微团的速度；i、j分别代表X、Y方向；Γk、Γω为和k、ω的扩散率；Gk为由速度梯度引起的湍流动能生成项；Gω为由比耗散率的生成项；Yk、Yω分别为湍动能、比耗散率的发散项；Sk和Sω分别为湍动能用户自定义量、比耗散率用户定义量。

颗粒随流体运动过程中，与管壁的碰撞伴随着由摩擦和非弹性效应产生的能量转移和损失，主要表现在入射粒子和反射粒子在速度分量的改变，通常以碰撞前后速度分量之比定义恢复系数。现采取Forder等[20]提出的粒子与壁的碰撞回弹模型来描述恢复系数与冲击角度之间的关系。

(12)

式(12)中:en为法向恢复系数；et为切向恢复系数；vn为法向粒子速度分量；vt为切向粒子速度分量；θ为粒子碰撞角度(入射速度与表面切线之间的角度)；下标1、2分别为碰撞前、后的状态。

1.3 粒子冲蚀模型

冲蚀效应的严重程度主要通过冲蚀磨损率表示，其定义为运动的颗粒物的冲击造成的靶材上单位时间单位面积上的质量损失。选取COMSOL中的E/CRC冲蚀模型来模拟弯管冲蚀过程中的冲蚀效应。该模型由E/CRC冲蚀与腐蚀联合研究中心(Erosion/ Corrosion Research Center)和Tulsa大学提出[21]，冲蚀磨损率EM的计算公式为

(13)

式(13)中:K为侵蚀模型系数，缺省值为2.17×10-9；BH为表面材料的布氏硬度，缺省值为200；V为入射粒子速度，m/s；Vref为入射粒子的参考速度，缺省值为1，m/s；FS为颗粒形状系数，缺省值为0.2，代表圆形颗粒;n为模型指数，缺省值为2.41；α为粒子入射角，rad。

2 模拟计算基础

2.1 建立模型

基于参考文献[15]所提出的超临界CO2输送管道的设计方案和运行参数来建立基准计算模型。建立如图1所示的几何模型，其结构由短管段-弯管段-长管段三部分组成，流体由短管段进入，长管段流出。管道直径D=200 mm，管壁厚度δ=10 mm，曲率半径R=200 mm，短管段的长度L1=400 mm，为了使弯管后段的湍流充分发展，出口直管段长度设置直管段长度L2=800 mm。选取研究区域弯管体积的1/2作为建模区域，减少所需的计算内存资源和求解时间。

图1 弯管模型的几何结构

当超临界CO2输送管道外界温度在10～25 ℃，水体积含量大于21 000 ppm并且压力大于约8 MPa时，CO2水合物易生成且对弯管冲蚀较为严重[15]。管内两相流中液相为超临界CO2，在40 ℃时密度为1 100 kg/m3，动力黏度为3×10-5Pa·s；离散相为CO2水合物颗粒，密度为900 kg/m3，粒径分布在5×10-6～200×10-6m。管道的管壁材料选取实际应用于CO2超临界与密相传输的管材X80管线钢，其材料属性数据如表1所示。

表1 材料物性参数

2.2 边界条件

参考李井洋[22]对大口径输送管道中粒子冲蚀因素的研究方法，使用COMSOL计算粒子冲蚀问题分两步进行，使用湍流模块求解连续相的流场分布，在此基础上使用流体流动颗粒跟踪模块求解颗粒运动和冲蚀效应。

湍流模块中，以图1中的面5、面7为速度入口边界，设置为充分发展的流动，平均流速为2 m/s；面42为压力10 MPa的出口边界；施加方向为Z轴负方向的重力场；设置YZ剖面为对称边界。

流体流动粒子跟踪模块中，为粒子指定密度和直径，其密度为900 kg/m3，平均粒径为80 μm，粒子体积浓度设为2.1%；流体与管壁接触面设定壁反弹边界条件，同时增设E/CRC冲蚀边界；设定入口面为面7，速度场继承自湍流流场计算的结果；出口边界为面42，壁条件设为消失；流域内的粒子受到斯托克斯曳力作用，湍流弥散模型设为离散随机游走；设置YZ剖面为对称边界。

求解器分步进行研究，使用稳态研究求解湍流流场及对管壁的应力，使用瞬态求解计算2 s内粒子跟随流体流动状态，时间步长为10-5s，容差均设为10-3以保证精度；采用内置的MUMPS求解器进行计算。

2.3 网格划分

采取非结构化单元进行网格绘制，网格模型整体概况如图2所示。边界层处细化为5层以求解边界处的湍流流动状况，生长因子为1.2。以求解得到流场的最大流速为参考标准，设置多组密度的网格进行计算进行网格无关性验证。计算结果如图3所示，网格单元数在到达18×104时求解结果基本稳定，最终的网格划分方案网格单元数为178 587，最小质量0.148，平均质量0.713。

图2 弯管模型的网格模型

图3 网格数量与对应的最大流速

3 结果分析

3.1 流场和冲蚀率分析

平均入口流动速度2 m/s，粒径为80 μm，弯曲角度为90°的弯管，中心对称面上的超临界CO2流体速度和压力分布分别如图4(a)、图4(b)所示。

由图4(a)所示的速度分布云图可知，上游短管段的湍流流场充分发展，管道中心区域的流速高于近壁区的流速，平均流速约为2 m/s；进入弯管区域后，弯管内侧流体获得较大流速，约为3 m/s，外侧区域流速约为1.5 m/s；流体由弯管段向直管段过渡时，发生了一定程度的流动分离，靠近内侧管壁附近的流动速度较小而靠近外侧区域流速较大。在下游长管段中，高流速区域贴近外侧管壁且截面逐渐缩小，管内整体流速降低到2 m/s。根据流场的分布可知，当水合物粒子随流体流动时，受曳力和粒子自身向心力作用，在长管段外侧的壁面容易遭受一次粒子的冲蚀破坏。

由图4(b)所示的压力分布云图可知，弯管附近存在压力变化，弯管内侧由于流体的流速较快，压力低于外侧区域，压差约为4 000 Pa；超临界CO2流经弯整个管段过后产生了约8 000 Pa的压降。

图4 对称面上超临界CO2流场分布

图5为直角弯管在前面所设条件下壁面受粒子冲蚀影响云图，由图5中可知，受到CO2水合物冲蚀较为严重的区域集中在弯管段外侧管壁与长管段交接处，对应图4中沿直管壁的高流速区域的起点位置，管道壁面所受到的最大冲蚀速率为1.37×10-9kg/m2/s。

图5 基准条件下弯管表面冲蚀云图

3.2 流速大小对冲蚀效应影响

为研究入口平均流速大小对弯管冲蚀磨损速度的影响，设置入口平均流速从1 m/s变化到5 m/s，速度每增加0.5 m/s进行一次计算，不改变其他条件。弯管壁面的最大冲蚀率与入口平均流速关系如图6所示。由图6可知，对于直角弯管模型，当入口平均流速增大时，弯管最大冲蚀率随之增大，基本呈现出二次函数增长的关系。

图6 最大冲蚀速率随入口速度变化曲线

另外，平均入口速度发生改变时，管壁受到冲蚀影响区域的位置也将改变。图7为不同的入口平均流速条件下弯管受到冲蚀影响区域的对比。从图7中可知，流速增大时，冲蚀效应增强，冲蚀中心更加集中并向弯管拐角处偏移。造成这种现象的原因是，水合物粒子在弯管中随流体运动时，除了受流体曳力影响外还受到自身惯性力作用。当流体的流动速度增加时粒子速度随之增加，此时惯性力对粒子运动轨迹的影响增强，流体流动方向突然改变时，粒子容易脱离流体运动方向，与管壁发生碰撞，导致冲蚀效应增大且冲蚀中心位置向下方移动。

图7 冲蚀区域受入口速度大小变化的影响

3.3 颗粒直径对冲蚀率影响

为研究颗粒直径大小对弯管冲蚀磨损速度的影响，设置CO2水合物的颗粒直径从20 μm变化到250 μm。颗粒的直径与弯管壁受到的最大冲蚀率之间关系曲线如图8所示。由图8可知，随着粒子粒径增加，最大冲蚀率随之增加；当微粒的直径在20～80 μm变化时，最大冲蚀速率受粒径变化的影响较小，基本呈二次函数关系，而颗粒直径在80～250 μm变化时，粒子直径的增长对冲蚀率的影响较小，呈线性增长关系。图9为不同粒子粒径对管壁造成冲蚀区域的对比。由图9可知，粒径逐渐增大时，冲蚀区域的位置和形状发生改变，由直管段外侧壁面的弥散分布转变为弯管外壁的集中椭圆形分布。这是由于，由于粒径较小时，Stokes数值小，粒子跟随性强，与管壁发生接触时位置通常在直管段，冲击角度与管壁夹角较小，对管壁冲蚀主要以切削作用为主。粒子与管壁碰撞时的动能与粒子速度和粒子质量的增加均相关，因此冲蚀率的增长曲线接近二次函数，但数值变化较小；当粒子直径增大到一定程度时，Stokes数值较大，粒子跟随性较弱，在弯管段受到向心力加速的程度弱，粒子容易脱离流体与管壁发生碰撞，弯管段的管壁发生正碰，因此冲蚀程度较为严重。此时粒子速度接近初始速度，粒子与管壁碰撞时的动能只受到粒子质量变化的影响，冲蚀增长率与粒径呈一次函数关系。

图8 最大冲蚀速率随粒径变化曲线

图9 不同粒径条件下的冲蚀区域分布

3.4 弯管角度对冲蚀率影响

为研究弯管角度(长短直管的夹角)对弯管受到冲蚀磨损速度的影响，设置弯管角度分别为30°、45°、60°、90°、120°、135°、150°。管壁受到的最大冲蚀率与弯管角度的关系如图10所示。由图10可知，管壁受最大冲蚀率在弯管角度小于90°时呈上升趋势，当弯管角度大于90°时呈下降趋势；角度为90°时，弯管受到的冲蚀效应最大。角度在90°前，最大冲蚀率将随角度增大而平缓增加，角度超过90°以后，冲蚀效应迅速减弱。

图10 不同弯管角度下的最大冲蚀率变化曲线

图11为不同角度弯管受到粒子冲蚀区域位置对比。从图11中可知，当弯管角度小于90°时，粒子对管壁造成的一次冲击区域较集中；当弯管角度大于90°时，冲蚀影响区域较分散，且随着夹角增大，冲蚀区域分散程度增加。弯管角度较小时，管内流体流动方向发生剧烈偏转，粒子受到流体曳力和向心力不足以抵消粒子自身的惯性力，粒子与弯管的冲蚀区域与直角弯管的冲蚀区域几乎相同；而流体在夹角大于90°的弯管流道中流动时，方向偏转小，粒子在与管壁碰撞前将随流体流动得更远；与管壁发生碰撞时，主要以切削为主，无法对管壁造成较大损害。

图11 不同弯管角度条件下冲蚀区域分布

4 结论

利用COMSOL湍流流动模块和颗粒跟踪模块，耦合求解超临界CO2输送管道受到CO2固体水合物冲蚀效应，得出以下结论。

(1)随着管道入口平均流速的增加，水合物颗粒的速度随之增加，CO2固体水合物对弯管的冲蚀程度将增大，冲蚀区域由弯管段外侧较远处的弥散分布转变为弯管外侧拐角处集中分布。采取较小的输送速度能够有效降低水合物对管壁的冲蚀。

(2)随着CO2水合物颗粒的粒径增大，管壁上受到的最大冲蚀率不断增大，粒径较小时，冲蚀不明显；粒径超过约80 μm时，冲蚀速率的增长与粒径呈一次线性关系。颗粒粒径增加将使冲蚀区域由远离拐角的直管外侧弥散分布向拐角靠拢聚集。

(3)固体水合物颗粒随超临界CO2流体流经不同角度的弯管时，对管壁的冲蚀程度不同，直角弯管更容易受到粒子冲蚀。当弯管的角度大于90°时，管壁受粒子冲蚀的区域相对分散，且角度增大，冲蚀影响区域会扩散但冲蚀程度会迅速减弱；当弯管角度小于90°时，粒子对管壁的冲蚀效应会随角度的减小而缓慢降低，但是冲蚀发生的主要区域位置不会发生较大改变。