引 言

水平管气液两相分层流是石油天然气输送工程和化工工程中最常见的两相流流型之一^[1]，其特点是除了气液两相与管壁间存在相互摩擦外，由于密度差，第一相(气相)流速明显大于第二相(液相)，还存在着气液界面相互作用^[2]。目前对气壁剪切和液壁剪切的预测已取得较好的结果^[2-4]，但对界面剪切应力的预测尚没有一致定论，也是这一领域的重点和难点^[5-7]。

目前已建立的机理模型通常是通过构建分层流封闭模型或互补的局部模型，预测分层流流动的压降、持液率和剪切作用等^{[3, 7]}。不同模型的区别主要体现在对界面波和界面形状不同的处理方法^[8-10]。例如Taitel Y和Dukler A E^[5-6]建立的一维两流体模型假设界面为一水平光滑壁面，其摩擦因子等于气壁摩擦因子。这种假设在湍流较强的流动中预测效果并不理想，为此，一些研究者区分了光滑界面和波动界面^[11-13]，并分析了界面波对气液界面摩擦因子本构关系的影响^[11]；还有一些学者将气液界面处理为圆弧凹面^[3]，或者将气液界面处理为水平粗糙壁面^{[1, 14]}，建立了不同的分层流模型。在前人工作基础上，Sidi-Ali K和Gatignol R^[15] 建立了将气液界面定义为一个无滑移的水平移动壁面的CFD模型，本文称之为移动壁面模型，并由移动壁面摩擦速度得到气液界面剪切应力，从而改进了Taitel Y和Dukler A E模型。然而该模型预测值仍然明显小于实验间接测量值，这是由于移动壁面模型假设的界面速度为液速平均值，远小于实际的界面速度，同时，这类模型的界面性质并不能引起足够的湍能损失。

本文主要目的是利用FLUENT软件建立基于气相流动的CFD模型，并尝试将气液界面处理为一个固定的、具有均匀剪切应力的水平滑移壁面，称之为剪切滑移壁面模型，以提高对气液界面剪切应力的预测。由于可以直接考虑壁面性质对气相流场分布的影响，使得预测方法更为简单有效。

1 模型建立

依据Taitel Y和Dukler A E的两流体模型^[5-6]，将气液两相单独讨论，分别对气液两相建立动量方程

$ - A_{\rm{G}} \left( {\dfrac{{{\rm{d}}p}}{{{\rm{d}}z}}} \right) - \tau _{{\rm{WG}}} S_{\rm{G}} - \tau _{\rm{I}} S_{\rm{I}} + {\rm{g}}\rho _{\rm{G}} A_{\rm{G}} \sin \delta = 0 $

(1)

$ - A_{\rm{L}} \left( {\dfrac{{{\rm{d}}p}}{{{\rm{d}}z}}} \right) - \tau _{{\rm{WL}}} S_{\rm{L}} + \tau _{\rm{I}} S_{\rm{I}} + {\rm{g}}\rho _{\rm{L}} A_{\rm{L}} \sin \delta = 0 $

(2)

式中：$A_{\rm{G}}$，$A_{\rm{L}}$-气相和液相截面面积，m²；

$\tau _{{\rm{WG}}}，\tau _{{\rm{WL}}}，\tau _{{\rm{I}}}$-气壁、液壁和气液界面剪切应力，Pa；

$S_{\rm{G}}，S_{\rm{L}}$-气相和液相湿周，m；

$S_{\rm{I}}$-管道横截面气液界面水平长度，m；

g-重力加速度，g=9.8 m/s²；

$\rho _{\rm{G}}，\rho _{\rm{L}}$-气相和液相密度，kg/m³；

$\delta $-圆管倾斜角度，rad；

dp/dz-轴向压降，Pa/m。

图 1为气液两相分层流示意图，各参数见图 1。

式(1)反映了气相在一封闭固壁管内流动情形。$\tau _{{\rm{I}}}$可由式(3)求得

$ \tau _{\rm{I}} = \dfrac{D}{4}\left( {\dfrac{{2h_{\rm{L}} }}{D} - 1} \right)\left( {\dfrac{{{\rm{d}}p}}{{{\rm{d}}z}}} \right) - \dfrac{{\cos ^{ - 1} \left( {2h_{\rm{L}} /D - 1} \right)}}{{\sqrt {1 - \left( {2h_{\rm{L}} /D - 1} \right)^2 } }}\left[{\dfrac{D}{4}\left( {\dfrac{{{\rm{d}}p}}{{{\rm{d}}z}}} \right) + \tau _{{\rm{WG}}} } \right] $

