天然气水合物垂直管内流动特性的数值模拟

刘宝玉, 侯朋朋, 商丽艳, 李 萍, 吕振波

(1. 辽宁石油化工大学 石油天然气工程学院, 辽宁 抚顺 113001; 2. 辽宁石油化工大学 化学化工与环境学部, 辽宁 抚顺 113001)

天然气水合物是水分子和甲烷等气体分子在低温(<10 ℃)、高压(3～5 MPa)、气体体积分数大于其溶液条件下形成的一种具有笼状结构的化合物,是一种清洁无污染的新型能源[1].深水天然气水合物主要以大块天然气水合物混合少量杂质和颗粒状分散于细颗粒的沉积物中两种形式存在,可以利用集矿机在海底对天然气水合物进行收集,并进行泥沙分离,将固态水合物提升至海平面[2].深水中天然气水合物的固态开采法结合了水力提升金属锰结核的技术及方法,水力输送系统的垂直提升管是核心设备之一,考虑到整个系统中的可靠性、安全性以及压力损失等方面,需对其进行研究.

在固体开采法中,天然气水合物是提升到海面分解收集平台后再进行相平衡的破坏.在此提升过程中可始终当成无相变的固液两相管道水力输送,关于管内流动特性分析,不少学者都做出了显著的贡献[3-6].文献[4-5]对海底水合物开采系统中用湍流模型对开采管内固液和固液气的流动介质进行了模拟分析,得到管径、速度、流量等参数对流场的影响规律.文献[6]对R11水合物浆液在复杂管内从压降、体积分数、流速、湍流耗散率等方面的流动特性进行了试验和模拟分析.目前,分析相关参数对流动特性的影响运用Fluent数值模拟方法已成为试验方法的有效代替手段.

与一般的固液两相流相比,水合物浆液数值模拟需要考虑黏度模型的选择,以管道中浆液体积分数变化,确定黏度值.文中采用传统的CFD/Fluent模拟技术,通过编制UDF对水合物浆液黏度进行自定义,着重研究水合物颗粒相的管输流变和管流特性,通过研究水合物颗粒粒径、体积分数、速度的变化对水合物浆液流动影响以及管道运输的水力坡度得到水合物浆管道安全运行最优工况.

1 输送参数选择

定义水合物浆液参数如表1所示.

表1 两相初始参数

改变流速v,分析浆液在垂直管道内速度变化规律(工况1-4),然后分别改变颗粒粒径ds、颗粒密度ρs、颗粒体积分数Cv等参数(工况5-12),见表2.分析计算不同工况条件下管道的体积分数分布规律和压降,对其管道流动特性进行研究.

表2 各工况的参数选择

2 模型和计算方法

2.1 垂直管几何模型

计算模型长度为4 m,管径为300 mm的垂直管,利用ICEM对模型计算域进行网格划分,采用六面体网格,横截面上网格划分采用“古钱币”形式将垂直管划分为128 250个单元的三维结构网格(进出口端面为15×15,周边4×10×15,轴向为150个单元),如图1所示.

图1 模型三维结构网格

2.2 网格无关性检验

对于圆管内湍流流动的模拟,需确定网格与计算结果之间的无关性.本工作主要研究浆液流动的变化,流动的变化与速度密切相关.初始参数计算收敛结束时,选取管道中心线上3点位置水合物浆液的速度分布作网格无关性检验.3点坐标分别为A(0, 1, 0),B(0, 2, 0),C(0, 4, 0).不同网格划分的网格数与各点速度分别如表3和图2所示.由图2可知,随网格数增加,1-3划分方式各点速度逐渐增加,第4种网格划分方式各点速度比第3种网格划分方式增加不明显.由于第 4 种网格划分方式过密,使模拟消耗时间过长,加大计算成本,因此取第3种网格划分方式作为计算网格.

表3 不同网格划分的网格数

图2 不同网格划分下的各点速度

2.3 基本假设

为了数值模拟的可行性与结果可靠性,进行以下假设: ① 流动模拟不受边界层影响; ② 固液两相均为不可压缩、连续的流体,其他物性参数均为常数,主次相分别为水和天然气水合物颗粒; ③ 天然气水合物颗粒为均一的球形颗粒,且不考虑相变以及相间的传热.

2.4 控制方程

连续性方程:

(1)

动量方程:

(2)

式中:ρ为液相密度,kg·m-3;t为时间,s;ui和uj为速度分量,m·s-1;Xi和Xj分别为x和y坐标,m;p为压力,Pa;μ为动力黏度,Pa·s;τij=-ρuiuj为Reynolds应力.

2.5 湍流模型

采用FLUENT软件中的标准k-ε模型,其中湍动能k和湍动耗散率ε是两个基本未知量,与之相对应的输送方程为

(3)

(4)

