水力机械的流动模拟与数值计算

十七水力机械的流动模拟与数值计算
1 调研背景概述
本报告主要查阅了2010年及以前的相关期刊和学术会议论文、有关专著和毕业论文等，选择能够代表本学科发展动态的科研成果作为调研的基本材料，包括：
国内期刊：《中国科学A辑》、《水动力学研究与进展A辑》、《水利学报》、《清华大学学报》、《水力发电学报》、《农业机械学报》、《机械工程学报》、《工程热物理学报》、《水科学进展》、《排灌机械》、《水力发电》、《流体机械》、《中国水利》、《大电机技术》、《东方电机》。

国际期刊：International Journal for Numerical Methods in Fluids
Transactions of the ASME—Journal of Fluids Engineering
International Journal of Heat and Fluid Flow
Computer & Fluids
Computers &Structures
Journal of Fluid Mechanics
Transactions of ASME: Journal of Turbomachinery,
International Journal of Computer Applications in Technology
专著：An introduction to computational fluid dynamics: the finite volume method 计算流体动力学分析——CFD软件的理论与应用
学术会议：Proceedings of the ASME Fluids Engineering Conference
24th IAHR Symposium on Hydraulic Machinery and Systems，2008
23rd IAHR Symposium on Hydraulic Machinery and Systems，2006
22nd IAHR Symposium on Hydraulic Machinery and Systems，2004
2 选择本专题进行调研的原因、必要性及意义
真机试验、模型试验和数值计算是水力机械性能特性分析的主要手段。

数值计算具有灵活性强、周期短、成本低、可预测性强以及可视化程度高等多方面的优势，为研究者提供了一个有效的研究手段。

近年来，随着计算机技术的高速发展以及计算流体动力学和数值求解
技术的不断完善，数值计算越来越受到人们的重视，并广泛应用于水力模型的性能预测及优化设计、水力机械非定常流动特性(湍流及湍流压力脉动)研究、空化两相流研究、固液两相流研究以及考虑流固耦合的结构动力学研究等领域。

也正是基于数值计算的以上优点，国内外很多学者基于CFD技术对水力机械内部流动进行数值模拟研究，其中涉及水力机械流动计算方法和计算模型的研究、采用数值方法对流动规律的研究以及采用数值方法解决工程实际问题的研究。

数值计算的方法和结果需要较强的理论分析来支撑，需要试验来验证、完善。

对当前水力机械内部流动的数值模拟和计算研究现状进行分析总结，有利用水力机械流动理论以及试验研究的发展，更有利于水力机械内部流动数值研究的发展。

3 近年来该专题发展新动向和值得关注点
3.1基于CFD的水力机械性能预测及优化设计
在20世纪50-70年代，水力机械中的流动计算只能更多的采用理想化假设的一元或二元理论，流动计算的目的也是仅仅估算出流道中平均流动的速度和压力，用于指导流道设计。

在这样的流动计算理论的基础上，所设计水力机械的整体能量性能较差，如效率不高、抗空蚀能力差，根据试验结果进行改型是主要的技术手段，经验起了决定性作用。

在计算机得到普遍使用以后，现代意义下的水力机械中的流动计算才开始出现。

吴仲华教授所提出的S1、S2两类流面的通用理论首先得到了发展，基于该理论开发了理想流体的准三维流动计算方法，满足了水力机械设计工况对流速场的要求。

为考虑粘性对性能计算的影响，还发展了边界层与内部理想流动的迭代计算等方法，它们在很大程度上推动了水力机械流动计算的发展。

假定粘性引起的摩擦作用只在很薄的边界层范围内，忽略粘性项，直接求解纳维-斯托克斯方程组，出现了全三维的欧拉法。

Moore把主流区看作非粘性流动，再结合边界层的粘性流动计算，可以考查回转通道内二次流对边界层发展的影响，计算结果与试验结果相当一致。

Nishi着重研究叶轮内部边界层和主流区的自律调整作用，多数情况下都存在尾流－射流结构，并被很多试验证实。

事实上水具有粘性，要想进一步提高流动计算的真实性，就必须全方位考虑粘性流动的效应。

Martelli以S2相对流面上的二维Reynolds时间平均N-S 方程和k-ε湍流模型，用有限差分法和时间推进法计算了离心泵叶轮内部二维粘性流动，以此来指导设计。

