液力透平内气液两相流动的数值模拟

毕智高 刘甜甜 谢 勰 相玉琳 王金玺

（1.榆林学院化学与化工学院；2.陕西省低变质煤洁净利用重点实验室；3.长庆油田分公司规划计划处）

泵作液力透平（PAT）回收液体富余压能因具有结构简单、安装方便、占地面积小、可批量生产及运行维护费用低等优点［1，2］，而广泛应用于油气储运、石油化工及煤化工等行业。 这些工业流程中作为能量回收的工质往往呈气、 液两相状态，入口含气率（IGVF）对液力透平的性能有着直接的影响［3］。

近年来，关于旋转可逆式水力机械中气液两相流动的研究主要集中于泵［4～12］工况，对PAT工况的研究多以纯液为工质，而气液两相条件下的研究报道相对较少。 文献［13～16］采用CFD软件分别分析了IGVF对PAT的外特性、内流规律和所受径向力的影响，文献［17］通过瞬态数值模拟，获得了气体在PAT内的分布规律和流场的压力脉动变化规律，但是还不够深入。 笔者在借鉴泵内气液两相流研究成果的基础上， 考虑入流流型，基于双流体模型， 以空气和水作能量回收工质，对一固定径向导叶式PAT进行定常数值计算， 旨在为气液两相条件下PAT的流动特性及其规律提供部分理论参考。

1 数值方法

1.1 Euler-Euler双流体模型

Euler-Euler多相流模型可分为均相流和非均相流模型。 均相流模型是基于均质平衡流理论，将气液两相流视为均匀混合物，不计相间速度滑移，并且忽略了客观存在的两相界面、相间相互作用等许多重要因素［18］，导致求解误差较大；非均相流模型将气液两相之间的界面视作移动的边界，考虑相间速度滑移、质量与动量传递等，每相流体都有各自的流场并通过相间传递单元进行传递， 即每相都有各自的速度场和温度场，但共用压力场，通过相间作用力和热量传递使得两相速度和温度得到平衡，因而较符合实际。

笔者采用非均相流模型，液相为不可压连续相，气相为不可压分散相，相传递单元采用Particle模型，动量传递方式采用Schiller Nauman模型，假定气液两相入流为泡状流型条件。

1.2 控制方程

连续性方程：

动量方程：

式中 fk——与叶轮旋转有关的质量力；

k——代表液相（l）或者气相（g）；

Mk——相间作用力；

pk——k相压力；

wk——k相相对（滑移）速度；

ρk——k相密度；

βk——k相体积分数；

μk——k相动力粘度。

气、液两相和两相混合工质满足：

式中 Qg、Ql、Q——气相、 液相和两相混合工质的体积流量，m3/s；

βl——体积含液率，%；

βg——体积含气率（气相浓度），%。

1.3 相间作用力

式中 CD——无量纲曳力系数；

dB——气泡直径；

ReB——气泡雷诺数；

wR——滑移速度；

νl——液相运动粘度。

虚拟质量力Mk，VM是由两相间的相对加速引起的，按文献［20］中的公式计算如下：

式中 aVM——虚拟质量加速度；

CVM——虚拟质量系数， 通常为气相含量的函数，球形气泡取0.5。

2 计算模型与边界条件

2.1 计算模型及网格划分

以一径向导叶式PAT为研究对象， 工况下的性能参数和主要结构参数依据文献［21］。 采用对复杂几何边界适应性强的非结构四面体网格对计算域进行网格划分，并对网格进行了无关性验证， 当水力效率波动小于0.45%时即满足网格无关性假设。 考虑到计算机的运算能力、计算精度和计算要求，最终确定网格数量在120万，且网格质量均在0.3以上。

2.2 假设条件

气相为粒径均匀的球形，其直径远小于流道特征尺寸；不考虑两相重力；忽略系统内部、系统与外界间的质/热量传输， 系统内无化学反应发生，温度恒定。

2.3 边界条件及数值求解

液相采用能更好处理高应变率和流线弯曲程度较大流动的RNG k-ε湍流模型， 中等湍流强度，气相采用湍流零方程模型；光滑壁面，固壁对液相为无滑移作用、气相为自由滑移，近壁区的壁面函数采用scalable； 热量传递采用等温模型，设定温度为25℃； 设定质量流量进口，IGVF分别为5%、15%、25%、40%，压力出口；叶轮流体域设置为旋转，其余计算域设置为静止，导叶-叶轮及叶轮-出水管间的动静交界面设置为Frozen Rotor模式。 控制方程的离散采用基于有限元的有限体积法，对流项和湍动能相均采用二阶迎风格式求解， 计算迭代步数设置为2 000， 求解残差类型RMS设为10-5。

3 结果分析

3.1 试验验证

