考虑自由液面的灯泡贯流式水轮机数值模拟研究

范志锋，林巧锋，黄叶雯，芦 月，陈齐灯，郑 源

（1．福建华电电力工程有限公司，福建省福州市 350001；2．中国华电集团有限公司福建分公司，福建省福州市 350001；3．河海大学，能源与电气学院，江苏省南京市 211100；4．河海大学，水利水电学院，江苏省南京市 210098；5．福建华电金湖电力有限公司，福建省将乐县 353399）

0 引言

贯流式机组在低水头、大流量及潮汐能开发利用等方面具有其他机组无法比拟的优势[1-3]，其中，灯泡贯流式机组因具有结构紧密、水力损失小及建设周期短的特点，近年来发展迅猛、应用广泛[4-6]，对机组性能和稳定性的研究得到越来越多的重视[7-10]。目前研究贯流式机组水力性能的方法主要是试验和数值模拟[11]，但由于贯流式机组流道较短、水体惯性占比较大，模型试验测得的运行参数与真机误差较大，而现场试验又费时费力，无法反映机组内部流态，因此，数值模拟是研究贯流式机组水力性能的主要手段。

此外，由于贯流式机组运行水头低，纵横尺度相当，水头与转轮直径之比较小，所以上下游水库中的自由液面以及水体重力因素对机组内部流动的影响较其他类型水轮机更明显，研究时要考虑自由液面对其水力性能带来的影响[12-13]。在自由液面方面，国内外学者开展了相关研究，周大庆[14]等采用VOF模型对抽水蓄能电站水轮机工况进行数值模拟，结果表明：基于VOF模型的数值计算具有较好的计算精度，能更真实地反映机组的实际运行条件；Benchikh等[15]基于VOF方法对Darrieus水平轴水轮机进行了多相数值计算，结果表明：多相模型模拟结果与试验相比误差更小，且能更真实地预测机组功率和转矩系数；赵亚萍[16]等对三叶片卧式贯流式水轮机进行数值模拟，在模拟时通过求解体积函数的输移扩散方程来确定自由水面的位置，得出需要考虑自由液面才能满足几何相似及流动相似的特性；邵方凯[17]对漂浮式垂直轴水轮机建立了在自由液面下周期性旋转模型，结果表明：自由液面的存在会在低速比下增加机组功率，高速比下则会降低。可见，自由液面对水力机械数值模拟的结果存在不同程度的影响，而此类研究中对灯泡贯流式机组的探讨较少，且大多没有考虑水体重力和上下游水池带来的影响，少数考虑到此因素的研究也只是把水池表面设置为无滑移壁面[18-21]。

基于上述研究现状，本文以福建某灯泡贯流式水轮机为研究对象，构建包含上下游水池在内的全过流系统三维模型，在考虑水体重力的基础上，针对水池表面与大气相通的特点，结合VOF模型对不考虑自由液面和考虑自由液面两种情况分别进行数值模拟，以更符合机组实际运行状态，并从能量特性、不同测点压力脉动及内部流场特性等方面对比单相工况和两相工况的差异，揭示上下游水池自由液面对灯泡贯流式机组水力性能的影响。

1 计算模型与计算方法

1．1 机组参数与模型试验

本文以福建某灯泡贯流式机组模型装置为研究对象，该机组水轮机型号为GZTF07B-WP-560，额定功率为21MW，额定转速为93.75r/min，额定水头为9.3m，最大水头11.9m，最小水头4.7m，额定出力为21.875MW，额定流量为255.86m3/s，叶片数为4，活动导叶数为16。该机组的模型水轮机试验是在哈尔滨大电机研究所的高水头试验I台上进行的。试验台模型效率综合测试误差小于±0.25%，模型效率重复测试误差小于±0.10%。模型的效率、空化、压力脉动等试验能在同一试验台同一模型上进行，试验台主要测量仪器设备见表1所示。

表1 模型试验台主要测量仪器设备Table 1 Main measuring instruments and equipment of model test bench

