1、水轮机转轮固液两相三维紊流计算及磨损预估高忠信，周先进，张世雄，陆 力(中国水利水电科学研究院机电所)摘 要：本文在两流体模型的基础上，使用贴体坐标的有限体积法，利用有科氏力修正的 -Ap 两相流紊流模型建立了水轮机转轮内部的三维泥沙固液两相紊流计算模型。模型中求解了考虑压力对固粒相影响的时均 Reynolds 动量方程，并利用泥沙相的速度和浓度分布，数值计算了转轮的泥沙磨损。为了验证模型的可行性，用该模型计算了刘家峡 HL001 模型水轮机转轮内部的三维泥沙固液两相紊流流动，数值预估了转轮的泥沙磨损并与磨损试验结果进行了比较。关键词：水轮机；固液两相紊流；数值模拟作者简介：高忠信(1961-
2、)，男，河北故城人，高级工程师，研究方向：计算流体力学及水力机械中的内部流动。我国有泥沙磨损的水轮机约占全部机组总容量的 40%，三峡建成后，预计汛期河水中也含有一定数量的泥沙，造成对水轮机的危害。所以研究水轮机内部的三维固液两相流动是十分必要的。本文在两流体模型的基础上，使用贴体坐标的有限体积法，建立了三维粘性固液两相流数值计算模型，利用 -A p两相流紊流模型计算了刘家峡 HL001 模型水轮机转轮内部的三维泥沙固液两相紊流。本文与以往的研究 1不同，考虑了压力对固粒相的影响。最后根据两相流动的计算结果预估了转轮的泥沙磨损并与试验结果进行比较。1 基本方程两流体模型中固粒相假定为拟流体。两
3、相的流动分别用一组流体运动守衡方程描述，其中包含了相间相互作用项。流体和拟流体共存于流动域中，各占一定的体积率。在与水轮机转轮一同旋转的直角坐标系中，流体相和固粒相的基本方程如下：1.1 流体相的基本方程 连续方程和动量方程为：(1)-2/3 c c ij+Gi c c+ c c(FCi+FIi)+Sci(2)式中：坐标 xj=(x,y,z)；流速 ucj=(uc,vc,wc)；重力加速度 Gi=(0,0,-g)；科氏力 Fci=(2v c,-2u c,0)；离心压力 FIi=(2 2x,2 2y,0)； c、 c、p 表示流体相的体积浓度、密度、压力。 为旋转角速度；S ci为相间作用力，假
4、定 Sci=；n 为固粒相粒子种类的个数。 F pki的矢量形式为 pk，表示第 k 种粒子所受阻力 DK、升力 LK和虚拟质量力 Vk的合力，即pk= Dk+ Lk+ Vk (3)阻力 Dk由下式计算：Dk= k(3/4Dpk) pCd| pk- c|( c-pk)(4)式中：D pk为颗粒的直径， pk、 c为颗粒相和流体相速度的矢量形式，C d为阻力系数，是颗粒 Reynolds,Rep的函数由 Drew 和 Lahey2的工作，升力 Lk和虚拟质量力 Vk的计算式如下：Lk= kCL c c- pk( c) (5)Vk= kCV cD c/Dt-D pk/Dt (6)式中：C L为升力
5、系数，取值 0.25。C V为虚拟质量系数，取值 0.5。另外 e=+ t。1.2 固粒相的基本方程 固粒相粒径为 Dpk的粒子所满足的连续方程和动量方程为：(7)(8)式中：速度 upki=(upk,vpk,wpk),重力加速度 Gi=(0,0,-g)，科氏力 Fci=(2v pki,-2u pki,0),离心力 FIi=(2 2x,2 2y,0)， pk, p表示固粒相的体积浓度、密度。1.3 紊流模型 连续相紊动能 及耗散率 方程为：(9)(10)其中： t=C c 2/，P k为紊动能 k 的生成项，由下式确定(11)Gc为 Howard3 等提出的科氏力的修正项，由下式确定(12)各
6、系数的意义和取值与清水紊流计算相同。对固粒相，使用 Hihze-Tchen 颗粒湍流粘性系数模型 4 ，称之为 Ap模型。1.4 两相间体积浓度的关系 流体相的体积浓度与固粒相体积浓度满足下式：,其中，n 为固粒相粒子种类的个数。2 数值方法和边界条件以上的基本方程可表示成统一的形式(13)式中： 为扩散系数，S 为源项， 为广义变量。引入变换函数 j= j(x,y,z)(j=1,2,3 ,)，将上述方程转换到计算空间 =(，)=( ， ， )可表示成：(14)这里，U j=(j=1,2,3 U,V,W)为反变速度，J 为雅可比矩阵，U j和 qij定义参见文献5。对方程(14)使用有限体积法
7、和交错网格进行离散，其中扩散项和源项采用二阶中心差分格式对流项采用混合格式(Hybrid Scheme)。基于 Simple 算法从流体和固粒两相的连续方程出发导出压力修正方程 6。边界条件的设定与清水紊流计算基本相同，只是固粒相的固壁处设定为滑移边界条件。3 计算结果为检验本文所采用两相流计算模型的合理性，对有磨损试验结果的刘家峡HL001 模型转轮进行数值模拟 7 。以下分别就试验方法和条件，计算条件的设定及计算结果和试验结果的比较进行说明讨论。3.1 试验方法和条件 该试验是在中国水科院机电所的浑水实验台上进行的，采用的方法是易损涂层法。试验中在模型转轮上涂上三层不同颜色的涂料，表层为灰
8、色，中间层为红色，底层为兰色。每层厚度平均约 5 道。试验时间为42min，试验时间定为 42min 是由于经 42min 试验后，各磨损部位及其相对强度都已便于区别和比较，试验时间太短磨损部位及强度不宜区分，太长磨损严重的部位涂层将被磨掉。由于涂层有一定的厚度，如果表层在某一区域变暗则表示该区域已呈现轻度的泥沙磨损。磨损的显示，由灰至红至兰色磨损的强度逐渐增强，因此从面积和颜色可以比较相对磨损强度。对HL001 模型转轮作了最优工况和限制工况的泥沙磨损试验，含沙量约为16.5kg/m3。在计量磨损量时，除了通过拍照反映各部位的磨损情况外，还用透明描图 1 HL001 模型转轮磨损试验结果图纸
9、和彩色铅笔将各磨损部位的位置及破坏面积进行描绘。试验结果如图 1 所示。各个部位的磨损面积如表 1 所示。按每层 5 道(100 道=1mm)计算，各工况的不同表面磨损的体积如表2 所示。表 1 转轮上各磨损区域的平均破坏面积区域及面积/mm 2层次A B C D E F最优工况 底层 365.3 223.8 0.0 5.46 0.0 79.62n 1=75 中层 527.2 293.0 0.15 9.46 15.4 145.02Q 1=916 表层 897.8 384.5 67.65 9.46 240.1 472.02限制工况 底层 411.2 238.2 7.3 0.0 68.3 116.
10、7n 1=68.4 中层 1028.5 398.2 37.3 5.0 254.5 303.4Q 1=840 表层 1798.5 564.0 158.1 5.0 2126.5 303.43.2 计算条件3.2.1 转轮的几何参数和计算网格 转轮型号为 HL001，转轮直径 D=0.25m，导叶高度 B=0.05m，叶片数 Z=14。计算区域为一个流道，进口向来流方向作了延伸，出口向下游方向作了延伸，整个计算区域的三维计算网格是 702127。表 2 转轮上各磨损区域的破坏体积正面/mm 3 背面/mm 3 下环面/mm 3 上冠面/mm 3 最优工况 89.52 49.68 34.83 0.0
11、174.03n 1=75Q 1=916 47.92% 28.55% 20.01% 0.0 100%限制工况n 1=68.4 161.81 70.66 36.18 0.0 268.65Q 1=840 60.23% 26.30% 13.47% 0.0 100% *仅对叶片间的磨损体积进行计算3.2.2 计算工况及泥沙条件 对表 3 所列最优工况和限制工况进行计算。其中进口含沙浓度由含沙量转换成了体积率，并假定为均匀分布。泥沙粒径设定为0.024mm，比重为 2650kg/m3。表 3 计算工况及泥沙条件工况 水头/m 流量/(m 3/s) 转速/(r/min) 进口含沙浓度 (体积率)/% 粒径
12、d50/mm最优工况 13.11 0.1893 986.48 0.6226 0.024限制工况 12.15 0.2064 1081.67 0.6226 0.0243.3 固液两相流的计算结果讨论 对表 3 所列工况实施清水单相流和固液两相流的数值模拟。为减少篇幅这里仅对固液两相流的计算结果进行讨论。图 2 给出了最优工况叶片正背面的压力分布，这里需要指出的是此压力分布不是实际的压力分布，而是在规定进口点(I=2,j=2,k=2)压力为零的情况下，算出的相对压力分布。计算结果给出的低压区位置与现场观察实际转轮的发生空化的区域基本吻合。图 3 给出了最优工况相对流速(V= cVc+ pkVpk)矢
13、量分布。清水与浑水情况下叶片背面最低压力和靠近叶片背面最大流速的比较见表 4，从中可知，叶片背面靠近出口的最低压力有所降低，这表明浑水情况下空化系数比清水情况要小。最优工况的泥沙体积率在叶片正，背面的分布如图 4 所示，由于泥沙的比重比水的比重大使得泥沙在科氏力的作用下偏向叶片的正面，造成叶片正面的泥沙浓度远高于叶片背面的浓度，浓度最高的区域发生在叶片正面的出口处，而浓度最低的区域出现在叶片背面，有关汽轮机的两相流也出现了类似现象。上冠与下环间，下环的泥沙浓度高于上冠的泥沙浓度，这是重力和离心力双重作用的结果。另外计算结果表明靠近叶片的地方泥沙浓度的变化梯度很大。表 4 最小压力与最大流速叶片
14、背面最小压力/m 靠近叶片背面的最大流速度/(s/m)工况两相流计算 清水计算两相流计算 清水计算 c水+ pk泥沙 水-13.71 -12.17 18.78 17.02最优工况I=16,j=1,k=6 I=16,j=1,k=6 I=31,j=2,k=4 I=39,j=2,k=4-11.31 -10.51 19.86 17.26限制工况I=38,j=1,k=16 I=28,j=1,k=17 I=30,j=2,k=4 I=39,j=2,k=4图 2 最优工况叶片背面的压力分布 图 3 最优工况相对流速矢量分布图 4 最优工况下叶片正、背面泥沙体积率分布图 5 最优工况叶片正面、背面及下环面的磨损
15、强度分布图 6 限制工况叶片正面、背面及下环面的磨损强度分布3.4 磨损强度的计算3.4.1 磨损强度的计算方法 与水轮机磨损有关的项目有，含沙量、粒径、形状、级配、硬度、水轮机材质以及水轮机相对速度和流态。野崎次男 8 根据秘鲁18 个水电站有关磨损的调查资料，给出估算水轮机磨损强度的公式如下：I=PEV Z (15)式中，I 为磨损强度(mm/h)； 为水轮机形式和部件的差异系数，对混流式的水轮机转轮，取 0.453210-7(mm/h);V 为泥沙与水轮机的相对速度(m/s)；z 为相对流速的修正系数；PE=P Xayk1k2k3，P 为水流中悬移质含沙量(g/l);a 为泥沙粒径，以
16、0.05mm 为基准(a=1.0)；k 1为泥沙粒径系数；k 2为泥沙硬度系数；k 3为水轮机材料的抗磨系数；x 为含沙量 P 的修正系数；y 为泥沙粒径的修正系数。3.4.2 有关系数的取值 与磨损有关的系数取值如表 5 所示。表 5 与转轮磨损有关的系数取值差异系数/(mm/h)泥沙粒径系数k1泥沙硬度系数k2水轮机材料的抗磨系数 k3含沙量 P 的修正系数 X泥沙粒径的修正系数 y相对流速修正系数z0.453210-7 1.0 1.0 1.0 1.0 1.0 3.03.4.3 磨损强度的计算结果 将通过转轮内两相流计算得到的泥沙浓度及速度分布结果代入公式(15)，计算出两个典型工况转轮的
17、磨损强度分布。图 5 为最优工况下，叶片正面，背面及下环面的磨损强度分布，从计算结果看，磨损部位发生在叶片正面靠近出口三角区域，叶片背面靠近下环的条形区域及下环面上靠近正面的条形区域。叶片正面、下环面的磨损部位与实验结果基本一致。其中正面的磨损部位与实验结果吻合的最好，计算获得叶片背面的重磨损部位较实验结果要大一些。图 6 给出限制工况叶片正面，背面及下环面的磨损强度分布，与最优工况下相比磨损部位基本相同，从两个工况的磨损强度来看，最优工况的磨损强度较限制工况要轻，这与实验的结果相一致。表 6 给出了转轮上各磨损区域的破坏体积及百分比。表 6 转轮上各磨损区域的破坏体积正面/(mm 3/h)
18、背面/(mm 3/h) 下环面/(mm 3/h) 上冠面/(mm 3/h) 最优工况 6.52 5.32 1.30 0.16 13.3n 1=75 Q 1=916 49.02% 40.00% 9.77% 1.20% 100%限制工况 8.38 7.14 1.44 0.24 18.09n 1=68.4 Q 1=840 48.65% 41.96% 8.13% 1.27% 100%*仅对叶片间的磨损体积进行计算由于磨损实验采用的是涂层方法，计算与实验的磨损强度无法进行数值上的比较，但就磨损部位和各各过流表面的磨损强弱的顺序而言，计算与实验结果基本吻合，计算表明叶道间个表面的磨损顺序是正面背面下环面上
19、冠面，实验也是这样的顺序，同时计算表明不同工况的磨损是限制工况最优工况，实验也为同样结果。图 7 给出了刘家峡转轮的实际磨损图形，从磨损部位看与计算结果非常一致。图 7 刘家峡水电厂 HL001 转轮真机磨损情况4 结论(1)对刘家峡 HL001 模型转轮进行了固液两相流及磨损强度的三维数值模拟，通过计算与实验结果的对比，表明采用的固液两流体模型是可行的，数值方法也是正确的。(2)通过对最优工况和限制工况固液两相流及磨损强度的数值模拟，表明开发的固液两相流计算软件能够反映不同工况水轮机转轮的磨损情况，在一定程度能够代替模型实验对水轮机转轮的磨损部位和磨损强度进行数值预估。参 考 文 献：1 吴
20、玉林，何燕雨，等.水轮机转轮内部的三维固液两相紊流计算J.水利学报，1998，(3)：17-212 Drew D A,Lahey R T Jr.The Virtual mass and Lift Force on a Sphere in Rotating and Inviscid Flow.Int.J.of Multiphase Flow,1987,13(1):113-121.3 Howard J H G,Patankar S V,Bordynuik R M.Flow Prediction in Rotating Ducts Using Coriolis-Modified Trubulence
21、 ModelsJ.J.Fluids Eng.ASME,1980,102:456-461.4 周力行.湍流两相流动与燃烧的数值模拟M.北京：清华大学出版社，1991.5 高忠信，周先进，等.水轮机转轮内部的三维紊流计算与性能预估J.水利学报，2001，(7).6 Moura L F M,Rezkallzh K S.Numerical Study on The Two-Phase Flow Distribution In a T-JunctionJ.International Jour for Numerical Methods in Fluids,1993,17:257-2707 姚启鹏，赵琨.刘家峡水轮机改造增容浑水模型试验研究报告R.中国水利水电科学研究院机电所，1991.8 野崎次男.水车磨耗考虑沉沙地容量断面决定J.电力土木，1989，(1)：143-152.
