CN108763785A - 一种水轮机调节系统实时仿真方法、装置和电子设备 - Google Patents

Abstract

本发明公开了一种水轮机调节系统实时仿真方法、装置、电子设备。所述方法包括：获取水电站基本参数、基本常数，进行变量定义，通过所述基本参数、基本常数、所述变量设定基本关系式；建立基于特征线法的所述水电站的水轮机引水系统的水击特征线数学模型；建立基本特征线法系统计算原理框图；建立去掉水头损失项、引入水头到流量的反馈的基于基本特征线法系统计算原理框图的改进的线型模型计算原理框图，依据所述改进的线型模型计算原理框图写出所述水轮机引水系统的z传递函数，获得所述改进的线型模型；建立所述改进的线型模型进行仿真。

Description

一种水轮机调节系统实时仿真方法、装置和电子设备

技术领域

本发明涉及水轮机仿真领域，特别是指一种基于改进后特征线法的水轮机调节系统实时仿真方法、装置和电子设备。

背景技术

特征线法是目前最常用的求解水轮机调节系统的管道系统水力瞬变过程的方法，其优点主要有：稳定性准则可以建立；边界条件容易编写程序；可适用于各种管道的水力瞬变分析。因此是所有离散算法中具有较好精度的一种。基本特征线法可以表现出弹性水击的性质，但弹性水击波衰减较慢，阻尼较小，收敛性差，尤其是在水流切断后更是如此，而实时仿真算法通常需要连续运行，因此对仿真迭代在各种工况下的收敛性有更严格的要求，否则可能导致仿真结果发散，严重失真。

发明内容

有鉴于此，为克服弹性水击波衰减较慢，阻尼较小的缺陷，在水轮机调速系统试验装置的非线性实时仿真计算中，本发明提出一种水轮机调节系统实时仿真方法、装置和电子设备，实现了对水击波的阻尼控制，取得了良好的效果。

第一方面，本发明提供了一种水轮机调节系统实时仿真方法，包括：

获取水电站基本技术数据，利用所述基本技术数据进行基本参数和变量定义、基本常数定义、计算用变量和常数定义、基本关系式定义；

利用所述基本参数和变量、所述基本常数、所述计算用变量和常数、所述基本关系式建立基于特征线法的所述水电站的水轮机引水系统的水击特征线数学模型；

利用所述水击特征线数学模型建立基本特征线法系统计算原理框图；

在所述基本特征线法系统计算原理框图基础上，建立去掉水头损失项、引入水头到流量的反馈的线型模型计算原理框图，依据所述线型模型计算原理框图得到所述水轮机引水系统的z传递函数，获得改进的线型模型；

利用所述改进的线型模型进行建模仿真。

在一些实施方式中，所述获取水电站基本技术数据的步骤中，通过实地测量、监控数据采集的方法获取所述水电站基本技术数据。

所述基本参数和变量为：

v—水击波速；L—等效管道长；D—等效管道直径；H_r—计算额定净水头；Q_r—计算额定流量；Q_A—水轮机入口流量；Q_B—水库入口流量；F—等效管道截面积；g—重力加速度；T_s—采样周期；y—接力器位置相对值；a—导叶开度相对值；j—轮叶开度相对值；H—水轮机实际作用水头；N₁₁—水轮机单位转速；Q₁₁—水轮机单位流量；M₁₁—水轮机单位力矩；—水轮机力矩相对值；D₁—水轮机转轮直径；

所述基本常数为：

T_r＝2L/v—水击相长；h_w＝vQ_r/(2gFH_r)—管道常数；T_w＝h_wT_r＝LQ_r/(gFH_r)—水流加速时间常数；

所述计算用变量和常数为：

—水轮机入口相对流量；—水库入口相对流量；—水轮机入口相对压力上升值；—水头损失系数；C_a＝1/(2h_w)—计算常数；

所述基本关系式为：

a＝f_a(y)—导叶开度相对值与接力器位置相对值间的非线性关系；

—水轮机单位转速计算式；

Q₁₁＝f_Q(N₁₁,a)—水轮机单位流量非线性特性；

—水轮机流量相对值；

M₁₁＝f_M(N₁₁,a)—水轮机单位力矩非线性特性；

—水轮机输出力矩相对值；

所述水电站的水轮机引水系统的水击特征线数学模型是单管单机简单情形下的；

所述水击特征线数学模型为：

其中，引水管道末端假定为A点，引水管道首端假定为B点，作用于所述水轮机的静水头选为基准水头H_r，忽略管道损失因素，h_B,k＝1，T_s＝L/v，C_a＝T_s/T_w，

在一些实施方式中，所述建立去掉水头损失项、引入水头到流量的反馈的线型模型计算原理框图的步骤中，为所述水轮机相对流量；所述水轮机引水系统的z传递函数为：

其中，F为水头到流量的反馈系数，K_l＝0。

在一些实施方式中，所述依据所述线型模型计算原理框图写出所述水轮机引水系统的z传递函数之后还能够包括，通过重采样技术，基于所述水轮机引水系统的z传递函数获得具有可调整计算步长的水轮机引水系统离散化数学模型：

其中，

a＝T_sn+T_so

b＝T_sn-T_so

c＝F+C_a

d＝F-C_a

T_so＝L/v

C_a＝T_so/T_w＝L/(vT_w)

A＝b²-a²＝(T_sn-T_so)²-(T_sn+T_so)²＝-4T_snT_so

B＝(a²-b2)＝-(b²-a²)＝4T_snT_so

T_sn为所述新的采样周期，T_so为所述z传递函数对应的采样周期。

在一些实施方式中，所述获得具有可调整计算步长的水轮机引水系统离散化数学模型之后还能够包括，对所述水轮机引水系统离散化数学模型进行修正，获得修正后引水系统数学模型：在所述具有可调整计算步长的水轮机引水系统离散化数学模型的分子常数项上补入一个水头稳态误差系数K_e；

修正后引水系统数学模型为：

所述水头稳态误差系数

其中，h_A(∞)是人为设定的值。

在一些实施方式中，所述获得具有可调整计算步长的水轮机引水系统离散化数学模型之后，能够建立所述具有可调整计算步长的水轮机引水系统数学模型进行仿真。

