一种预测水平及近水平管道中气液分层波动流界面形态、持液率和压力梯度的新模型

《International Journal of Heat and Fluid Flow》：A new model to predict the interface configuration, liquid holdup, and pressure gradient of gas-liquid stratified-wavy flow in horizontal and near horizontal pipes

　　

# 水平管道中气液两相流新型模型

- 提出了一种水平管道气液两相流新型模型。
- 发展了基于能量最小的气液构型模型。
- 考虑了流型结构和界面波对剪应力的影响。
- 提出了界面摩擦因子的经验关联式。
- 液相含气率和压力梯度的预测误差分别为14.42%和17.95%。

## 引言

气液分层波状流是工业过程中常见的流型，特别是在油气行业。在天然气生产和输送过程中，管道中轻质烃类和水的积聚可能导致额外的压力损失、内部腐蚀和冰堵（Zhang和Sarica, 2011, Paolinelli和Nesic, 2021, Meng等, 2001）。因此，准确预测液相含气率和压力梯度对于湿气管道的设计和运行至关重要。过去几十年，人们在开发经验关联式和机理模型以预测流场关键参数方面付出了巨大努力（Badie等, 2000）。

早期的经验关联式，如Lockhart-Martinelli（1945）和Beggs-Brill（1973），由于其简单性和便利性而广泛应用于现场。然而，这些第一代模型缺乏严格的理论框架，且在不同流态条件下的预测结果往往难以解释。为解决这一问题，Taitel和Dukler（1976）提出了一种基于各相动量平衡方程的数学模型，假设气液界面形状为光滑平坦的形态（图1(a)），且界面摩擦因子与气壁摩擦因子相等（$f_i = f_g$）。该平面模型为两相流的物理机制提供了基础性见解，但后续研究表明，其预测性能对所假设的界面几何形状和剪应力闭合关系高度敏感（Wang等, 2021）。在实际中，当气相流速较高时，气液界面通常变得弯曲并覆盖有波浪，这会显著增加界面摩擦。为克服这些局限性，一系列第二代机理模型被开发出来，根据其各自的改进可大致分为两类：(i) 界面几何形状的改进，和 (ii) 界面剪应力闭合关系的改进。

**关于界面几何形状：** Hamersma和Hart（1987）以及Hart等（1989）提出的ARS（表观粗糙表面）模型引入了一个凹形界面，其切向方向具有均匀厚度的液膜，如图1(b)所示。虽然这一改进优于平面界面假设，但该膜分布本质上是由重力和界面剪切驱动的特殊情况，可能不具有普适性。双圆模型（Chen等, 1997）和MARS（改进的表观粗糙表面）模型（Grolman和Fortuin, 1997）进一步完善了几何闭合关系，但它们仍然依赖于预设的界面形状，分别如图1(c)和(d)所示。

Brauner等（1998）首次尝试基于能量考虑确定界面构型，但其模型未考虑气液两相的动能。随后，Du（1999）将该模型扩展到包含液相动能，用于小径管中的凝结流动。在我们之前的研究中，在零净液流（ZNLF）条件下观察到凸形界面形状（Wang等, 2021），并基于最小能量原理成功预测了凸形界面形状。然而，在一般的共流气液分层流中，界面形状由湍流气相与液层之间的复杂相互作用决定，这意味着必须考虑从气相到液相的机械能传递。现有模型尚未完全纳入界面曲率对两相动能的影响，这导致界面形状预测可能存在不准确，特别是在界面曲率显著时。

**关于界面摩擦因子：** 普遍认为波浪是增强界面阻力的主要原因，摩擦因子与波浪特性（如波长和波幅）密切相关（Chen等, 1997, Grolman和Fortuin, 1997）。ARS模型采用表观粗糙度方法估算界面摩擦因子。对于湍流气体流动，摩擦因子由Eck（1966）获得的公式(1)计算：

$$f_i = 0.0625 \log_{10} \frac{1}{5\text{Re}_g} + \frac{k/D}{3.715^2}$$