(3)

式中：D-圆管直径，m；

$h_{\rm{L}}$-液面高度，m。

可见只有准确预测压降和气壁剪切应力才能更好地预测界面剪切应力。由于CFD预测$\tau _{{\rm{WG}}}$ 时受到网格和湍流模型的影响(同时影响压降预测)，使用上式很难得到满意结果^[16]。现在假设$\tau _{{\rm{I}}}$为已知变量(作为数值计算初始值)，应用滑移壁面条件建立模型，通过对比模拟和实验的气相剖面速度分布直接预测$\tau _{{\rm{I}}}$，从而避开应用$\tau _{{\rm{WG}}}$和dp/dz预测$\tau _{{\rm{I}}}$。为便于研究计算，壁面的剪切应力做均匀处理，只考虑轴向z方向正向剪切应力，忽略水平方向x方向的剪切应力，即所谓的剪切滑移壁面模型，见图 2。单向均匀化处理的结果忽略二次流的影响，壁面边界处理的结果忽略了表面张力的影响，而二次流和表面张力对界面剪切预测精度的提高很有限^[17]。

实际的气液界面剪切应力分布并不均匀，靠近管壁的剪切应力小，界面中心线位置最大^{[3-4, 14, 17-18]}。剪切滑移壁面模型假设剪切均匀分布，在气液两相的壁面交界处的不连续性很强，这一特征与实际流动非常相似^[4]。此外，移动壁面模型的均匀界面速度分布与实际两相界面速度分布也不相同^[18]，两种类比模型对界面的均匀简化处理的目的之一是为了使用简单边界条件来获得气速的收敛解，同时，简化条件对气相主流区气速分布影响并不大^{[4, 14-15]}，这是本文模拟方法的主要依据和条件。

2 模型方程 2.1 控制方程

建立气相连续性方程和动量方程

$ \dfrac{{\partial \rho }}{{\partial t}} + \nabla \cdot \left( {\rho U_{\rm{i}}} \right) = 0 $

(4)

$ \dfrac{\partial }{{\partial t}}\left( {\rho U_{\rm{i}}} \right) + \nabla \cdot \left( {\rho U_{\rm{ij}}} \right) = - \nabla p + \rho {\rm{g}}_{\rm{i}} + \\ \hspace{2em}\nabla \cdot \left\{ {\mu \left[{\nabla U_{\rm{i}} + \left( {\nabla U_{\rm{i}}} \right)^T } \right]} \right\} $

(5)

式中：

$\rho$-密度，kg/m³；

t-时间，s；

$U_{\rm{i}}$-流速张量，m/s；

$g_{\rm{i}}$-重力加速度张量，m/s²；

$-\nabla p$-压力梯度，Pa/m；

$\mu$-动力黏度，Pa·s。

使用$k-\varepsilon$两方程湍流模型(对于本问题，$k-\omega$模型的预测效果没有$k-\varepsilon$模型好)，湍流脉动动能k(m²/s²)和湍流耗散率$\varepsilon$(m²/s³)的传递方程为

$ \dfrac{\partial }{{\partial t}}\left( {\rho k} \right) + \dfrac{\partial }{{\partial x_{\rm{i}} }}\left( {\rho kU_{\rm{i}} } \right) = \dfrac{\partial }{{\partial x_{\rm{j}} }}\left[{\left( {\mu + \dfrac{{\mu _{\rm{t}} }}{{\sigma _{\rm{k}} }}} \right)\dfrac{{\partial k}}{{\partial x_{\rm{j}}}}}\right] +\\ \hspace{3em} \rho P_{\rm{k}} - \rho \varepsilon $

(6)

$ \dfrac{\partial }{{\partial t}}\left( {\rho \varepsilon } \right) + \dfrac{\partial }{{\partial x_{\rm{i}}}}\left( {\rho \varepsilon U_{\rm{i}} } \right) = \dfrac{\partial }{{\partial x_{\rm{j}} }}\left[{\left( {\mu + \dfrac{{\mu _{\rm{t}} }}{{\sigma _\varepsilon}}} \right)\dfrac{{\partial \varepsilon }}{{\partial x_{\rm{j}}}}} \right] +\\ \hspace{3em} C_{1\varepsilon } \rho P_{\rm{k}} \dfrac{\varepsilon}{k} - C_{2\varepsilon } \rho \dfrac{{\varepsilon ^2 }}{k} $