Shi et al以二维Reynolds时均N-S 方程和考虑了旋转、曲率效应的湍流模
型k-ε，用SIMPLE法计算了圆柱形叶片的离心泵叶轮内部的二维湍流流动。

将湍流理论应用于水力机械三维流动计算始于计算机技术得到迅猛发展的20世纪80年代。

1986年，Shyy和braaten首次应用k-ε湍流模型对水轮机尾水管的稳态流动进行计算研究，从技术上验证了k-ε模型在尾水管流动计算的可行性。

1987年，瑞士流体机械协会组织全世界范围内的水轮机内部流动计算，不限计算方法，进行蜗壳、转轮和尾水管的流动计算比较。

通过计算结果与试验比较得出：基于k-ε模型的湍流计算方法可用于水轮机稳定流动计算和性能预估，误差在4%之内；对单部件的计算可得到各部件的主要流动现象，基本上可满足设计高性能水轮机的需要。

之后，水力机械的三维流动计算成为大家研究的热点。

Albert（2002）领导的GAMM工作组就是较早应用CFD技术计算水轮机内部流动比较成功的例子。

CAD与CFD技术的联合使用，将CFD分析结果用于指导水力设计，形成了满足工业应用的反问题设计方法。

现在，国外一些大的水轮机制造公司采用无转轮投标，就是基于CAD－CFD技术来完成的。

在转轮设计阶段，采用先进的流场计算软件，通过对不同几何参数组合的转轮进行流动数值模拟，可以得到性能较优良的转轮。

基于CAD－CFD技术设计的转轮最高效率可以达到94％以上。

CFD技术能够较准确地预测水轮机在较大运行范围内的能量特性，最大限度地减少试验和模型加工费用。

Yang等（1998）利用CAD和CFD联合研究三峡X形水轮机叶片，研究表明，X形叶片比传统转轮叶片在非设计工况下更稳定、可靠，这对流量和水头变化大的三峡水电站尤其重要。

Kvaerner公司和加拿大GE公司针对我国鲁布格、大朝山、三峡等电站，采用现代CFD技术，并进行广泛的模型试验，设计出了高性能的X型叶片水轮机转轮。

1999年，VGS和ALSTOM为三峡左岸水电站研发设计的水轮机模型效率分别达到95.26%和94.54%； 2003年，ALSTOM为三峡右岸水电站研发设计的水轮机模型效率达到95.06%。

Gehrer等（2006）使用CFD方法对轴流转桨式水轮机转轮进行了优化。

此外，很多研究者采用CFD进行水力机械内部流动的定常计算，分析流动细节，指导水力设计。

Muntean等（2004）对蜗壳和活动导叶内的流动进行了模拟，得出不同工况下叶道涡的信息。

Tomas等（2004）使用流动模拟方法对轴流转桨式水轮机的水力设计提出了改进意见。

国内在这方面的研究也取得了很大的进步，在理论和实践方面做了很多有益地研究和探索。

1990年，中国水利水电科学研究院研发设计的水轮机模型效率率先突破93%；1999年水轮机模型效率率先突破94%。

通过三峡技术引进，东方电机厂和哈尔滨电机厂分别引进了CFD软件CFX-TASCFlow及配套计算机硬件，用于转轮的水力设计，使得转轮的设计水平有
了很大的提高，迅速缩短了与国外的差距，提高了在国际市场上的竞争力。

2003年，哈尔滨电机厂和东方电机厂为三峡右岸水电站研发设计的水轮机模型效率分别达到94.64%和94.59%。

CFD还应用于鱼友水轮机的优化设计。

ARL/NREC小组采用准三维水流模型来优化鱼友水轮机转轮的几何尺寸，利用三维CFD技术分析鱼友水轮机的水流特性，进一步完成鱼友水轮机的优化（George等，2005）。

Voith小组利用CFD工具分析水轮机水流速度和压力分布，进行鱼友水轮机设计，提出了转桨式和混流式鱼友水轮机的设计理念（Odeh等，2000）。

美国IT电力公司关于鱼友水轮机的概念设计也是基于CFD技术完成的（Tests validate design of small fish-friendly turbine, 2007）。

CFD技术在水泵水力设计方面的应用也达到了工程水平（Hornsby，2002）。