其中相对粗糙度$k/D$与管壁湿润部分$\theta$上的液膜平均厚度相关联：

$$k/D = 2.3\frac{\delta_l}{D}, \quad \delta_l = \frac{\varepsilon}{4\theta}$$

这为界面粗糙度提供了一种可行的工程近似，但未将$f_i$与实际的波浪特性直接联系起来，且粗糙度关联式的经验性质限制了其普适性。Grolman和Fortuin（1997）指出，使用光滑管摩擦因子计算界面摩擦因子不够准确，在考虑波浪对界面摩擦因子的影响后，采用了改进的Eck（1966）模型来预测界面摩擦因子。与ARS模型不同，为了将振荡波视为相对粗糙度，采用迭代程序计算摩擦相对粗糙度$k/D$。然而，波浪的形成和生长过程，以及波浪特性（包括波幅、波长等）与界面摩擦因子的关系无法通过MARS模型预测。

Chen等（1997）开发了修正关联式，将$f_i$与气速和波状流的转换准则联系起来，摩擦因子表示为公式(2)：

$$\frac{f_i}{f_g} = 1 + 3.75\varepsilon\Theta\left[\frac{u_{sg}}{u_{sg,t}} - 1\right]^{0.20} - 0.08$$

其中

$$u_{sg,t} = 4\varepsilon\frac{\mu_l}{\sqrt{\rho_l(\rho_l - \rho_g)g\sigma}}$$

此处$u_{sg,t}$为转换到分层波状流的最小表观气速。正如Mascarenhas等（2015）所指出的，波浪同时贡献了形体阻力和粘性阻力，且$f_i$对波浪特性的依赖关系是多维和复杂的。上述方法的一个关键局限在于准确预测$u_{sg,t}$至关重要。在工况条件下，若$u_{sg}$低于$u_{sg,t}$，这些关联式可能不成立。

Banafi和Talaie（2014）观察到了一种新的气液界面形状，如图1(e)所示。由于气相速度头向静压头的能量转换，靠近管壁的液膜上升并形成一层。基于管道中两种流体相几何参数的适当闭合关系，以及考虑了涟漪影响的界面摩擦因子拟合经验关联式，建立了一种预测界面形状、压力降和液相含气率的机理模型。在该模型中，$f_i$可通过公式(3)预测。迄今为止，一种不依赖于经验阈值且明确考虑波浪特性的通用且物理合理的关联式仍然是一个挑战。

当$\text{Re}_g < 17758$时：

$$f_i = 0.02066 \log_{10} \frac{1}{5\text{Re}_g} + \frac{k/D}{3.715^2}$$

其中$k/D = 0.01249\varepsilon\Theta$

当$\text{Re}_g \geqslant 17758$时：

$$f_i = 0.07452 \log_{10} \frac{1}{5\text{Re}_g} + \frac{k/D}{3.715^2}$$

其中$k/D = 3.5196\varepsilon\Theta$

综上所述，尽管现有模型已取得进展，但在气液分层波状流的机理建模中仍存在两个关键缺口：(1) 缺乏一种在一致的能量最小化框架下同时考虑两相动能贡献的通用界面形状预测方法；(2) 缺乏一种直接从不稳定理论导出的波浪特性中建立的、而非依赖经验粗糙度或速度阈值近似的稳健界面摩擦因子闭合关系。

因此，本文旨在开发一种新的机理模型，同时预测水平和近水平管道中的界面构型、液相含气率和压力梯度。本工作的具体创新体现在以下方面：

(1) 提出了一种基于能量最小化原理的新界面形状预测方法，明确纳入了气液两相的动能。与先前的基于能量的模型不同，本文方法不预先设定特定的几何形状（如平面、凹面或凸面）。而是通过最小化系统总能量来确定最可能的界面构型，从而在各种流动条件下为界面曲率提供更加通用和物理一致的描述。

