双流体模型中人工扩散的影响数值分析

韩朋飞，郭烈锦*，李清平，姚海元，程 兵

(1. 西安交通大学动力工程多相流国家重点实验室，陕西 西安 710049；2. 中海油研究总院，北京 100028)

双流体模型中人工扩散的影响数值分析

韩朋飞1，郭烈锦1*，李清平2，姚海元2，程 兵2

(1. 西安交通大学动力工程多相流国家重点实验室，陕西 西安 710049；2. 中海油研究总院，北京 100028)

研究了双流体模型中加入人工扩散项后不同扩散项系数对模拟的影响。对于竖直管的Water Faucet问题，逐步增加人工扩散系数后非物理振荡完全消除。对于水平管流动病态区域的工况，不采用扩散系数在中等精细度网格上计算的结果同实验比较吻合，但无法得到网格无关解。采用Water Faucet问题确定的人工扩散系数进行计算，水动力学不稳定性被抑制，液塞频率被低估。因此，人工扩散的量随不同流动条件发生变化，需要具体确定。

双流体模型；人工扩散；段塞流；数值模拟

0 引 言

管道内多相流流动模拟对石油工业的正常生产及安全保障有非常重要的意义。目前，石油开采正在向海洋尤其是深海进发，因此超长距离的管道多相混输对管内流动的预测有着更高和更迫切的需求。

根据以往国内外的研究，一维的漂移流模型和双流体模型是模拟预测超长距离管道流动的有力工具。漂移流模型本质上属于混合物模型，无法对相界面上波动的发展进行预测，因此精确度有限。Lin等[1]很早就从理论上证明了双流体模型可以预测界面上波动的发展，并且Issa等[2]也使用双流体模型在数值模拟中捕捉到了界面波动的实时发展。因此，双流体模型是对长距离管线流动进行准确预测的最好手段。但是，传统的标准双流体模型在某些区域不能保证双曲性，界面上小的波动会被无限放大。这种特征表现在数值模拟当中，就是当网格步长不断减小时，不断有更小的波动被计算分辨，而计算对这些波动的放大系数不是有界的，在波长极小时趋于无穷大，从而无法得到唯一正确的收敛结果。

针对这一问题，很多学者采用各种方法进行了研究。一方面，Zanotti等[3]从纯数学的角度上考虑，通过矩阵预处理实现了无条件的双曲性，尽管对算例的模拟显示瞬态结果也是准确的，但是理论上预处理影响了模型的瞬态项，只适用于稳态问题。另一方面，很多学者尝试引入其他物理机制来保证双曲性，得到了一些有用的结论。Chung等[4]发现虚拟质量力的加入并不是完全有益的，有时它会影响计算的稳定性及精确度。Holmas等[5]则指出表面张力系数作用范围有限，人为增大表面张力系数虽然可以抑制小波长不稳定性的发展，但是这部分非物理的能量会转移到长波里，而通过人工增加二阶扩散项耗散掉不稳定的小波长波动的能量是一种合理的方法。Song等[6]通过考虑截面速度不均匀引入动量通量参数，也在广泛的区域内得到了良态的系统。但是，Issa等[7]发现动量通量参数的预测结果同实验偏离。他们经过同实验的对比发现Holmas的人工扩散方法能够在网格不断加密的情况下得到同实验吻合的预测结果[8]。但是，他们只对水平管的一个工况进行了研究，提出的系数对竖直管计算的作用还不清楚，需要进一步研究。

本文针对目前保持双流体模型双曲性的研究，对双流体模型模拟竖直管内的流动进行探索，寻求获得正确预测结果所需人工扩散项系数，并验证此组系数计算水平管段塞流时的精确度。

1 控制方程

加入人工扩散项后的一维双流体模型如下：

(1)

(2)

(3)

(4)

αl+αg=1，

(5)

式中：下标l和g分别表示液相和气相；下标w和i分别表示在壁面处和相界面处；α表示相含率；ρ表示密度；u表示速度；p表示压力；g表示重力加速度；hl表示液膜厚度；τ表示切应力；θ表示管道倾角；S表示湿周；A表示管道横截面积。气相为理想可压缩气体，液相为不可压缩流体。

式(1)～(4)左边第三项二阶偏导即为加入的人工扩散项，φal和φag分别表示液相和气相连续方程的人工扩散项系数，μal和μag分别表示液相和气相动量方程的人工扩散项系数。剪切应力的表达式如下：

(6)

式中：f为Fanning摩擦系数；两相的相对速度ur=ug-ul。

对于水平管和竖直管，两相的湿周计算分布如图1所示。

图1 湿周及两相分布示意图Fig.1 Schematic representation of wetted perimeters and two-phase distributions

水平管相界面上的摩擦系数采用Taitel-Dukler模型[9]：

(7)

竖直管相界面上的摩擦系数采用Wallis[10]的模型：

(8)

