水泵水轮机增减负荷过程三维流动特性大涡模拟分析

胡秀成，张立翔

（昆明理工大学 建筑工程学院，云南 昆明 650500）

1 研究背景

水泵水轮机内部是一种复杂的湍流流动，具有普遍性和特殊性，而“S”形逆变流、动静干涉流、工况变迁流和叶道旋转切变流等是当前研究的重点。李君等［1］采用全流道CFD数值模拟的方法对“S”形区域流场特性进行了分析。游光华等［2］采用导叶不同步预开装置（MGV装置）来改善“S”形流动特性，增强了机组的运行稳定性。但改善不稳定的MGV装置也存在一些弊端，例如，Liu等［3］发现采用MGV装置时压力脉动幅值是采用同步导叶时压力脉动幅值的2倍，流动失稳现象也显著增加；Xiao等［4］发现MGV破坏了流场的对称性，流动更加复杂，在一定程度上增大了压力脉动幅值。陶然等［5］和王焕茂等［6］分析了泵工况的驼峰特性。Hasmatuchi等［7］在导叶与转轮之间注入气泡采用特殊的图像叠加技术观察到制动工况导叶流道内的失速形态，测得失速转速约为额定转速的70%。夏林生等［8］选用SST-SAS湍流模型采用三维非恒定数值模拟方法对某低比转速水泵水轮机同一开度下四象限内不同工况进行了数值分析。

一般情况下，电力系统调度有静止、发电、发电方向调相、抽水和抽水方向调相等5种基本运行模式，而每种模式间的工况变换排列组合多达20余种，工况常一日变换数次，甚至一小时变换数次。在工况变换的增减负荷过程中，由于流态的巨变可能会诱发抽水蓄能系统发生剧烈的耦合振动，因此相关的分析研究对于抽水蓄能系统的安全稳定运行有着重要意义。在当前的技术条件下，对于诸如水力机械这类流道几何空间复杂且具有高雷诺数的湍流数值模拟问题，大涡模拟（large-ed⁃dy simulation，LES）无疑是最具有潜力和可实现的湍流数值模拟方法。本文采用商业软件ANSYS FLUENT17.1，选用LES湍流模型［9-13］对某水泵水轮机由额定水轮机工况转为额定水泵工况过渡过程轨迹线上的10个工况点对应的流动状态进行三维非定常数值模拟，对小流量工况的流动结构变化特性进行分析，探索工况变迁与不稳定流动结构变化特性间的关系。

2 计算方法

2.1 计算区域及网格设计 本文以某水泵水轮机为模拟对象，其主要参数为：转轮直径4678.5 mm，转轮进口高度656.25 mm，额定转速300 r/min，额定水头200 m，额定水轮机工况和水泵工况导叶开度分别为23度和22度。计算区域如图1所示，包括蜗壳区域、导水机构区域（19个固定导叶和20个负曲率活动导叶）、转轮区域（7个叶片）和尾水管区域。

根据不同区域的流动特性及大规模并行计算的特点进行网格的设计和分区划分，流固界面为非滑移界面，转轮旋转区域形成的动静界面上设置滑移界面，采用HyperMesh实现各区域网格划分。蜗壳区域（除鼻端）、导水机构区域、尾水管区域采取结构网格划分，其余区域采取适应性强的四面体非结构化网格划分。导水机构区域、转轮区域网格最大尺度控制不超过15 mm，其中叶道及近壁区域的网格尺度控制为不超过6 mm，第一层网格离壁面的距离y+≤25；蜗壳进口区域与尾水管区域网格尺度控制分别为不超过40 mm和50 mm。按此原则划分各计算区域形成基本网格方案A，其总单元数为19 685 714个，局部区域网格如图2所示。

为检验网格设计的合理性，在网格方案A的基础上，分别采用1倍和2倍加密方式，形成网格方案B和C。以额定水轮机工况作为网格检验计算工况，比较叶片压力面压力分布及转轮水力矩的数值结果，判断网格的合理性。通过比较，网格方案B和C的计算结果十分接近，转轮水力矩相差仅为0.49%，故选用网格方案B作为本文数值计算网格。

图1 计算区域

图2 导叶中心处网格剖面图

2.2 计算方法及初始条件 选用大涡模拟Smagoringsky-Lilly亚格子应力动力模式，该模式克服了经典Smagoringsky常规模式的缺点。在网格方案B网格规模的原则下，采用贴体坐标下的有限体积法和非交错网格技术进行空间离散。采用二阶全隐式格式进行时间离散，其中对流项采用二阶迎风格式，源项和扩散项采用二阶中心格式。

