基于改进内特性并与外特性相融合的水轮机建模方法与流程

本发明属于水轮机调节系统技术领域，具体涉及一种基于改进内特性并与外特性相融合的水轮机建模方法。

背景技术：

目前，应用于水轮机调节系统仿真的水轮机建模方法，其原理均是根据水轮机综合特性曲线中各参数的变化趋势并结合飞逸特性曲线等约束后而进行的拓展得到，没能充分考虑水轮机内在规律，模型过度依赖个人经验，具有较大的随意性。

技术实现要素：

本发明的目的是提供一种基于改进内特性并与外特性相融合的水轮机建模方法，解决了现有方法没能充分考虑水轮机内在规律，拓展结果过度依赖个人经验，具有较大的随意性的问题。

本发明所采用的技术方案是，基于改进内特性并与外特性相融合的水轮机建模方法，对水轮机各部件的能量损失进行了详细划分，建立了水轮机各部件的能量损失表达式，进而得到水轮机能量平衡方程，联立水轮机流量调节方程得到改进的水轮机内特性方程，然后将水轮机外特性作为样本数据，采用最小二乘法与遗传算法相结合的方法，辨识改进的水轮机内特性方程的相关参数，得到水轮机过渡过程仿真模型。

本发明的特点还在于：

建立了水轮机各部件的能量损失表达式，进而得到水轮机能量平衡方程具体为：

将水轮机中的能量损失分为容积损失、导叶入口撞击损失、叶片入口撞击损失、叶片出口旋转损失和圆盘摩擦损失：

容积损失功率表示为：

P_c＝λ_cρgQ_tH_t

其中λ_c为容积损失系数，表示漏水流量与总流量的比值；ρ为水密度，kg/m³；g为重力加速度，m/s²；Q_t为水轮机流量，m³/s；H_t为水轮机水头，m；

导叶入口撞击损失功率表示为：

其中，α_g为导叶入流角，rad；α_fg为固定导叶出流角，rad；D₀为导叶分布圆直径，m；b₀为导叶高度，m；ρ为水密度，kg/m³；Q_t为水轮机流量，m³/s；

叶片入口撞击损失功率表示为：

其中，α₀为导叶出流角，rad；β₁为叶片入口安放角，rad；D₁为水轮机标

称直径，m；n_t为水轮机转速，r/min；

叶片出口旋转损失功率表示为：

其中，r₂为水轮机出口半径，m；A₂为水轮机出流面积，m²；β₂为水轮机叶片出口安放角，rad；

圆盘摩擦损失功率表示为：

其中，k为圆盘损失系数；

水轮机中能量损失总和表示为：

P_L＝P_c+P_g+P₁+P₂+P_f

将水轮机各部件的能量损失代入能量守恒方程：

P_w＝P_t+P_L

其中，P_w为流过水轮机的水体失去的能量，W，P_w＝ρgQ_tH_t；P_t为水轮机获得的能量，W，P_t＝πM_tn_t/30，M_t为水轮机力矩，Nm；

利用相似原理将n_t、Q_t和M_t用单位转速n₁₁、单位流量Q₁₁和单位力矩M₁₁表示，经过整理可得基于能量平衡的表达式：

其中，a₇＝-1.2×10^-5kρ。

联立水轮机能量平衡方程和水轮机流量调节方程得到改进的水轮机内特性方程，具体为：

其中，α为导叶角度，Δα₀为导叶角度α与导叶出流角α₀之间的差角，Δα_g为导叶角度α与导叶入流角α_g之间的差角，参数a₂、a₅、a₆、a₇、b₁～b₈及Δα₀、Δα_g根据水轮机高效率区域的实测数据辨识得到。

