肖业祥，郭 豹，张 瑾，梁权伟，刘 洁

（1．清华大学水沙科学与水利水电工程国家重点实验室 & 能源与动力工程系，北京市 100084；2．东方电气东方电机有限公司研究试验中心，四川省德阳市 618000；3．中国长江三峡集团有限公司机电技术中心，四川省成都市 610095）

0 引言

我国西南地区水电技术可开发装机容量为44982万kW，约占全国的68%，该地区水力资源主要以高水头资源为主，开发潜力巨大[1]，但受诸多不利条件影响，开发难度大。冲击式水轮机的应用水头范围非常宽，在30～3000m内均可适用，也是一种适合中高水头电站的重要机型，有望在我国西南地区的高水头水力资源开发中发挥重要作用[2]。此外，在抽水蓄能机组的应用中，也有着采用冲击式水轮机三机式方案的电站——西藏羊卓雍湖抽水蓄能电站[3]。采用冲击式水轮机三机式方案的优点较为明显：①可适用更高水头条件；②可以“水力短路”运行，即用水轮机的出力对水泵入力进行补偿，使水泵工况实现入力可调；③工况转换和负荷响应速度比常规可逆式机组更高。近年来，冲击式水轮机也被应用于反渗透式海水淡化装置的余压能回收一体机中，节能效果明显[4]。

冲击式水轮机中，水流经喷射机构后加速出射，在大气环境中短时间经过后，喷入曲率急剧变化的水斗表面，并快速从水斗流出。在上述流动过程内，存在着复杂的水气两相流动；当流动含沙时，其内部流动特性会变得更为复杂，泥沙颗粒在水气两相流作用下会与喷射机构和水斗壁面发生剧烈的碰撞，从而导致相应过流部件的泥沙磨损[5]。

在冲击式水轮机的含沙三相流动特性方面，宋文武等[6]对泥沙颗粒在冲击式水轮机水斗的流动过程进行了数值研究，定性分析了泥沙粒径、浓度等参数对水斗表面磨损程度的影响。近些年来，葛新峰等[7]对冲击式水轮机喷射机构的泥沙磨损特性进行了数值和试验研究，探讨了粒径和浓度对磨损特性的影响。肖业祥等[8-9]对冲击式机组的喷射机构及水斗内部的含沙流动特性进行研究，分析了泥沙颗粒在其内部的运动过程，而磨损的预测结果也与实测数据取得了定性的一致。Leguizamon[10]基于拉格朗日方法发展了冲击式水轮机的含沙三相流动模拟方法，并通过实测的磨损值对计算结果进行了间接的验证。Bajracharya等[11]对Chilime电站机组进行了长期泥沙磨损跟踪测量分析，给出了过流部件磨损量与机组效率下降间的量化关系。Padhy等[12]在冲击式水轮机模型试验装置中，对水斗磨损情况与效率间的关系也进行了分析。Felix等[13]对HPP Fieschertal电站磨损的水斗进行跟踪测量，记录了在不同磨损程度下机组的水力性能参数。Rai等[14]对水斗的磨损特性进行了模型试验和真机跟踪实测分析，系统地研究了水斗表面的泥沙磨损特性。

本文采用数值方法，预测了冲击式水轮机内的水气沙三相非定常流动过程，分析了泥沙颗粒的运动特性及其与流体的分离机制；通过综合分析冲击角度、次数以及相对速度的变化，对旋转水斗的泥沙磨损特性进行了定性的预测，与实测磨损结果进行了对比分析。开展冲击式水轮机的泥沙磨损特性预估研究，对于提升装置运行稳定性具有重要意义。

1 计算模型与数值方法

1.1 机组模型与运行参数

本文以转轮水斗数为19枚的冲击式水轮机为研究对象，设计运行水头456m，设计转速500r/min，额定出力3.73MW，单位转速n11为41r/min。根据从喷射机构到转轮整体的几何结构特点和参数，建立了如图1所示的三维流场计算模型。计算域整体被划分为4个区域，分别为喷射机构、射流域、静止域和转轮域。

图1 数值计算三维模型Figure 1 Computational simulation model

1.2 计算网格

网格的划分和布局对计算结果影响较大，首先对网格的无关性进行了验证分析，本文最终确定的整体模型网格节点数约为1500万，如图2所示[9]。由于几何模型较为复杂，计算域的大部分区域采用非结构化网格，少部分区域如射流自由出流段采用了结构化网格。除此之外，还重点对射流的水气交界面附近以及水斗表面网格进行了局部加密，用以准确地捕捉射流轮廓以及模拟旋转水斗内表面上的流动形态。