D为管道直径。水平管气相与壁面间的摩擦系数同样采用Taitel-Dukler的模型[9]：

(9)

水平管液相与壁面间的摩擦系数采用Spedding-Hand模型[11]：

(10)

竖直管液相与壁面间的摩擦系数采用Wallis模型[10]：

(11)

雷诺数的表达式为

(12)

2 数值方法

气液相的连续和动量方程式(1)～(4)通过有限体积法在同位网格上进行求解，采用Rhie-Chow插值以消除压力速度的失耦。压力修正方程采用两相混合的质量守恒方程导出，其形式如下：

(13)

式中：D/Dt为所给量的物质导数。

为了提高计算的稳定性，时间离散格式为一阶欧拉隐式后向差分，空间离散格式采用二阶迎风。为了保证计算收敛性，压力与速度的耦合采用SIMPLE算法，并设定全局最大Courant数为0.1，根据Courant数自动选取时间步长，计算收敛的判据为1×10-6。

3 结果及讨论

3.1 Water Faucet问题

Water Faucet问题是检验数值计算结果精确度的常用标准方法，其过程如图2所示。其结构为一段12 m长、内径为1 m的竖直管道，上部入口含气率为0.2，气相和液相的速度分别为0和10 m/s，出口压强固定为105Pa。气相的密度为1.16 kg/m3，液相的密度为1 000 kg/m3，整个系统的温度固定为50 ℃。此问题的解析解为

(14)

.

(15)

当采取不同的扩散项系数和网格控制容积(cv)数量时，模拟得到的t=0.5 s和t=2.5 s时的结果和解析解的结果对比分别如图3和图4所示。从图3可以看出，当未引入人工扩散时，双流体模型是不稳定的病态问题，网格密度为500 cv时，含气率曲线在梯度剧变处发生了大幅振荡。而如图4显示，未引入人工扩散的计算并不能得到稳定的稳态结果。

图3 不同扩散项系数及网格密度在t=0.5 s时模拟与解析解比较Fig.3 Comparison of simulation and analytical results for different diffusion coefficients and mesh densities at t=0.5 s

图4 不同扩散项系数及网格密度在t=2.5 s时模拟与解析解比较Fig.4 Comparison of simulation and analytical results for different diffusion coefficients and mesh densities at t=2.5 s

当扩散系数采用Issa等[8]的系数时，含气率的非物理振荡已被减弱但尚未消失。若进一步增大动量扩散系数，计算没有明显的改善，因此逐步增大质量扩散系数。质量扩散系数从0.1开始逐步增大时的计算结果与解析解的比较如图5所示。从图中可从看到随着质量扩散系数的逐步增大，非物理波动的幅度逐渐减小；当质量扩散系数达到0.7时，计算结果已经不存在非物理振荡。虽然由于扩散的影响，含气率的曲线变得更加平滑，但是能够得到符合物理过程的无振荡解，损失的精度是可以接受的。当网格密度进一步加密到1 500 cv和3 000 cv时，仍旧没有非物理振荡出现，可以认为此时能够得到所求解问题的网格无关解。

图5 不同质量扩散项系数在t=0.75 s时模拟与解析解比较Fig.5 Comparison of simulation and analytical results for different diffusion coefficients at t=0.75 s

3.2 水平管段塞流

Woods等[12]列举的内径为0.095 m水平管内各种表观气速uSG和表观液速下uSL的段塞流液塞频率如图6所示。本文选取了气液相入口表观速度分别为15.9 m/s和1.2 m/s远离稳定边界的病态区域工况，分别使用较粗网格不加入人工扩散和人工扩散系数φ=0.7，μ=1×10-5进行计算。

图6 Woods等[2]列举的液塞频率实验数据Fig.6 Experimental data for slug frequency presented by Woods et al[12]

图7 模拟得到的沿管道无量纲液位高度随时间的变化(φ=μ=0)Fig.7 Time trace of simulated non-dimensional liquid film height along the pipe (φ=μ=0)

在网格密度为1 500 cv时，计算得到的无量纲液位高度(hl/D)的瞬态曲线随时间的变化如图7所示，图中相邻曲线间的时间间隔为0.4 s。计算初始状态管道内无量纲液位高度为均匀的0.5，随着气液间相互作用的发展，约2 s时相界面的波动开始出现并在约3 s时发展为液塞。由于初始条件设置，液塞前面的液膜厚度较大，液塞长度能够不断增长，因此需等第一个长液塞流出后才能反映真实的流动。从后续的发展曲线来看，并未出现此种超长的液塞。在第一个超长液塞流出管道后，统计得到的液塞频率为0.88 Hz，与0.79 Hz的实验结果相比，误差为11.4%。虽然计算结果同实验相比很接近，但是当加密网格后，进一步被解析的短波波动使得计算很快发散，无法得到网格独立解。