不少水泵制造商自行开发了满足特定需求的CFD软件，将流体动力学分析应用到产品的设计过程中，如NREC公司开发的pbCFD软件（Andersen，2003），德国慕尼黑大学开发的CNS3D软件等（Frobenius，2002）。

Martelli等（1990）以S2相对流面上的二维Reynolds时均N-S 方程和k-ε湍流模型，用有限差分法和时间推进法计算了离心泵叶轮内部二维粘性流动，以此来指导设计。

Shi等（1992）以二维Reynolds时均N-S 方程和k-ε湍流模型，用SIMPLE法计算了圆柱形叶片离心泵叶轮内部的二维湍流流动。

Manish等（2000）对带有叶片扩散流道蜗壳的离心泵叶轮流道中的流动做了定量的PIV可视化试验研究，并采用雷诺平均和大涡模拟方法数值模拟其内部流动规律，得到水力机械中引起机组效率下降和不稳定运行的脱流、回流、断面二次流等流动现象，以此指导水力机械的优化设计，使水力机械的能量特性和抗空蚀性能进一步提高。

Goto等（2002）基于数值模拟和三维反问题方法，提出了水泵导叶体的优化设计方法，根据给定的环量分布和轴面几何形状，进行导叶设计。

汤方平等（2005）采用径向平衡流动模型和二维叶栅面元法叶片造型，设计了高比转速轴流泵水力模型，并对设计结果进行了流体动力学分析。

曹树良等（2005）耦合反问题求解和数值模拟，成功设计了气液混输泵。

清华大学、华中科技大学、中国农业大学以及中国水利水电科学研究院等国内的一些科研机构也在应用CFD技术对水力机械（水轮机和泵）进行优化设计和性能预测等研究，如韦彩新教授等对湖南省东江水力发电厂4号机组进行了整体数值分析，并对转轮的性能改善提出了新的方案。

通过三维流动计算比较不同设计方案，可以得出一些指导性的设计经验，进一步提高转轮的能量特性和空化性能，开发优秀的水力模型。

CFD技术广泛应用于一些新建电站和泵站
的前期论证工作，为电站开发合适的水力模型；还应用于一批老电站和泵站的改造项目，经过数值仿真和优化，效率、空化性能和稳定性均取得了不同程度的改进。

CFD技术已成为水力模型开发和性能预测的重要工具。

在传统设计方法中引入全三维数值模拟技术来研究水力机械中的流动规律，全面预测水力机械的运行特性并反过来指导、优化设计，即水力机械设计的正问题和反问题，已成为水力机械研究发展的一个重要方向。

CFD技术的应用使水力模型的研制周期和成本明显下降, 设计准确度大大提高。

不足之处在于，现在的计算软件和计算方法还有一定的局限性, 需要与传统设计方法配合使用, 但随着流场计算方法的完善, 它将成为今后水力机械主要的设计手段之一。

3.2基于CFD的水力稳定性研究
3.2.1 水力机械内部压力脉动分析
（1）计算方法和湍流模型
近年来，应用最广泛的湍流解法基于雷诺时均N-S 方程（简称RANS方法）。

常用的湍流模型包括k-ε双方程模型、k-ω双方程模型、Reynolds 应力模型和代数应力模型等。

RANS 方法的最大好处是其数值计算量小，计算效率高。

但是，由于这种方法要在时间域上对N-S 方程中的瞬态物理量做平均处理，因此在求解非定常流动问题时遇到一定困难。

解决瞬时湍流物理量计算的最直接方法是直接数值模拟（DNS）方法，即直接用瞬时的N-S方程对湍流进行计算。

但是，试验测试表明（Versteeg，1995），在高Reynolds 数的湍流中可能包含尺度为10 ~ 100μm的涡，在一个0.1× 0.1m2 的流动区域内，要描述所有尺度的涡，计算的网格节点数将高达109到1012量级。

同时，湍流脉动的频率约为10kHz，这就要求将时间的离散步长取为100μs以下。

对于这样的计算要求，现有的计算机能力还是比较困难的。

此外，即使能真正得到这些细节，对于解决实际问题也不一定有太大的意义。

这是因为，从工程应用的观点上看，重要的是湍流所引起的动量、质量、能量及其它物理量的输运，是整体的效果。