图2 数值计算整体网格Figure 2 Mesh for the computational domains

1.3 边界条件与计算设置

计算域的入口设置在喷射机构入口处，根据设计水头设置总压入流条件，泥沙颗粒也从该入口流入；静止域表面和靠近水斗根部区域的表面设置为静压出口；所有固定壁面给定无滑移的边界条件。采用交界面的设置来连接不同的计算域，动静交界面间采用滑移网格技术连接。采用非定常数值计算，即流体和颗粒的模拟都随时间步长的增加而向前推进。非定常计算的时间步长设置为7.89×10-5s。颗粒非定常运动轨迹的计算则在每个流体时间步完成后进行，即对应流体时间计算完成后，基于完成后的流场进行颗粒运动的耦合计算。

整个流动过程中包含水—气—沙三相，数值计算时将水—气两相视为连续相，采用基于欧拉方法的VOF法模拟；由于沙粒的体积分数较低，故将其视为离散相，采用拉格朗日方法对它的运动进行追踪。在连续相模拟中，采用不可压缩流体假设和雷诺时均的方法模拟水气两相的流动，湍流模型选择SST-CC模型，数值计算方法采用Couple算法，选取二阶迎风格式。

1.4 磨损预估模型

根据电站在汛期3个月的泥沙实际监测数据，过机泥沙的平均浓度为1.31 kg/m3，过机泥沙颗粒粒径范围在0.001～0.23mm之间，粒径分布通过Rosin-Rammler方法进行拟合[8]，并植入Fluent进行计算。磨损速率可以用质量损失率描述，即给定时间内由于颗粒撞击导致的靶材体积或质量减少量。本文采用Mansouri[15]基于沙粒—不锈钢冲击试验所确定的磨损预估模型，该模型借鉴了Oka建立的冲击角度函数（磨损速率与冲击角度的关系），以确定不同冲击角度下磨损速率的相对值[16]。

2 喷射机构的预测结果与分析

2.1 分析截面的选取

为分析泥沙颗粒在不同过流部件内的运动规律，在沿流动方向的四个位置上选取了四个分析截面，具体位置如图3（a）所示，图中坐标0点为喷射机构出口位置；另选取了喷射机构0°和45°位置上的两个轴面作为分析截面，如图3（b）所示，并采用截面附近±5mm的空间用以表示对应截面的颗粒分布，筛选范围如图3（b）中虚线所示。

图3 流动分析截面的选取Figure 3 Location of flow behavior sections

2.2 含沙三相流动特性

图4（a）、（b）中给出了0°和45°两个截面的流体速率分布图，两个截面上的喷针尖端附近都出现了一个明显的低速区，这是由于喷针表面持续发展的边界层导致的；针尖低速区也影响了自由射流区的速度分布，即射流中心的速度也相对较低，但沿着流动方向，这种影响在逐渐减弱，射流在远离喷针尖部区域的速度均匀性有所提高。还可以发现，在尖端低速区的作用下，喷针尖端也形成了局部的高压区，如图4（c）、（d）所示。图4（e）、（f）为粒径1.1e-4m的颗粒在两个不同截面上的分布情况，而两图中颗粒分布情况近似一致，因此可认为，颗粒在喷射机构和射流不同角度轴截面中的分布均无明显差异，沿周向的分布较为对称。

图4 喷射机构的射流与泥沙分布Figure 4 Flow and sediment distribution of the injector and jet

喷射机构及射流中主流的运动特性在很大程度上决定了颗粒的运动轨迹，而垂直主流方向的速度则会影响颗粒相对壁面的运动趋势，是影响磨损特性的重要因素。将流体和颗粒速度投影到如图3（a）所示的径向截面上，得到喷射机构收缩段和自由射流段不同截面上的含沙三相流动状态，如图5所示。每个截面的流态分布图中，左侧为颗粒的分布，颜色刻度表征截面上速度分量的大小；右侧为截面上流速投影分量的流线图，描述垂直主流方向上的二次流。