模型试验中的能量试验是在水头为常数、无空化条件下进行的，每个工况点采集时间为60s，试验包括整个水轮机运行范围：模型试验水头为4～2.5m，模型活动导叶开度为24°～90°，桨叶开度为5°～45°，根据模型与原型换算公式，可得到原型水轮机能量试验结果，其中，在额定工况下，即水头9.3m、导叶开度为68.5°、桨叶开度为36.1°时，稳定运行时对应的机组流量为260.19m3/s，转轮力矩为2086.2 kN·m，整机效率为86.2%。

1．2 计算域及网格划分

对该机组进行额定工况下三维数值模拟计算，利用UG软件对水轮机全流道进行同等尺寸建模，三维模型包括上游进水池、进水流道、导叶、转轮、出水流道、下游出水池等六部分构成，上下游水池高度差为额定水头，进出水池均为三面过水，图1为灯泡贯流式水轮机全流道结构示意图。在考虑自由液面数值模拟时，仅在上下游水池的顶端加入一定高度的空气域，空气层高度与液体高度接近1：1。

图1 灯泡贯流式水轮机全流道结构示意图Figure 1 Structural diagram of the whole passage of bulb tubular turbine

本文利用ICEM软件对模型进行网格划分，在不考虑自由液面模拟时，采用非结构化四面体网格对进出水池、进出水流道进行划分，对流态复杂的导叶和转轮区域进行网格加密。在考虑自由液面模拟时，由于气体和水体在密度上有较大的差异，为了能较为清晰地捕捉两相介质在自由液面处的相界面，需要对上下游水池自由液面交界面进行网格加密，因此在上下游水池部分采取结构化网格进行划分，图2为考虑自由液面水轮机各部分网格划分示意图。

图2 考虑自由液面水轮机各部分网格划分示意图Figure 2 Schematic diagram of grid division of each part of hydraulic turbine considering free surface

在不考虑自由液面的额定工况模拟中进行网格无关性验证，所以现采取5种不同数量的网格划分情况，总数依次为550万、625万、790万、845万以及910万，对机组进行额定工况下的稳态数值模拟计算，具体计算网格数及其效率如表2所示，从表中可以发现网格数超过845万之后机组的效率变化已经很小，出于对时间成本和计算精度的考虑，最终选取方案4的网格数进行后续计算。

表2 不同网格划分方案及机组额定工况效率模拟结果Table 2 Simulation results of different grid division schemes and unit efficiency under rated working conditions

1．3 考虑自由液面VOF模型

本文基于Fluent19.2软件平台VOF多相流模型进行三维湍流数值模拟计算，VOF模型是建立在固定的欧拉网格下的表面跟踪法，其特点在于通过建立体积函数来表征流体所占单元的体积分数，以此创建和追踪自由表面。在单元中，若q相流体体积分数为αq，那么可能存在三种情况，一是αq为0即单元里不存在第q相流体，二是αq为1即单元里充满了第q相流体，三是0＜αq＜1即单元里包含了第q相流体和其他单相或者多相流体的界面。在本文中，分别定义αa和αw为贯流式机组上游水池和下游水池内气体和液体所占的体积分数，则气体的体积分数可用液体的体积分数表示，即：

其中，液体的体积率函数αw的控制方程表达式为：

引入VOF方法的k-ε湍流模型与单相流的k-ε湍流模型形式完全相同，只有密度ρ和动力黏度μ的具体表达式不同，在此系统中的物性参数均由体积率来控制，其表达式如下：

式中：ρw——液体相的密度；

ρa——空气相的密度；

μw——液体相的分子黏性；

它以特殊铝材为材料，可应用于几乎所有行业的输送机结构中，运输绝大多数种类的货物。它被组装在OTU换向台上，并配有电源、B L D C电机和控制器，可根据特定的应用程序进行设置。

μa——空气相的分子黏性。

通过对液体的体积分数αw的求解，密度ρ和动力黏度μ即可通过式（3）和式（4）得出。

1．4 计算工况及边界条件

