刘松舰¹，张志昌¹，潘海荣²，王裕征³

（1．西安理工大学，陕西西安710048；2．青海省万立水电有限责任公司，
青海西宁810001；3．青海省水利水电设计研究院，青海西宁810001）

    摘　要：本文从调压井的基本方程出发，采用Fortran语言和龙格库塔方法编制了阻抗室调压井水位波动的数值计算程序，用此程序计算了青海地盘子水电站阻抗式调压井两种阻抗孔面积下在增荷或甩荷时的水位波动过程，并与模型试验结果相比较，结果表明计算结果与实测结果相符。这种计算方法对其它实际工程有参考价值。
    关键词：调压室；阻抗式；水位波动；数值计算 　　

1　引　言
　　在有压引水式和混合式水电站中，当压力引水道较长，往往设置调压井来减小引水道中的水击压强，调压井中的水位波动为非恒定流问题。对这种非恒定流问题传统的解决方法是图解法和解析法。解析法快捷、简单，但只能计算最高涌波水位和第二振幅，图解法虽然能计算水位波动过程，但计算繁琐而且精度差。随着计算机软、硬件性能的大幅提高，利用数值模拟计算调压井水位波动已成为调压井设计计算的主要方法。
　　目前，利用数值计算调压井水位波动过程的程序主要为BASIC语言^{［1，2，3］}，计算方法有欧拉法、改进欧拉法和龙格库塔法。这些计算方法和程序在目前的WINDOWS系列操作系统下使用不方便，数据输出、处理方式单一，不便于其它数据处理软件调用。故本文采用精度较高的龙格库塔方法，利用最适用于数值计算领域的FORTRAN语言在VISUALFORTRAN编译平台下编制了阻抗式调压井水力瞬变过程的数值计算程序。该程序使用方便、灵活、便于掌握，计算速度快、精度
高。本文将这一程序用于青海省地盘子水电站调压系统的设计计算，结果表明计算结果和实验结果吻合。
2　调压井水位波动的基本方程
　　调压井水位波动的基本微分方程为：
    能量方程：

式中：L——引水
     Z——调压井波动水位（以库水位为基准面，向下为正）；
     V——引水道流速；
     f——引水道断面面积；
     F——调压井断面面积；
 　　h_w——水头损失；
     Q——水轮机引用流量，见图1。(XMW)

    由公式（1）可得阻抗式调压井的能量方程为：
  

式中　K——引水道的沿程水头损失和局部水头损失系数之和；
     K₀——为调压井阻抗孔的阻抗系数；
     Q_f——进出调压井的流量（以进入调压井为正）。
　　由公式（2）可得阻抗式调压井的连续方程为：

3　阻抗式调压井数学模型差分方程
　　对上述方程求解可采用精度较高的四阶龙格库塔方法。其基本公式为：

根据以上公式，可得阻抗式调压井数学模型的差分议程为：

式中　Q_e（n）——为第n时段引水道中的流量；
　　Z_（n）——第n时段调压井中的水位。
　　式中的K₁、K₂、K₃、K₄以及S₁、S₂、S₃、S₄可以根据式（7）、（8）、（9）、（10）求得。
4　边界条件
　　水轮机导叶按线性规律关闭或开启，当水轮机突然甩负荷时，导叶关闭时流量变化可用下式［4］表示：

式中　Q_（n）——电站引用流量；
      Q₀——水轮机甩荷前正常运行时的最大引用流量；
      ΔT——计算时间步长；
      T₀——导叶关闭时间；
　　  n——计算步数。
　　当水轮机突然增荷时，导叶开启时流量变化可用下式［4］表示：

式中　Q₁——水轮机增荷前电站的引用流量；
      T₁——电站引用流量增至Q₁时导叶开启所需的时间；
     ——电站引用流量增至Q₀时导叶开启所需的时间；
     －T₁——电站由Q₁突增至Q₀时导叶开启所需的时间；
