陈宏举 王靖怡 康琦 吴海浩 宫敬
（1.中国石油大学（北京）北京 102249；
2.中海油研究总院有限责任公司 北京 100028；
3.中国石油天然气管道工程有限公司 河北廊坊 065000；
4.清华大学 北京 100084；
5.深圳清华大学研究院 广东深圳 518057）
油气水多相混输技术能够降低分输管道的建设费用，但多相流动的不稳定性和气液流动的交替性可能导致管道内压力、持液率和气液相流速发生周期性波动［1-2］，对精确预测多相流动现象提出了更高要求。目前国内多相流仿真模拟大多采用多相管流模拟软件OLGA和LedaFlow，二者都是基于流体动力学机理的瞬态模拟工具［3-6］。中国使用多相管流模拟软件一直受国外的制约，每年要缴纳高额的使用授权费用。
为了减少对国外商业软件的依赖，推动多相管流模拟软件国产化发展，中国石油大学（北京）宫敬教授团队自主研发了动态多相管流模拟软件MPF［7-9］，目前处于测试阶段。本文对MPF、OLGA和LedaFlow软件的理论基础进行介绍，以实际工程管线稳态工况为例，比较3个软件在预测压降、温降、持液率和气液相速度等方面的差异；
并研究瞬态工况下3个软件对严重段塞流的处理情况，从而对比分析MPF与成熟商业软件OLGA和LedaFlow在稳态、瞬态工况方面的预测能力，为国产化软件MPF的进一步完善提供参考。
1.1 OLGA
OLGA采用基于6个守恒方程进行求解的拓展双流体模型，包括气相、液相或液膜、液滴3个质量守恒方程，液滴的气相、液相或液膜2个动量方程及1个混合体系能量方程［10］。该模型采用欧拉法和欧拉形式的拉格朗日边界跟踪形式，通过人工假定条件或依据实验关联式反映本构关系。质量守恒方程中考虑了控制液滴夹带和沉积的质量传递项，由此得出气相、液相或液膜、液滴的质量守恒方程，计算气相、液相之间的质量分布情况。此外，结合流体流动状态的影响及气液两相之间因相变引起的质量交换，计算出多相流各项参数［11］。拓展的双流体模型包括质量守恒方程、动量守恒方程、能量守恒方程［4］。
1.2 LedaFlow
LedaFlow的核心计算模型可分别求解油气水三相的质量、动量及能量的守恒方程［5］。通过研究9个独立物理场的特点，建立每个独立的物理场所对应的质量方程以及每个连续区域的多混合模型所对应的动量守恒方程，LedaFlow内核模型会求解每一相的能量守恒方程，计算同一截面上三相的温度，故可以获得更多更准确的关于流体温度的信息。LedaFlow的每个计算单元拥有16个方程，包括9个独立场的质量守恒方程，连续相的3个动量方程，连续相的3个能量方程及1个体积守恒方程［5］。与OLGA相比，LedaFlow拥有更多的质量守恒方程，包含3个连续相方程，6个分散相方程，以及针对每个场的质量守恒方程、针对连续相中每一相的动量方程和针对连续相中每一相的能量方程。
关于动量方程的处理，LedaFlow认为在多相流混输管道中移动介质的界面张力相对较小，因此可以忽略不同界面的压差。采用交错网格贡献边际法进行求解，把速度分配给网格边界，把其余参数分配到网格内部。方程组内函数解析使用逆向差分法，从而提高计算稳定性。通过修正后的SIMPLE算法进行求解［12-13］。
1.3 MPF
MPF采用双流体模型，以变量在管道截面上的平均值为基础，分别将两相流视为连续介质而建立控制方程，分别描述各相的质量、动量、能量等性质［14］。MPF采用2个质量方程、2个动量方程、2个能量方程及1个混合能量方程组成的复杂7方程模型来研究两相流动［15］。
各相质量方程：
式（1）～（5）中：k为g、l时分别表示气相与液相；
αk表示截面相含率；
ρk表示密度；
uk表示速度；
t表示时间；
x表示沿线管道的坐标；
˙mk表示质量源变化率；
p表示压力；
Fk表示除压力外的其他表面力与质量力的合力；
Ek表示总能；
Hk表示焓；
Wk表示能量源（包括外力做功和流体与外界换热）。
关于能量方程的处理，MPF采用各相能量方程相加的方式来计算混合温度，保证了流体总焓守恒。与OLGA、LedaFlow相比，MPF对摩阻系数的选取（表1）更适用于中高持液率工况。
表1 MPF气-液两相流动摩阻封闭关系式Table 1 MPF closures for friction in gas-liquid flow
2.1 管线运行基础参数
某实际气田新建一条管道输送天然气至处理厂，输气规模64.8×104m3／d，设计压力4.0 MPa，出站压力3.5 MPa，进站压力2.5 MPa，出站温度45℃，进站温度3℃，总传热系数5.57 W／（m2·℃），全长约20.96 km。管道沿线纵断面如图1所示，流体组成如表2所示。
图1 某气田湿气管线纵断面图Fig.1 Profile diagram of a wet gas pipeline in a gas field
表2 某气田湿气管线流体组成数据Table 2 Fluid composition data of a wet gas pipeline in a gas field
2.2 3个软件稳态模拟计算结果对比
2.2.1 压降
3个软件水力计算的内嵌模型较为相似，且都是给定终点压力、起点温度和起点输量反算出起点压力。为了验证软件计算结果的准确性，筛选出湿气管线的10组运行数据（表3），模拟得到不同工况下3个软件起点压力及压降计算结果（图2）：①设计工况10中，OLGA和LedaFlow起点压力计算值稍低于实际设计工况的起点压力，相对误差约为-3%；
MPF起点压力计算值稍高于实际设计工况的起点压力，相对误差为3%，3个软件计算结果均在误差允许范围内。②工况1～9，3个软件的起点压力计算值均高于实际运行的起点压力，但3个软件之间的计算结果差距很小；
OLGA与Leda-Flow起点压力和压降的计算结果更接近，MPF起点压力计算结果误差高于另外2个软件。结合设计工况10的计算结果来看，可能是现场气质组分与给定设计工况下的组分存在差异，现场输送的原料气重组分占比更低，管道持液率偏小，因此压降小；
所以3个软件反算出的起点压力均大于实际运行情况。③与OLGA、LedaFlow相比，MPF的摩阻计算模型为BeggsBrill，水力算法为BeggsBrill-Modified，对摩阻系数的选取更适用于中高持液率工况，MPF压力计算值偏大可能与软件内置水力模型和摩阻系数的选取不同有关。但整体上，3个软件的压降计算偏差均在可接受范围内，可用于工程模拟计算。
图2 某气田湿气管线不同工况下3个软件起点压力及压降计算结果Fig.2 Starting point pressure and pressure drop of a wet gas pipeline in a gas field under different conditions caculated by three software
表3 某气田湿气管线运行情况Table 3 Operating conditions of a wet gas pipeline in a gas field
2.2.2 温降
基于现场运行数据，得到不同工况下3个软件温降计算结果（图3）：MPF温降变化与LedaFlow较为吻合，OLGA温降值偏低，但3个软件计算结果偏差不大。MPF热力算法为Deng，传热模型为K＿U，采用3个能量方程分别计算各相的温度和混合温度。LedaFlow能够计算同一截面上三相的温度，从而达到求解每一相的能量守恒方程的要求。可见，MPF和LedaFlow对于传热模型的处理更加完善，认为其具备更细致的传热计算功能。
图3 某气田湿气管线不同工况下3个软件温降计算结果Fig.3 Temperature drop of a wet gas pipeline in a gas field under different conditions caculated by three software
2.2.3 其他参数
以设计工况为例，比较稳态运行下3个软件的全线持液率计算结果（图4）：3个软件的持液率预测存在微小差异，OLGA在管线低点存在少量积液，LedaFlow持液率保持在0.001左右，MPF全线基本没有积液产生。
图4 某气田湿气管线设计工况下3个软件全线持液率计算结果Fig.4 Holdup of a wet gas pipeline in a gas field under design condition caculated by three software
3个软件对全线气液相速度的计算结果（图5）表明：MPF计算出的液相速度振荡幅值与气相速度基本重合，相较OLGA和LedaFlow的计算结果，MPF计算出的液相速度振荡幅值范围较大，但整体速度的变化趋势与另2个软件的结果相同。从水力段塞形成的角度分析，在水力段塞形成过程中，部分时刻会出现等效于均相的流动效果，段塞流动即液塞、分层流间歇出现，段塞体中会携带一些微小气泡，其速度与段塞体基本一致，但OLGA与LedaFlow假设段塞内部为纯液相，因此图5a、b中气液相流速变化波动较小，且相差较大，而MPF充分考虑段塞体内气相流动的影响，更加贴合实际。
图5 某气田湿气管线设计工况下3个软件预测的全线气液相速度Fig.5 Gas-liquid velocity of a wet gas pipeline in a gas field under design condition caculated by three software
3.1 卧底-悬链线立管运行基础参数
以国外某深海油田一条实际运行的卧底-悬链线立管为例［16］对3款软件的瞬态计算功能进行对比。立管高度约1 720 m，入口温度90℃，出口压力1.5 MPa，气液比50 Sm3／Sm3；
管径203.2 mm，壁厚0.009 m；
保温材料为聚丙烯，厚度0.03 m。对该管线模型进行简化，得到其高程和环境温度变化情况（图6）。
图6 某卧底-悬链线立管系统高程和环境温度变化Fig.6 Variation of elevation and ambient temperature of a undercover-catenary riser system
3.2 3个软件瞬态模拟计算结果对比
3.2.1 沿线持液率
3个软件均采用黑油模型，模拟得到该系统在同一时刻不同流量对应的沿线持液率变化（图7）：3个软件的沿线持液率变化总体趋势较为一致；
管线下倾管内液体积聚，持液率迅速升高直至为1，上倾立管内持液率基本稳定在1左右，出口处持液率迅速下降；
随着入口流量的增大，沿线持液率走向更加稳定。
图7 3个软件预测的某卧底-悬链线立管系统沿线持液率Fig.7 Holdup of a undercover-catenary riser system caculated by three software
3.2.2 立管底部压力
3个软件计算得到的卧底-悬链线立管系统立管底部压力（图8）均表现出良好的周期性变化，随着入口流量的增加，立管底部压力的波动周期逐渐缩短，符合真实物理规律［17-18］。MPF预测立管底部压力的变化呈现出周期性波动趋势，符合有回落的严重段塞流循环，具有稳定的平台期。LedaFlow预测结果为无回落的严重段塞流循环，与有回落循环类似，但液塞喷发过程表现不是很剧烈，有回落循环过程的液塞长度会高于无回落循环。OLGA预测结果为有回落循环过渡至无回落循环，立管底部压力值也更接近LedaFlow。由于各软件的严重段塞流预测模型存在差异，判别下倾管内气体进入立管的划分方式不同，导致三者的预测结果存在差异。MPF立管底部最大压力值高于另外2个软件。这是因为：MPF严重段塞流模型在液体回落过程中进行了简化处理，MPF预测结果均为有回落循环，模型默认气液喷发后的回落阶段仅有液相，而立管中气相全部喷出，在下一个液塞积累过程，认定立管内液位以下部分为纯液相，而实际管内流动非常复杂。邱伟伟等［19-20］通过实验观察也证实了在气液喷发过程中多相流动的复杂性，液体回落阶段立管内也并非纯液体。因此，MPF对于立管底部最大压力的计算值高于OLGA和LedaFlow的计算结果，且具有稳定的平台期。整体上看，在严重段塞流的模拟中，MPF的表现稍逊于OLGA和LedaFlow。
图8 3个软件预测的某卧底-悬链线立管系统立管底部压力Fig.8 Pressure at the bottom of the riser in a undercover-catenary riser system caculated by three software
3.2.3 段塞周期和液塞长度
3个软件模拟得到不同流量下卧底-悬链线立管系统段塞周期和液塞长度的变化情况（表4）：①3个软件的液塞长度计算结果均超过了立管的高度，这是严重段塞流形成的显著标志；
②随着入口流量的增加，液塞长度逐渐缩短，段塞周期也不断缩短。
表4 3个软件预测的某卧底-悬链线系统液塞长度和段塞周期Table 4 Liquid plug length and slug fluctuation period of a undercover-catenary riser system caculated by three software
以OLGA的计算结果为基准，对比不同流量下LedaFlow和MPF计算结果的偏差（表5）。可见：①MPF的段塞周期与OLGA的段塞周期偏差较小；
②随着入口流量的增大，LedaFlow的段塞周期变化幅度较为保守，不同流量下均为无回落的循环过程，因此最大液塞长度小于另2个软件的计算结果。
表5 LedaFlow和MPF在某卧底-悬链线系统液塞长度和段塞周期的计算偏差（以OLGA为基准）Table 5 Error calculation of liquid plug length and slug fluctuation period for a undercover-catenary riser system caculated by LedaFlow and MPF software（based OLGA）
3个软件对液塞长度的预测存在差异。MPF对于最大液塞长度的计算值高于LedaFlow和OLGA，主要原因在于：液体积累阶段立管内充满液体，最大液塞长度的计算差异主要是下倾管内液塞长度的不同导致的；
MPF立管底部压力计算值高于另2个软件，因此基于公式计算下倾管内液塞长度相应也会偏大。此外，由于MPF对回落过程进行简化处理，默认立管内全部为液相，但事实上油气水三相流体在管内的流动过程十分复杂，实际液塞积累不一定全部充满液相直至气液喷发，且严重段塞流又会引发强烈的周期性波动，从而使MPF对回落过程的简化处理导致了最大液塞长度的计算结果高于OLGA和LedaFlow的计算结果。因此，在后续开发中建议进一步研究完善MPF内嵌段塞流的计算模型。
3.2.4 立管底部气相速度
3个软件模拟得到的不同流量下卧底-悬链线立管系统立管底部气相速度随时间的变化规律（图9）：随着入口流量的增加，立管底部气相速度总体呈现上升趋势，气相速度的波动周期在减小。在液塞积累过程中，OLGA和LedaFlow的气相速度计算结果并不完全为0，说明两款软件在计算严重段塞流时，上升立管内并非纯液相，这2个软件在计算液塞积累时将立管内气体流动情况考虑在内。MPF软件在此方面需要进一步考虑液塞形成中气速的精确描述。
图9 3个软件预测的某卧底-悬链线系统立管底部气相速度Fig.9 Gas velocity at the bottom of the riser in a undercover-catenary riser system caculated by three software
MPF与OLGA、LedaFlow在计算管线的压降、温降、沿线持液率、立管底部压力、段塞周期与液塞长度及立管底部气相速度等方面存在差异，表现在：①稳态工况下，3个软件水力计算与现场实际数据相对误差绝对值在3%以内，MPF和LedaFlow热力计算模型能够精确到每一相，二者计算结果更为接近；
气液相速度方面，MPF对水力段塞流形成过程的处理更符合实际物理过程。②严重段塞流瞬态工况下，3个软件的持液率波动位置及频率变化十分吻合，OLGA与MPF的波动周期相差很小，呈现了良好的预测结果，但由于MPF对严重段塞流模型进行了简化处理，故其对立管底部压力、最大液塞长度的计算结果均高于OLGA和LedaFlow的计算结果，立管底部气相速度波动预测不够准确。建议在后续研究中，进一步完善MPF中的段塞流等模型，并结合更丰富、完整的实验数据和现场运行数据进行测试更新，全方位提升MPF的功能及可靠性，从而使其能准确应用于实际工程的各类工况计算。