不考虑自由液面：选取RNGk-ε湍流模型对控制方程进行封闭，在计算中方程组压力项为二阶中心差分格式，速度项、湍动能和湍动能耗散率均为二阶迎风格式，求解方法选定为SIMPLEC算法，转轮区域设置为滑移网格，固定旋转速度为额定转速93.75r/min，固体壁面设置为无滑移边界，临近固体壁面处采取标准壁面函数，不考虑壁面粗糙度对流场的影响。最大迭代步数为20步，收敛残差为10-5，设置叶片旋转每三度所经过的时间作为步长，约为0.005s。上游水池进口设置压力进口，下游水池出口设置压力出口，考虑到重力场对水头较低的贯流式机组内部流态的影响，通过Fluent中的用户自定义函数设置上下游水池进出口压力沿水池高度变化。

考虑自由液面：打开VOF模型，添加明渠，上下游水池进出口均给定液面高度，水池顶部空气域设为开放面，空气层区域水相分数为0，即充满空气；其余各部分水的体积分数设置为1，其他条件与不考虑自由液面时相同，考虑自由液面后的流体初始状态如图3所示。

图3 考虑自由液面水轮机初始流场Figure 3 Initial flow field of hydraulic turbine considering free surface

2 计算结果对比与分析

2．1 能量特性

对于不考虑自由液面和考虑自由液面的工况分别进行三维数值模拟，从计算开始至16s后计算收敛，各监测数据逐渐趋于稳定值，选取稳定后相同计算时间内的结果进行对比。表3所示为不考虑和考虑自由液面水轮机的能量特性对比，从表中可以发现，在机组水头和转速相同的条件下，不考虑自由液面时，机组的流量较小、转轮力矩较大，效率达到91.1%；考虑自由液面时，机组的流量升高了1.79%，转轮力矩减小了6.65%，效率降低了7.5%，且与模型试验的结果更接近，误差更小。这说明贯流式机组在实际运行时，受到重力场和自由液面的影响，机组过流量有少许增加，而转轮出力有大幅下降，机组效率有明显降低。这种差异会影响贯流式机组的初期设计，导致机组出力不足与设计值误差较大的情况发生。

表3 不考虑和考虑自由液面水轮机的能量特性对比Table 3 Comparison of energy characteristics of hydraulic turbines without and without free surface

2．2 压力脉动规律

为探究贯流式机组流道内各位置压力变化情况，分别在导叶前、转轮前及尾水管进口、出口四个截面上布置若干压力监测点，考虑重力场影响，在每个截面外缘位置以及截面上同一纵向位置分别布置三个监测点，监测点布置图如图4所示。为了使计算结果更准确，采样时间为12个周期，选取最后两个周期进行压力脉动分析。压力脉动的时域值是通过数值计算结果直接获得的，而频域值则需要进行傅里叶变换才能获取。本文取无量纲的压力脉动系数Cp来反映各过流部位的压力脉动特性，具体公式为：

图4 水轮机压力脉动监测点布置图Figure 4 Layout of pressure fluctuation monitoring points of hydraulic turbine

式中：pi——监测点在某一时刻的压力值；

通过对考虑自由液面模拟下各监测点压力脉动的初步分析，发现各监测截面上，外缘位置监测点的波动情况无明显差异，因此，这里只对比截面上同一纵向位置上的三个点，图5为考虑自由液面的两相工况各截面纵向分布的三个点的压力脉动频域图，从频域特性上看，导叶进口处的压力主频均出现在四倍转频fn处（4fn），且压力脉动系数Cp值都很小，纵向三个点中位于中间位置的A6’点Cp值最高为0.0013，其次频Cp值为主频幅值的38%，除此以外出现了很多幅值频率不定的干扰信号，说明导叶进口在一定程度上出现了低频压力脉动，这可能是由于导叶进口处靠近进水池，水位的波动带来的影响，但整体数值较小，流态相对平稳；转轮进口处的压力主频也在4fn处，在8fn、12fn、16fn处出现了不同程度的谐波，其压力波动情况在四个截面中最强烈，三个点中外缘位置的B4′点Cp值最高为0.258，这是因为此位置距离转轮最近，受到导叶和转轮的动静干扰；尾水管进口处的压力主频也在4fn处，外缘位置的C4′点Cp值最高为0.023，同时其高频区域起伏也最大，倍频信号幅值呈减小趋势，说明水流在进入尾水管时对壁面造成了一定的冲击，三个点中越靠近中心位置，谐波也逐渐消失，表明在尾水管进口处轴心区域流态较稳定；尾水管出口处的压力频率主要在fn处、0.78fn、0.44fn处，这是由于此位置距离转轮区域较远，受到转轮旋转的影响较小，此外，该截面低频压力脉动占主导地位，其中D4′点主频Cp值最高为0.056。