　　 Z₀——库水位。
　　为了真实反映阻抗孔的阻抗系数，本文采用的阻抗系数为实测值。
　　以上计算方法也适用于简单式调压井，对于简单式调压井，其阻抗系数为零。
5　试验与验证
　　为了验证数值计算成果，现以青海省地盘子水电站调压井水位波动为例进行验证。该水电站为典型的引水式电站，发电系统由压力引水隧洞、调压井、压力钢管、厂房、尾水渠、升压站及管理所组成。压力隧洞从隧洞进口至调压井轴线总长4 771．74 m。比降i＝1／500，采用圆形断面，直径为4．4 m，糙率为0．015，设计引用流量33 m³／s。调压井拟采用圆孔阻抗式，调压井底部高程2 670．986 m，水库设计水位2 682．0 m，校核水位2 683．4 m，水轮机关闭时间4．0 s，水轮机开启时间20 s。在调压井中设有闸门孔，孔宽0．9 m、长3．4 m，它实际上增加了阻抗孔的面积，故在计算时将其与阻抗孔面积统一考虑，本文还对取消闸门孔的工况进行了计算。电站压力引水系统布置如图2所示。

　　在考虑闸门槽的情况下，已知调压井直径为12 m，阻抗孔直径为2．81 m（闸门孔面积已换算成当量直径），则采用程序计算得在校核水位四台水轮机同时弃荷时，调压井的最高、最低水位分别为2 697．408 m、2 672．979 m；在设计水位三台水轮机正常运行而第四台增荷时，调压井最高、最低水位分别为2 678．328 m、2 674．511 m。以上数值计算和试验结果见图3和图4。

　　本文还计算了闸门槽去掉后调压井的水位波动，在校核水位四台水轮机同时弃荷时，计算得调压井最高水位与最低水位分别为2 692．86 m和2 678．076 m；在设计水位三台水轮机正常运行而第四台增荷时，调压井最高水位与最低水位分别为2 678．1 m和2 674．883 m。如图5、图6所示。

　　由以上计算可知，数值计算与实测结果相符，如阻抗孔直径为2．81 m时，最高涌波水位计算值2 697．408 m，实测值2 697．586 m，最低涌波水位计算值2 672．979 m，实测值2 673．573 m。而且由波形图也可以看出实测波形图与计算波形图一致。由图中还可以看出阻抗孔面积对水位波动影响较大，阻抗孔面积越大，则阻抗系数越小，波动越明显，反之，波动减弱。本文还计算了四台机组同时开启这种最不利的开启方式，结果表明，调压井最低水位为2 664．5 m，已空底，和采用文献［5］的计算方法结果一致，实测确为空底。
6　结　论
　　（1）本文从调压井的基本方程出发，采用龙格库塔方法对阻抗式调压井的水位波动过程进行了研究，在WINDOWS环境下，使用VISUALFORTRAN编译平台，编制了FORTRAN90计算程序。该程序数据处理灵活、易于调用、便于掌握，计算速度快、精度高。
　　（2）用该程序计算青海省地盘子水电站调压井的水位波动过程，计算中阻抗系数采用实测值。计算的最高涌波水位为2 697．408 m，实测值为2697．586 m；最低涌波水位计算值为2672．979m，实测值为2 673．573 m，计算结果与实测资料相符。
　　（3）本文还计算了四台机组同时开启这种最不利的开启方式，结果表明，调压井已空底，与实验结果一致，所以在机组开启方式上要加以限制，避免出现调压井空底现象。
　　（4）本文对不同阻抗孔口尺寸的调压井的水位波动过程进行了计算，由计算可知，随着阻抗系数的增大，调压井的水位波动减弱。
　　（5）由以上计算可以看出，用本文编制的计算程序完全可以预测调压井水位波动的过程。计算精度高，满足工程实际应用。
　　（6）该程序亦可用于简单调压井的设计计算。

参考文献 　
［1］　王世泽．水电站建筑物［M］．北京：水利电力出版社，1986．
［2］　王树人．调压井水力计算理论与方法［M］．北京：清华大学出版社，1983．
［3］　刘启钊，彭守拙．水电站建筑物［M］．北京：水利电力出版社1993．
［4］　秋元得三著，支培法，等译．水击与脉动压力［M］．北京：电力工业出版社．1981．
［5］　中华人民共和国电力行业标准．DL／T5058－1996水电站调压室设计规范［S］．北京：中国电力出版社．1997．