混流式水轮机特性曲线在多重边界条件下的分区方法

马伟超，杨桀彬，赵志高，杨威嘉，杨建东

（武汉大学水资源与水电工程科学国家重点实验室，武汉 430072）

0 引 言

随着农业现代化的发展，通过修建引调水工程解决农业和生态环境用水问题是一种必然的发展趋势[1]。以重力流输水的长距离引调水工程首末段水头差往往较大，常建设消能水电站，以合理运用其剩余水头[2]，在满足引调水工程安全的前提下，充分发挥其灌溉、发电等水资源综合利用效益。水电站及水电机组的安全、稳定、高效运行也因此越来越受到水利部门的重视，对水电机组建模仿真精度提出了更高的要求[3-5]。水轮机综合特性曲线是进行水电站运行仿真与调节保证分析的基本数据来源，特性曲线的完整程度、疏密程度及准确性直接影响仿真结果的可靠性。然而，水轮机厂家所提供的混流式水轮机综合特性曲线试验结果范围相对较小，较难满足水轮机各种工况过渡过程仿真需求[6]。因此，在计算前需要对水轮机特性数据进行扩展和补充，外延内插水轮机特性曲线[7-8]，目前主要处理方法有外特性法和内特性法。

外特性法充分利用水轮机综合特性曲线和飞逸特性曲线的模型试验数据，通过纯数学方法合理预估、推测各物理量的变化规律。阮文山[9]假定力矩特性曲线将交于一点，并采用二次曲线模型扩展，此方法计算简单有效，但在低转速区域时误差较大；张培等[7,10]和王珊等[11-13]利用神经网络搭建了水轮机数学模型，有效减小了数值仿真误差，但曲线外延精度仍存在较大的提升空间，刘冬等[14]在水轮机神经网络模型中引入评价和修正函数，有效提高特性曲线的预测精度；张蓉生等[15]通过Delaunay三角剖分及分片线性插值的方式扩展了水轮机效率曲线；郑源等[16-17]将水轮机特性曲线进行简单分区处理，其区域特征描述更精确，但该方法计算量大，各分区边界点存在不连续现象。上述外特性方法由于缺乏水轮机水力特性理论的支撑，所得结果与机组运行的真实情况存在一定差异。

内特性法通过理论推导研究了水轮机结构参数之间的内在联系，充分反映了水轮机内部特征与水流流动规律，可以有效减小外特性法对个人经验的依赖。常近时[18-19]推导了水轮机广义基本方程与全特性曲线的解析表达式，方程充分反映了水轮机内部特征，但未利用外特性试验数据，计算精度较低；赵林明等[20]基于水轮机流量调节方程，提出了水轮机线性模型；杨建东等[21-22]推导了水轮机零转速和交点边界及可逆式水泵水轮机的特征交点公式；门闯社等[23-24]改进了水轮机内特性模型，并结合试验数据辨识内特性参数，但方程求解易陷入病态，参数辨识结果可能与其物理意义不符，且转速较高时，内特性方程可能无解。

针对上述局限性，本研究综合考虑水轮机内外特性，结合多重边界条件、外特性试验数据与内特性模型，提出了一种混流式水轮机特性曲线分区处理方法；对某水轮机进行实例分析，辨识其内特性参数，外延内插其特性曲线，并与典型的外特性法与内特性法对比；最后将完整的特性曲线应用于过渡过程试验反演及对比，以期获得较高的仿真精度，更好地应用于工程实际。

1 混流式水轮机特性曲线的边界条件

水轮机特性曲线通常采用单位流量、单位力矩与导叶开度、单位转速间的关系式表达。基于流量调节方程和能量平衡方程，可建立水轮机内特性模型[24]，如式（1）所示。

式中M11为单位力矩，N·m；Q11为单位流量，m3/s；n11为单位转速，r/min；a1~a7为结构参数，kg/m3，具体表达式见文献[24]，其中，a1i、a3i、a4i与开度α有关，下标i表示第i条开度线对应的参数，其他参数与开度无关。

将边界条件作为特性曲线分区的特征点，约束各分区内特性曲线的延拓范围。边界条件由特性曲线外延内插的范围确定，主要包括以下5类：开度线与n11=0轴、Q11=0轴、M11=0轴的交点、零开度线和单位力矩的交点，分别简称为零转速条件、零流量条件、飞逸条件、零开度线和交点条件。

1.1 边界条件表达式

1.1.1 零转速条件