(7)

式中：$\mu _{\rm{t}}$-湍动黏度，Pa·s，$\mathop \mu \nolimits_{\rm{t}} = \rho \mathop C\nolimits_{\rm{\mu }} \mathop k\nolimits^2 /\varepsilon$；

$P_{\rm{k}}$-湍流生成项，m²/s³；

$x_{\rm{i}}，x_{\rm{j}}$-空间张量，m；

经验常量：

$C_{\rm{\mu}}$=0.09；

$C_{1\varepsilon}$=1.44；

$C_{2\varepsilon}$=1.92；

$\sigma _{\rm{k}}$=1.0；

$\sigma _\varepsilon$=1.3。

对于稳态问题，上述各方程非稳态项为零。

2.2 边界条件

入口为平均速度入口

$ U = U_0 $

(8)

式中：$U_0$-速度，m/s。

出口边界条件选择压力出口，各量在出口法向n的扩散率为零

$ \dfrac{\partial }{{\partial n}}\left( {U, k, \varepsilon } \right) = 0 $

(9)

管壁处，速度满足无滑移条件

$ U = 0 $

(10)

滑移壁面处，FLUENT在边界上计算切向速度，滑移壁面不使用壁函数，同时，给出壁面剪切应力边界条件

$ \tau _{\rm{I}} = \tau _0 $

(11)

式中：$\tau _0$-剪切应力，Pa。

无滑移管壁面使用低雷诺数两层湍流边界层模型，即黏性支层和对数律层^[19-20]，则

$ {U_\tau } = \sqrt {{\tau _{\rm{I}}}/\rho } = C_\mu ^{1{\rm{ /}}4}{k^{1{\rm{ /}}2}} $

(12)