参数a₂、a₅、a₆、a₇、b₁～b₈及Δα₀、Δα_g根据水轮机外特性数据辨识得到，水轮机特性曲线中的每一个工况点都可以用一个4维向量表示，记为(α,n₁₁,Q₁₁,M₁₁)，那么水轮机综合特性可以看作有限个工况点的集合，记为{(α_{_i},n_{11_i},Q_{11_i},M_{11_i})|i＝1,2,…,n}，其中i表示第i个工况点，n为工况点总数量，遗传算法中的种群记为{(Δα_{0_j},Δα_{g_j})|j＝1,2,…,m}，其中j表示第j个个体，m为种群中个体总数量；

具体为：

步骤a：首先产生随机的{(Δα_{0_j},Δα_{g_j})|j＝1,2,…,m}，做为遗传算法的初始种群；

步骤b：针对式采用最小二乘法得到个体j的系数b_{1_j},b_{2_j},a_{2_j}的表达式为：

其中，B₁＝[b_{1_j} b_{2_j} a_{2_j}]^T，C₁＝[c_{1_1} c_{1_2} c_{1_3}]^T，

对应的拟合误差为：

步骤c：针对式：

采用最小二乘法辨识其余参数b₃～b₈及a₅～a₇得到表达式为：

其中，B₂＝[b₃ b₄ b₅ b₆ b₇ b₈ a₅ a₆ a₇]^T，C₂＝[c_{2_1}c_{2_2} … c_{2_9}]^T，

对应的拟合误差为：

采用误差E_{1_j}与E_{2_j}的平均值作为整体误差，表示为：

步骤d：将E_max-E_{1_j}-E_{2_j}作为个体j的适应度进行筛选操作，E_max为种群中{E_{1_j}+E_{2_j}|j＝1,2……m}的最大值，对筛选后的个体进行交叉和变异来操作，产生新一代的种群，然后重新从步骤b开始执行直至所有个体足够接近为止；

此时所得个体(Δα_{0_j},Δα_{g_j})及对应的b_{1_j}～b_{8_j}，a_{2_j}，a_{5_j}，a₆__j，a_{7_j}为模型参数，从而得到水轮机过渡过程仿真模型。

本发明的有益效果是：本发明基于改进内特性并与外特性相融合的水轮机建模方法，通过建立改进水轮机内特性方程，采用水轮机外特性曲线和辨识改进内特性方程相关参数实现水轮机内特性与外特性相融合的建模方法，该方法能够提高水轮机调节系统的仿真精度。

附图说明

图1是本发明方法中导叶入口速度三角形示意图；

图2是本发明方法中叶片入口速度三角形示意图；

图3是本发明方法中水轮机出口速度三角形示意图；

图4是本发明方法中参数辨识流程图；

图5是水轮机HLN574综合特性曲线图；

图6是水轮机HLN574飞逸特性曲线图；

图7是绘制的仿真与实测的水轮机流量特性曲线对比图；

图8是绘制的仿真与实测的水轮机力矩特性曲线对比图；

图9是由水轮机HLN574绘制各导叶角度下的x与y关系图；

图10是不同模型对甩负荷过渡过程的影响仿真图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明基于改进内特性并与外特性相融合的水轮机建模方法，对水轮机各部件的能量损失进行了详细划分，建立了水轮机各部件的能量损失表达式，进而得到水轮机能量平衡方程，联立水轮机流量调节方程得到改进的水轮机内特性方程，然后将水轮机外特性作为样本数据，采用最小二乘法与遗传算法相结合的方法，辨识改进的水轮机内特性方程的相关参数，最终得到水轮机过渡过程仿真模型。

为便于分析，在不失一般性的前提下，做以下4条假设：假设1：水为理想液体，即不考虑水体的可压缩性及粘性；假设2：水轮机及流道壁面为刚性壁面，其物理形状不随受力的变化而变化；假设3：流道内水体的流动为有压流，不考虑水体的空化影响；假设4：水轮机内的流动为轴对称流动。

水轮机流量调节方程为：

