ABSTRACT

Simulation program based on the theory of drift-flux and hydraulic model was finished in consideration of the influence of pipe angle on distribution parameter C₀ and drift velocity V_gj within thermodynamic system two-phase flow. A research was carried out using a simplified hydraulic model in which the void fraction was calculated by iteration; according to the computing result of both adiabatic and uniform heated two-phase pipelines, the comparison was given, and the details showed that theoretical analysis was well agreed with experimental results and provided a reference for engineering application.

Keywords:Pipeline, Two-Phase Flow, Drift-Flux, Voidfraction

一种基于漂移模型的管内气液两相流分析方法

吴  巍，唐胜利

重庆大学动力工程学院，重庆

Email: 20111002075@cqu.edu.cn, tytsl@cqu.edu.cn

收稿日期：2014年5月7日；修回日期：2014年5月23日；录用日期：2014年6月5日

摘  要

针对热力系统管道内气液两相流动，在漂移通量理论和水力模型基础上，综合考虑管道倾角对空泡份额分布参数C₀及漂移速度V_gj的影响，提出了一种迭代求解空泡份额的简化水力模型，该模拟的计算结果与绝热管及均匀加热管内两相流动试验数据相比有很好的吻合度，证明了该模型的可用性，为工程实际应用提供了参考。

关键词

管道，两相流，漂移通量，空泡份额

1. 引言

目前气液两相流的研究和应用大多为管内流动，近十年来随着社会生活及工业生产对石油及天然气的依赖，越来越多的国内学者将精力集中在如何解决与油气运输相关的两相流问题，而对热力系统中的汽水管道两相流动问题研究却减少，这在很大程度上是由于该领域已经有了较多的研究成果并基本能够满足一些实际的工程问题需要。但在实际的生产中仍然存在较多需要解决的汽水两相流动问题，如在常规火电厂中经常会遇到锅炉水冷壁壁内水汽流动传热特性、加热器疏水管等管道设计，核电站二回路蒸汽系统基本都处于饱和湿蒸汽区，其工作过程及相关参数同样需要使用两相流的分析方法进行描述和计算。大多数情况下，当气相含量低于5%或大于96%时则按照单相流进行计算。

漂移通量模型(DFM)介于分相模型与均相模型之间，既考虑了各相平均速度之间的差异，克服了均相模型的缺点，同时也避免了分相模型中需要考虑相界面之间复杂关系所增加的计算量。本文结合目前国内外对汽水管道的研究成果，在绝热管道内两相流动计算模型的基础上，提出有热量输入的均匀受热管道气液两相流计算方法，并采用了迭代的方式来求解空泡份额，将计算值与文献中试验值对比，两者有较好吻合，可以为工程实际应用提供参考。

2. 理论模型

解决两相流问题的漂移通量模型是由Zuber和Findlay[1] 于1965年首次提出，他们认为，在分析两相流截面含气率时，应当同时考虑气液相之间的滑动以及各自的密度变化，在模型中建立截面平均以及权重平均值的概念。其控制方程中使用平均速度代替了分相模型中需要求解的两相速度，减少了未知变量个数，与此同时在离散相质量方程及混合动量方程中包含了能够反映气液两相之间滑移情况的漂移速度，避免了二流体模型的繁难复杂并能够较好反映管内各相的流动情况。

2.1. 两流体模型

两流体模型通用方程[2] 如下：

(1)

质量、动量和能量平衡式中各物理量如表1所示：

若求解(1)则还需要其他本构方程即：相与壁面间的摩擦与传热方程以及相界面上的质量、动量、能量转换关系。而到目前为止，如何准确确定两相流相界面的物理过程以及几何特性变化仍然是一个需要深入研究的问题，同时计算该模型也将导致较为繁复的计算。

2.2. 简化漂移模型

漂移模型实际上是建立在二流体模型的基础上，结合大量的实验研究结论以解决工程问题。由上述

通用方程以及漂移通量模型[3] 可得到本文简化模型控制方程组：

连续性方程

(2)

动量方程

(3)

模型中变量定义：

(4)

(5)

(6)

(7)

由以上各式可以得到：

(8)

由漂移速度可以计算得到气液两相的速度分别为：

(9)

(10)

其中：—管道倾角˚

—截面含气率

—加权平均速度m/s

—截面压强Pa

—混合工质密度kg/m³

,—气相、液相密度kg/m³

—平均漂移速度m/s

—摩擦压降Pa

为了求解上述方程组，还需要知道与气液相间滑移相关的漂移参数以及工质的物性函数表达式，国内文献中大多将漂移模型参数即分散系数与漂移速度记作干度与压力的函数，而在我国电站锅炉水动力计算中亦是将两相滑动比表达为以压力、平均速度与容积含气率作为自变量的关系式，虽考虑到管道倾斜角因素，但计算较为繁琐，需要分情况讨论计算。近几十年来国内外学者通过不断试验、总结，得出了较多经验关系式，其中以Ishii等为代表，在近段时间发表的文章中，Swanand M. Bhagwat[4] 等人统计了60多种不同试验条件下的8255个试验数据，针对不同工质将不同漂移参数表达式计算结果作了对比，发现考虑了管道倾角与内径等因素的函数关系式能够更准确地计算漂移参数，进而得到精度较高的截面空泡份额值。本文采用上述文献中所提出的漂移系数表达式，并采用Smith[5] 提出的关系式作为求解模型空泡份额的迭代初始值。

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