包括零转速时的单位流量和单位力矩。单位流量可根据单位出力、效率和单位转速间的关系，通过洛必达法则计算[21]，如式（2）所示。由式（1）流量调节方程推导可得单位力矩，如式（3）所示。

式中f1(n11)与f2(n11)分别为某条开度线上单位转速关于单位出力P11(W)和水轮机效率η的函数。

1.1.2 零流量条件

零流量条件为水轮机绝对流量=0时对应的单位转速和单位力矩。由于零流量时速度三角形的特殊性，很难直接计算对应的单位转速，因此先推导转轮进口处相对流速W1在进口切线方向（图1中虚线方向）的速度分量Wb1=0时的单位转速，并对流量系数进行修正，近似得到绝对流量=0对应的单位转速。分析转轮进口处的速度，可列方程组如式（4）所示。

式中下标0表示导叶出口，下标1表示转轮进口；U为圆周速度，m/s；D为断面直径，m；n为水轮机转速，r/min；V为绝对速度，m/s；W为相对速度，m/s；μ为孔口出流的流量系数，此时情况与大孔口出流近似，流量系数μ可取0.85～0.9，考虑到叶片为流线型，可取为0.9[25]；H为工况水头，m。

将式（4）整理，零流量时的单位转速可表达为

式中μ′为相对零流量与绝对零流量转化的修正系数，可取 0～0.4，对于常规混流式水轮机，零流量时的单位转速一般略大于飞逸单位转速，以此为判断依据，通过试算确定修正系数。

在零流量处，转轮转速较大，高速的水流使得叶片进口处存在低压区，会诱导流体在低压区形成涡流产生水力损失，故很难求出零流量处的单位力矩理论解[22]。

1.1.3 飞逸条件

由于水轮机厂家只提供某一范围开度的飞逸工况散点数据，需要对飞逸特性曲线进行扩展。当机组飞逸且导叶完全关闭时，单位转速和单位流量均为0，以此作为小开度区域延拓的依据。飞逸曲线表达如式（6）所示[18]。

式中nc为飞逸状态下机组的瞬时转速，r/min；Qc为飞逸状态下机组的瞬时流量，m3/s；b0为导叶高度，m；r2为转轮出口半径，m；A2为转轮出口断面面积，m2；β2为水流出口角，(°)。

由于水轮机参数信息往往难以获得，且当水轮机发生飞逸时，机组振荡剧烈，流量与水头等参数出现大幅振荡[26-28]，此时远离最优工况，方程中内特性参数与实际参数也存在偏差，故很难直接使用该式指导飞逸曲线的扩展。

补充零开度的飞逸流量与飞逸转速后，可通过式（7）拟合飞逸条件。飞逸转速拟合系数可利用 MATLAB的nlinfit函数求解，求解时需保证至少已知5个飞逸试验点；飞逸流量拟合系数可通过最小二乘法求解。

式中α为导叶开度，(°)；C1Q,C2Q,B0n,B1n,B2n,B3n,B4n为拟合系数；n11c为飞逸状态下机组的瞬时单位转速，r/min；Q11c为飞逸状态下机组的瞬时单位流量，m3/s。

1.1.4 零开度线

零开度时，流量为0；将Q11|α=0=0代入式（1）可知，此时力矩特性曲线为关于单位转速的二次曲线，如式（8）所示。

1.1.5 交点条件

假定在中高转速区域，单位力矩也为关于单位转速的二次曲线[9]，结合零开度线的单位力矩表达式M11=a7n112，各条力矩特性曲线的交点可写作（n11p，a7n11p2），通过最小二乘法可求出各条开度线的拟合系数和交点边界，称为力矩二次曲线模型，如式（9）所示。

式中BM1i，BM2i为第i条开度线的拟合系数。

1.2 简化内特性模型与边界条件计算

上述边界条件中，零转速时的流量和飞逸条件可通过外特性的试验数据直接计算，而零转速时的单位力矩、零流量时的单位转速、零开度时的力矩特性曲线和交点边界条件表达式中包含内特性模型中的参数a1i、a7、βb1、Δα0，由于原内特性模型方程求解易陷入病态，精确解难以求取[24]，不宜直接使用，故对内特性模型进行改进，选取水轮机高效率区试验数据辨识时，叶片水流不发生脱流，可取叶片出口安放角与叶片出流角相同，即βb2≈β2，对式(1)简化，对参数a1～a7进行合理消元、替换，统一方程量纲，如式（10）所示。