由图5（a）、（b）可见，沿流动方向，喷射机构内流体垂直于流向的速度分量在逐渐增大。这主要是喷射机构内收缩流道使压力能转换为动能导致的，说明喷射机构流道的收缩不仅会引起主流方向速度的提高，也会增加垂直主流方向的流动速度。还可发现，由于导流板的存在，在肋板的正下游处形成了明显的二次流现象，这些区域的流线存在不同程度的弯曲，截面内的速度分量也出现一定程度的增加。这种二次流现象在-2d0截面处最为明显，而在-d0截面处，由于轴向速度较大，该现象虽表现得不明显，但其强度仍然很高。

图5（c）、（d）给出了自由射流区不同截面上的投影速度分布。图中可见，在四个导流板的下游方向，卡门涡导致的二次流仍然在发展；由于空气持续地被卷入射流，这种二次流效应表现得更为明显，在每个导流板对应区域都出现了较大尺度的旋涡团，且这些旋涡会随时间变化而交替发展变化。流速在二次流发展的区域仍然较高，但速度的最大值沿着流动方向（d0→2d0）在下降，这说明二次流的强度有所减弱。在导流板引起二次流的下游区域，颗粒垂直于流向的分量速度仍然较大，说明此时颗粒仍受二次流的影响。

图5 沿流向不同截面的流态与颗粒分布（Dp=1.1e-4m）Figure 5 Flow and sediment distribution at different sections along the flow direction

2.3 喷射机构的磨损特性

根据数值计算结果，通过综合分析喷射机构中喷嘴与喷针表面的泥沙冲击角度、冲击次数以及相对速度等变量的分布，得到喷射机构表面的磨损形态分布，并与文献实测的磨损结果进行了对比分析，如图6所示。其中，图6（a）、（b）为预测的磨损分布，图6（c）、（d）为文献中的磨损形态照片。数值预测的喷针和喷嘴磨损结果显示，喷针表面的磨损集中在出口附近，且受上游导流板的影响，磨损量沿周向分布不均，存在四个沿流向的条状严重磨损区域；喷嘴的磨损现象也集中出现于出口附近的颈部，且其磨损程度显著高于喷针表面，其分布沿周向相对较为均匀。通过对比数值预测结果和实际情况，预测的喷针与喷嘴表面的磨损形态与实测的部位较为吻合。

图6 喷射机构磨损特性的对比分析Figure 6 Comparative analysis of the injector sediment abrasion

3 旋转水斗的预测结果与分析

冲击式水轮机水流经喷射机构加速后，形成集束性较强的自由射流并流入旋转的水斗而做功。在进入水斗前，流体主要受到静止的喷射机构作用，可近似认为是定常流动；而进入水斗后由于水斗的旋转，射流在水斗内部为带有移动自由边界的水膜流，表现出强烈的非定常特性。数值计算分析时，选择水斗分水刃刚接触到射流时为起始时间步，以此时水斗的位置为初始位置，以此位置为基准定义水斗旋转角θB，如图7（a）所示。图7（b）给出了水斗在接收射流的过程中，其相对扭矩随水斗旋转角的变化。随着水斗旋转角的增大，射流逐渐喷入水斗做功后全部流出，其扭矩先增大后减小为0，呈类抛物线形状，本计算工况下，单枚水斗在接收射流做功的一个周期内，转轮约旋转80°。

图7 水斗旋转角θB与相对扭矩Figure 7 Bucket rotation angle θB and relative torque

3.1 非定常流动特性

射流进入旋转水斗形成具有自由表面和移动边界的光滑水膜流。随着转轮的旋转，这种水膜流覆盖的区域将发生明显的变化。当该流动挟带泥沙颗粒时，在水膜流的作用下，沙粒的运动轨迹也将发生明显的变化，同样表现出强烈的非定常特征。图8中分别给出了在不同的水斗旋转角/位置下，水斗内表面的流态分布。射流喷入水斗后的流动大致经历三个阶段，图8（a）～（c）分别选取了这三个阶段的典型位置进行分析。A位置的开始阶段，射流经分水刃流入表面，形成水膜并改变运动方向流向水斗侧边，这一阶段射流能不受阻挡地进入水斗，且流量逐渐增大。B位置的中间阶段，射流与分水刃近乎垂直，水膜几乎充满水斗表面，完整的射流持续喷入水斗且流量近乎不变，但进入水斗射流的尾部已开始被下一水斗切割。C位置的后期阶段，部分不完整的射流喷入水斗，进入的流量逐渐减小，水膜在水斗表面的覆盖面积也逐渐减小，直至最后一滴水从水斗完全流出。

