基于行波法的段塞流瞬态捕捉模型建立与验证

2020-08-17王冬旭胡其会李玉星李爽王权

化工进展 2020年8期

王冬旭，胡其会，李玉星，李爽，王权

（1 中国石油大学（华东），山东省油气储运安全省级重点实验室，山东青岛266580；2 中国市政工程中南设计研究总院有限公司，湖北武汉430060）

段塞流是油气混输管线中的常见流型，当管线处于启动、停输、清管、输量变化甚至正常生产时，都有可能产生[1]。其流动的间歇性将引起管道中持液率和压力的急剧波动，并使得下游管道及设备承受间歇性应力冲击[2-4]。同时，管道末端产生的大段塞会引起下游处理设备中液位剧烈波动。为保证管线和下游处理设备的安全生产，必须掌握段塞流的特征规律。为此，自20 世纪40 年代起，国内外学者对段塞流开展了相关研究[5-7]。

Dukler 等[8]最早建立了气泡和段塞在时间和空间内呈均匀分布的段塞流稳态模型。此模型可以预测段塞流重要参数的平均值，如段塞长度、气泡平移速度和压降。Ishii[9]在20世纪70年代建立的欧拉双流体模型广泛应用于核工业。Bendiksen 等[10]建立的OLGA模型在石油工业中被广泛应用，是最早以双流体模型为基础开发并使用的瞬态计算代码之一。Zheng 等[11]通过研究简单起伏管路内的段塞流建立了一个SINK&SOURCE 模型。 Nydal 和Banerjee[12]，Taitel 和Barnea[13]也建立了段塞模型，以上模型都需要额外的段塞流参数，如段塞形成条件、段塞频率、段塞长度等。

应用于双流体模型的线性稳定性理论经常被用来确定分层流的稳定性。Taitel 和Dukler[14]将该理论应用于不考虑剪切应力的理想流体，得到非黏性Kelvin-Helmholtz（IKH）稳定性准则。该准则可以预测伯努利效应产生的吸力大于流体重力时气液界面波的生长。Barnea 和Taitel[15]随后通过对黏性流场的线性稳定性分析，得出黏性Kelvin-Helmholtz（VKH）稳定性准则。Issa和Kempf[16]论证了当气液流速处于IKH 和VKH 准则之间时，使用交错网格的双流体模型能够捕获分层流过渡到段塞流过程中气液界面产生的扰动。气液两相流在初始条件下自由发展，减少了对物理现象模型（流型转变、段塞形成等）的依赖。Issa 等[16]利用此模型模拟了V 型两相流起伏管道内的段塞频率、段塞平均长度和持液率，得到与实验数据一致的波动趋势。

Renault[17]以Issa 的研究为基础，基于IKH 和VKH准则实现了分层流过渡到段塞流的模拟，捕捉段塞前后界面运动过程，取得了较好的预测效果。尽管Renault 模型基于IKH和VKH 准则捕捉段塞前后液面变化，然而其求解过程中将液相方程转化为浅水方程，并在单元格之间采用黎曼（Riemann）精确解，使得求解速度较慢。因此，需对模型进一步改进，使其保持精度的同时提高计算效率。

许仁义等[18-20]将行波法与Riemann 精确解相结合，对稀疏波采用单波近似，实现了浅水方程的快速求解，实现了高分辨率、高计算精度的效果。因此本文尝试将行波法引入Renault 模型，建立新的段塞捕捉模型及求解方法。同时为保证行波法适用于所有计算单元，对可能出现的干区采用薄液膜进行简化处理，在保证模型精度的条件下进一步提高求解速度。

1 段塞捕捉模型建立

1.1 模型基本方程

Renault 建立的段塞捕捉模型假设条件为：①液相不可压缩；②理想气相；③相对于液相动量，气相动量可以忽略；④求解液相动量方程时，气相为局部不可压缩。计算域内网格被划分成段塞网格和气泡网格两类。段塞网格存储平均液相速度（Uslug）、气泡网格存储含液率（β）、平均液速（Ul）、气体表观流速（USg）、压力（p）存储在交错网格上，计算网格如图1 所示，模型方程如式(1)～式(5)。

1.2 闭合关系式

图1 计算网格

液相和壁面、气相和壁面以及气液界面间的剪切应力关系如式(6)～式(8)。液相与壁面的摩阻系数、气相于壁面的摩阻系数如式(9)和式(10)，气液相界面摩擦因子采用由Cohen 与Hanratty[21]提出的常数值如式(11)，气泡平移速度采用Bendiksen[22]关系式，如式(12)。

紊流情况下，根据弗鲁德数（Fr）计算参数C0、Ut0。

1.3 模型求解

Renault 模型求解过程中气相和液相方程的求解是交替进行的。气相方程采用一阶迎风格式化简为式(15)。

由于气相方程符合三对角矩阵形式，因此采用追赶法即可求解每个网格中n+1 时刻的气相状态。

液相方程可化简为浅水方程形式，如式(16)和式(17)。