其中，H_t为水轮机水头，m；Q_t为水轮机流量，m³/s；η_t为水轮机效率；n_t为水轮机转速，r/min；r₂为水轮机出口半径，m；A₂为水轮机出流面积，m²；β₂为水轮机出口安放角，rad；b₀为导叶高度，m；α₀为导叶出流角，rad。

水轮机力矩为：

其中ρ为水密度，kg/m³；g为重力加速度，m/s²。

将式(2)带入式(1)并将n_t、Q_t及M_t分别用单位转速n₁₁、单位流量Q₁₁及单位力矩M₁₁表示，经过整理可以得到如下表达式：

其中，D₁为水轮机标称直径。

水轮机的能量平衡关系可以表述为：流过水轮机的水体失去的能量等于水轮机获得的能量与损失能量之和。表示为：

P_w＝P_t+P_L (4)

其中，P_w为流过水轮机的水体失去的能量，W，P_w＝ρgQ_tH_t；P_t为水轮机获得的能量，W，P_t＝πM_tn_t/30；P_L为损失能量，W。

其中能量损失项最为复杂，随水轮机工况变化水轮机不同区域损失的能量具有较大变化，为此将水轮机中的能量损失分为容积损失、导叶入口撞击损失、叶片入口撞击损失、叶片出口旋转损失和圆盘摩擦损失，并逐一进行讨论。

容积损失

从水轮机缝隙中流失的水体所具有的能量即为容积损失，其损失功率可以表示为：

P_c＝λ_cρgQ_tH_t (5)

其中λ_c为容积损失系数，表示漏水流量与总流量的比值。

导叶入口撞击损失

水体在导叶入口处与导叶叶片发生撞击，损失部分能量。导叶入口速度三角形如图1所示，水体流入导叶时流速为v₀，流入导叶后沿导叶骨线方向流动，流速为w₀，水体流入导叶前、后流量保持不变，即v₀和w₀沿轴面分速度v_m相同。由速度三角形可知导叶入口处的损失为v_g分量，即损失水头可以表示为：

其中，v_g为导叶入口损失速度分量，m/s；α_g为当前导叶入流角，rad；α_fg为固定导叶出流角，rad；D₀为导叶分布圆直径，m。

导叶入口撞击损失功率可表示为：

叶片入口撞击损失

水体在水轮机入口处与叶片撞击损失部分能量。叶片入口速度三角形如图2所示，水体绝对速度为v₁，水轮机切向速度为u₁，水体相对流速为w₁，将导叶与叶片之间空隙的流动按照等势流考虑，那么叶片入流角为导叶出流角α₀。水体与叶片撞击前以相对速度w₁流动，与叶片撞击后相对流速方向沿叶片骨线方向流动，损失为w_c分量，叶片入口损失水头可以表示为：

由速度三角形几何关系可知：

水轮机切向流速为：

叶片入口绝对速度为：

联立(8)、(9)、(10)和(11)可得：

入口撞击损失功率为：

叶片出口旋转损失

水体在水轮机出口处因旋转而在蜗壳中损失部分能量，水轮机出口速度三角形如图3所示，水体相对流速为w₂，水轮机圆周速度为u₂，水体绝对流速为v₂，w₂和v₂在轴面的分量相同为v_m2，损失则为v_u2分量，叶片出口旋转损失水头可表示为：

圆盘摩擦损失

水轮机在水中旋转与周围的水体摩擦所损失掉的能量为圆盘摩擦损失。损失功率可根据经验公式进行计算，其表达形式为：

水轮机中能量损失总和表示为：

P_L＝P_c+P_g+P₁+P₂+P_f (17)

将水轮机各部件的能量损失代入能量守恒方程，并利用相似原理将n_t、Q_t和M_t化为单位转速n₁₁、Q₁₁和M₁₁，经过整理可得基于能量平衡的表达式：

其中，a₇＝-1.2×10^-5kρ。

改进的水轮机内特性方程建立