图1为PAT试验台，通过试验测试来验证数值计算的准确性。

图1 PAT试验台示意图

用一台比转数ns=84.5的单级单吸离心泵作PAT在IGVF为20%时， 性能参数数值计算与试验测量的结果对比如图2所示。 由图2可见，两者虽有误差（误差的主要因素是数值计算过程中忽略了前、后腔内部的流体），但在最高效率时的数值计算与试验测量的结果误差小于5%，可见笔者所采用的数值计算方法可以较好地对PAT性能进行预测［17］。

图2 泵性能参数数值计算与试验测量的结果对比

3.2 外特性曲线

气液两相条件下，PAT的外特性参数透平压头H（m）、透平效率η（%）和混合工质流动密度ρf（kg/m3）的计算公式分别为：

式中 g——重力加速度，m/s2；

G——混合工质的质量流量，kg/s；

M——叶轮扭矩，N·m；

pin、pout——进、出口总压，Pa；

ρg、ρl——气、液相密度；

ω——叶轮旋转角速度，s-1。

图3为不同IGVF下PAT的外特性曲线，表1列出了PAT最高效率时的性能参数。

结合图3、表1可知，不同IGVF下PAT的流量-效率、流量-压头和流量-功率曲线的总体趋势基本一致，即随着两相混合工质流量的增加，PAT的效率先较快增加后缓慢减小，压头和功率随着流量的增加而增加。 随着IGVF的增加，相同流量下PAT的效率和对外输出的功率降低，压头升高，小流量下的效率降低和压头升高较大流量下的明显，而ρl远大于ρg导致了大流量下的功率降低较小流量下的明显。PAT最高效率时的流量发生偏移，与IGVF为5%时相比，最高效率分别下降了1.2%、2.2%、4.5%。

图3 不同IGVF下PAT外特性曲线

表1 不同IGVF下PAT最高效率时的性能参数

3.3 蜗壳流道气相分布

图4为蜗壳流道各截面位置，图5为不同IGVF下PAT最高效率时蜗壳流道各截面气相分布情况。 由图5可见，同一IGVF下，工质从进口沿蜗壳流道（从截面1～8）气相分布和气相浓度均匀程度均降低。

图4 蜗壳流道各截面位置

图5 蜗壳流道各截面气相分布

截面1靠近PAT进口端，气相均匀分布，数值接近初始IGVF， 与泡状流型的入流假定条件相符，经过截面1后，气相向流道半径较小的蜗壳截面中下方汇聚堆积，形成高含气区，在此可能发生流型转换和相态分离，这是因为两相混合工质沿蜗壳流道在做减压增速运动，而液相的密度较大，所受到的惯性离心力大于气相，导致大量液相裹挟少量气泡绕蜗壳外缘壁面运动，而在液相的排挤作用下，大量气泡被迫向蜗壳内缘出口处转移，从而形成高含气区。

另外，截面8处气相浓度最低，且气相均布于整个截面；同一截面处的气相浓度随IGVF的升高而增大。

3.4 导叶叶轮内气相分布

图6给出了不同IGVF下PAT最高效率时导叶叶轮内截面上的气相分布情况。 由图6可见，PAT各过流部件中蜗壳内气相分布相对最均匀，各过流部件的气相浓度均随着IGVF的升高而增大，气相分布的不均匀度增加，流动紊乱度增强，这也是随IGVF的增加PAT效率下降的原因。 汇聚堆积在蜗壳截面下方的气相因惯性随液流顺势就近进入附近的导叶和叶轮流道，因此蜗壳流道内高含气率截面对应中心角所包含范围及其附近区域的气相浓度也相对较高，高低含气区域呈现出明显的分界（图6b）。 高含气率叶轮流道内气泡聚向叶片吸力面，且随着IGVF的增加，气相聚集程度增强，气泡聚并成为大的气团而滞留于叶轮流道，改变了有效过流面积，甚至可能造成“气堵”。 这是由于在有限叶片数的叶轮旋转效应和叶形曲率的影响下，液相在进入叶轮后受到较大的惯性离心力和科氏力作用改变了运动状态而倾向叶片压力面运动，而气相受到的惯性离心力和科氏力则相对较小，但在叶片压力面和吸力面之间压力梯度的驱动和液相的排挤下，导致了叶轮流道内的气相向叶片吸力面聚集。 对比图6c中的液相速度流线，发现在滞留气团的上游，存在对应漩涡，表明叶轮流道内漩涡形成与气相的聚并相关。 因此，在对含气工况下应用的固定径向导叶式PAT进行优化设计时， 可考虑通过改变导叶开度、 导叶叶轮的叶片数及其出口面积等方案，控制流动状态和调节压力梯度以削弱气相在叶轮流道内聚并与滞留的程度。

