垂直管内水合物浆液流动特性的CFD-PBM数值模拟

2019-01-17侯朋朋王春华商丽艳韦雪蕾刘志明

石油化工高等学校学报 2019年1期

侯朋朋，王春华，潘振，商丽艳，韦雪蕾，刘志明

（辽宁石油化工大学石油天然气工程学院，辽宁抚顺113001）

在中国天然气水合物市场深度分析及“十三五”发展战略报告中重点指出，要加快天然气水合物试开采进度，推动产业化进展。目前，国际上对天然气水合物的勘探研究一直处于探索阶段，近期我国已在南海“神狐海域”试开采成功，但现在对于天然气水合物开发基础和技术方法的研究仍需进一步探索[1]。

在“地层流体抽取”试采法中，将天然气水合物从海底输送至海面的水力输送垂直硬管是核心设备之一。近年来，国内外许多学者对天然气水合物浆液的多相流动进行了数值模拟[2-6]和实验研究[6-11]。宋光春等[12]利用水合物颗粒数量密度的群体平衡模型（PBM），分析了该模型下的水合物颗粒变化，但未考虑气相的直径变化。丁麟等[13]采用双流体模型和欧拉方法对气体水合物浆液系统的多相流流动进行了模拟，结果表明，水合物的形成与压降和流体的相变有关，然而并没有对相变后的流动进行分析。Xu H L等[14]以海底天然气水合物的固态开采法系统中的开采管为研究对象，对管道内水合物浆液的固液两相流动进行了CFD模拟，得到管径、流量、颗粒粒径等参数对开采系统的影响规律，由于是基于固液两相不发生相变的情况下进行的模拟，因此对水合物浆液的流动还需考虑气液固三相流动。

目前，国内外对海底天然气水合物管道输送的相关研究文献较少。由于海底天然气水合物开采输送过程中不断有气体分解出，因此本文以水合物浆液的气相为主导系统，流动模型中代入PBM方程，考虑了水合物分解后气相的聚集与破碎对管道内流动特性的影响，并对浆液的湍动能耗散率、气相的初始浓度与速度影响气泡分布情况进行了分析，从而对深海水合物开采管道的浆液流变特性提供了一定的基础数据和技术支持。

1 数值模拟模型的建立

1.1 几何模型及网格划分

计算模型为长度3 m，管径300 mm的垂直管，利用ICEM对模型计算域进行网格划分，采用六面体网格，横截面上网格划分采用“古钱币”形式将垂直管划分为128 250个单元的三维结构网格（见图 1）。

图1 模型三维结构网格Fig.1 3D mesh of model

1.2 边界条件

运用Fluent软件建立PBM模型，定义进口边界为速度入口，速度为浆体输送速度，因为湍动能和湍动耗散率难以估计，用计算获得的湍流强度和水力直径代替，湍流强度取4%，水力直径为0.3 m。定义出口边界条件为充分发展的边界条件。固液两相参数采取文献[15]中工况参数，气相用PBM模型对直径为0.5～2.0 mm的气泡分五组进行模拟。在壁面处，流体相采用壁面函数法和无滑移边界条件，颗粒相在壁面设置有滑移条件。颗粒间碰撞归还系数为0.95，并加以重力加速度等环境因素。

1.3 数值计算模型

1.3.1 控制方程将垂直管内气液固三相流动过程简化为恒温体系，即不考虑能量方程，同时采用Euler多相流模型模拟垂直管内气液固三相间的相互作用，水合物浆液的三相看作充满计算区域的连续介质进行处理，通过连续方程和动量方程求解每一相，得到控制方程分别为连续性方程式(1)和动量方程式(2):

式中，ρ 为液相密度，kg/m3；t为时间，s；ui、uj为速度分量，m/s；xi、xj分别为 x和 y坐标，m；p为压力，Pa；μ为动力黏度，Pa·s；τij=-ρ为 Reynolds应力，、为对应的脉动量。式（1）、（2）实质上是质量守恒方程和动量守恒方程在流体力学中的具体体现。1.3.2 湍流模型流体在垂直管道水利提升时速度变化较大，同时产生较大的湍动能以及湍动耗散率，本文采用realizableεk-ε湍流模型。运输方程为式（3）：

式中，k为湍动能；ε为湍动耗散率；Gk为由于平均速度梯度引起的湍动能k的产生项；模型常数C1ε、C2ε、Cμ、σk和 σε的取值为 1.44、1.92、0.09、1.00、1.30；μt为湍动黏度，可表示成 k和 ε的函数关系式（4）:

1.3.3 相群平衡模型气泡直径大小的变化与聚并及破碎过程有关，需要添加一个平衡方程来描述气泡在体系中的平衡，该方程称为相群平衡模型[16]。研究中气泡间会发生连续不断的反应，并对气泡大小进行分组，这样可对浆体中气泡流动进行模拟研究，从而对水合物浆液气液固多相流动机理进行更直观准确的模拟研究，气泡群体平衡方程可用式(5)表示：

式中，n(v,t)为气泡分布函数；a(v,v’)为气泡聚并函数；b(v)为气泡破碎函数。

群体平衡模型的求解方法有离散法、标准动量法和动量积分法。本文采用离散法，此模型第k组气泡的输运方程为式（6）:

式中，ρg为气相密度，kg/m3；Bk,c为气泡聚并生成函数；Dk,c为气泡聚并消亡函数；Bk,b为气泡破碎生成函数；Dk,b为气泡破碎消亡函数。

天然气水合物分解出的气泡破碎聚并的原因主要考虑为湍流碰撞，文献[17]中Luo与Lehr的气泡聚并破碎模型被广泛应用。本文采用Luo等基于破碎机制建立的气泡破碎模型对水合物分解出的气泡分布大小进行模拟。

1.4 模拟方法可靠性

1.4.1 网格无关性验证为保证CFD数值模拟结果的准确性并尽量减少计算量，首先进行网格无关性验证。选取内径为300 mm的垂直管，网格尺寸为 4.0、5.0、6.0、6.5 mm。图 2为在入口气泡平均直径为0.45 mm、气含率为0.3条件下所得气泡平均直径沿垂直管轴向的分布。4组模拟结果具有相同的气泡大小分布趋势，网格尺寸为6.0 mm与6.5 mm的轴向平均气泡大小基本相同，此时所得到的气泡大小与另外两组存在一定误差，因此为方便计算，采取6.0 mm的网格尺寸。

图2 网格无关性检验Fig.2 Grid independence test

1.4.2 实验验证依据相关实验[18]，选取气泡直径分布和流动压降对数值模型的可行性进行验证。表1为不同气相体积分数与不同流速的7种工况。表2为气泡从直井2.500 mm递减到0.200 mm的6组气泡直径分布。

表1 模拟工况Table 1 Simulate condition table

表2 气泡初始直径分布Table 2 Initial bubble diameter distribution

表3为单位管长压降模拟值与实验值的对比。由表3可知，工况5-7单位管长模拟值与实验值的变化趋势相同且两者间的相对误差均小于12%。因此该模型可较好地模拟水合物在管道中流动的压降流动变化规律。

表3 实验和模拟条件下单位管长压降对比Table 3 Comparison of pressure drop of unit tube length under test and simulation conditions

图3为工况1—4气泡直径分布模拟值与实验值的对比。由图3可知，模拟得到的气泡直径分布与实验值气泡直径分布规律相同，均呈现正态分布。其中在0.5 mm以下的气泡占比较大，大于1.5 mm的气泡占比较少，因此该数值模型可较好的模拟水合物颗粒在分解出气相后的气液固三相流动。

图3 气泡直径分布实验值与模拟值对比Fig.3 Comparison of experimental and simulated values of bubble diameter distribution

2 结果与分析

2.1 湍动能耗散率对流动特性的影响

图4为气含率0.1时，不同湍动耗散率下轴向距离气泡个数与气泡体积分数分布。由图4(a)可知，气泡个数随轴向距离呈现先增后减的趋势，这是由于气泡在入口处发生明显的破碎现象，由于浆液在入口处受到湍动能耗散率的影响较大，使气泡碰撞破碎强度增强，气泡个数随之增多；在管道轴向距离2 m处气泡的聚并与破碎作用趋于平衡，气泡个数保持不变。当湍动能耗散率较低，即浆液湍动较弱时，气泡的破碎作用相对较弱，气泡的数量增加速率较慢，随着湍动能耗散率增大，气泡的破碎作用增强。图4(b)为气含率0.1时改变浆液湍动耗散率得到的气泡体积分数分布。由图4(b)可知，由于气含率0.1时水合物颗粒刚开始分解，因此浆液的湍动能耗散率的变化只是使气泡大小及其分布宽窄有所变化，气泡的最高体积分数在0.2～0.3，变化较为平稳。