式中:k为湍动能；ε为湍动耗散率;Gk是平均速度变化引起的湍动能k的产生项,模型常数C1ε,C2ε,Cμ,σk和σε的取值为1.44,1.92,0.09,1.0,1.3;μt为湍动黏度,可表示成k和ε的函数:

μt=ρCμk2/ε.

(5)

2.6 混合相黏度模型

两相流双流体模型是把颗粒相与液相作为相互渗透的拟连续介质的流动,为了求解动量方程组,即要获得液体相的黏度物性参数,还要建立固体颗粒拟流体的黏度计算模型.因此不再采用常数作为黏度值,计算公式[7]:

μs=μl(1+2.5C+10.05C2+0.002 73e16.6C),

(6)

变换后可得到固相的黏度计算公式:

(7)

式中:μs为浆液黏度,Pa·s;μl为水的黏度,Pa·s;Cs为浆液的阻力系数;Cl为水的阻力系数.在此基础上,编制UDF,实现流体黏度随着浆液体积分数的变化而变化.

2.7 边界条件

运用Fluent软件建立三维欧拉模型,定义进口边界为速度入口,速度为浆体输送速度,因为湍动能和湍动耗散率难以估计,用容易通过计算来获得的湍流强度和水力直径代替,湍流强度取4%,水力直径为0.3 m.定义出口边界条件为充分发展的边界条件outflow.在壁面处,液相设置为无滑移边界条件和采用壁面函数法,颗粒相设置有滑移条件,即与壁面碰撞时的弹性碰撞反射.颗粒间碰撞归还系数为0.95,并加以重力加速度等环境因素.

3 计算算例及结果分析

3.1 模型验证

颗粒的体积分数是体现浆体在管道水力提升中颗粒分布的主要特性，由于颗粒在管道内分布的对称性，将数值模拟得到的沿半径方向颗粒体积分数Cv除以管道中心的颗粒体积分数值Cvc与试验值[8]做对比，如图3所示，当体积分数分别为小于5%，5%～20%和40%～45%时，试验值与模拟值最大误差分别为18%，14%和19%，模拟数值与试验数值的相对误差均控制在20%以内.可见，该模型较好地反映了天然气水合物颗粒在管内流动的浓度特性.

图3 径向相对位置颗粒体积分数试验值与模拟值对比

3.2 速度分布

垂直管道模拟的基本参数是ds=0.02 mm,Cv=20%,ρs=1 190 kg·m-3,在不同流速下,浆液在垂直管道内的速度分布云图,如图4所示.

图4 不同入口速度,z=0轴向浆液速度和出口径向速度分布

从图4中可知,水合物浆液在垂直管道中的流速梯度分布随着入口流速的增大,变化较为明显.当入口速度为1.02 m·s-1时,流速不足以使管道底部水合物颗粒全部运动,管道底部出现不动或滑动的水合物颗粒淤积层[9],整个浆体在出口处的流动并未出现明显的梯度递减,易造成管路的堵塞;而当流速增大时,底部不动的固体颗粒开始运动,此时管道轴向的速度梯度沿管径向外呈明显的梯度递减,壁面处的流速最小,可使浆体平缓流出管道.

固体颗粒在垂直管道中的轴向速度沿径向分布可分为3种:第1种是抛物线型的分布;第2种是受限塞状流行分布,颗粒向管道中部集中;第3种是混合型分布,该管道中部速度分布类似一型,管壁附近类似二型.通过模拟判断出水合物浆液在垂直管道水力提升时,速度分布呈现混合型,由图4d竖直管速度云图知,颗粒平均流速约等于浆液最大流速的0.85倍,在中心附近,速度分布更为平缓,形成以恒定速度运动的颗粒核.

3.3 体积分数分布

水合物浆液体积分数能够较好的反映水合物颗粒在垂直管中的流动情况,以水合物颗粒的体积分数来表征浆液的体积分数,以初始工况的中轴线上的流速与体积分数的变化曲线来分析管道内浆液流动变化,如图5所示,从图中可知,当水合物浆液以1.43 m·s-1的流速开始流动时,浆液体积分数从20%快速升至22.53%;流速从1.45 m·s-1到1.47 m·s-1时,浆液体积分数保持基本不变;流速从1.47 m·s-1提升到1.49 m·s-1时,浆液体积分数下降至20%;随着浆液的流速继续增大,浆液体积分数则保持在20%不变.

图5 体积分数分布与流速的关系曲线

图6中分别表示不同浆液初始速度、颗粒粒径、颗粒密度和颗粒体积分数在垂直管中中轴线上的浆液体积分数和动力黏度的变化情况.基于Thomas动力黏度方程,对模拟水合物浆液流动时的动力黏度进行了UDF编制,采用文献[10]提出的适用于固液两相流动的双层流动模型,模拟得到的水合物浆液黏度,再用Thomas方程计算浆液黏度,由于模拟和计算方法都采用Thomas黏度方程求解,因此两者在一定范围内结果相近.从图6可得水合物浆液的动力黏度在管道输送初始阶段都有明显的上升,在输送至一定管道距离后,动力黏度达到最大值,随着管输距离的增大,水合物浆液的动力黏度能够达到平稳值.其波动的原因：一是浆液输送过程中颗粒相间的碰撞造成颗粒分布不均匀;二是颗粒的曳力作用对水合物浆液的体积分数影响较大,与液相有一定的速度差.较高的黏度易造成管路中颗粒聚集或管道堵塞的情况发生.