在一些实施方式中，所述对所述水轮机引水系统离散化数学模型进行修正之后，能够建立所述修正后引水系统数学模型进行仿真。

第二方面，本发明提供了一种水轮机调节系统实时仿真装置，包括：

定义模块，用于获取水电站基本技术数据，利用所述基本技术数据进行基本参数和变量定义、基本常数定义、计算用变量和常数定义、基本关系式定义；

特征线建模模块，用于利用所述基本参数和变量、所述基本常数、所述计算用变量和常数、所述基本关系式建立基于特征线法的所述水电站的水轮机引水系统的水击特征线数学模型；

第一绘图模块，用于利用所述水击特征线数学模型建立基本特征线法系统计算原理框图；

第二绘图模块，用于在所述基本特征线法系统计算原理框图基础上，建立去掉水头损失项、引入水头到流量的反馈的线型模型计算原理框图，依据所述线型模型计算原理框图写出所述水轮机引水系统的z传递函数，获得改进的线型模型；

仿真模块，用于利用所述改进的线型模型进行建模仿真。

在一些实施方式中，所述定义模块还用于，

通过实地测量、监控数据采集的方法获取所述水电站基本技术数据。

所述基本参数和变量为：

v—水击波速；L—等效管道长；D—等效管道直径；H_r—计算额定净水头；Q_r—计算额定流量；Q_A—水轮机入口流量；Q_B—水库入口流量；F—等效管道截面积；g—重力加速度；T_s—采样周期；y—接力器位置相对值；a—导叶开度相对值；j—轮叶开度相对值；H—水轮机实际作用水头；N₁₁—水轮机单位转速；Q₁₁—水轮机单位流量；M₁₁—水轮机单位力矩；—水轮机力矩相对值；D₁—水轮机转轮直径；

所述基本常数为：

T_r＝2L/v—水击相长；h_w＝vQ_r/(2gFH_r)—管道常数；T_w＝h_wT_r＝LQ_r/(gFH_r)—水流加速时间常数；

所述计算用变量和常数为：

—水轮机入口相对流量；—水库入口相对流量；—水轮机入口相对压力上升值；—水头损失系数；C_a＝1/(2h_w)—计算常数；

所述基本关系式为：

a＝f_a(y)—导叶开度相对值与接力器位置相对值间的非线性关系；

—水轮机单位转速计算式；

Q₁₁＝f_Q(N₁₁,a)—水轮机单位流量非线性特性；

—水轮机流量相对值；

M₁₁＝f_M(N₁₁,a)—水轮机单位力矩非线性特性；

—水轮机输出力矩相对值；

在一些实施方式中，所述特征线建模模块还用于，所述水电站的水轮机引水系统的水击特征线数学模型是单管单机简单情形下的；

所述水击特征线数学模型为：

其中，引水管道末端假定为A点，引水管道首端假定为B点，作用于所述水轮机的静水头选为基准水头H_r，忽略管道损失因素，h_B,k＝1，T_s＝L/v，C_a＝T_s/T_w，

在一些实施方式中，所述第二绘图模块还用于，为所述水轮机相对流量；所述水轮机引水系统的z传递函数为：

其中，F—水头到流量的反馈系数，K_l＝0。

第三方面，本发明提供了一种水轮机调节系统实时仿真电子设备，包括：存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现如上述任意一项所述的方法。

从上面所述可以看出，本发明提供的基于改进后特征线法的水轮机调节系统实时仿真方法、装置和电子设备，在水电机组调节系统仿真算法中针对流量与水头的关系，引入水头到流量的反馈系数F，通过调整F可以改变系统的阻尼特性，实现了对水击波的阻尼控制，可以模拟真实引水管道中水击波的衰减过程。

进一步的，在具有水头到流量反馈的模型的基础上，提出重采样技术，建立引水系统计算采样周期可以任意设定的模型，实现了引水系统和发电机-负载模型仿真时采样周期的协调统一。

更进一步的，在计算采样周期可以任意设定的模型的基础上提出了引入水头稳态误差系数K_e，造成一个与实际情况更为相符的稳态水头损失，实现可以人为地调整其在稳态时的水头损失值，使仿真模型更接近实际情况。

以上所述采用基于改进后特征线法在水轮机调节系统实时仿真方法，使得仿真过程比较贴近于机组实际过渡过程，具有良好的仿真效果。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例的水轮机调节系统实时仿真方法流程图；

图2为本发明实施例的水轮机调节系统实时仿真装置结构示意图；

图3为本发明实施例的典型水电站引水系统示意图；

图4为本发明实施例的典型水电站引水系统的特征线网格图；

图5为本发明实施例的单管单机引水系统的特征线计算示意图；

图6为本发明实施例的单管单机引水管道单段计算的基本特征线法计算原理框图；

图7为本发明实施例的单管单机引水系统的基于特征线法的改进线性模型计算原理框图；

图8为本发明实施例的F＝0反馈作用下的系统响应示意图；

图9为本发明实施例的F＝0.001反馈作用下的系统响应示意图；

图10为本发明实施例的F＝0.01反馈作用下的系统响应示意图；

图11为本发明实施例的F＝0.1反馈作用下的系统响应示意图；

图12为本发明实施例的重采样后引水系统响应示意图；

图13为本发明实施例的引入水头损失后系统响应示意图；

图14为本发明实施例的甩25％负荷过程真机与模型比较示意图；

图15为本发明实施例的甩50％负荷过程真机与模型比较示意图；

图16为本发明实施例的甩75％负荷过程真机与模型比较示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

需要说明的是，本发明实施例中所有使用“第一”和“第二”的表述均是为了区分两个相同名称非相同的实体或者非相同的参量，可见“第一”“第二”仅为了表述的方便，不应理解为对本发明实施例的限定，后续实施例对此不再一一说明。

本发明实施例提供了一种水轮机调节系统实时仿真方法、装置、电子设备，下面结合附图对本发明所提供的方法的技术方案做进一步的详尽的阐释。

作为本发明的一个实施例，参考图1，为本发明实施例的水轮机调节系统实时仿真方法流程图。

所述基于改进后特征线法的水轮机调节系统实时仿真方法，包括以下步骤：

步骤101、获取水电站基本技术数据，利用所述基本技术数据进行基本参数和变量定义、基本常数定义、计算用变量和常数定义、基本关系式定义。