在确定系数时用到的拉普拉斯变量表示无量纲管道内径的倒数，由Kataoka and Ishii (1987)首先提出并由式计算，其中分别、表示液相表面张力以及两相密度差。另外，在式(17)中只有当其后条件同时满足才取相应的值，其中表示基于气相折算速度定义的Froude数，表达为式(18)。为了在迭代计算时能够提高收敛速度，模型中采用式(19)作为迭代初值，通过运算发现该方法能够较准确快速地确定空泡份额初始值。

3. 模型求解

3.1. 物性参数计算

水和水蒸汽作为典型的常规工质被广泛地应用于动力工程及其相关学科领域，本文依据1997年蒸汽性质国际协会公布的IAPWS-IF97热力学性质新标准[6] ，利用Fortran编制了物性计算程序，实现水和蒸汽性质的计算机求解。

程序采用模块化方式编制，对每个区域的自由能函数表达式中需要使用的常数数组用module进行封装，而将各物性参数的计算式编制为子函数，在主程序计算中通过调用的方式求解所需参数，分区情况如图1所示。

3.2. 摩擦阻力系数

在确定摩擦阻力系数时，通常采用尼古拉兹粗糙管经验公式，而工程实践中使用莫迪图更符合实际工业管道特性。通过计算发现尼古拉兹的部分判定公式计算值与实际取值之间存在较大误差，如光滑管与粗糙管过渡区之间使用的Colebrook计算公式，该公式为隐式表达，在计算摩擦系数时需要通过迭代求解，而部分简化式与实际取值偏差较大。综上所述，本文对经验公式[7] 进行适当的修正，通过线性插值法求取过渡区域摩阻系数。

(20)

(21)

通过对比常用经验式计算值与莫迪图各曲线取值，在公式(20)基础上提出摩阻系数修正计算式(21)，式中为常数修正项，对图中未给定曲线采取插值方式进行计算。以相对粗糙度为例，简化公式计算值与实际值之间存在较大误差，并且平均相对误差约为19%，而对于经过修正之后的计算式，其计算结果与实际取值之间的最大相对误差不超过2%，相关系数。另外，光滑管区使用布拉休斯公式以及阻力平方区采用尼古拉兹粗糙管公式均能达到较高精度，对于过渡区与阻力平方区的分界线，本文采用皮勾(B.J.S.Pigott)推荐式。

图1. IAPWS-IF97公式及分区

3.3.求解原理及步骤

本文研究稳态情况下管内流动，因此略去公式(2)、(3)中时间项，并使用Fortran编制了计算程序。计算中以管内两相干度即热平衡含气率沿管道长度方向的变化量作为反映热量交换的参数，另外计入管内压降引起的工质参数变化，通过分段积分求解管内气液两相的流动参数。由于存在流动压降，因而管段节点与其相邻节点的物性参数有所不同，文中假设各节点间参数为线性变化由于绝热管道内因压降而产生的新汽量很少，故管内干度变化可以忽略不计，通过给定干度沿管长的变化来反映热量输入，对于均匀加热管道则可认为沿管长方向干度为线性变化。

根据上述原理，在对每个管段建立守恒式后，由初始的工况条件得到工质物性参数，同时根据式(19)计算出对应的作为管道入口边界条件，利用考虑管道倾角等参数影响的漂移参数计算式可以迭代求解到较准确的初始空泡份额。在求解第节点处离散方程组时，已经得到第节点的参数、、、、和，此时假设代表管段节点处压强值为，在程序中以压差形式给定即，由性质程序计算得到在条件下的工质物性参数，结合此处热平衡含气率 (若为绝热管道则)，可以计算出各参数、，、、然后将其带入动量离散方程中得到该管段节点处计算压强值,通过比较与之间的差异进行迭代修正直到满足设定的精度。需要说明的是，在迭代计算某节点处的空泡份额时，气相漂移速度中的空泡份额使用的是前一节点处值，通过式(11)~(17)对其进行修正，迭代至满足精度要求时得到该节点参数值。程序计算流程如图2所示。

4. 计算结果与试验值比较

本文对比所用实测参数均取自文献[8] [9] ，针对非加热管以及均匀加热管内两相流动情况作了详细计算，并对计算结果与实测值进行了对比及分析。另外两条对比曲线为锅炉水动力计算方法以及均相模型计算方法结果。

图2. 程序计算流程图

4.1. 非加热管内流动参数计算值与实测值对比

图3所示为15.1 Mpa时，质量流速在370 kg/m²∙s~1900 kg/m²∙s范围条件下，垂直上升绝热管中两相流空泡份额分布情况。图4所示为3 Mpa时，质量流速在552.98 kg/m²∙s~656.71 kg/m²∙s条件下，垂直下降绝热管中汽水混合物的空泡份数值额模拟结果。