采用前文验证的扩散系数在网格密度为3 000 cv的细网格上进行模拟，得到的无量纲液位高度的瞬态曲线随时间的变化如图8所示。从图8可以看出，液面高度的波动被抑制，曲线变得平滑，液塞形成的频率降低，统计得到的液塞频率降低为0.58 Hz，与实验相比减小了26.5%，显然此时人工扩散量是过大的。但是，进一步加密网格未能显著改善预测结果，证明人工扩散量在精细网格上仍旧过大，得到的收敛结果对液塞频率的预测偏低。可以预测当扩散量很大时，双流体模型的性能类似于混合物模型。

图8 模拟得到的沿管道无量纲液位高度随时间的变化(φ=0.7,μ=1×10-5)Fig.8 Time trace of simulated non-dimensional liquid film height along the pipe(φ=0.7,μ=1×10-5)

4 结 语

通过加入人工扩散可以使双流体模型在原本病态的区域内转变为良态问题。经过一系列的数值实验，得出在人工扩散项系数为φ=0.7，μ=1×10-5时，竖直管Water Faucet问题非物理振荡完全消失。对于水平管内流动病态区域的工况，φ=0.7，μ=1×10-5时水动力学不稳定性被抑制，液塞频率降低。水平管人工扩散项的系数需要根据不同工况决定，有待进一步研究。

[1] Lin P Y, Hanratty T J. Prediction of the initiation of slugs with linear stability theory[J]. International Journal of Multiphase Flow, 1986, 12(1): 79.

[2] Issa R I, Kempf M H W. Simulation of slug flow in horizontal and nearly horizontal pipes with the two-fluid model[J]. International Journal of Multiphase Flow, 2003, 29(1): 69.

[3] Zanotti A L, Méndez C G, Nigro N M, et al. A preconditioning mass matrix to avoid the ill-posed two-fluid model[J]. Journal of Applied Mechanics, 2007, 74(4): 732.

[4] Chung M-S, Lee S-J. A modified semi-implicit method for a hyperbolic two-fluid model[J]. Appl Numer Math, 2009, 59(10): 2475.

[5] Holmas H, Sira T, Nordsveen M, et al. Analysis of a 1D incompressible two-fluid model including artificial diffusion[J]. IMA Journal of Applied Mathematics, 2008, 73(4): 651.

[6] Song J H, Ishii M. The well-posedness of incompressible one-dimensional two-fluid model[J]. International Journal of Heat and Mass Transfer, 2000, 43(12): 2221.

[7] Issa R I, Montini M. Applicability of the momentum-flux-parameter closure for the two-fluid model to slug flow[C]. 6th International Symposium on Multiphase Flow, Heat Mass Transfer and Energy Conversion, 2009: 712.

[8] Issa R I, Montini M. The effect of surface tension and diffusion on one-dimensional modelling of slug flow instabilities[C]. 7th International Conference on Multiphase Flow, 2010.

[9] Taitel Y, Dukler A E. A model for predicting flow regime transitions in horizontal and near horizontal gas-liquid flow[J]. AIChE Journal, 1976, 22(1): 47.

[10] Wallis G B. One-Dimensional Two-Phase Flow[M]. New York: McGraw-Hill,1969.

[11] Spedding P L, Hand N P. Prediction in stratified gas-liquid co-current flow in horizontal pipelines[J]. International Journal of Heat and Mass Transfer, 1997, 40(8): 1923.

[12] Woods B, Fan Z, Hanratty T. Frequency and development of slugs in a horizontal pipe at large liquid flows[J]. International Journal of Multiphase Flow, 2006, 32(8): 902.

NumericalAnalysisoftheEffectofArtificialDiffusionTermsonTwo-FluidModel

HAN Peng-fei1, GUO Lie-jin1, LI Qing-ping2, YAO Hai-yuan2, CHENG Bing2

(1.StateKeyLaboratoryofMultiphaseFlowinPowerEngineering,Xi’anJiaotongUniversity,Xi’an,Shaanxi710049,China; 2.CNOOCResearchInstitute,Beijing100028,China)

When artificial diffusion terms are introduced in two-fluid model, the effects of different artificial diffusion coefficients are studied. For vertical pipe Water Faucet problem, when the artificial diffusion coefficients are gradually increased, nonphysical oscillations can be eliminated. For horizontal pipe flow in ill-posed region， the calculated results without diffusion terms agree with experimental data but the mesh independent solution cannot be obtained. When the coefficients calibrated in Water Faucet problem are adopted, hydrodynamic instability is suppressed to some extent and slug frequency is underestimated. Therefore, the amount of artificial diffusion varies with different flow conditions and should be determined specially.

two-fluid model; artificial diffusion; slug flow; numerical simulation

O359

A

2095-7297(2015)01-0023-05

2014-12-29

国家科技重大专项(2011ZX05026-004-02)、国家自然科学基金重点项目(51236007)

韩朋飞(1983—),男，博士研究生，主要从事气液两相流方面的研究。

*通信作者