在本步骤中，参考图1，所述获取水电站基本技术数据的步骤中，通过实地测量、监控数据采集的方法获取所述水电站基本技术数据，进行基本参数和变量定义：

v—水击波速；L—等效管道长；D—等效管道直径；H_r—计算额定净水头；Q_r—计算额定流量；Q_A—水轮机入口流量；Q_B—水库入口流量；F—等效管道截面积；g—重力加速度；T_s—采样周期；y—接力器位置相对值；a—导叶开度相对值；j—轮叶开度相对值；H—水轮机实际作用水头；N₁₁—水轮机单位转速；Q₁₁—水轮机单位流量；M₁₁—水轮机单位力矩；—水轮机力矩相对值；D₁—水轮机转轮直径；

定义基本常数：

T_r＝2L/v—水击相长；h_w＝vQ_r/(2gFH_r)—管道常数；T_w＝h_wT_r＝LQ_r/(gFH_r)—水流加速时间常数。

定义计算用变量和常数：

—水轮机入口相对流量；—水库入口相对流量；—水轮机入口相对压力上升值；—水头损失系数；C_a＝1/(2h_w)—计算常数。

定义基本关系式：

a＝f_a(y)—导叶开度相对值与接力器位置相对值间的非线性关系；

—水轮机单位转速计算式；

Q₁₁＝f_Q(N₁₁,a)—水轮机单位流量非线性特性；

—水轮机流量相对值；

M₁₁＝f_M(N₁₁,a)—水轮机单位力矩非线性特性；

—水轮机输出力矩相对值。

步骤102、利用所述基本参数和变量、所述基本常数、所述计算用变量和常数、所述基本关系式建立基于特征线法的所述水电站的水轮机引水系统的水击特征线数学模型。

所述水电站的水轮机引水系统的水击特征线数学模型是单管单机简单情形下的；

所述水击特征线数学模型为：

其中，引水管道末端假定为A点，引水管道首端假定为B点，作用于所述水轮机的静水头选为基准水头H_r，忽略管道损失因素，h_B,k＝1，T_s＝L/v，C_a＝T_s/T_w，

步骤103、利用所述水击特征线数学模型建立基本特征线法系统计算原理框图。

本步骤中，所述基本特征线法系统计算原理框图中的符号“|·|”表示绝对值运算。

步骤104、在所述基本特征线法系统计算原理框图基础上，建立去掉水头损失项、引入水头到流量的反馈的线型模型计算原理框图，依据所述线型模型计算原理框图得到所述水轮机引水系统的z传递函数，获得改进的线型模型。

在本步骤中，所述线型模型计算原理框图中，为所述水轮机相对流量；所述水轮机引水系统的z传递函数为：

其中，F为水头到流量的反馈系数，K_l＝0。

步骤105、利用所述改进的线型模型进行建模仿真。

本步骤中，进一步的，所述仿真软件可以是MATLAB的simulink。

在一些实施例中，进一步的，在步骤104之后，通过重采样技术，基于所述水轮机引水系统的z传递函数还获得具有可调整计算步长的水轮机引水系统离散化数学模型：

其中，

a＝T_sn+T_so

b＝T_sn-T_so

c＝F+C_a

d＝F-C_a

T_so＝L/v

C_a＝T_so/T_w＝L/(vT_w)

A＝b²-a²＝(T_sn-T_so)²-(T_sn+T_so)²＝-4T_snT_so

B＝(a²-b²)＝-(b²-a²)＝4T_snT_so

T_sn为所述新的采样周期，T_so为所述z传递函数对应的采样周期。

在一些实施例中，更进一步的，所述获得具有可调整计算步长的水轮机引水系统离散化数学模型之后还可以包括，对所述水轮机引水系统离散化数学模型进行修正：在所述具有可调整计算步长的水轮机引水系统离散化数学模型的分子常数项上补入一个水头稳态误差系数K_e；修正后引水系统数学模型为：

所述水头稳态误差系数

其中，h_A(∞)是人为设定的值。

作为另一个实施例，参考图1，为本发明实施例的水轮机调节系统实时仿真方法流程图。

步骤101、获取水电站基本技术数据，利用所述基本技术数据进行基本参数和变量定义、基本常数定义、计算用变量和常数定义、基本关系式定义。

在本步骤中，所述获取水电站基本技术数据的步骤中，通过实地测量、监控数据采集的方法获取所述水电站基本技术数据，进行基本参数和变量定义：

v—水击波速；L—等效管道长；D—等效管道直径；H_r—计算额定净水头；Q_r—计算额定流量；Q_A—水轮机入口流量；Q_B—水库入口流量；F—等效管道截面积；g—重力加速度；T_s—采样周期；y—接力器位置相对值；a—导叶开度相对值；j—轮叶开度相对值；H—水轮机实际作用水头；N₁₁—水轮机单位转速；Q₁₁—水轮机单位流量；M₁₁—水轮机单位力矩；—水轮机力矩相对值；D₁—水轮机转轮直径；

定义基本常数：

T_r＝2L/v—水击相长；h_w＝vQ_r/(2gFH_r)—管道常数；T_w＝h_wT_r＝LQ_r/(gFH_r)—水流加速时间常数。

定义计算用变量和常数：

—水轮机入口相对流量；—水库入口相对流量；—水轮机入口相对压力上升值；—水头损失系数；C_a＝1/(2h_w)—计算常数。

定义基本关系式：

a＝f_a(y)—导叶开度相对值与接力器位置相对值间的非线性关系；

—水轮机单位转速计算式；

Q₁₁＝f_Q(N₁₁,a)—水轮机单位流量非线性特性；

—水轮机流量相对值；

M₁₁＝f_M(N₁₁,a)—水轮机单位力矩非线性特性；

—水轮机输出力矩相对值。

所述获取水电站基本参数和基本常数的步骤中，通过实地测量、参数收集、参数设定的方法获取所述水电站基本参数、基本常数，进行变量定义通过如下方式得到：

参考图3，对本发明所用到的数学模型和计算方法中使用的基本参数、基本常数、变量和基本计算公式进行定义：

第一步，基本参数和变量定义

v—水击波速，米/秒；

L—等效管道长，米；

D—等效管道直径，米；