(2) 建立了一种将$f_i$与波浪特性联系起来的物理基础界面摩擦因子模型。我们基于能量守恒进行Kelvin-Helmholtz不稳定性分析，在波动过程中考虑重力对表面波的影响以量化波浪特性。然后直接从这些波浪特性建立新的界面摩擦因子关联式，从而克服了经验粗糙度基础或阈值依赖关联式的局限性。

(3) 建立了完全集成的机理模型，将上述两个子模型在统一的动量和能量平衡框架内相结合。该模型能够在不依赖特定几何假设或波浪状态经验转换准则的情况下，同时预测界面形状、液相含气率和压力梯度。新模型经过大量实验数据验证，包括本研究数据和文献中综合的分层波状流数据库。预测性能与现有的第二代模型进行了比较，证明了所提出方法在广泛工况条件下的优越性和通用适用性。

## 实验

实验装置设计用于测量水平管道中气液分层波状流的液相含气率和压力梯度，如图2所示。该装置主要由螺杆压缩机、离心泵、储液罐、气液混合器、测试段、气液分离器和数据采集系统组成。测试段为一根长4米、内径40毫米的透明亚克力管。管道使用支撑装置仔细调平。

## 基于能量考虑的界面构型

分层流是气液两相流的常见流型，流动过程始终伴随界面波现象，这被称为波状分层流。根据界面波的强度和形式，波状分层流可分为二维波流、三维波流和滚动波流。当界面波幅相对较小时，波浪可以忽略不计，流动可简化为理想光滑分层流。

## 求解方法

如前所述，控制方程可表示为$f(\varepsilon) = 0$，这是液相含气率的隐函数。注意，给定$\varepsilon_0$值，$\phi_0^P$被确定，对于每个指定的$\phi^*$值可获得$\phi_0$值。在特定的流体性质参数（$\rho$, $\mu$, $\sigma$等）和流动条件（$D$, $\beta$, $u_{sg}$, $u_{sl}$等）下，液相含气率应满足控制方程，即$f(\varepsilon_0) = 0$，然后，压力梯度通过稳态气相或液相动量方程计算得出。

## 结果与讨论

为进一步分析气液构型，推导了无量纲形式。将几何和剪应力闭合关系（分别为公式(15)和(17)）代入组合动量平衡方程（公式(14)），得到以下无量纲方程：

$$We_{sg} f_g \frac{\pi - \phi_0}{1-\varepsilon^3} + f_i \frac{\pi - \phi^*}{\sin\phi_0} \frac{\varepsilon^3}{1-\varepsilon^3} \sin\phi^* - We_{sl} f_l \frac{\phi_0}{\varepsilon^3} = 0$$

其中

$$We_{sg} = \frac{\rho_g u_{sg}^2 D}{\sigma}, \quad We_{sl} = \frac{\rho_l u_{sl}^2 D}{\sigma}$$

上述公式可表示为液相含气率的函数$f(\varepsilon) = 0$，这表明气液两相流中的液相含气率可由该方程确定。

## 结论

通过考虑界面曲率对气相和液相动能的影响、表面自由能和势能，并结合最小能量原理，阐明了气液构型的机制。通过无量纲控制方程的分析，理论结果表明，气液共流的构型与表观气相Weber数、表观液相Weber数以及液相含气率有关。

## CRediT作者贡献声明

王俊杰：资源提供，资金获取，数据整理。赵敏：撰写——初稿，数据整理。柴集峰：撰写——审阅与编辑。曹驿奇：调查研究。潘杰：资源提供。

## 利益冲突声明

作者声明，他们不存在任何已知的可能影响本文所报告工作的竞争性财务利益或个人关系。

## 致谢

作者衷心感谢国家自然科学基金（项目编号：52304038，52274064）和西安石油大学研究生创新基金（YCX2512024）的资助。

作者还感谢审稿人提出的宝贵意见，帮助改进了本文。

王俊杰 | 赵敏 | 柴集峰 | 曹驿奇 | 潘杰