正是基于这一特性，数值计算可以考虑放弃对极小尺度涡的模拟，只将比临界尺度大的湍流运动通过瞬时N-S方程直接计算出来，而小尺度涡对大尺度涡运动的影响则通过一定的模型在针对大尺度涡的瞬时N-S方程中体现出来（王福军，2004）。

这就是目前备受关注的大涡模拟法（Large Eddy Simulation，简称LES方法）。

近几年内，由于计算机硬件水平有了很大提高，大涡模拟法开始在非定常计算中发挥作用。

在分析水压力脉动时，LES 方法直接求解瞬态N-S 方程，从而得到较广泛的应用（Wang，
2006；Jiang，2007）。

采用大涡模拟理论进行水轮机的流动计算始于美国学者Song（1996），他采用Smagorinsky的涡粘性公式模化亚格子雷诺应力，首次对水轮机主要的过流部件进行了计算。

1999年，杨建明对涡粘性系数中的湍动能和耗散率采用湍流时均处理和k-ε模型中的概念和模化方法处理，建立了基于大涡模拟思想而方程结构上类似于时均k-ε模型的大涡模拟-双方程模型，对水轮机尾水管和转轮内的流动进行了计算，计算结果与试验结果接近。

Chen等（2000）使用大涡模拟方法计算了具有蜗壳、固定导叶、活动导叶和转轮的水泵-水轮机内的三维流动。

Thomas和Albert等（2002）利用改进的k-ε模型和大涡模拟方法，对尾水管中涡带进行模拟，并对两种模型得到的结果进行比较。

Byskov等（2003）采用大涡模拟方法对六叶片离心泵叶轮的内部流场进行了研究，分别计算了设计工况和几种偏离设计工况的流场，结果表明，设计工况下的流场预测并没有发现明显的分离流动，但在偏离设计工况下，叶轮流道内部的流场与设计工况差别很大，发现了流道内的延迟现象，并将计算结果与RANS时均的计算结果和试验值进行了比较，发现RANS时均方法并没有预测出这个叶轮流道内的延迟现象。

Nagahara等（2005）对多级离心泵内湍流进行了大涡模拟，并将计算结果与RANS计算结果和试验值进行了比较，指出大涡模拟方法对非定常湍流及压力脉动的预测结果与试验值更加接近。

Yamanishi（2007）、Chisachi（2006）和Liu（2005）等也采用LES 方法对水力机械进行了流场计算，取得了不错的结果。

LES方法在模拟含有较大尺度的涡旋流动时，优势非常明显，但如果将其应用到含有较小尺度的涡旋流动时，其数值计算的网格必须足够小才能得到所需要的计算精度（Davidson，2006），因此很难将它用于近壁区的低雷诺数流动。

这样，近壁区流动的数值模拟均不采用LES 方法，多数情况下采用形式简单的壁面函数来近似处理。

而在分析水力机械压力脉动所导致的结构动态特性时，壁面上的压力场分布又极其重要，因此，在LES方法应用的问题上遇到了困难。

URANS（Unsteady RANS，简称URANS）方法的成熟改变了雷诺平均等同于时间平均的概念（Wegner，2004）。

出现了将LES 方法与URANS 方法耦合求解核心区流动和近壁区流动的新方法（Davidson，2006）。

由于LES方法与URANS 方法对湍流的描述模式不同，如何建立一种界面强制模型将LES方法和URANS方法结合在一起来解决水力机械流场压力脉动的问题也是该方法的应用难点之一。

（2）并行算法
并行算法作为一种数值计算的工具而非流动数值模拟的模型，在计算量庞大的流动计算问题上必然会得到越来越广泛的应用，在求解复杂而且巨大的水力机械内部流场时，通过把
流动区域分成若干块子区域，子区域间通过公共边界上节点信息的耦合条件进行相互约束和交换，以实现相邻区域解的光滑过渡，从而实现复杂区域整体流场的并行计算，这将提高计算的效率。

目前在国内很少有人能够对水力机械过流部件内的流场做整体计算，一般都是各部件单独计算，单独给出边界条件，而这些边界条件很难从物理上准确提出。

在采用LES方法的前提下，仅就一个工况下的水轮机/水泵压力脉动非定常计算而言，使用目前配置最高的微机也要数天才能完成。