H_r—计算额定净水头，米；

Q_r—计算额定流量，米³/秒；

Q_A—水轮机入口流量，米³/秒；

Q_B—水库入口流量，米³/秒；

F—等效管道截面积，米²；

g—重力加速度，米/秒²；

T_s—采样周期，秒；

y—接力器位置相对值，无量纲；

a—导叶开度相对值，无量纲；

j—轮叶开度相对值，无量纲；

H—水轮机实际作用水头；

N₁₁—水轮机单位转速；

Q₁₁—水轮机单位流量；

M₁₁—水轮机单位力矩；

—水轮机力矩相对值；

D₁—水轮机转轮直径，米。

第二步，定义基本常数：

T_r＝2L/v—水击相长；

h_w＝vQ_r/(2gFH_r)—管道常数；

T_w＝h_wT_r＝LQ_r/(gFH_r)—水流加速时间常数，秒；

第三步，定义计算用变量或常数：

—水轮机入口相对流量；

—水库入口相对流量；

—水轮机入口相对压力上升值；

—水头损失系数；

C_a＝1/(2h_w)—计算常数。

第四步，定义基本关系式：

a)导叶开度相对值与接力器位置相对值间的非线性关系：

a＝f_a(y) (1)

b)水轮机单位转速计算式：

c)水轮机单位流量非线性特性：

Q₁₁＝f_Q(N₁₁,a) (3)

d)水轮机流量相对值：

e)水轮机单位力矩非线性特性：

M₁₁＝f_M(N₁₁,a) (5)

f)水轮机输出力矩相对值：

如果为转浆式或贯流式机组，水轮机单位流量非线性特性和水轮机单位力矩非线性特性又可表示为：

Q₁₁＝f_Q(N₁₁,a,j) (7)

M₁₁＝f_M(N₁₁,a,j) (8)

步骤102，利用所述基本参数和变量、所述基本常数、所述计算用变量和常数、所述基本关系式建立基于特征线法的所述水电站的水轮机引水系统的水击特征线数学模型。

所述水击特征线数学模型为：

所述水击特征线数学模型通过如下方式得到：

为合理利用水能并兼顾水电站开发的可行性和经济性，水电站的引水系统的布置方式和结构通常是多种多样的。一般而言，主要由上下游水库、压力引水管道、水轮机流道、尾水管等共同组成。当压力引水管很长时，为了减少水力动态过程中的水击压力和改善水轮机调节系统的动态品质，有时还需设置上游调压井或尾水调压井。引水管道的布置通常可分为单管单机和单管多机的方式。对于较为复杂的引水管道还可能包括各种管道阀门和岔管等部件。

本步骤中，采用单管单机简单情形下的引水系统，示意图如图3所示。

压力引水系统最显著的特点是其固有的水击现象。水击现象是指当压力水管末端的流量发生变化时，水管内将出现非恒定流现象，其特点是随着流速的改变，压强有显著的变化。通俗地讲，当水轮机的导叶开度发生变化时，水轮机的流量和力矩都跟着发生变化，而流量的变化会在引水系统中引起水击，反过来又会引起水轮机力矩的变化。因此，在系统的动态过程中，引水系统的各个元件彼此相互依存、相互影响，构成非常复杂的有着内部联系的统一体。显然，要准确描述其完整的特性，需建立十分详细的数学模型，这将是一个非常困难的工作，考虑到计算的目的和应用范围，在所述引水系统的建模中，假定：

1)引水系统仅按简单的单管对单机情形考虑，即忽略各引水管道的相互作用；

2)流速在整个管道截面上均匀分布；

3)不考虑水中含气量对水击波速的影响，不考虑水的汽化状态；

4)近似采用稳态阻力损失公式。

考虑弹性水击的影响，使用特征线法进行引水系统的建模。列出描述管道内某流体隔离体的连续性方程和运动方程：

连续方程：

运动方程：

其中，Q为管道中某一过流断面在时间t的流量(m³/s)；H为管道中某一过流断面在时间t的水头(m)；x为相应过流断面到规定点的距离(m)；D为管道直径(m)；F为管道截面积(m²)；v为水击波速(m/S)；f为管道水力摩擦阻力系数。在摩擦损失项内引入绝对值符号，以便使水体摩擦力与流速方向相反。

利用特征线法求解式(9)和式(10)描述的一阶双曲型偏微分方程组。

特征线法是目前最常用的求解管道系统水力瞬变过程的方法，其优点主要有：稳定性准则可以建立；边界条件很容易编写程序；可适用于各种管道的水力瞬变分析。因此是所有离散算法中具有较好精度的一种。然而，将其直接用于仿真装置中水轮机引水系统的水击计算，存在一个问题：特征线法的基本要求是，计算步长要严格等于计算管段长/水击波速。当管道较长时，为避免计算本身可能引起的计算不稳定性，一般需要分段计算，这又不可避免地增加了计算量，在计算资源有限的实时仿真装置中，一般是不希望的。此外，也不利于参与仿真的各环节在计算步长上的协调统一。

为求解式(9)和式(10)描述的偏微分方程组，令

并用因子λ进行如下线性组合

L＝L₁+λL₂＝0

得

整理可得

如果H＝H(x,t)和Q＝Q(x,t)，则他们的全微分可写为

定义

即

将式(12)、(13)、(14)、(15)代入到式(11)中得下述两组方程：

把偏微分方程式(9)和式(10)转化为独立变量t的常微分方程式。式(16)、式(17)分别为所述单管单机引水系统的正、负特征线方程，其中dx/dt＝±v表示的直线称为特征线(方程的约束)的路线，如图4所示。

在(16)式中，用dt乘以第二个方程，然后沿C⁺积分得：

对运用Streeter(斯特瑞特)近似积分公式，得

因而有

同理可得

整理(18)和(19)式可得另一种形式

设在R和S点对应的时刻为t＝(k-1)T_s，T_s为时间计算步长(或称采样周期)，为简单起见，后续过程中表示为(k-1)。P点对应的时刻为k，则上两式写成其对应的递推形式，即

所述引水系统的管道长L，分成N段，每一段长度为取时间步长则(22)和(23)变为：

令水轮机静水头H_r为水头基准值，Q_r为流量基准值，将式(24)和式(25)方程两端同时除以并以表示流量的相对值，以表示水头的相对值，同时再进一步令