由图中曲线可看出，垂直上升和下降绝热管管内两相流动的计算结果能够较好地吻合试验数据。具体地，垂直上升流动时，均相模型中空泡份额等于容积含气率，因此均相模型计算结果始终大于采用滑动比模型的锅炉水动力方法计算结果，本文模型计算得到的两相速度差值略大于滑动比模型结果，因此最终得到的值也略小于后者。垂直下降流动时，低含气率区域中干度的微小增量导致空泡份额的急剧上升，如图4所示，由于气相流量远小于液相流量，气相以泡状形式在管内流动，两相之间作用力较大，故此时漂移速度较小甚至可能出现负值，因此使得在低含气率区域空泡份额值接近容积含气率值，在高含气率区则与上升管内的流动情况一致。

图5所示为管内工作压力在18.1 Mpa左右时，不同流量条件下压力降随干度的变化情况。本文在计算摩擦压降时，摩擦系数处理为均相条件下的取值，因此在流型由泡状流转变为弹状或环状流时，会出

图3. 垂直上升绝热管内空泡份额分布

图4. 垂直下降绝热管内空泡份额分布

现一定的偏差，在合理的范围之内。

4.2. 均匀加热管内流动参数计算值与实测值对比

如图6所示，工作压力为18.1 Mpa时，均匀加热时垂直上升管内汽水两相流空泡份额沿管长方向在不同干度条件下的实际分布情况与各模型的计算结果对比，试验数据流量范围为378 kg/m²∙s~1042 kg/m²∙s，通过最小二乘法分别拟合得到在不同流量条件下的各模型计算值曲线。可以看出，均匀加热管内空泡份额的分布与前述绝热管内情况类似，即整体的变化趋势相同，但由于压力条件的区别而使得在相同干度情况时，因此使得计算得到的空泡份额值与绝热管内有所差异。综上所述，在输入恒定均匀热量条件下，本文模型计算结果同样能够较好吻合实测数据。

图7所示为管内压力分布计算结果对比，管内两相流量范围为378 kg/m²∙s~1892 kg/m²∙s，对图中数据偏差较大的点，由于加热使流体的密度减小，对加速压降有影响，同时加热还使壁面产生气泡，并与流体进行动量交换，从而相对减小了摩擦压降，因此其误差在合理范围之内。

图5. 垂直绝热管内压降计算结果

图6. 均匀加热上升管内空泡份额分布

4.3. 综合对比

图8所示为压降误差分析结果，本文模型计算得到的管内压力分布能够较好吻合实验数据，与实测值的平均相对误差集中5%左右，最大相对误差为17.7%，对于锅炉水动力标准计算结果，其相对误差集中在5%~10%之间，但其中有较多数据点与实际值偏差较大，最大相对误差为25.2%；另外，工程中常用到的基于分相模型的Friedel方法计算结果在高质量流量条件下偏离实际情况较大，并随着含气率的增加而增大。

图9所示为空泡份额误差分析结果，本文计算值与滑速比模型计算值基本集中在5%的误差范围内，而由均相模型得到的空泡份额值与实测值的相对误差主要集中在5%~10%，但当空泡份额处于0.3至0.7之间时，其结果偏离实际值较大，部分数据处于10%误差范围之外，出现这种偏差的原因主要是由于有热量输入的两相流管内，沿管长方向流体的流型处于时刻变化中，加热管中随着截面含气率的增加，弹状流及环状流为主要流型，因此实际流动过程偏离均相模型程度则增大。

图7. 均匀加热垂直管内压降计算结果

图8. 管内压降计算误差对比结果

图9. 管内空泡份额计算误差对比结果

5. 结论

本文结合已有的研究成果，比较了简化漂移模型预测结果与文献试验结果间的异同，由各对比图可看出，对管内汽水两相流动的空泡份额及压力分布进行了较准确地预测，能够较好地吻合试验数据。具体的结论如下：

1) 本文模拟结果与滑动比模型及均相模型计算结果之间有微小差异且精度稍高于后两者；另外，由于文中模型在计算空泡份额时采用的是考虑了管道倾角影响的统一表达式，而滑动比模型在计算截面含气率时，随着管道倾斜角度的改变需要选择不同的计算式，而且在非垂直情况时还需对公式进行修正，计算相应的修正系数，因此本文对截面含气率的计算在表达式使用的方便性上也较优于滑动比模型。

2) 无论在垂直上升管还是下降管内，相同参数条件下，压力降预测值与其他模型计算结果同实测数据相比集中在较小误差范围，通过分析可知两相流型为导致误差产生的主要因素。由目前工程中常使用的几种计算方法结果可以看出，本文计算方法能够较为方便地预测管内流动情况，同时能够保证在较高的精度范围内吻合实际情况，为工程应用提供了参考。

致  谢

作者诚挚地感谢审稿专家给予本文宝贵的修改意见和建议以及期刊编辑对于本文的付出。

参考文献 (References)