式中Y=M11/n11Q11，kg/m2；X=Q11/n11，m3·min/(r·s)；b1,b2为结构参数，kg/m3，具体表达式见文献[24]；Δα0为导叶出流角α0与导叶开度角α的差角，(°)；ρ为水密度，kg/m3；g为重力加速度，m/s2；容积效率η0近似取为0.995[29]。

于是需要辨识的参数有：b1，b2，Δα0，a2，a3i，βb1，a7。式（10）可通过最小二乘法求解，应当注意的是，针对每一条开度线i，其参数a1i和a3i需视作独立的变量。

参数辨识精度直接影响边界条件计算结果的准确性。本研究给出了参数βb1的验证方法：叶片入口安放角βb1可通过速度三角形近似计算验证：考虑到水轮机在偏离最优工况较远的高转速区域的水力性能，取叶片安放角βb1比水流进口角β1略小3°～10°，而水流进口角β1表达式可根据速度三角形关系和式（4）进行推导，结果如式（11）所示，代入工况参数即可求解。

2 多重边界条件下分区拟合特性曲线

2.1 分区处理策略

根据水轮机厂家提供的模型综合特性曲线和飞逸特性曲线，对特性曲线进行外延和内插，其中外延是对已有试验数据的开度线沿单位转速轴向两端延拓；内插是根据零开度条件和外延的开度线，插值没有试验数据的中间开度线。

特性曲线在小开度区域、低转速区域和中高转速区域的性质差异显著，为了精确描述区域特征，以模型综合特性曲线外围和 5个边界条件为界，将水轮机特性曲线分作5个区域，如表1所示。

表1 水轮机特性曲线分区Table 1 Partitioning of turbine characteristic curves

由于流量调节方程不适用于高转速区域，力矩二次曲线模型在低转速区域误差明显，能量平衡方程仅在制动工况前描述较为准确，综合考虑，引入不同拟合方法拟合各区的特性曲线，区与区之间，用三次B样条曲线[30]将各区的曲线段重新拟合，平滑连接，以解决分区处理在分界点处出现的插值误差。

流量特性曲线外延内插的次序为：

1）结合零转速条件，通过流量调节方程与能量平衡方程，将曲线外延至S2区域；

2）结合飞逸条件，通过能量平衡方程与力矩二次曲线模型，将曲线外延至S3区域；

3）结合飞逸和零开度条件，采用三次多项式分片拟合的方法，对飞逸工况前已外延的开度线进行插值；

4）结合零流量条件，通过三次B样条曲线，同时实现S5区域的外延和全区域的平滑连接。

力矩特性曲线外延内插的次序为：

1）结合零转速条件，通过流量调节方程与能量平衡方程，将曲线外延至S2区域；

2）结合飞逸和交点条件，通过力矩二次曲线模型，将曲线外延至S3区域和S5区域；

3）结合零转速、飞逸和交点条件，通过三次B样条曲线，对全区域内已外延的开度线平滑连接；

结合零开度条件，采用三次多项式分片拟合的方法对已外延的开度线进行内插。

2.2 分区处理

2.2.1 S2区域外延

在S2区域，水轮机效率较低，内特性模型与水体实际流动情况存在偏差，故在S2区域的外延中只使用内特性方程的形式，忽略部分参数的物理意义，舍弃距离S2区域较远的飞逸试验数据，仅使用高效率区试验数据及零转速条件求解。经过大量试验与方法比选，最终仅保留原内特性方程组中容积效率项a5，将其余内特性参数替换为拟合参数，此时拟合精度较高。将式（1）简化为式（12）：

式中A1i~A4i为第i条开度线的拟合系数；A6与A7为拟合系数；a5=30ρgη0/π，η0≈0.995。

式（12）可通过最小二乘法求解。方程求解后均可整理为M11=f(Q11,n11)的形式，联立消去M11后，可看作是以n11为参数，Q11为变量的一元三次方程，利用盛金公式[31]进行求解。

2.2.2 基于力矩二次曲线模型的外延

基于此模型，流量特性曲线可外延至S3区域，力矩特性曲线可外延至S3与S5区域。相较于S2区域，本部分外延的不同点在于流量调节方程已不再满足，需将式（12）替换为式（9）的力矩二次曲线模型，求解方法与S2区域类似。

2.2.3 力矩特性曲线段的平滑连接

除小开度区域外，各区域内的力矩特性曲线延拓完成后，需要连接各区域的曲线段，在全区域内使用三次B样条曲线将力矩特性曲线段重新拟合连接。

2.2.4 特性曲线的内插

此时，需要对飞逸工况前的流量特性曲线和全区域内的力矩特性曲线进行内插。可使用三次多项式分片插值的方式实现：将单位转速和开度归一化处理，以单位流量和单位力矩为变量，补充零开度条件与插值开度的零流量条件，插值计算。

