水轮机部分水轮机动静干涉的非定常流动数值模拟肖惠民1,2, 于波1,2（1. 武汉大学动力系，武汉 430072；2. 水力机械过渡过程教育部重点实验室，武汉 430072）[摘 要]在水轮机领域，流动分析技术已经成为设计开发的常规工具。
随着水轮机单机容量和尺寸的加大，其运行稳定性日益受到重视，这使得水轮机内部非定常流动的数值计算成为必要。
本文讨论了动静干涉、移动网格的不同数值计算方法，并应用“滑移网格技术”对混流式水轮机进行了非定常流动计算，分析了部分负荷工况下的流态和压力脉动。
[关键词]水轮机；非定常流动；动静干涉；移动网格；数值模拟[中图分类号] TK730.2 [文献标识码] A [文章编号] 1000-3983（2009）05-0033-05 Numerical Simulation of Unsteady-flow in Hydraulic Machinery Induced by Rotor-stator InteractionXIAO Hui-min, YU Bo( 1. Department of Power Engineering, Wuhan University, Wuhan 430072, China;2. Key Laboratory of Transients in Hydraulic Machinery, Ministry of Education, Wuhan 430072, China )Abstract: In the field of water turbine computational fluid dynamics is routinely used today in research and development as well as in design. A special attention is paid to the water turbine stability with the increasing of unit power and size. In order to get solutions to the stability problem an unsteady flow analysis is necessary. Different numerical schemes are discussed for the rotor-stator interaction and moving grid. The rotor-stator coupling by application of sliding mesh is shown on the example of a complete Francis turbine. Computational results are analyzed for part load point of operation.Key words: water turbine; unsteady flow; rotor-stator interaction; moving grid; numerical simulation1 引言目前，计算流体动力学（CFD）技术已广泛用于水轮机的设计和研究中。
水轮机单一部件的定常（稳态）求解，如蜗壳、串联叶栅、转轮和尾水管等，已用于日常的工作中。
而且还经常进行多部件或整个水轮机的定常耦合模拟。
然而，在水轮机中还存在着由各部件相互作用引起的非定常流动问题，为解决这些问题就必须采用非定常流动分析，即需要进行与时间相关的数值计算。
非定常流动可分为两大类[1]。
第一类是由于外界引起的非定常流，即由于非定常边界条件或流动域几何形状随时间变化引起的非定常流，如导叶、阀门的开启和关闭，动静干涉等。
第二类是由于自激引起的非定常流，如紊流运动、旋涡脱落（卡门涡街）、不稳定涡（如尾水涡带）等。
此时边界条件、流动域的几何形状并没有变化，但是产生了非定常流。
以上两类流动还可能联合作用引发非定常流，如流动诱发振动，并引起流动域几何形状的变化。
所以这些现象都有可能在水轮机和水泵中出现，并需要采用不同的算法进行求解。
本文讨论了第一类非定常流动中动静干涉的数值计算方法，并应用“滑移网格技术”对某混流式水轮机进行了非定常流动计算。
2 动静干涉动静干涉指旋转部件和固定部件间流动的相互干扰。
以水轮机为例，在活动导叶和转轮之间、转轮和尾水管之间就存在着动静干涉。
转轮旋转切割导叶出口的尾迹，导致转轮内的流动为非定常流；转轮出口的非轴向水流在尾水管内形成旋转水流，严重时形成螺旋状涡带。
水轮机中这种相互干扰引起的不稳定流可能诱发水力振动，引起疲劳破坏，甚至在导叶、转轮叶片或尾水管管壁产生裂纹，影响机组的安全运行。
2.1 数学模型目前，模拟动静干涉的数学模型从简单到复杂有以下四种[2]：（1）简单的“顺序”计算。
每个部件在各自的参考系中求解，将上游部件的结果作为下游进口边界条件。
但部件间不存在干扰，计算结果完全取决于进口边界条件。
（2）级平均法。
在转动部件和非转动部件间的交界面上，传递所有变量沿圆周方向上的平均值。
（3）冻结转动部件。
计算时，转轮位置相对导叶固定，因此必须计算转轮相对导叶的几个不同位置，从而获得有用的结果。
（4）非稳定计算。
直接模拟各个部件之间的相对运动，动静部件间的界面随着时间而更新。
通过交界面传递通量，实现每一时间步长内部件间的流场耦合。
非稳定计算最接近水轮机的真实流动。
但由于蜗壳内水流的非对称性，以及随后流到串联叶栅和转轮内的水流也呈非对称性，再加上导叶和转轮节距不等，就形成这样一个事实：不存在周期性的水流，因此计算时必须考虑水轮机转轮和串联叶栅的所有流道。
以混流式水轮机为例，导叶数一般在18～24之间，转轮叶片数在7～15之间。
若要对整个流道进行三维模拟，就需要上百万的计算网格，这将带来相当大的计算量，需在高效、高性能大容量计算机或并行计算机上完成，而且要耗费相当长的时间。
与定常流相比，非定常流计算工作量为定常流的5～10倍，倍数因问题和所采用紊流模型的不同而异。
随着计算机技术及CFD 技术的发展，现在已可以进行非定常流动分析[1]。
2.2 移动网格 在非稳定计算的过程中，旋转部件（或移动部件）的网格随时间步长的递增，按转轮的旋转速度（或固壁的移动速度）同向移动，从而实现旋转部件网格与固定部件网格的相互移动。
动静干涉移动网格的模拟可用不同的方法解决，图1为其中的3种方法[3]。
（1）捕捉单元。
交界面上的网格单元因为网格移动而变形。
变形增大直到与下一个单元处于同一水平位置，此时连接单元变形最小，这意味着捕捉到了这个单元（图1a ）。
这种方法的优点是不需要插值，从而保持求解的精度。
缺点是要求界面两侧采用相同的网格，以获得连接单元，并避免单元变形太大，从而使网格生成比较复杂，特别是对于复杂几何体网格生成几乎是不可能的。
（2）动态边界条件。
独立计算旋转和固定部件，由边界条件实现流场耦合，边界条件在迭代过程中不断更新（图1b ）。
对上游节点值插值得到下游节点值（速度和紊流量）。
对下游边界单元的压力、动量通量和紊流通量积分，并传递到上游边界单元。
为获得精确的插值和积分结果，采用重叠网格，以确保下游部件的所有节点一直在上游网格内，整个上游边界完全在下游网格内。
例如，转轮网格延伸到靠近活动导叶出水边，双列叶栅网格延伸到靠近转轮进水边。
转轮与尾水管的边界也用类似方法处理。
（3）滑移网格。
滑移网格法允许交界面两侧的网格相互滑动，且交界面两侧的网格不重叠、无须对齐。
但需计算交界面两侧的通量，并使其相等。
计算时，在每一个新的时间步长内，根据交界面两侧网格的相对位置求取新的重合面，再通过重合面传递通量，从而实现每一时间步内流场的耦合。
实现过程如图1c 所示。
图1c 中A-B 、B-C 和D-E 、E-F 构成交界面，d-b 、b-e 和e-c 构成交界面上的重合面，而剩余的a-d 和c-f 形成周期或壁面边界。
若要计算通过重合面从单元I 和III 到单元IV 的通量，不使用D-E 面，而使用d-b 和b-e 两个重合面。
本文采用“滑移网格”法实现通量的传递，导水机构－转轮、转轮－尾水管间的交界面分别位于转轮进口前和转轮出口后。
且因进行的是全三维的流动计算，交界面为完全重合的面，因此没有剩余边界面生成。
图1 移动网格示意图捕捉单元时间步长1时间步长2时间步长3流动方向（a ）捕捉单元（b ）动态边界条件部件1流动方向部件2移动网格单元区域1交界面区域1交界面区域2单元区域2（c ）滑移网格流动方向IIIIII IVVVIab d ec f ABCDE F3 数值方法3.1 物理模型本文研究是一个混流式模型水轮机，由蜗壳、23个固定导叶、24个活动导叶、带16个X 型叶片的转轮及窄高型尾水管组成，如图2所示。
水轮机每个过流部件的网格单独生成，采用非结构化网格，网格节点数共约40万个，网格单元数共约156万个。
图2 计算区域3.2 紊流模型和离散方法计算基于不可压缩雷诺时均Navier-Stokes 方程，紊流计算采用RNG ε−k 模型（基于重整化群理论的ε−k 模型）[4]：M b k jeff k j Y G G x kx Dt Dk −−++∂∂∂∂=ρεμαρ)( （1） 1322()() eff k b j j D C G C G Dt x x kC R kεεεεεεεεραμερ∂∂=++∂∂−−（2） 采用贴体坐标下的有限体积法和非交错网格对上述方程进行空间离散，时间离散采用二阶完全隐式格式。
压力项采用二阶中心差分格式，其他项采用二阶迎风差分格式。
采用SIMPLEC 算法，实现压力和速度变量的分离求解。
根据模型试验水头设置进出口边界条件，即给定蜗壳进口压力和尾水管出口压力。
在壁面处采用无滑移边界条件。
初始条件（计算初值）由级平均法的稳态计算获得[5]。
4 结果分析对不同开度、不同水头下的九个运行工况进行了非定常计算。
下面就高部分负荷工况（额定水头、56 %max a 、78 %r P ）的计算结果做详细分析。
4.1 流态分析图3为串联叶栅和转轮在某一时间步的瞬态压力分布。
由图可见，从蜗壳到转轮出口压力逐渐下降；在转轮进水边可以清楚地看到滞流点；由于受蜗壳内不均匀水流分布的影响，固定导叶进口周向压力分布也不对称；从活动导叶上的压力分布可看出，沿导叶高度方向压力变化很大，顶部的压力非常高，而底部的压力较低，其原因是由于水流从径向向轴向的突然转变造成的。
图4为串联叶栅中间剖面在某一时间步的瞬态速度分布。
从图中再次观察到由蜗壳引起的串联叶栅不对称入流。