液相求解过程较为复杂，采用戈杜诺夫（Godunov）方法[22]，在每个单元边界求解浅水方程Riemann问题的精确解，进而构成整个流场的数值解[式(18)和式(19)]，求解过程如图2（图2中xj与xj+1左右两侧均为稀疏波）。

图2 戈杜诺夫方法求解过程

1.4 模型简化

由于Renault 模型中液相求解需要对网格边界产生的稀疏波区域持液率（βLM、βRM）进行核算，计算量较大。为此，本文提出采用行波法求解单元边界Riemann问题，对稀疏波采用单波近似，同时用薄液膜代替干区，使行波法适用于干区计算，从而简化求解过程，提高计算效率。

图3 Riemann问题精确解结构

图4 行波法求解过程

对比图2、图4 求解过程发现，对稀疏波采用单波近似，并用行波法求解液相方程极大地简化了求解过程。浅水方程求解过程会出现干区情况，为使行波法适用于所有计算单元，使用薄液膜代替干区，薄液膜的液量来自相邻单元。为了解行波法处理干湿边界问题时的性能，采用一种长100m、宽1m、初速度为0、水面不连续的一维渠道浅水方程作为研究对象。用0.1m 的液膜代替干区，初始条件如式(22)。

计算范围x∈[0,100m]被分成100 个网格，空间步长dx=1m。溃坝初始时刻如图5所示。初始时刻水面间断处（x=50m）产生向左传播的稀疏波和向右传播的激波。图6 为CFL=0.8、t=5s 时的计算结果。

从图6可以看出，用薄液膜代替干区的简化方法得到的数值解与精确解相一致，可以反映液相在干湿边界的流动状态。将Renault 模型单元边界存在的干区均用薄液膜代替。此简化处理使得行波法适用于所有网格计算，极大地提高了运算速度，简化后的全新计算流程如图7所示。

图5 初始液面状态

图6 溃坝模拟

图7 计算流程

实验发现段塞运动过程中头部呈竖直断面、尾部呈曲面（图8）。为使计算结果与实际段塞运动过程相一致，需对生成的段塞进行特殊处理。段塞头部速度Ufront根据液体质量守恒核算，气泡头部速度UBendiksen由Bendiksen 关系式核算，计算过程如图9。

图8 段塞运动

图9 段塞计算过程

2 模型验证

2.1 实验系统

本实验的介质采用空气和水，实验装置示意图如图10 所示。实验管道内径D=40mm，总长45m，其中水平测试段管长10m。测试管段采用透明有机玻璃材质，以便观察管内气液流型变化。管道混合器出口处布置有压力传感器P0，水平段内布置有间距均为0.8m 的双平行电导探针（CP1～CP3）与压力传感器（P1～P3）。采用Iotech6220 采集卡进行持液率与压力信号的采集，采样频率为50Hz，采样时间为60s。实验中的气体表观流速范围为0～4m/s，液体表观流速范围为0～1.6m/s。

2.2 实验数据验证

为验证模型的段塞捕捉能力，以实验管道为基准进行模拟计算。管道剖面尺寸如图11，模型计算域为[0,10]m，网格宽度Δx=0.02m，管径D=40mm，时间步长Δt=0.01s，模拟时间10min。气体表观流速范围为0～4m/s，液体表观流速范围为0～1.6m/s。入口采用质量流量边界条件，出口采用压力边界条件，计算结果如下。

图10 实验装置示意图

图11 管道剖面示意图

图12 为气体表观流速0.9m/s、液体表观流速0.45m/s、时间t=30s时的模拟结果。图13为CP1处1min内实验与计算持液率对比。从图中可以看出，随着管道长度的增加，液面渐渐产生波动，流型由分层流转变为波浪流；随着波动幅度的增大，流型进一步发展成为段塞流；计算CP1 处1min 内产生的段塞数略大于实测值，持液率变化与实测结果相一致，且段塞头部呈竖直面，气泡头部呈曲面；本模型可以捕捉段塞流生成与发展过程，实现段塞流压降精确计算。

图12 t=30s持液率计算结果

图14、图15 分别为模拟实验工况30min 内的平均压降、平均段塞长度值与实验值对比。图16为两者计算值与实测值相对误差概率密度分布。从图中可以看出，压降、段塞长度的相对误差分别在25%、30%以内，且主要分布在20%以内，压降计算准确度要高于段塞长度。

图16 30min内压降与段塞长度的相对误差概率分布

为进一步了解本文模型简化方法对计算效率的影响，对比Renault 模型与本文简化模型在50、100、200、300、400 网格数下的计算耗时。计算硬件条件CPU 为Intel i7-6700 3.4GHz、内存为16GB，操作系统为64 位Windows 7，运算软件为Matlab 2017，对比结果如图17所示。

由图17 可以看出，两种模型的计算耗时与网格数均成正相关，随着网格数的增加，计算量增大，耗时逐渐增加；当网格数较低时，两种模型计算耗时差距较小，随着网格数的增加，本文采用薄液膜代替干区与行波法求解液相方程组的简化方法优势逐渐明显，运算时间相较于原模型平均减少28%。