2.2.5 流量特性曲线在S5区域的外延与平滑连接

S5区域偏离高效率区较远，水力性能差，每条开度线仅包含飞逸条件和零流量条件2组数据，数据过少，很难结合理论公式外延，故结合飞逸工况前的计算结果，采用三次B样条曲线对全区域内的流量特性曲线拟合外延。这样将拟合范围从 S5区域扩展到全区域，既通过补充特性曲线数据，将流量特性曲线合理外延至S5区域，又实现了飞逸工况前各区流量特性曲线的平滑连接。

3 实例计算

3.1 实例模型及方法

水轮机HLD563-F13模型转轮进口直径0.365 m，模型几何比尺为19.44。其模型综合特性曲线及飞逸特性曲线由水轮机厂家开展模型试验得到，如图2所示。采用Delaunay三角剖分-分片三次多项式插值[15]加密等效率线，共提取效率线与开度线的交点 215个，飞逸特性曲线数据16个，即试验数据共231个。按本文方法、典型外特性法[9]和内特性法[24]对该特性曲线进行外延内插，并对比拟合结果以及过渡过程实测反演结果。各方法的计算流程见图3。

3.2 内特性方程参数辨识与边界条件

3.2.1 内特性方程参数辨识

内特性参数辨识结果及对比见表2和表3，其中典型内特性法未直接计算出叶片入口安放角βb1和参数a3i，可通过最小二乘法进一步计算得到。

参数b1和b2用于计算零转速条件，βb1用于计算零流量条件，a7用于计算零开度线和交点条件，其余辨识的参数只包含在简化内特性模型中，并未参与边界条件计算。由表2和表3对比参数辨识结果可知，参数b1、b2、a3i和βb1的偏差较大，其中b1结果相差112.07 kg/m3，b2结果相差82.78 kg/m3，βb1结果相差36.45°，a3i随着开度增加，结果偏差逐渐变小，最大偏差为5.39×105kg/m3。上述偏差较大的参数中，βb1在边界条件计算中最为重要，对其进行验证：以额定工况为计算工况，将相关参数代入式（11）：额定开度αr=34.74°，额定转速nr=100 r/min，额定水头Hr=80 m，原型转轮进口直径D1=7.1 m，假定水轮机导叶为标准化导叶，取D0/D1=1/1.1，可近似计算得β1=60.12°，即βb1的合理取值范围为 50.12°～57.12°。由此可知本文方法的参数辨识结果精确性更高。

表2 不同方法下部分水轮机参数辨识结果及对比Table 2 Identification results and comparison of turbine characteristic parameters by different methods

表3 不同开度及方法的参数a3i辨识结果及对比Table 3 Identification results and comparison of a3i in different GVO by different methods(×105kg·m-3)

3.2.2 边界条件

根据对应公式和参数辨识的结果，取流量系数为0.9，在满足零流量时的单位转速略大于飞逸时的单位转速前提下，近似取修正系数为0.3，各开度的零转速和零流量时单位参数计算结果如表4所示。

表4 零转速条件和零流量条件下单位参数计算结果Table 4 Results of unit parameters in zero discharge and zero speed conditions

对于飞逸条件：补充零开度的飞逸数据后，代入外特性试验结果，根据式（7）进行拟合，沿开度方向扩展的结果如图4所示，两曲线的拟合决定系数均为0.9997，拟合精度高，扩展结果可靠。

对于零开度条件：代入参数辨识的结果a7=−0.04 kg/m3，可得零开度时的流量特性曲线为Q11|α=0=0,力矩特性曲线为M11|α=0=−0.04n112；对于交点条件，代入模型综合特性曲线和飞逸特性曲线试验数据后，根据式（9），采用最小二乘法求解，可得交点条件(n11p，a7n11p2)结果为（175.31, −1215.13）。

3.3 水轮机特性曲线外延内插结果

水轮机厂家提供了从开度10°开始，每隔2°，直到最大开度40°，共16条开度线的外特性试验数据，图5a与图5b为此16条开度线外延的结果；图5c与图5d为从0°起，每隔1°，直到最大开度40°的开度线内插结果，其中开度线均从下至上依此增大。