图6a中,改变入口速度的大小,得到的体积分数与黏度值,在工况1,4中,水合物浆液输送至入口1.5 m处,体积分数已经达到最大值,而工况3的最大值发生在管道后半段,此工况下水合物颗粒在管道初始阶段是可以相对稳定输送,由于滑脱作用,颗粒相的速度远远滞后于液相速度,在管道3 m处颗粒聚集较明显,使得管路出现堵塞的不安全流动,出于效率考虑不建议速度较低的输送速度.图6b中,改变水合物颗粒粒径大小,工况5,6中管道会发生多次聚集现象,由于粒径的增大,颗粒与管道壁面产生的摩阻损失也不断增大,颗粒相流动受阻,固液两相间速度差不能快速趋于稳定;细小颗粒在管道中的体积分数基本保持不变,颗粒粒径越小,其剪切稀释性变化较明显,使得浆液能够维持在初始体积分数.当颗粒相出现多次聚集时,管道中浆液均质性变差,不能够满足高效率提升水合物浆液.

图6 各工况下对体积分数与动力黏度影响

图6c中,改变水合物浆液颗粒的密度,对管道的稳定流动影响较大,在工况9中,管道入口1 m处颗粒相会出现频繁的聚集,此时在较小的速度下,难以让水合物浆液输送出管道;颗粒密度越大,颗粒容易形成胶团,聚集在管道输送初始阶段,此时浆液的动力黏度增大,造成管路堵塞等不安全流动.图6d中,浆液初始体积分数的变化对管道中体积分数最大值变化影响较明显,在工况10中,管道最大体积分数为11.57%,而工况12最大体积分数和黏度分为45.02%,0.007 5 Pa·s,黏度较高使得管道的压力损失较大,且高体积分数的水合物浆液管道中,在管道出口处才能达到初始体积分数,极易在出口前发生管路堵塞.

3.4 管道压降

管道压降用水力坡度来表示,选取入口和出口2个断面的压力进行计算,计算公式如下:

(8)

式中:I为浆体的阻力损失,Pa·m-1;Δp为出入口的压差,Pa.

根据式(8)计算各种工况条件下的水力坡降,得出各个因素对管道压降的影响情况,如表4所示.

表4 水合物浆液流动参数模拟结果与计算结果对比

对比3组水力坡度[11],说明数值模拟可以在一定程度上体现压降变化规律.由表中可以看出,垂直管道中压降随流速、粒径、颗粒密度、浆液体积分数成正比关系.当粒径在0.01～0.02 m时,水力坡度变化缓慢,当颗粒粒径大于0.03 m时,水力坡降变化较大;颗粒密度越大表明水合物含量较小,所以建议颗粒密度在1 020～1 190 kg·m-3;随着体积体积分数的增大,管道水力坡度先增大后减小,主要因为体积分数变化导致浆液密度变化,有利于颗粒的运动,促使管道运输能力增强,建议颗粒体积分数选择20%上下;为了使垂直管更有效率的进行水力提升,在水力坡度影响不大的情况下,建议速度控制在1.4～1.8 m·s-1.在垂直管输送水合物浆液时,由于实际工况影响因素较多,与理论值有一定偏差,为减小压降,以及考虑经济因素,降低设备能耗,可选取以上工况进行弹性调节.

4 结 论

1) 水合物浆液在垂直管道中的速度沿轴向不断增大,管道流速从中心沿径向呈梯度递减,径向流速趋于圆形分布,固体颗粒完全悬浮,且颗粒平均流速约为浆液最大流速的0.85倍.

2) 水合物浆液的体积分数变化与动力黏度趋势相似,其中水合物颗粒粒径与密度对浆液体积分数变化较大;初始输送体积分数对浆液动力黏度影响最大,在速度为1.4～1.8 m·s-1、颗粒粒径为0.01～0.02 m、颗粒密度为1 020～1 190 kg·m-3、初始体积分数控制在20%时管路能够稳定流动,不易发生颗粒聚集和管路堵塞.

3) 结合浆液中固相的体积分数变化的特点,在开采水合物垂直管中,对浆液编制UDF,使颗粒相的黏度随着体积分数的变化而变化,模拟试验的结果较为相似.对管道中不同工况下的水力坡度进行计算,基于能量理论建立的公式则比较全面的衡量了各个因素的影响,可对开采管路的压降进行简单的计算与流动特性预测.