水轮机工况转为水泵工况的增减负荷过渡过程计算按给定转速，改变导叶开度的方式设定。水轮机工况采取流速进口，压力出口；水泵工况采取质量进口，自由出流。选取过渡过程轨迹线上水轮机工况点及水泵工况点各5个进行非定常计算，参数如表1所示，表中Pn为额定工况。

表1 计算工况及参数

每个工况先进行定常计算，并以定常计算结果作为非定常计算的初始条件。考虑转轮转速以及计算规模，确定转轮每旋转1.5°（相当于0.833 ms）作为1个时间步长［14］，即转轮每旋转一周计算240个时间步，每个时间步收敛残差目标值0.001。计算时长按7个旋转周期考虑，选择已经稳定的第6个周期（T）计算结果6T时刻进行分析。由于计算规模较大，使用昆明理工大学PowerCube-S01云立方汉柏高性能计算系统中一定数量的CPU核进行并行计算。由于不同规模的并行计算选择核数过多计算效率反而下降，本文选取了节点1中的60个CPU核进行并行计算。每个工况计算需要约72～96 h，计算总计耗时约720 h。

3 计算结果及分析

3.1 水轮机工况减负荷过渡过程 图3—图8为水轮机各个工况流道内涡核图，各涡核图均为涡量大小0.1时的涡核。水轮机额定工况100%Pn时，水流在进入转轮前流体的压力和速度分布较均匀，与导叶端部略有轻微撞击，如图3所示，导致导叶端部附近有少量涡结构发展。水流进入转轮，流速梯度急剧变化，湍动逐步剧烈。理论上水轮机额定工况时转轮进口水流速度与转轮进口叶片骨线切线方向基本一致，即接近满足最优工况时无撞击进口条件，但由于动静干涉等影响，存在较小攻角。水流经叶片前缘的阻碍后沿叶片内弧和背弧流入叶道，在三维扭曲叶片曲面的黏性作用和由突起物引起的逆压梯度的共同影响下，流动在曲面上产生流动分离，在叶片吸力面形成脱流，流态失稳，脱流产生的旋涡在叶道内形成复杂的涡结构，即叶道涡。如图3所示，叶道涡由一系列不同尺度的涡组成，约占据了叶道1/3空间。叶道涡沿着吸力面流向下游，汇聚一起，相互干扰，又不断串级分叉，演化为许多小尺度涡。水流流出转轮进入尾水管，在尾水管直锥段和弯肘段内流线相互缠绕，流动状态快速变化，最终在扩散段趋于均匀。

图3 水轮机工况100%Pn流道内涡核

图4 水轮机工况75%Pn流道内涡核

图5 水轮机工况50%Pn流道内涡核

图6 水轮机工况25%Pn流道内涡核

图7 水轮机工况1%Pn流道内涡核

图8 水轮机工况1%Pn叶道涡核

随着导叶开度减小，水流流量减少，水流对叶片端部攻角增大，与导叶端部撞击加强，叶道内流动分离越来越显著，叶道涡尺度增大，叶道中管状涡和马蹄涡的分布已十分明显，表明流态呈失稳状态，如图4、图5所示。随着导叶开度进一步的减小，活动导叶间流道受限，转轮内流态紊乱加剧，流线相互缠绕且杂乱无章，叶道涡尺度相对较小但空间范围广，各叶道间的流态分布呈明显不对称状态。如图7、图8所示，水轮机工况1%Pn时，导水机构区域甚至出现了相对较多的涡结构，转轮内叶道涡结构尺度小但较为均匀，且几乎占满了整个转轮叶道空间，在叶片出水边靠近下环侧有众多涡带，转轮出口水流不满足法向出口条件，在水流绝对速度的圆周分量作用下在尾水管中心位置出现低压并形成偏心尾水涡带，此外尾水管扩散段存在回流等现象。