管道常数─h_w＝vQ_r/(2gFH_r)

计算常数─C_a＝1/(2h_w)

基值流量水头损失相对值─

得工程上常用的水击特征线方程：

其中，C'₊—正特征线，C'_-—负特征线，对于每一时段，C_a＝1/(2h_w)是常量，C'₊、C'_-取决于管道特性和前一时刻k-1的状态，是己知的。式(28)、(29)中的两个未知量和可通过联立求解以下两个方程来确定，即

利用上式仿真水力过渡过程时，从k＝0时的定常流状态开始，计算出所述管道每一个计算截面上的和的起始值利用上两式，计算出求取k＝1时和的值，求取k＝2时和的值，以此类推，逐时步进行计算，求出管道上所有网格节点处的在所需仿真时间内的全部过渡过程参数。由于在边界点缺少正特征线或负特征线，需要补充边界条件方程，才能求解边界节点的动态参数。

对于本发明所述单管单机的情形，可把有压管道简化为一根等效管，不必求中间点的压力和流量，如图5所示单管单机引水系统的特征线计算示意图，对引水管道末端(假定为A点)，仅需应用方程(28)，即

而对引水管道首端(假定为B点)，仅需应用方程(29)，即

因此，为了求出B与A点的压力、流量还需寻求边界条件。A点为水轮机的入水口，其边界条件包含在水轮机特性中，需与水轮机特性联立；B点为水库，并认为水库水位在过渡过程中不变。如果假定作用于该水轮机的静水头也被选为基准水头H_r，并忽略管道损失等因素，易知h_B,k＝1，相应地，式(26)～(29)就变为

其中，h_w＝vQ_r/(2gFH_r)，C_a＝1/(2h_w)，Δx＝L；

如果以增量相对值表示水头(小写字母h，下同)，即

式(30)～(33)又可以进一步写成更加简洁的形式

其中，C”₊是C'₊的简洁表示形式的字母代号。由特征线法的基本前提条件，计算步长还应满足条件T_s＝L/v。式(35)～(37)是在所述单管单机情形下，所述水轮机引水系统的水击特征线数学模型。

考虑到工程人员的实际情况和计算上的简洁，计算前，如果给定等效管道长度L、水击波速v及水流加速时间常数T_w，则：

式(35)～(37)中使用的参数汇总为：

T_s＝L/v (38)

C_a＝T_s/T_w (39)

计算的初始条件为：

a₀＝f_a(y₀)

Q₁₁₀＝f_Q(N₁₁₀,a₀)

步骤103，利用所述水击特征线数学模型建立基本特征线法系统计算原理框图。

本步骤中，基于式(35)～(37)描述的所述水轮机引水系统的水击特征线数学模型，用如图6的系统计算原理框图描述，图6中符号“|·|”表示绝对值运算。

所述图6的系统计算原理框图通过如下方法得到：

基本特征线法可以表现出弹性水击的性质，但弹性水击波衰减较慢，阻尼较小，尤其是在水流切断后更是如此。这种情况通常与实际不符。人们一般认为弹性水击波的衰减与水头损失系数K_l关系，水头损失系数越大，衰减越快。但比较两图可知，这种欠阻尼情况与水头损失系数有一定的关系，但相关性并不大，即使K_l由0.02增加到0.1，当水流切断后，阻尼情况也没有改善。这也进一步说明了基本特征线法存在的两个问题：

1)水头波的阻尼问题，尤其是在水流突然切断后。

2)计算采样周期问题。当管道较长，水击波速度一定，则引水系统计算采样周期通常较长T_s＝L/v。如管道长1000米，水击波速1000米/秒，则为满足特征线法要求，计算采样周期应为1秒，这对于系统其它部分，如机发电机-负载模型的仿真来说，显然是太长了。当然可以考虑采用对不同的部分，采用不同的采样周期，但这样一来，对于实时仿真的实现来说，造成了不必要的麻烦。理想的情况是：将基本特征线法通过一定的改进，使得采样周期可以任意设定，这样就可实现引水系统和发电机-负载模型仿真时采样周期的统一。

步骤104，在所述基本特征线法系统计算原理框图基础上，建立去掉水头损失项、引入水头到流量的反馈的线型模型计算原理框图，依据所述线型模型计算原理框图得到所述水轮机引水系统的z传递函数，获得改进的线型模型。

所述线型模型计算原理框图如图7所示；其中，为所述水轮机相对流量；

所述水轮机引水系统的z传递函数为：

其中，F—水头到流量的反馈系数，K_l＝0。

所述线型模型计算原理框图、所述水轮机引水系统的z传递函数通过如下方式得到：

第一步，建立去掉水头损失项、引入水头到流量的反馈的基于基本特征线法系统方框图的线型模型计算原理框图。

大量的仿真研究表明，要实现对水击波的阻尼控制，一个很好的方法是引入水头到流量的反馈。事实上，当对正常运行的水轮机进行仿真时，通常水头对流量有较大的影响。在非线性仿真时，这种影响体现在：

Q₁₁＝f_Q(N₁₁,a)

影响的作用过程为水轮机工作在大流量时，这种反馈作用就强，水轮机工作在0开度时，Q₁₁为零，这种反馈作用消失。这是水轮机在甩负荷时，当导叶开度关到0开度时，水击波阻尼最小的原因。

本步骤中，为避免上述情况发生，引入一个固定的水头到流量的反馈，使得在零流量时，这个固定反馈仍能起到一定的作用。

在图6所示单管单机引水管道单段计算的基本特征线法系统计算原理框图上，去掉水头损失项、引入水头到流量的反馈后，计算原理框图变为图7所示情形，F为水头到流量的反馈系数，为水轮机相对流量。仅讨论由到h_A,k间的线性模型。

第二步，列写并整理基于所述线型模型计算原理框图的所述水轮机引水系统的z传递函数。

首先，依据图7列写所述水轮机引水系统的离散形式(z传递函数)方程

Q_B(z)＝Q_A(z)z^-1-C_aH_A(z)z^-1 (41)

Q_A(z)＝FH_A(z)+Q(z) (43)

其次，所有变量的z变换，均用变量的大写字母表示，则Q_A(z)变为Q_A。将式(43)代入式(42)和式(41)得