水轮机流量调节方程与水轮机能量平衡方程共同描述了水轮机特性，将式(3)与式(18)联立可得：

其中，a₂、a₅、a₆和a₇仅与水轮机结构及尺寸有关，a₁、a₃、a₄与水轮机结构与尺寸有关外还与水轮机当前导叶角度有关。将a₁、a₃和a₄的表达式可以写成：

其中，

其中b₁～b₈仅与水轮机结构及尺寸有关，与导叶角度无关。将水轮机导叶角度α与导叶出流角α₀之间的差角记为Δα₀，将导叶角度α与导叶入流角α_g之间的差角记为Δα_g。将式(20)带入式(19)，并考虑到导叶角度与导叶入流角及导叶出流角之间的关系即可得到：

式(21)即为改进后的水轮机内特性方程，将式中的M₁₁消掉后即可得到水轮机流量特性表达式f₁(α,n₁₁,Q₁₁)＝0，同理消掉Q₁₁即可得到水轮机力矩特性表达式f₂(α,n₁₁,M₁₁)＝0。若能根据水轮机高效率区域的实测数据辨识得到式中的a₂、a₅、a₆、a₇、Δα₀、Δα_g以及b₁～b₈，则能够得到水轮机仿真模型。

参数辨识

由于式(21)为超越方程组难以采用传统的辨识方法辨识，遗传算法虽能够对非线性方程组中的参数进行辨识，但当系统维度较高、计算规模较大时容易陷入“早熟”，辨识结果难以得到保障。表达式(21)中需要辨识的参数多达14个，采用遗传算法时种群维度高达14维，很难避免因算法“早熟”而陷入局部最优。但式(21)中仅有Δα₀和Δα_g两个参数参与了三角函数运算，其余参数均为多项式系数，因此可以采用遗传算法给出Δα₀和Δα_g的种群，然后采用最小二乘法辨识其余参数，通过寻找误差最小值获得最优Δα₀和Δα_g个体，这样就可以将遗传算法的维度降低至2维，可以有效避免算法陷入“早熟”。

水轮机特性曲线中的每一个工况点都可以用一个4维向量表示，记为(α,n₁₁,Q₁₁,M₁₁)，那么水轮机特性可以看作有限个工况点的集合，记为{(α_{_i},n_{11_i},Q_{11_i},M_{11_i})|i＝1,2,…,n}，其中i表示第i个工况点，n为工况点总数量。遗传算法中的种群可以记为{(Δα_{0_j},Δα_{g_j})|j＝1,2,…,m}，其中j表示第j个个体，m为种群中个体总数量。辨识流程如图4所示。

步骤a：首先产生随机的{(Δα_{0_j},Δα_{g_j})|j＝1,2,…,m}，做为遗传算法的初始种群；

步骤b：针对式21-1采用最小二乘法得到个体j的系数b_{1_j},b_{2_j},a_{2_j}的表达式为：

其中，B₁＝[b_{1_j} b_{2_j} a_{2_j}]^T，C₁＝[c_{1_1} c_{1_2} c_{1_3}]^T，

对应的拟合误差为：

步骤c：针对式21-2，由于式(22)已经得到b₁，采用最小二乘法辨识其余参数b₃～b₈及a₅～a₇得到表达式为：

其中，B₂＝[b₃ b₄ b₅ b₆ b₇ b₈ a₅ a₆ a₇]^T，C₂＝[c_{2_1}c_{2_2} … c_{2_9}]^T，

矩阵中系数较多，各系数的表达式不再逐一罗列。

对应的拟合误差为：

采用误差E_{1_j}与E_{2_j}的平均值作为整体误差，表示为：

步骤d：将E_max-E_{1_j}-E_{2_j}作为个体j的适应度进行筛选操作，E_max为种群中{E_{1_j}+E_{2_j}|j＝1,2……m}的最大值，对筛选后的个体进行交叉和变异来操作，产生新一代的种群，然后重新从步骤b开始执行直至所有个体足够接近为止；