如图9所示，在其中一叶道内设置5个与叶面形状大小相同的虚拟面，叶片压力面近壁面S1、叶道空间1/4处面S2、叶道中面S3、叶道空间3/4处面S4、叶片吸力面近壁面S5，分别观察这5个面上涡的分布，以便揭示叶道中涡的演化规律。为便于观察，一般情况可将涡结构分解为流向涡和展向涡两个分量。但由于转轮内叶道为三维扭曲空间，其流场流向复杂，将涡分解成流向涡和展向涡较困难，为此本文采用平均涡量值进行比较分析。如图10—图15所示为水轮机工况50%Pn时各虚拟面附近涡核分布情况，如图16所示为水轮机偏工况虚拟面平均涡量。可看出：除1%Pn水轮机工况外其它各偏工况虚拟面涡量线均为管状的鱼钩型，叶片吸力面近壁面S5平均涡量值最大，叶片压力面近壁面S1次之，其余依次为S4、S3、S2。平均涡量值越大的区域涡结构越多，吸力面附近涡结构最多，表明吸力面附近流动分离相对显著。涡结构的演化伴随着能量的耗散，可定义涡能耗散系数We代表涡能，Q代表流量，H代表水头。计算可得：同一水头下，1%Pn水轮机工况涡能耗散系数最大，其余依次为25%Pn水轮机工况、50%Pn水轮机工况、75%Pn水轮机工况。由此可知：开度越小，涡能耗散系数越大，涡结构能量损耗越大。

根据上述分析以及如图3—图8所示结构分布，展示了水轮机工况不同开度下的叶道内涡系结构及分布特征。由图可知，不同攻角下，转轮叶道内涡系结构不同。攻角较小时，稳定螺旋点在叶道内发展为较小的涡结构，类似于羊角涡，促使主涡非对称发展。攻角较大时，发展为马蹄涡，通过牵制主涡涡核，抑制主涡向非对称发展。并随着攻角继续增大而演变为一个复杂的多涡系（称之为主体涡系），如图7所示，不仅包含了主涡还包含了复杂的二次涡结构，对于非定常马蹄涡系还包含了多种非定常运动模态的旋涡结构。攻角中等时，叶道内既形成羊角涡，又形成马蹄涡，还存在管状涡。

图9 叶道等分截面示意图

图10 水轮机工况50%Pn时面S1附近涡核

图12 水轮机工况50%Pn时面S3附近涡核

图11 水轮机工况50%Pn时面S2附近涡核

图13 水轮机工况50%Pn时面S4附近涡核

图14 水轮机工况50%Pn时面S5附近涡核

图15 水轮机工况50%Pn下某一叶道涡量涡核区

图16 水轮机偏工况虚拟面平均涡量

受测量手段和分析方法的制约，对涡系的研究还较为有限。各类涡结构的形态和功能不尽相同，例如羊角涡也叫直立涡、螺旋涡、龙卷涡等，是流动中三维分离的结果。羊角涡物面形态为分离螺旋点，呈龙卷形状。羊角涡的出现频率较高，结构尺度较小，涡强度较弱，但它通过主涡来体现它的水动力作用和对主涡系发展的影响。

在研究湍流时，人们常注重主涡系的发展，一般很少观察或注意到羊角涡的存在及其特性。管状涡是流动结构在叶道空间中发展的结果，由于管状涡属于拧旋涡结构，具有周向分量，外形经相邻叶片壁面的拉展形如水管状，在一般的数值计算中较难捕捉到，故在文献中极少提及，本文捕捉的管状涡位于叶道中部位置，在叶道出口位置消失。马蹄涡涡量强度较大，涡的形成机制相对管状涡而言较为单一，数值计算中容易捕捉，一些学者对平面二维流场马蹄涡进行大量研究，如陈启刚等［15］圆柱绕流产生的马蹄涡进行了研究，但三维扭曲空间内的马蹄涡仍鲜有报道。Adilozturk等［16］认为“螺旋形流线是判断马蹄涡是否存在，确定马蹄涡位置和形状的主要依据，涡结构具有螺旋形旋转特征”。据此根据流线分布特性判断出本文计算中马蹄涡出现的位置如图17、图18所示，这与如图5、图7所示中捕捉到的马蹄涡位置相符。但Jeong等［17］认为“流线的外形随着观察坐标系而改变，即流线不满足伽利略不变性，由此流线方法识别出的涡结构不具有普遍意义”。参考二维空间研究成果［15］，为保证识别出的涡结构既具有螺旋形流线特性，其几何形态又满足伽利略不变性，定义旋转强度λci为马蹄涡识别变量。

设三维流场中任一点的速度梯度矩阵为：

式中：U、V、W分别为该点x、y、z方向上的速度分量。

该速度梯度矩阵特征方程可表示为λ3+a1λ2+a2λ+a3=0，其中当判别式（其中为正时，有一实根和两个共轭复根λ=λcr±λcii，根的虚部λci即为旋转强度，其大小反映旋涡结构的强度。但需证明“λci满足伽利略不变形，并当流场任意一点的速度梯度矩阵具有复特征根时其周围流场呈螺旋形”。