$ U^ + = \dfrac{U}{{U_\tau }} = \left\{ {\begin{array}{*{20}c} {\dfrac{{\rho U_\tau }}{\mu }l = l^ +, \quad l^ + \leqslant 11.225} \\ {\dfrac{1}{\kappa }\ln \left( {El^ + } \right)\, , \quad l^ + > 11.225} \\ \end{array}} \right. $

(13)

式中：$U_\tau$-摩擦速度，m/s；

l-点到壁面距离，m；

$\kappa$-von Kármán常数，无因次，$\kappa$=0.4187^[19-20]；

E=9.793^[19-20]。

上标：

+-壁面无因次参数，黏性支层的无因次厚度为11.225^[19-20]。

3 模拟方法

Strand Ø^[21]在长35.0 m、直径0.1 m的水平圆管中进行了空气-水两相分层流系统实验，液相表观速度为0.1 m/s。利用FLUENT对前7组(R1~R7)进行模拟。管长10.0 m，直径0.1 m，仅模拟气相流动，考虑重力项(y方向-9.81 m/s²)。管壁内表面建立第一层0.000 1 m、增长率1.2的4层边界层；对于滑移壁面，考虑到较大的速度梯度，对网格进行了相应的细化，见图 3。

模拟的基本过程是通过动量方程更新初值气速和剪切应力，经过压力修正使更新值满足连续方程，在湍流方程的控制下，流场和源项得到更新，将模拟得到的流场与实验观测数据(气相流动分布)进行对比，然后进一步更新初值，重复这个过程，直到模拟流场与实验流场吻合得最好。对于速度分布特征的对比，选择管道轴向5 m处(z=5 m，认为流动已充分发展)的两条线速度分布进行对比：中心垂线速度(z=5 m，x=0)和高0.065 m的水平线速度(z=5 m，y=0.065 m)。为做统一对比，引入无因次高度Y和无因次宽度X

$ Y = \dfrac{{y - h_{\rm{L}} }}{{D - h_{\rm{L}} }}, \;X = \dfrac{{x + x_{\rm{0}} }}{{2x_{\rm{0}} }} $

(14)

式中：$x_{\rm{0}} $-y=0.065 m时中心垂线到管壁的水平距离，m；$x_{\rm{0}} $≈0.047 7 m。对垂线速度的最大速度$U_{\rm{p}}$和最大速度的无因次高度$Y_{\rm{p}}$进行重点对比观察。

4 模拟结果与对比分析 4.1 模拟结果

图 4和图 5对比了Strand Ø的实验、Sidi-Ali K和Gatignol R与本文的剪切滑移壁面模型得到的速度剖面分布图。图 4为垂线速度分布图；图 5为水平线速度分布图，流场关于x=0近似对称，故仅比较0.5≤X≤1.0的速度，实验R6组测量数据有明显误差^{[15, 21]}，不做对比。

除R6组外，剪切滑移壁面模型模拟的$U_{\rm{p}}$与实验相差不超过1.8%，$Y_{\rm{p}}$与实验值相差0~4.8%，模拟结果基本上反映了实验的流场特征。

表 1给出了实验值以及移动壁面模型和剪切滑移壁面模型的气相速度及剪切应力的模拟结果。其中，移动壁面模型通过摩擦速度求得$\tau_{\rm{WG}}$和$\tau_{\rm{I}}$；本文的$\tau_{\rm{I}}$为模拟最后更新的收敛值，$\tau_{\rm{WG}}$通过摩擦速度求得。

剪切滑移壁面模型的滑移壁面滑移速度分布不均匀，利用动量方程预测液相速度$U_{\rm{L}}$必须提供经验关联式来封闭方程，例如依据实验得到的持液率和含液率的关系可以得到$U_{\rm{L}}$值。移动壁面模型的移动壁面速度就是液相平均速度^[15]。

4.2 摩擦因子表达式

剪切应力和摩擦因子通常使用下式建立计算关系^{[5-6, 13]}

$ \tau _{{\rm{WG}}} = f_{\rm{WG}} \dfrac{{\rho _{\rm{G}} U_{\rm{G}}^2 }}{2} $

(15)

$ \tau _{\rm{I}} = f_{\rm{I}} \dfrac{{\rho _{\rm{G}} \left( {U_{\rm{G}} - U_{\rm{L}} } \right)^2 }}{2} $

(16)

式中：$f_{\rm{WG}}$-气壁摩擦，无因次；

$f_{\rm{I}}$-界面摩擦因子，无因次；

$U_{\rm{G}}$-气相截面平均流速，m/s；

$U_{\rm{L}}$-液相截面平均流速，m/s。

由于移动壁面模型和剪切滑移壁面模型将气液界面类比为壁面，即只考虑气相固壁单相流动，故使用与式(15)相似的方法计算$\tau_{\rm{I}}$^[11]

$ \tau _{\rm{I}} = f_{\rm{I}} \dfrac{{\rho _{\rm{G}} U_{\rm{G}}^2 }}{2} $

(17)

Taitel Y和Dukler A E认为气液界面与气壁具有相同的摩擦因子^[5-6]，使用Blasuis公式计算$f_{\rm{WG}}$和$f_{\rm{I}}$^[5]

$ f_{\rm{I}} = f_{{\rm{WG}}} = 0.046Re_{\rm{G}}^{ - 0.2} $

(18)

移动壁面模型得到的$f_{\rm{WG}}$和$f_{\rm{I}}$分别为^[15]

$ f_{{\rm{WG}}} = 1.14Re_{\rm{G}}^{ - 0.45} $

(19)

$ f_{{\rm{I}}} = 0.94Re_{\rm{G}}^{ - 0.427} $

(20)

对于剪切滑移壁面模型，根据模拟收敛结果，计算出气壁和界面摩擦因子，可以拟合得到$f_{\rm{WG}}$和$f_{\rm{I}}$的计算公式(图 6)

$ f_{{\rm{WG}}} = 0.266Re_{\rm{G}}^{ - 0.317} $

(21)

$ f_{{\rm{I}}} =0.3965Re_{\rm{G}}^{ - 0.336} $

(22)

气相雷诺数$Re_{\rm{G}}$和水力直径$D_{\rm{IG}}$为^[13-15]

$ Re_{\rm{G}} = \dfrac{{\rho _{\rm{G}} U_{\rm{G}} D_{{\rm{IG}}} }}{{\mu _{\rm{G}} }}, \;D_{{\rm{IG}}} = \dfrac{{4A_{\rm{G}} }}{{S_{\rm{G}} + S_{\rm{I}} }} $

(23)

液相雷诺数$Re_{\rm{L}}$和水力直径$D_{\rm{IL}}$为^[13-15]

$ Re_{\rm{L}} = \dfrac{{\rho _{\rm{L}} U_{\rm{L}} D_{{\rm{IL}}} }}{{\mu _{\rm{L}} }}, \;D_{{\rm{IL}}} = \dfrac{{4A_{\rm{L}} }}{{S_{\rm{L}} }} $

(24)

式(21)和式(22)的适用范围在9 400$\leqslant Re_{\rm{G}}\leqslant$ 50 000，21 000 $\leqslant Re_{\rm{L}} \leqslant$ 30 000。

4.3 预测结果对比分析

使用Strand Ø实验的气速值，维持液高、管径和介质，对比Taitel Y和Dukler A E、移动壁面模型和剪切滑移壁面模型的预测公式计算效果。表 2给出了Taitel Y和Dukler A E(下标“-TD”)、移动壁面模型和剪切滑移壁面模型预测值，预测值分别用$\tau_{\rm{WG}}$和$\tau_{\rm{I}}$代表气壁剪切应力和界面剪切应力，图 7对比了$\tau_{\rm{I}}$的预测结果。

在9 400$\leqslant Re_{\rm{G}}\leqslant 50 000$和21 000 $\leqslant Re_{\rm{L}}\leqslant 30 000$内，由于在计算气壁剪切应力时，受到$k-\varepsilon$湍流模型和网格划分方案的影响，移动壁面模型和剪切滑移壁面模型在预测气壁剪切应力时存在较大误差，相对误差超过60%。对于界面剪切应力预测误差，移动壁面模型和剪切滑移壁面模型的平均相对误差分别为-20.12%和3.16%，平均绝对误差分别为20.12%和16.30%，方根平均误差百分比分别为27.35%和19.71%，平均实际误差分别为-0.073 9和0.009 6，平均绝对实际误差分别为0.073 9和0.048 4，方根平均实际误差分别为0.124 9和0.067 0。可见，剪切滑移壁面模型的预测值最接近实验间接测量值。 Taitel Y和Dukler A E按照单相流动对气壁剪切应力进行预测(Blasuis公式)，这种方法被证明是行之有效的，但在预测界面剪切应力时，该模型只适用于雷诺数更低的气液两相湍流流动情形^{[8, 22]}。

预测效果证明，气相流场的观测在预测气液界面剪切应力时具有重要意义。由于剪切滑移壁面模型的近壁面速度为不均匀的滑移速度，速度的重新分布一方面引起了更大的湍能耗散，另一方面导致模拟得到的气速稍小(见表 1)，所以在相同气相雷诺数条件下，剪切滑移壁面模型的界面剪切应力预测值大于移动壁面模型，甚至大于实验值(R1，R3~R5)，也正是这个原因，使得剪切滑移壁面模型的整体预测效果更接近实际情况。 但是，提高预测精度的代价是减弱了模型的封闭性。此外，移动壁面模型和剪切滑移壁面模型都是将气液界面类比为水平壁面进行简化计算，因而无法充分考虑界面波或界面形状对摩擦阻力增强的影响^{[2, 7, 12, 23]}，特别是当雷诺数较大时，水平壁面模型的预测值必然小于实际值(R6，R7)。

从模拟方法角度分析，移动壁面模型的移动壁面性质取决于界面剪切应力的计算结果，而不是决定界面剪切应力大小；而剪切滑移壁面模型直接考虑了界面性质对气相流场分布的影响，因此方法更为简单有效。

5 结 论

(1) 建立了将气液分层流界面处理为具有均匀剪切应力的滑移壁面模型，即剪切滑移壁面模型，利用CFD软件模拟气相流动，得到新的界面摩擦因子表达式(式(22))，对水平管气液两相分层流界面剪切应力进行了预测。

(2) 与Taitel Y和Dukler A E及Sidi-Ali K和Gatignol R的计算模型相比，剪切滑移壁面模型得到的界面剪切应力预测值与Strand的实验间接测量值吻合得最好；但是，水平壁面模型无法考虑界面波与界面形状对界面作用的影响，建议当气相雷诺数(至少)大于50 000时，必须考虑界面形状特性。

(3) 对于模型封闭性，剪切滑移壁面模型需要借助更多的实验经验关系式来封闭模型，如持液率(或液高)和液相平均速度等。

(4) 剪切滑移壁面模型考虑了水平壁面性质对气相流场分布的影响，从而可以直接控制关键参数，为CFD方法预测水平管气液两相分层流界面剪切应力提供了一种简单而可行的思路。