图6 导叶叶轮内气相分布

3.5 叶片表面气相分布

图7给出了不同IGVF下PAT最高效率时叶片表面气相分布情况。 由图7可见，不同IGVF下，叶片表面气相分布并不均匀，IGVF较低时压力面气相主要分布在叶片高压边侧，而叶片低压边侧有少量气相聚集。 随着IGVF的增加，气相向叶片表面扩散，整个压力面上叶片高压边附近的气相浓度明显高于其他区域，且前盖板侧气相浓度高于后盖板侧。 整个吸力面上叶片低压边附近的气相浓度明显高于其他区域，这也与气团在叶轮流道内滞留的位置相符。

图7 叶片表面气相分布

对比图8的气相速度矢量分布可见， 两相混合工质进入叶轮流道后，气相在叶片压力面的速度方向指向前盖板侧，表明气相速度场的分布与气相聚集位置相对应，而叶片吸力面中部的气相则沿叶片分别向高压边和低压边的两侧运动，这是因为叶片吸力面中部区域流体的总压 （机械能）高于两侧高压边和低压边。 叶片吸力面逆流而上的气相与从导叶进入的混合工质在叶轮进口区域碰撞掺混，同时还受叶轮、导叶和动/静叶栅的干涉扰动，导致该区域流动紊乱，产生漩涡；顺流而下的气相汇入气团而滞留于叶轮出口区域附近。

图8 叶片表面气相速度矢量分布（IGVF=25%）

3.6 气液两相滑移速度分布

图9为IGVF=25%时PAT中截面气液两相滑移速度分布情况。 由图9可见，在导叶、叶片头部及尾部等局部滑移速度为负值，即气相速度小于液相速度。 总体而言，PAT内气相速度大于液相速度，这种差异在可能发生“气堵”的叶轮流道进口区域中更为明显。 这是由于旋转叶轮与固定导叶之间剧烈的干涉扰动， 使该区域的流动紊乱，相间作用增强，由式（6）可知，气液两相相互作用越强，两相间的速度差越大，从而导致了气液两相滑移速度的增大。 与图6c对比分析后可知，漩涡区上游的两相滑移速度大于下游的，而漩涡下游区域也是气相聚集处，这表明气相聚结为气团而发生滞留， 降低了该处气液两相的滑移速度，导致了相态分离，同时也说明漩涡上游区域的流动较为紊乱。

图9 PAT滑移速度分布（IGVF=25%）

3.7 不同流量下PAT气相分布

图10为不同流量下IGVF=25%时PAT中截面上的气相浓度与气相流线分布情况。 由图10可见， 不同流量下PAT流道内的气相分布规律较为相似，蜗壳内气相分布相对均匀，高、低含气区域分界明显，部分叶轮流道靠近叶片吸力面存在气相聚结区，气相聚结区存在相态分离。 随着流量的增加， 叶轮流道内气相聚结程度和范围增强（2.2Qd共1个流道，3.0Qd共3个流道，3.4Qd共4个流道）。 这是因为流量越大， 混合工质的动能也越大， 同时气液两相各自的体积含量成比例增加，相间的相互作用增强。气相流线在2.2Qd时较为紊乱，存在较为明显的漩涡和脱流，加剧了能量耗散，随着流量的增大，这种情况有所改观。

图10 不同工况PAT气相浓度及气相流线分布（IGVF=25%）

4 结论

4.1 随着IGVF的增加，相同流量下PAT的水力效率和输出功率降低而压头升高，最高效率时对应的流量发生偏移，与IGVF为5%时相比，最高效率分别下降了1.2%、2.2%和4.5%。

4.2 相同IGVF下， 沿蜗壳流道截面气相浓度和气相分布均匀度均降低，气相聚集于蜗壳截面中下方；同一蜗壳截面处的气相浓度随IGVF的升高而增大。

4.3 PAT各过流部件中，蜗壳内气相分布相对最均匀，随IGVF和流量的升高，各过流部件的气相浓度和分布的不均匀度均增加， 存在明显的分界。

4.4 叶轮流道内气相聚向叶片吸力面， 随着IGVF和流量的增加，聚集程度和范围增强，气相聚并成为气团而滞留于叶轮流道， 发生相态分离，可能造成“气堵”；气团附近上游存在对应漩涡， 表明叶轮流道内漩涡的形成与气相聚并相关；IGVF较低时叶片压力面气相主要聚集在高压边侧，随着IGVF的增加，气相在叶片表面扩散，叶片压力面高压边附近靠近前盖板侧的气相浓度较高，叶片吸力面低压边附近的气相浓度较高。

4.5 PAT内的气相速度总体上大于液相速度，且在气团滞留的叶轮流道进口区域较明显。