图17 水轮机工况50%Pn叶道内流线

图18 水轮机工况1%Pn叶道内流线

强烈的非定常多涡系在叶道内产生很大的剪切应力和冲刷碰撞作用，增加了能量的耗散，增加流体机械的阻力、引起较大的流动能量损失、减小流体机械效率、增强流动噪声及机械振动。当涡系非对称时，在转轮上会诱导出很大的侧向力和偏心力矩。这也是不稳定的复杂流动结构变化常导致机组失稳，轻则厂房噪音增大、机组振摆增大、有功突变、调速器接力器微动频繁、并网不成功等，重则引发机组事故的重要水力因素之一，如云南某梯级电站4号发电机（600 MW）注入式定子一点接地保护电源触点因长期空载振动以致松脱导致保护动作机组跳闸。因此，在工程实践中应注意防范。

3.2 水泵工况增流量过渡过程 水泵工况1%Pn流道内流线见图19，各水泵工况流道内涡核见图20-图24。水泵水轮机在水泵工况小流量区运行时，水流在尾水管内流动平稳，进入转轮后与叶片进水边发生撞击产生漩涡结构，并在叶片出水边附近产生涡带（见图19—图20）。由于拥塞叶道内流线呈螺旋状上升，水流流出转轮进入导水机构；由于活动导叶间流道狭小，水流无法快速通过，导致了回流现象，且在活动导叶尾端脱流产生一系列涡结构。涡结构随着水流流向蜗壳，在固定导叶间大量溃灭。相对于水泵工况其它流量区工况，水泵工况小流量区运行时转轮进口水流与叶片进水边攻角较大，产生撞击，脱流较为严重。

图19 水泵工况1%Pn流道内流线

随着活动导叶开度增大，水流增大，流态越趋于平稳，流线平滑且均匀，在转轮流道内流线不再呈螺旋状，但仍可观察到少量的活动导叶尾迹。伴随着流量进一步加大，由于水流拥塞，流线越趋不均匀与不平顺。尾水管内流线相互缠绕，且水流存在旋转分量。转轮流道内流线与小流量工况时相似，呈螺旋状上升。且在尾水管弯肘段略有少量涡结构。此外，在转轮叶片压力面进口边靠近下环侧附近因撞击产生有大量旋涡结构，部分流体在漩涡结构作用下沿尾水管边壁反向流入尾水管，下环侧周向外侧流道的反向流形成分离涡。

图21 水泵工况25%Pn流道内涡核

图22 水泵工况50%Pn流道内涡核

图23 水泵工况75%Pn流道内涡核

图24 水泵工况100%Pn流道内涡核

水泵工况转轮内涡结构相对较少，导水机构区域内涡结构变化较大。如图25—图29所示，为各水泵工况导水机构中心剖面涡量（绝对值），如图30所示为导水机构中心剖面平均涡量。由此可见，水泵工况100%Pn下导水机构涡量绝对值最大，水泵工况1%Pn下最小，其它工况居中。水泵工况100%Pn时导水机构涡量绝对值受流量拥塞影响较大，而水泵工况1%Pn时最主要影响因素来自壁面黏性。

图25 水泵工况1%Pn导水机构中心剖面涡量

图26 水泵工况25%Pn导水机构中心剖面涡量

相对于水泵工况，水轮机工况转轮内涡量最高，涡系发育更为明显，涡结构较多，大部分分布在叶道前半位置，且随导叶开度变化明显。如图30所示，水泵工况导水机构平均涡量高于对应水轮机工况。且水泵工况额定工况下平均涡量最大，空载时最小。而水轮机工况反之，空载时平均涡量最大。在导水机构区域，水轮机工况涡结构变化较小，且不易观察到涡结构，而水泵工况涡结构较多且随导叶开度变化明显。

上述分析表明，水泵工况和水轮机工况在转轮区域及导水机构区域涡结构特征基本相反，但存在共性，即：水流从大空间流向小空间，流道受阻，涡结构增多且变化显著；反之，流道较顺畅，涡结构变化不明显且相对较少。

图27 水泵工况50%Pn导水机构中心剖面涡量

图28 水泵工况75%Pn导水机构中心剖面涡量

图29 水泵工况100%Pn导水机构中心剖面涡量

图30 导水机构中心剖面平均涡量

4 结论

本文采用基于Smagorinsky-Lilly亚格子应力模型的大涡模拟方法，对水泵水轮机额定水轮机工况转为额定水泵工况增减负荷过渡过程进行了三维非定常湍流数值计算，并重点分析水轮机工况及水泵工况小流量区的不稳定流动结构，探寻流态发生巨变主要因素。根据分析研究得出如下结论：