因此，在引入其他一些界面模型（如LES方法和URANS 方法结合），采用并行计算模式实属必然。

结合水力机械内部流场的计算方法和湍流模型，建立并行计算网，提高计算效率，是推动水力机械非定常研究的必要手段。

（3）压力脉动的研究
水力机械内部压力脉动特性的数值研究集中在水轮机不同过流部件压力脉动特性、关系以及尾水管内压力脉动与涡带特性分析。

Vatech和Sulzer等于2002年对混流式水轮机和可逆式水轮机的尾水管进行非定常解析，并将其结果与试验对比，发现计算得到的涡带所引起的压力脉动幅值和频率与试验结果非常相似。

Fan和Kubota等（2006）利用雷诺应力模型对混流式水轮机弯肘型尾水管内流动进行了长时间非定常计算，研究尾水管内死水区域与涡带的不规则运动，预测尾水管的压力脉动并分析其频率特性。

Thi等（2006）对水轮机内部动静翼型的干涉、部分负荷下尾水管内涡带、固定导叶后的卡门涡列等进行了非定常流动模拟，并结合模型、原型试验数据评估压力脉动数值预测的精度。

Braun（2005）、Zhou（2007）和Wang （2008）等也模拟了水轮机内部的流动，取得了水轮机计算领域很普遍的计算结果。

国内清华大学吴玉林等（2000，2004）预测了混流式模型水轮机尾水管内压力脉动，分析引水部件内压力脉动的特征。

肖若富等（2004）对尾水管内的低压涡带进行数值模拟，分析低压涡带的形成和运动规律。

王正伟和周凌九等对水轮机进行了非定常流动计算，转轮进口压力脉动主频为叶片转频，而蜗壳进口主频为叶片转频的2倍（2001）；分析了叶片表面和尾水管内压力脉动，结果表明：在部分负荷工况下，由尾水涡带引起的转轮内压力脉动主频是转频与涡带压力脉动频率之差；在同一工况下，叶片背面的压力脉动幅值一般大于叶片正面的值（2006）。

刘树红等（2004，2005）对水轮机整机压力脉动进行模拟发现，转轮前压力脉动主频为叶片转频，而在其他引水部件压力脉动主频为尾水管内涡带频率。

杨昌明等（2008）对混流式水轮机进行非定常计算，发现转轮进口处的压力脉动幅值达到水轮机水头的16.5％。

朱嵩等（2006）计算获得了分离流、二次流和多尺度旋涡结构。

从计算方法、湍流模型和物理计算域出发的研究都是为了提高CFD在水力机械领域的应
用，提高计算精度和计算效率。

但事实上，从近几年国内外发表的水力机械非定常流动数值模拟的文献来看，数值模拟在一定程度上可以预测压力脉动的频率与幅值，但预测幅值与试验存在较大的差异，离工程实用还有一定的距离，水力机械压力脉动的数值预测精度还有待提高。

3.2.2 水力机械流固耦合问题的研究
对复杂的流固耦合系统进行动力学分析，大致有二类方法：一类解析－数值法，即对结构采用有限元离散，而对流体则采用某种近似的解析关系式来式来描述；另一类是纯数值法，即对结构和流体都采用有限元离散，或采用有限元和边界元离散。

由于流固耦合问题的复杂性，求解主要立足于数值分析，尤其是有限元法。

研究者们在应用有限元法时通常采用三种途径研究流固耦合问题：①选择欧拉坐标系，在结构中以位移为未知量，对流体以压力或速度势为未知量。

由于两种介质采用不同的变量，不便应用通用程序；②选择拉格朗日坐标系，对结构和流体均以位移为未知量，沿流、固交界面上的每个结点的相容和平衡条件可自行满足，且可得到对称的总刚矩阵，易利用如SAP5、SSAP、ADINA等有限元分析通用程序；③选择任意拉格朗日—欧拉坐标法 （ALE坐标系），该坐标系可以任意速度在空间运动，若其速度为零，则为欧拉坐标系；若其速度等于质点速度，则为拉格朗日坐标系，可处理两相界面的协调及自由面问题。

早期的拉格朗日流体位移元存在零能频率及伪振型，从而引起数值解的不稳定性，影响该方法的广泛应用。

在近期研究中，采用了一些措施，如单元数值积分的降阶，改进流体元的插值函数等，克服了上述困难。