Q_B＝(FH_A+Q)z^-1-C_aH_Az^-1 (44)

将(45)带入(44)得

整理后得(F+C_a-(F-C_a)Z^-2)H_A＝Q(-1+Z^-2)

最后，对应的所述水轮机引水系统的z传递函数为

步骤105，确定反馈系数F，利用所述改进的线型模型进行建模仿真。

在一些实施例中，进一步的，可以利用MATLAB的Simulink进行建模仿真。

结合图8-11对本实施例作进一步的描述。

首先、定义参数：

等值引水管道长：L＝200米/秒

水击波速：v＝1000米/秒

计算步长：T_s＝L/v＝200/1000＝0.2

水流加速时间常数：T_w＝2s

水头到流量的反馈系数分别取为F＝0，F＝0.001，F＝0.01和F＝0.1，并假定，水轮机流量与接力器开度成正比，并在5秒内线性地由全关到全开或由全开到全关。所有变量初始条件均为0。

其次，在MATLAB的Simulink中，建立图7所示模型，代入定义参数，进行仿真计算，结果如图8、图9、图10和图11所示，其中图8(a)、图9(a)、图10(a)、图11(a)为反馈作用下的所述引水系统的时域响应过程，图8(b)、图9(b)、图10(b)、图11(b)为所述反馈作用下的所述引水系统的零极点分布图。

由图8、图9、图10和图11可以看出：

当水头到流量的反馈系数F＝0时，如图8所示,801为导叶开度，802为水压相对偏差值。系统有两个极点落在单位圆上。这时，系统是不稳定的，其时域响应表现为明显的等幅振荡；

当水头到流量的反馈系数F＝0.001时，如图9所示,901为导叶开度，902为水压相对偏差值；

当水头到流量的反馈系数F＝0.01时，如图10所示,1001为导叶开度，1002为水压相对偏差值；

随着人为引进的反馈系数的增加，系统原来在单位圆上的两个极点(箭头所指处)向圆心移动，从而使得系统响应输出的阻尼加大，衰减加快；

当水头到流量的反馈系数F＝0.1时，如图11所示,1101为导叶开度，1102为水压相对偏差值；系统的这两个极点移动到圆心。这时，系统呈现过阻尼特性，已经不具有振荡特性。

调整水头到流量的反馈系数F可以改变系统的阻尼特性，从而可以模拟真实引水管道中水击波的衰减过程。即使引入一个很小的反馈系数，对系统增加阻尼的作用也十分明显，但对于过渡过程后的第一个水击波的幅值影响不大。因此，在水轮机引水系统中引入水头到流量的反馈系数F是可行的，也是合理的。

在一些实施例中，进一步，获得所述水轮机引水系统的z传递函数之后还包括，通过重采样技术，将所述z传递函数转换为连续系统的传递函数，再将所述连续系统的传递函数按新的采样周期离散化，达到引水系统计算采样周期可以任意设定的目的，实现引水系统和发电机-负载模型仿真时采样周期的统一。基于所述水轮机引水系统的z传递函数获得具有可调整计算步长的水轮机引水系统离散化数学模型：

其中，

a＝T_sn+T_so

b＝T_sn-T_so

c＝F+C_a

d＝F-C_a

T_so＝L/v

C_a＝T_so/T_w＝L/(vT_w)

A＝b²-a²＝(T_sn-T_so)²-(T_sn+T_so)²＝-4T_snT_so

B＝(a²-b²)＝-(b²-a²)＝4T_snT_so

T_sn为所述新的采样周期，T_so为所述z传递函数对应的采样周期。

所述水轮机引水系统离散化数学模型通过如下方式得到：

为使得采样周期转化前后的系统频率特性尽可能保持一致，转化方法采用tustin法(双梯形变换法)。

所述z传递函数对应的采样周期为T_so，重采样后系统对应的所述新的采样周期为T_sn，对

的重采样算法进行如下推导：

替换离散系统传递函数中的变量z。然后用新的采样周期T_sn及

替换连续传递函数中的变量s。将式(49)代入式(48)，可得

整理得

令a＝T_sn+T_so，b＝T_sn-T_so，则有

将其代入到下式

得

进一步化简整理，得

为书写方便并不致引起异议，将直接改写成z，有

其中：

及

将式(50)简洁化，进一步令：

则式(50)又可写成

式(54)是所述具有可调整计算步长(采样周期)的水轮机引水系统离散化数学模型，其中的参数由式(52)和式(53)确定。

在一个实施例中，所述z传递函数对应的采样周期T_so＝L/v＝200/1000＝0.2s，指定的所述新的采样周期T_sn＝0.01s，所述水头到流量的反馈系数取为F＝0.01；利用式(54)得所述具有可调整计算步长(采样周期)的水轮机引水系统离散化数学模型为：

在MATLAB的Simulink中，建立上式模型，得系统的时域响应如图12所示，1201为导叶开度，1202为水压相对偏差值。

可见，引水系统模型经重采样后，仍然可以体现弹性水击的基本特性，并且可以任意调整计算步长(采样周期)。

在一些实施例中，进一步，所述获得具有可调整计算步长的水轮机引水系统离散化数学模型之后还包括，对所述水轮机引水系统离散化数学模型进行修正：在所述具有可调整计算步长的水轮机引水系统离散化数学模型的分子常数项上补入一个水头稳态误差系数K_e，修正后引水系统数学模型为：

所述水头稳态误差系数

其中，h_A(∞)是人为设定的值。

所述水头稳态误差系数K_e具体通过如下方式得到：

为能够用离散的线性模型描述水轮机的引水系统，取消运动方程中的摩擦损失项，即令水头损失相对值这与实际情况不符。为克服这一缺陷，对所述具有可调整计算步长(采样周期)的水轮机引水系统离散化数学模型(54)进行适当的修正。由离散系统终值定理

由于A+B＝0，故因此，在式(54)的分子常数项上补入一个水头稳态误差系数K_e，形成一定的稳态误差，即

由终值定理

因而

由于

故

式(55)即为所述补入水头稳态误差系数的离散化数学模型。

显然，系统引入了水头稳态误差系数K_e后，可以人为地调整其在稳态时的水头损失值。

在一个实施例中，在稳态时，使得水轮机在额定流量时的水头损失为h_A(∞)＝-0.02，则