（1）由于转轮进口水流速度与进口叶片骨线切线方向不一致，存在攻角，水流经叶片前缘的阻碍后沿叶片内弧和背弧流入叶道，在三维扭曲叶片曲面的黏性作用和由突起物引起的逆压梯度的共同影响下，边界层流动在曲面上发展产生三维流动分离，形成脱流，脱流的旋涡在叶道内形成复杂的涡系，即叶道涡。计算捕捉到了叶道中发展的管状涡。

（2）水流进口攻角不同，诱发的叶道涡涡系结构也不同。攻角较小时，管状涡的破裂将导致螺旋涡在叶道内发展为羊角涡；攻角较大时，发展为马蹄涡；攻角中等时，叶道内既形成羊角涡，又形成马蹄涡，攻角越大形成的涡系越复杂。

（3）水轮机工况时，叶片近壁面平均涡量大且吸力面高于压力面，流动分离相对显著，流道中间次之。导叶开度越小，涡结构能量损耗越大。

（4）水泵工况下导水机构区域内涡结构随开度变化显著，其平均涡量高于对应水轮机工况。且水泵工况额定工况下平均涡量最大，空载时最小。而水轮机工况反之，空载时平均涡量最大。水流从大空间流向小空间，流道受阻，涡结构增多且变化显著；反之，流道较顺畅，涡结构变化不明显且相对较少。

参 考 文 献：

［1］李君，王磊，廖伟丽.可逆式水泵水轮机“S”形区域内部流场特性分析［J］.农业工程学报，2014，30（15）：106-113.

［2］游光华，孔令华，刘德有.天荒坪抽水蓄能电站水泵水轮机"S"形特性及其对策［J］.水力发电学报，2006，25（6）：136-139.

［3］LIU J T，LIU S H，SUN Y K，et al.Numerical simulation of pressure fluctuation of a pump-turbine with MGV at no-load condition［C］//IOP Conference Series：Earth and Environmental Science.［S.L.］：IOP Publishing，2012 ．

［4］XIAO Y X，XIAO R F.Transient simulation of a pump-turbine with misaligned guide vanes during turbine mod⁃el start-up［J］.Acta Mechanica Sinica，2014，30（5）：646-655.

［5］陶然，肖若富，杨魏，等.可逆式水泵水轮机泵工况的驼峰特性［J］.排灌机械工程学报，2014，32（11）：927-930.

［6］王焕茂，吴钢，吴伟章，等.混流式水泵水轮机驼峰区数值模拟及分析［J］.水力发电学报，2012，31（6）：253-258.

［7］HASMATUCHI V，FARHAT M，ROTH S，et al.Experimental evidence of rotating stall in a pump-turbine at off-design conditions in generating mode［J］.Journal of Fluids Engineering，2011，133（5）：051104.

［8］夏林生，程永光，蔡芳，等.水泵水轮机四象限工作区流动特性数值分析［J］.水利学报，2015，46（7）：859-868.

［9］张耀良，朱卫兵.张量分析及其在连续介质力学中的应用［M］.哈尔滨：哈尔滨工程大学出版社.2004：198-219.

［10］吴玉林，刘树红，钱忠东.水力机械计算流体动力学［M］.北京：中国水利水电出版社，2007.

［11］王福军.计算流体动力学分析-CFD软件原理与应用［M］.北京：清华大学出版社，2004.

［12］MENEVEAU C，KATZ J.Scale-invariance and turbulence models for large-eddy simulation［J］.Annual Review of Fluid Mechanics，2000，32：1-32.

［13］REBOLLO T C，LEWANDOWSK R.A variational finite element model for large-eddy simulations of turbulent flows［J］.Chinese Annals of Mathematics，（Series B），2013，34（5）：667-682.

［14］郭涛，张立翔.混流式水轮机小开度下导水机构内湍流特性和叶道涡结构研究［J］.工程力学，2015，32（6）：222-230.

［15］陈启刚，齐梅兰，李金钊.明渠柱体上游马蹄涡的运动学特征研究［J］.水利学报，2016，47（2）：158-164.

［16］ADILOZTURK N，AKKOCA A，SAHIN B.Flow details of a circular cylinder mounted on a flat plate［J］.Jour⁃nal of Hydraulic Research，2008，46（3）：344-355.

［17］JEONG J，HUSSAIN F.On the identification of a vortex［J］.Journal of Fluid Mechanics，1995，285：69-94.