图5 两相工况各截面纵向点压力脉动频域图Figure 5 Frequency domain diagram of pressure fluctuation at longitudinal points of each section under two-phase working condition

将不考虑自由液面和考虑自由液面的压力脉动频域进行比对，选取两种工况下各截面上波动幅值最高的对应点进行对比，如图6所示为不考虑和考虑自由液面各截面压力脉动频域对比图，从图中可以发现四个截面处单相和多相下的压力脉动频域变化规律基本相同，且主频均在同一位置，不同的是，气液两相工况整体在频域图中出现了很多幅值频率不等的干扰信号，杂频较多，特别是靠近进水池的导叶进口处和靠近出水池的尾水管出口处最为明显，这是因为考虑自由液面后机组运行受到空气的干扰，液面的变动需要经历一段时间后才能逐渐趋于稳定。此外，考虑自由液面后四个截面压力脉动的主频幅值均增大：导叶进口处的两相工况主频Cp值为单相时的1.34倍；转轮进口处的两相工况主频Cp值为单相时的1.27倍，谐波幅值无明显差异；尾水管进口处的两相工况主频Cp值为单相时的1.23倍，且两相工况在4fn之前区域出现了部分杂频；尾水管出口处的两相两相工况主频Cp值为单相时的1.43倍，且均为低频压力脉动。综合对比图来说，考虑自由液面时，机组上下游水池液面受到空气层的干扰，水相和气相交互引起水池内流态波动，这种波动随后进入流道内影响转轮旋转，继而影响整个机组的压力脉动状况。

图6 不考虑和考虑自由液面各截面压力脉动频域对比图Figure 6 Frequency domain comparison of pressure pulsation of each section without and without free surface

2．3 机组内部流态

为了进一步分析加入自由液面后内部流场水流分布情况，对进出水池自由液面的水相进行了对比分析。图7为两相工况进出水池自由液面处的水相分布。从图中可以看出，进出水池自由液面水相分布均为对称分布，不同的是进口处水相体积范围明显比出口多，最高值分布在边缘处，一定程度上呈现梯度变化，而出水池则没有了规律，最高值分布在中心区域。

图7 进出水池自由液面处水相分布图Figure 7 Distribution of water phase at the free surface of inlet and outlet pool

如图8所示为不考虑自由液面和考虑自由液面的机组纵向剖面压力分布云图。可以发现，在重力作用的影响下，进出水池的水体部分压力均呈现沿水深增加的梯度变化规律，且转轮前区域的压力整体都比转轮后区域的压力数值大，在转轮区域两个工况压力沿轴向都有一定的变化；考虑自由液面后，进出水池的底部高压区域变大，且与空气层交界处位置的压力与不考虑自由液面时相比增大2.4%，说明液相明显受到了气相的影响，此外，空气层部分压力不沿梯度变化，与实际情况相符。这是因为不考虑自由液面时进出水池被假设为密闭容器，水池顶部设置为固体壁面，而考虑自由液面后进出水池直接与大气相通，水池顶部设置为开放面，自由液面的存在增加了水池顶部的初始压力，因而造成水池内压力随水深的增加而增大，整体大于不考虑自由液面时的水池压力。

图8 纵向剖面压力分布对比云图Figure 8 Comparison cloud chart of pressure distribution in longitudinal section