T_sn＝0.01

带入到式(56)

这时，系统传递函数为：

在MATLAB的Simulink中，建立上式模型，得系统的时域响应如图13所示，1301为导叶开度，1302为水压相对偏差值，1303为稳态误差。

可见，在(54)式人为引入水头稳态误差系数K_e后，形成的(55)式可以造成一个与实际情况更为相符的稳态水头损失。

另一个实施例，获取某个电站基本技术数据，进行基本参数和变量定义、基本常数定义、计算用变量和常数定义、基本关系式定义：

等值引水管道长：L＝200米

水击波速：v＝1000米/秒

计算重采样步长：T_sn＝0.01秒

水流加速时间常数：T_w＝1.9秒

引入水头到流量的反馈系数F＝0.1

加入水头损失修正为h_A(∞)＝-0.02

经过步骤101～104之后，再进行所述具有可调整计算步长(采样周期)的水轮机引水系统离散化数学模型的求取，进一步求取所述补入水头稳态误差系数的离散化数学模型，在MATLAB的Simulink下，建立所述某个电站的所述补入水头稳态误差系数的修正后引水系统数学模型，经过计算，对结果进行处理后，得到图14至图16，分别示出了真实机组与实时仿真机组甩25％、50％、75％负荷过程中相同导叶开度变化过程对应的蜗壳压力变化趋势。其中，如图14所示，1401为仿真蜗壳压力，1402为实测蜗壳压力，1403为导叶开度；如图15所示，1501为仿真蜗壳压力，1502为实测蜗壳压力，1503为导叶开度；如图16所示，1601为仿真蜗壳压力，1602为实测蜗壳压力，1603为导叶开度。

比较图14至图16真机和仿真机输出波形图易知，在甩负荷后，两者机组频率的上升曲线和蜗壳压力变化曲线趋势基本一致。且对于不同水平的甩负荷，真机与模型间的过渡过程差异表现是不同的，这也体现了非线性动态系统的重要特征。总体来说，甩负荷过程中的机组转速的模型逼近程度更好一些。甩25％负荷时，模型的转速输出与真实机组的转速输出两者过渡过程几乎重合。甩50％负荷时，两者转速最大差异发生在转速上升最高点(约为4％)，但其后的过渡过程几乎重合。甩75％负荷时，两者的转速上升最高点的数值几乎相同，最大差异仅约为1％，但其后过渡过程又出现明显差异，最大差异约为2～4％。在甩较小负荷时(<50％负荷)，两者的压力上升最大值几乎相同，且过程趋势是一致的。但在甩75％负荷时，模型蜗壳水压略微大于真机的蜗壳水压，单过程趋势基本一致。因此可以看出采用基于改进后特征算法在水轮机调节系统实时仿真中的应用结果使得仿真过程比较贴近于机组实际过渡过程。

另一方面，本发明实施例还提供了水轮机调节系统实时仿真装置。参考图2，为本发明实施例的仿真装置结构示意图。

所述水轮机调节系统实时仿真装置，包括：

定义模块201，用于获取水电站基本技术数据，利用所述基本技术数据进行基本参数和变量定义、基本常数定义、计算用变量和常数定义、基本关系式定义；

特征线建模模块202，用于利用所述基本参数和变量、所述基本常数、所述计算用变量和常数、所述基本关系式建立基于特征线法的所述水电站的水轮机引水系统的水击特征线数学模型；

第一绘图模块203，用于利用所述水击特征线数学模型建立基本特征线法系统计算原理框图；

第二绘图模块204，用于在所述基本特征线法系统计算原理框图基础上，建立去掉水头损失项、引入水头到流量的反馈的线型模型计算原理框图，依据所述线型模型计算原理框图写出所述水轮机引水系统的z传递函数，获得改进的线型模型；

仿真模块205，用于利用所述改进的线型模型进行建模仿真。

在一些实施例中，进一步，在所述第二绘图模块204之后还可以包括：

调整计算步长模块，用于通过重采样技术，基于所述水轮机引水系统的z传递函数获得具有可调整计算步长的水轮机引水系统离散化数学模型。

具体的方法为：将所述z传递函数转换为连续系统的传递函数，再将所述连续系统的传递函数按新的采样周期离散化，获得具有可调整计算步长的水轮机引水系统离散化数学模型。

在一些实施例中，进一步，在所述调整计算步长模块之后还可以包括：

修正模块：用于对所述水轮机引水系统离散化数学模型进行修正，获得修正后引水系统数学模型，方法是，在所述具有可调整计算步长的水轮机引水系统离散化数学模型的分子常数项上补入一个水头稳态误差系数K_e。

上述实施例的装置用于实现前述实施例中相应的方法，并且具有相应的方法实施例的有益效果，在此不再赘述。

又一方面，本发明实施例还提供了水轮机调节系统实时仿真电子设备。

所述水轮机调节系统实时仿真电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现上述实施例中所述的基于改进后特征线法的水轮机调节系统实时仿真方法。

所属领域的普通技术人员应当理解：以上任何实施例的讨论仅为示例性的，并非旨在暗示本公开的范围(包括权利要求)被限于这些例子；在本发明的思路下，以上实施例或者不同实施例中的技术特征之间也可以进行组合，步骤可以以任意顺序实现，并存在如上所述的本发明的不同方面的许多其它变化，为了简明它们没有在细节中提供。

另外，为简化说明和讨论，并且为了不会使本发明难以理解，在所提供的附图中可以示出或可以不示出与集成电路(IC)芯片和其它部件的公知的电源/接地连接。此外，可以以框图的形式示出装置，以便避免使本发明难以理解，并且这也考虑了以下事实，即关于这些框图装置的实施方式的细节是高度取决于将要实施本发明的平台的(即，这些细节应当完全处于本领域技术人员的理解范围内)。在阐述了具体细节(例如，电路)以描述本发明的示例性实施例的情况下，对本领域技术人员来说显而易见的是，可以在没有这些具体细节的情况下或者这些具体细节有变化的情况下实施本发明。因此，这些描述应被认为是说明性的而不是限制性的。