特性曲线外延内插结果中开度线间均无交叉，同一开度线在分区边界位置未出现不连续现象；且在边界条件与试验数据的共同约束下，拟合精度较高，各区特性曲线演化规律能较好反映水轮机过流特性[32-33]：在低转速区域，单位流量较为平稳，单位转速对单位流量的影响不大，单位流量主要由导叶开度决定，进入中高转速区域后，开度线开始有下降趋势；在制动工况区域，各条开度线上的单位流量随着单位转速的增加而逐渐靠近，然后急剧降低，通过零流量边界后改变符号。

将3种方法对特性曲线的处理结果进行对比。由于经内插处理的开度线排列过于紧密，不便对比，仅节选开度为 2°、6°、10°、16°、22°、28°、34°和 40°的特性曲线对比,其中10°～40°的开度线为外延结果，有试验数据；2°和6°的开度线为内插结果，无试验数据。试验数据见图5a与5b。试验数据处的拟合决定系数为见表5。各方法在零转速和飞逸条件处计算的单位流量与单位力矩的误差对比如图6所示，其中零转速条件参考值由式（2）与（3）确定，飞逸条件参考值由式（7）确定，对比表明：

表5 单位流量与单位力矩试验数据处的决定系数Table 5 Determination coefficient of unit discharge and unit torque in experimental data

1）在试验数据处，本文方法既考虑了内特性模型，又有效提升了拟合精度：与典型内特性法相比，在试验数据处单位流量拟合决定系数提高了2.69%，单位力矩拟合决定系数提高了11.02%。

2）本文方法考虑了零转速条件，有效降低了零转速处的误差：相对于典型内特性法和典型外特性法，本文方法在零转速条件处单位流量的平均误差分别降低0.12、0.09 m3/s，在零转速条件处单位力矩的平均误差分别降低192.44、215.94 N·m；

3）通过分析制动工况区域的特点和外特性法计算方法可知，典型外特性法单位力矩的飞逸位置与交点位置单位转速明显偏高，而典型内特性法未计算到飞逸位置与交点条件；与典型外特性法相比，本文方法在大开度飞逸条件处单位力矩的平均误差降低89.68 N·m；

综上分析，本文方法综合考虑多重边界条件、外特性试验数据、简化内特性模型，采取分区处理方法实现，优势显著。

3.4 过渡过程试验反演

针对某电站实际情况，采用本文方法与典型外特性法计算的特性曲线进行过渡过程试验反演与对比，典型内特性法延拓范围不足，较难满足过渡过程仿真计算要求，故未比较。电站为“单管单机”引水布置，引水发电系统管道总长 1 087.7 m，水轮机组飞轮力矩为127 500 t·m2，工作水头为77 m，水轮机转轮进口直径为7.1 m，额定出力为367 MW，额定转速100 r/min，甩50%负荷过渡过程，比较 2种特性曲线计算过渡过程的蜗壳压力，其时域变化过程及 60～70 s内的精度分析与对比如图7所示。两种方法计算的蜗壳压力变化趋势与实测值基本一致；实测结果最大蜗壳压力为1.182 MPa，本文方法计算结果为1.202 MPa，与实测结果基本吻合；典型外特性法为 1.206 MPa；最大蜗壳压力最大误差从0.024 MPa降低至0.020 MPa，相对误差从2.03%降低至1.69%；随导叶关闭，水轮机进入小开度工况区域，61.4 s时，导叶完全关死，产生较大水击压强，蜗壳压力振荡明显。在此振荡区域内，本文方法计算结果振荡幅值更接近实测结果，反演效果更佳：其绝对平均误差为0.010 MPa，典型外特性法绝对平均误差为 0.023 MPa，降低0.013 MPa，相对误差从3.48%降低至1.47%。

综上分析，本研究提出的特性曲线处理方法在过渡过程试验反演的最大蜗壳压力及小开度区域动态过程时域响应均更加贴近实测结果。

4 结 论

本研究综合考虑多重边界条件、外特性试验数据与内特性模型，提出了一种特性曲线分区处理方法。该方法有效结合了内外特性的优点，既保证了外特性试验数据处的拟合效果，又能反映水轮机的水力特性。

特性曲线计算结果相比典型外特性法，本文结果在飞逸条件处的单位力矩平均误差降低89.68 N·m，有效提升制动工况区域拟合精度；相比典型内特性法，在外特性试验数据处的单位流量拟合决定系数提高了2.69%，单位力矩拟合决定系数提高了11.02%。，并可延拓特性曲线至制动工况区域。工程实例计算结果表明，本文方法与典型外特性法相比，蜗壳压力极值最大误差从2.03%降低至1.69%，小开度区域蜗壳压力反演相对误差从3.48%降低至1.47%，动态过程时域响应更接近实测结果。