此时所得个体(Δα_{0_j},Δα_{g_j})及对应的b_{1_j}～b_{8_j}，a_{2_j}，a_{5_j}～a_{7_j}为模型参数，从而得到水轮机过渡过程仿真模型。

实例验证

以水轮机HLN574为例对基于改进内特性并与外特性相融合的水轮机建模方法进行验证，水轮机HLN574综合特性曲线如图5所示，飞逸特性曲线如图6所示。

参数辨识过程

为了减少读取水轮机综合特性曲线时因插值带来的误差，读取等导叶开度线与等效率线的交叉点，综合特性曲线上共有146个可读取工况点，同时飞逸特性曲线上有9个导叶角度下的飞逸工况点，即工况点总数量n为155。

Δα₀和Δα_g的取值一般在15°内，采用浮点编码方式对0～15范围数进行编码，编码长度选11位，即编码精度为0.007，每一个个体中具有两条染色体，即分别代表Δα₀和Δα_g的取值，种群数量m选取为6，为了防止随机产生的初始种群较为接近而导致辨识过程异常退出，将以Δα₀和Δα_g所构成的二维空间均为划分为2×3个区域，每个区域内随机产生一个初始个体，组成初始种群。

采用式(26)计算出每一个个体所对应的误差E_{_j}，并求出种群中误差最大值E_max，然后采用如下适应度函数计算每一个个体适应度。

P_{_j}＝E_max-E_{_j} (27)

采用轮盘赌法复制出6个个体，随机进行两两匹配，并随机选取1～11的整数值作为交叉点进行交叉操作，变异概率选取为0.0001，经过变异操作后得到新的6个个体。

辨识结果及误差分析

通过辨识得到水轮机特征参数Δα₀,Δα_g，a₂,a₅,a₆,a₇,b₁～b₈如表1所示，采用该组参数绘制的水轮机流量特性曲线如图7中的虚线所示，实测水轮机特性曲线如图7中的实线所示，采用该组参数绘制的水轮机力矩特性曲线如图8中的虚线所示，实测水轮机特性曲线如图8中的实线所示。

表1水轮机HLN574特征参数辨识结果

由图可以看出拓展得到的小开度区域、低单位转速区域的特性曲线均能够与实测曲线较好吻合，但当单位转速较大的区域误差较大。为便于分析引起该区域误差的原因将式(3)整理为如下形式：

y＝a₁x+a₂ (28)

其中

可见当导叶角度给定时x与y满足线性关系，由水轮机HLN574绘制各导叶角度下的x与y关系如图9所示，在图9中x与y在较大范围内能够满足线性关系，但当x取较小值时x与y严重偏离了线性关系，说明n₁₁取较大值Q₁₁取较小值时水轮机流量调节方程已难以得到满足，这是由于水轮机工作在该区域时流体流态较为复杂，常伴随有局部旋涡、空化甚至旋转失速等现象，此时模型中的假设条件难以得到满足，进而引起误差增大。

过渡过程仿真分析

为了分析模型误差对过渡过程的影响大小，分别采用内特性模型和实验测量水轮机特性曲线进行过渡过程仿真，仿真电站引水管道长为1200m，水轮机组飞轮力矩为4500t·m²，水头为430m，水轮机直径为2.238m，额定出力为300MW，额定转速428.6r/min，导叶采用17s直线关闭规律，甩100％负荷过渡过程，为便于表示采用实测水轮机特性曲线的仿真称为案例1，采用内特性模型的仿真称为案例2，仿真结果如图10所示。采用本发明提出的内特性仿真模型的最高转速上升率为48.67％，最高水压上升率为12.19％；采用试验所得水轮机特性曲线仿真的最高转速上升率为49.45％，最高水压上升率为12.35％。可见本发明提出的水轮机内特性模型与实验曲线进行的仿真结果较为接近，能够达到工程应用的要求。