采用Wilson非协调元中的位移插值函数，能更好地模拟流场的真实变形，消除早期拉格朗日流体元存在的缺陷。

数值方法在研究流固耦合问题时，首先遇到的是建模问题。

国内外学者提出了多种数值计算模型的格式和方法，比如特征值算法、模态叠加法、零剪切格式、Us-P格式、附加质量法、Us-U f格式、Us-Ψ格式、Us-Φ-Po格式、Us-V f等格式。

在对水力机械进行考虑流固耦合的自振特性分析时，一般假设水为无粘、无旋的介质，且流体与结构分别作小扰动和小变形运动。

场变量的选择，对结构通常采用位移法，对流体则采用质点位移、压力或速度势。

国内外学者在这些方面作过一些研究，Ohashi等（1991）曾介绍了液固耦合的一种模型，但在考虑液体运动方程时，对流体对流项不加考虑，这种近似对转轮叶片内的液流是不合适的。

Liu等（1984）也提出一种计算附加质量的方法，在总质量矩阵中组装附加水质量的方法来近似考虑流体对结构的作用。

这种耦合的方法只增加了系
统的质量而不改变系统的刚度，同时没有考虑结构本身的外形，也具有一定的局限性（陶敏,1997）。

瞿伦富等（1998）以采用位移一速度势作为连续场变量的计算格式，首次在方程中考虑流动的三维复杂性，使方程更能符合液固耦合运动的本质，假设：1）液体无粘性，绝对运动有势，不可压缩及小的流场干扰；2）叶片在弹性变形范围内，不计叶片本身阻尼。

肖若富等（2001）通过振动相似分析获得结构在空气中和液体中的固有频率关系，计算转轮在空气中的固有频率和主振型，进而计算转轮在水中的固有频率。

这与真实情况有较大的误差，因为结构在水中固有频率的降低程度与许多因素有关，不能简单地乘以一个系数来计算得到。

此外，通过这种方法不能得到结构在水中的主振型。

谷朝红等（2001）把有限元法与边界元法相结合计算水轮机部件流固耦合振动问题，编制了结构在水下流固耦合振动的通用计算程序。

梁权伟等（2003,2004）采用全耦合的有限元法对某电站混流式转轮进行了模态分析，比较了转轮在空气和水两种介质中的前7阶固有模态，发现水介质对转轮不同阶次振动频率的影响各不相同，但转轮在两种介质中的同阶振型相似。

郑晓波等（2005）运用采用流固耦合技术，用有限元法求解流体与固体的耦合方程,对轴流式叶片进行了振动特性分析,并比较了轴流式叶片在空气中与水中固有频率的变化。

Rodriguez等（2006）研究发现，对于结构振幅较小时，结构振动对于流体域的影响一般较小，这时可将流体一结构的耦合效应简化为系统的附加质量、附加阻尼及附加刚度加以考虑。

曹良，张立翔（2008）针对混流式水轮机提出了一种新的采用附加质量法的简化计算其部件流固耦合效应的附件质量模型，并计对某机组固定部件的动力特性进行了计算，结果证明提出的附加质量模型不仅可以大幅度地节省计算量，而且能得到高精度的计算结果。

基于计算理论和计算机技术的迅速发展，水力机械流固耦合动态响应特性的数值研究向全三维和粘性方向发展。

水力机械流固耦合计算的研究集中在流场压力脉动的计算、流固耦合界面模型的研究以及固体域的计算三方面。

水力机械流固耦合的研究有两种模式，分别是强耦合模式和松耦合模式。

水力机械松耦合模式的流固耦合计算与强耦合计算相比，计算大大简化，计算量大大减低，但从理论上来看，强耦合计算模式考虑更周全，计算更精确。

强耦合流固耦合数值计算要考虑流场载荷在结构表面的加载和结构振动引起的流场边界变化。

Zhang等（2007）通过流-固系统的广义变分原理建立了描述强耦合流激振动的控制方程，采用Wilson-θ方法分析结构振动，采用LES方法计算流场，依次单独迭代计算流场和结构振动，同时记录水轮机叶片表面的压力和振动加速度。

计算结果表明，叶片振动、叶片外形的弯曲以及流动经过活动导叶产生的尾流对转轮流场特性随时间和空间的变化影响很大，。