图4 气含率0.1、不同湍动耗散率下轴向距离气泡个数与气泡体积分数分布Fig.4 Gas holdup 0.1,the number of bubbles and the number of bubbles at axial distance under different turbulent dissipation rates

图5 为气含率0.3时改变湍动能耗散率计算的气泡大小和气泡体积分数分布。由图5可知，在该气含率条件下，如果不外加机械能量，即湍动能在0.5 m2/s3时气泡体积分数表示的气泡大小分布呈现气泡不均的二次聚集，气泡分布很宽；当增大湍动能耗散率时，气泡破碎作用增强，气泡大小分布变为较小的单峰分布。当湍动能增大到1.5 m2/s3以上时，浆液流动由气泡不均匀分布转变为均匀分布的湍流。由此可见，气含率相同时，由于耗散率的不同，可导致流动处于不同的流动区域。

图5 气含率0.3、不同湍动耗散率下轴向距离气泡个数与气泡体积分数分布Fig.5 Gas holdup 0.3,the number of bubbles and the number of bubbles at axial distance under differ ent turbulent dissipation rates

2.2 气含率对流动特性的影响

对气泡尺寸以2.500～0.200 mm递减的方式分成5组，设置Ratio Exponent值为3.25。将水合物浆中初始气含率分别取0.10、0.25、0.30和0.35四种工况，5组气泡沿轴面流态分布所占气相体积分数分布如图6所示。

图6 不同气含率不同尺寸的气泡沿轴向的体积分数分布Fig.6 Axial distribution of bubbles with different gas holdup at different sizes

在气含率为0.10和0.25的条件下，在管道入口2 m前，气泡处于不断聚并的状态，气泡Sauter直径不断增大，但增加趋势变缓。在气含率为0.35的条件下，气泡在管道入口1 m左右即达到聚并和破碎平衡，之后Sauter直径不再发生变化；在气含率为0.10和0.25时气泡破碎速率较低，由于湍流涡体的动能很小，不足以克服小气泡的表面张力，因此在气含率较小时气泡破碎率较低；气含率大于0.25时，直径在2.5 mm的气泡聚并速率增加较快，气泡聚并作用增强。

图7为在不同气含率的条件下，管道中轴面上气泡平均直径分布云图。由图7可知，在垂直入口处的初始气泡直径0.55 mm左右，气泡聚并作用不明显，由于气泡聚并在管道中逐渐占主导作用，因此气泡会迅速增大，增至1.50 mm左右，在气泡聚并和破碎达到平衡时，气泡直径趋于稳定，在管道出口处气泡直径变化不大。在气含率较低时，气泡间距离较大，以小气泡为主，气泡聚并作用很弱，气泡的湍动速度不足以使气泡相互碰撞，因此由湍动涡体引起的气泡聚并作用基本可以忽略；当气含率在0.25到0.30时，气泡个数增多，气泡碰撞频率和聚并频率增大，导致气泡直径变大且分布较宽；当气含率进一步增大至0.35时，由于气泡聚并作用形成大气泡的体积分数明显增大，浆液的湍动耗散率随之增大，因此气泡直径大小分布沿轴向变化更快，图7(a)中，在管道入口1/4处，气泡直径接近2 mm，气相结构与低气含率条件下相比发生明显转变。

图7 不同气含率的管道轴向位置气泡平均直径Fig.7 Average diameter of bubbles in axial positions of pipes with different gas holdup

2.3 进气速度对气泡变化的影响

图8 、9分别为改变表观气相、表观液相和固相速度，其他两相不变时气泡和浆液上升速度径向分布的模拟计算结果。由图8（a）和图9（a）可以看出，表观气速为0.25 m/s时，气泡和浆液上升速度沿径向分布较为均匀，管道中心气泡三相之间的滑移速度与单气泡上升速度接近；当流动进入气体分解量较大后，浆液速度和气泡速度上升径向分布不均匀性增加，三相之间的滑移速度明显大于单气泡流动速度，与进气速度成线性关系，这主要是分解出的气泡聚并成大气泡的加速作用造成的。由图8（b）、（c）和图 9（b)、（c）可以看出，相同进气速度下，表观液速较高时，浆液速度和气泡上升速度径向分布较均匀，三相间的滑移速度与单气泡上升速度接近，表观液速降到1.02 m/s时，浆液速度和气泡上升速度沿径向分布均匀下降，三相滑移速度增大。