图8 旋转水斗内的流态发展过程Figure 8 The flow process in the rotating bucket

图9中分别为水斗在三个不同位置表面颗粒与流态的分布，图左侧为泥沙颗粒的分布，右侧为水膜流态。与水膜流的演变相似，颗粒分布的形状也经历了先扩大，然后覆盖大部分水斗表面并保持一段时间，最后因水膜的逐渐流出而减小的过程。在进入水斗前，颗粒随射流高速流动，当射流接触水斗并以水膜的形式扩散时，这些颗粒也会逐渐分散开来。图9（a）显示，在水膜的扩散过程中，泥沙颗粒有向水斗的根部迅速扩散的趋势。第二阶段从图9（b）可见，水膜覆盖了水斗表面，在靠近水斗出水边，颗粒出现了富集的现象，这主要是由于水斗曲率较大，局部离心力过大导致的。此外，在此过程中，泥沙颗粒在水斗出水边靠近根部区域较尖部更密集一些。在C位置的水膜萎缩阶段，如图9（c）所示，颗粒较均匀地分布在水斗表面，与B位置相比，颗粒向水斗出水边尖部区域扩散的趋势更为显著。

图9 水斗内部颗粒分布Figure 9 Particle distribution in the bucket

3.2 泥沙磨损特性预测

影响水斗表面磨损特性的因素主要有颗粒的冲击次数、冲击速度和冲击角度等。在数值模拟的过程中，这些变量均可采用自定义函数来记录相关数据。图10（a）～（c）中分别为预测的颗粒冲击次数、平均冲击角度和平均冲击速度。由于水斗磨损是一个比较短暂的过程，因此将颗粒的冲击次数取为一个累积值，并随着水斗的旋转而单调增加。图10中可见，在靠近根部的水斗出水边区域，颗粒的累计冲击次数最大，沿水斗出水边从根部到尖部冲击次数逐渐减小，靠近尖部的水斗内面的冲击次数较小。水斗出水边区域的平均冲击角度比较大，尤其是靠近水斗尖部的区域，沿水斗出水边向分水刃的方向，平均冲击角度会逐渐变小；颗粒平均冲击速度在水斗中间区域较大，尤其是靠近水斗尖部的区域，但相对于累计冲击次数和平均冲击角度，平均冲击速度分布相对较为均匀。

图10 预测的水斗表面泥沙冲击特性分布Figure 10 Predicted distribution of impact characteristics of sand on the surface of bucket

综合考虑泥沙冲击个参数的影响，预测了旋转水斗表面的泥沙磨损分布，如图11（a）所示。总体上来说，水斗的主要磨损区域，分布在分水刃处、水斗切口边以及水斗靠近根部的出水边，这与文献中现场实测的结果较为相似，如图11（b）所示。水斗表面的磨损分布与冲击次数分布最为接近，这表明颗粒冲击次数对水斗磨损影响最大，在大多数区域颗粒冲击次数在水斗磨损模式中起着更大的作用。

图11 实测与预测的水斗磨损特性对比Figure 11 Comparison of measured and predicted bucket abrasion characteristics

4 结论

本文采用数值计算方法，研究了冲击式水轮机从喷射机构至旋转水斗内部全过程的水—气—沙三相流动特性，并预测了喷射机构和水斗表面的泥沙磨损情况。主要结论如下：

（1）发展了一种适合冲击式水轮机内水—气—沙多相流动计算的数值方法，该方法结合VOF和DPM模型，将水气处理为连续相，沙粒处理为离散相，能较为准确预测静止的喷射机构以及旋转的水斗内部的三相流动特性。

（2）喷射机构中的四个导流板对射流形态及泥沙的运动影响较大，从而影响泥沙磨损部位的分布；喷射机构磨损严重的部位主要集中在喷嘴颈部，其次是喷针表面；受上游导流板的影响，喷针表面有四道沿流向的条状严重磨损区域。

（3）射流进入旋转水斗形成具有自由表面和移动边界的光滑水膜流，当挟带泥沙颗粒时，在水膜流的作用下，沙粒的运动轨迹也随水膜流动而扩散。考虑泥沙颗粒的冲击次数、冲击速度和冲击角度等影响因素，预测了水斗表面的泥沙磨损分布，与文献中的实测分布较为相似；更为准确的磨损预估还需考虑更充分的影响因素，包括沙粒物性、实际流态与材质等。