如图9所示为不考虑自由液面和考虑自由液面的导叶进口湍动能分布对比图。从图中可以发现两种工况下的湍动能分布规律基本一致，均为对称分布，且数值相差很小，两相工况略高于单相工况；高湍动能区域都在导叶区域上部靠近灯泡体壁面的位置，低值区域均在流道的中间部位；整个截面上湍动能值从外壁面到灯泡体壁面呈现先减小后增大的变化规律，不考虑自由液面的湍动能等值线较为流畅、过渡均匀，而考虑自由液面后的湍动能等值线有突变情况，且在靠近外壁面处突变更明显，这可能与考虑自由液面后，气相和水相发生碰撞，有能量交换的情况相关。

图9 导叶进口湍动能分布对比图Figure 9 Comparison of turbulent kinetic energy distribution at guide vane inlet

如图10所示分别为考虑自由液面和不考虑自由液面的转轮叶片吸力面和压力面的涡核区域速度分布对比图。从图中可以看出单相工况和两相工况转轮涡核区速度分布规律基本一致，压力面几乎没有涡核分布，吸力面进水侧分布较密集，从轮毂处向轮缘处逐渐增大，在叶片外缘侧尖端速度达到最大值；考虑自由液面后，两相工况叶片外缘侧的速度值高于单相工况，且高速区域范围明显大于单相工况，说明加入气相后使转轮吸力面的涡核速度增加，这是因为自由液面部分气相和液相的相互作用引起了液面的波动，这种波动会随液相的流动作用在转轮部分，而波动产生的能量增加了叶片区域涡核的动量，从而干扰了转轮的旋转，使转轮部分的水力损失增加，这也是造成转轮出力下降的原因。

图10 转轮叶片进出水面涡核区域速度对比图Figure 10 Velocity comparison diagram of runner blade in and out of water surface vortex core area

如图11所示为不考虑自由液面和考虑自由液面的尾水管内速度分布对比云图。从速度云图中可以看出两种工况速度都是呈现对称分布的，在尾水管入口处速度最高，且最高流速相差不大，同时在尾水管中间部分形成低速区域，说明此处有同程度的旋涡产生，不考虑自由液面时低速区域主要集中在中间部分，而加入气体相后低速区域范围扩大且主要分布在尾水管的中部及靠近出口部位，蓝色区域分布更为广泛，这是由于尾水管靠近下游出水池，受到出水池空气层的干扰，两相工况尾水管内流态更为混乱，容易形成扰流和回流，在实际运行时，尾水管部分区域还可能会出现真空情况。

图11 尾水管内速度分布对比云图Figure 11 Comparison cloud chart of velocity distribution in draft tube

3 结语

本文以福建某灯泡贯流式水轮机为研究对象，在考虑水体重力的基础上，采用VOF模型对不考虑自由液面和考虑自由液面两种情况分别对机组设计工况下的稳态运行过程进行数值模拟，探究加入气相后对机组能量特性、压力脉动及内部流场的影响，对比结果表明，考虑自由液面后：

（1）能量特性方面，水轮机的流量比单相工况时提高了1.79%，转轮力矩减小了6.65%，效率降低了7.5%，与模型试验的结果更接近，这种差异会影响贯流式机组的初期设计，导致机组出力不足与设计值误差较大的情况发生。

（2）压力脉动方面，各截面监测点主频没有发生变化，但压力脉动幅值均大于单相工况，靠近水池的导叶进口面和尾水管出口面幅值变化较其他面更大，且出现了很多幅值频率不等的干扰信号。

（3）内部流态方面，进出水池的底部高压区域变大，与空气层交界处位置的压力比单相时增大2.4%，导叶进口处湍动能等值线有明显突变，转轮吸力面涡核高速区域范围明显增大，尾水管内低速区域范围扩大且向靠近出口部位分布，内部流态更混乱。

综上所述，当贯流式机组数值模拟考虑重力项和自由液面后，更符合机组实际运行状况，受到气相和水相的交互影响，机组性能也会产生一定的差异，在设计机组和对机组运行状态预测时需要考虑自由液面的影响。