尽管已经结合了本发明的具体实施例对本发明进行了描述，但是根据前面的描述，这些实施例的很多替换、修改和变型对本领域普通技术人员来说将是显而易见的。例如，其它存储器架构(例如，动态RAM(DRAM))可以使用所讨论的实施例。

本发明的实施例旨在涵盖落入所附权利要求的宽泛范围之内的所有这样的替换、修改和变型。因此，凡在本发明的精神和原则之内，所做的任何省略、修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种水轮机调节系统实时仿真方法，其特征在于，包括：

获取水电站基本技术数据，利用所述基本技术数据进行基本参数和变量定义、基本常数定义、计算用变量和常数定义、基本关系式定义；

利用所述基本参数和变量、所述基本常数、所述计算用变量和常数、所述基本关系式建立基于特征线法的所述水电站的水轮机引水系统的水击特征线数学模型；

利用所述水击特征线数学模型建立基本特征线法系统计算原理框图；

在所述基本特征线法系统计算原理框图基础上，建立去掉水头损失项、引入水头到流量的反馈的线型模型计算原理框图，依据所述线型模型计算原理框图得到所述水轮机引水系统的z传递函数，获得改进的线型模型；

利用所述改进的线型模型进行建模仿真。

2.根据权利要求1所述的水轮机调节系统实时仿真方法，其特征在于，所述获取水电站基本技术数据的步骤中，通过实地测量、监控数据采集的方法获取所述水电站基本技术数据；

所述基本参数和变量为：

v—水击波速；L—等效管道长；D—等效管道直径；H_r—计算额定净水头；Q_r—计算额定流量；Q_A—水轮机入口流量；Q_B—水库入口流量；F—等效管道截面积；g—重力加速度；T_s—采样周期；y—接力器位置相对值；a—导叶开度相对值；j—轮叶开度相对值；H—水轮机实际作用水头；N₁₁—水轮机单位转速；Q₁₁—水轮机单位流量；M₁₁—水轮机单位力矩；—水轮机力矩相对值；D₁—水轮机转轮直径；

所述基本常数为：

T_r＝2L/v—水击相长；h_w＝vQ_r/(2gFH_r)—管道常数；T_w＝h_wT_r＝LQ_r/(gFH_r)—水流加速时间常数；

所述计算用变量和常数为：

—水轮机入口相对流量；—水库入口相对流量；—水轮机入口相对压力上升值；—水头损失系数；C_a＝1/(2h_w)—计算常数；

所述基本关系式为：

a＝f_a(y)—导叶开度相对值与接力器位置相对值间的非线性关系；

—水轮机单位转速计算式；

Q₁₁＝f_Q(N₁₁,a)—水轮机单位流量非线性特性；

—水轮机流量相对值；

M₁₁＝f_M(N₁₁,a)—水轮机单位力矩非线性特性；

—水轮机输出力矩相对值；

所述利用所述基本参数和变量、所述基本常数、所述计算用变量和常数、所述基本关系式建立基于特征线法的所述水电站的水轮机引水系统的水击特征线数学模型的步骤中，所述水电站的水轮机引水系统的水击特征线数学模型是单管单机简单情形下的；

所述水击特征线数学模型为：

其中，引水管道末端假定为A点，引水管道首端假定为B点，作用于所述水轮机的静水头选为基准水头H_r，忽略管道损失因素，h_B,k＝1，T_s＝L/v，C_a＝T_s/T_w，

3.根据权利要求1所述的水轮机调节系统实时仿真方法，其特征在于，所述建立去掉水头损失项、引入水头到流量的反馈的线型模型计算原理框图的步骤中，为所述水轮机相对流量；所述水轮机引水系统的z传递函数为：

其中，F为水头到流量的反馈系数，K_l＝0。

4.根据权利要求3所述的方法，其特征在于，所述依据所述线型模型计算原理框图写出所述水轮机引水系统的z传递函数之后还能够包括，通过重采样技术，基于所述水轮机引水系统的z传递函数获得具有可调整计算步长的水轮机引水系统离散化数学模型：

其中，

a＝T_sn+T_so

b＝T_sn-T_so

c＝F+C_a

d＝F-C_a

T_so＝L/v

C_a＝T_so/T_w＝L/(vT_w)

A＝b²-a²＝(T_sn-T_so)²-(T_sn+T_so)²＝-4T_snT_so

B＝(a²-b²)＝-(b²-a²)＝4T_snT_so

T_sn为所述新的采样周期，T_so为所述z传递函数对应的采样周期。

5.根据权利要求4所述的方法，其特征在于，所述获得具有可调整计算步长的水轮机引水系统离散化数学模型之后还能够包括，对所述水轮机引水系统离散化数学模型进行修正，获得修正后引水系统数学模型：在所述具有可调整计算步长的水轮机引水系统离散化数学模型的分子常数项上补入一个水头稳态误差系数K_e；

修正后引水系统数学模型为：

所述水头稳态误差系数

其中，h_A(∞)是人为设定的值。

6.根据权利要求4所述的方法，其特征在于，所述获得具有可调整计算步长的水轮机引水系统离散化数学模型之后，能够建立所述具有可调整计算步长的水轮机引水系统数学模型进行仿真。

7.根据权利要求5所述的方法，其特征在于，所述对所述水轮机引水系统离散化数学模型进行修正之后，能够建立所述修正后引水系统数学模型进行仿真。

8.一种水轮机调节系统实时仿真装置，其特征在于，包括：

定义模块，用于获取水电站基本技术数据，利用所述基本技术数据进行基本参数和变量定义、基本常数定义、计算用变量和常数定义、基本关系式定义；

特征线建模模块，用于利用所述基本参数和变量、所述基本常数、所述计算用变量和常数、所述基本关系式建立基于特征线法的所述水电站的水轮机引水系统的水击特征线数学模型；

第一绘图模块，用于利用所述水击特征线数学模型建立基本特征线法系统计算原理框图；

第二绘图模块，用于在所述基本特征线法系统计算原理框图基础上，建立去掉水头损失项、引入水头到流量的反馈的线型模型计算原理框图，依据所述线型模型计算原理框图写出所述水轮机引水系统的z传递函数，获得改进的线型模型；

仿真模块，用于利用所述改进的线型模型进行建模仿真。

9.一种水轮机调节系统实时仿真电子设备，其特征在于，包括：存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现如权利要求1至7任意一项所述的方法。