调压室的水力计算1. 调压室断面计算当上游死水位，下游为最低水位，最小水位H min=188.9m，三台机满发，引水道糙率取最小值，压力管道糙率取最大值，通过水轮机的流量为57m3s⁄，则此时的引水隧洞水头损失的计算如表格1，压力钢管水头损失的计算如表格2。

引水道应选可能的最小糙率0.012，压力管道应选择可能的最大糙率0.013。

表格1引水隧洞水头损失表表格2压力钢管水头损失表F Tℎ>Lf2αgH1=Lf2αg(H0−ℎw0−3ℎwT)=45.548m2其中H0——最小水头损失，H0=188.9m；ℎw0——引水隧洞损失，ℎw0=17.802+0.296=18.098；ℎwT0——压力管道水头损失，ℎwT=3.110+2.805=5.915m；L——引水隧洞长度，12662m；g——重力加速度，g=9.81m/s2 f——引水隧洞面积，16.62m2。

α——引水道阻力系数v0=Qf=5716.619=3.43m s⁄α=ℎw0v02=18.0983.4302=1.5385为了保证大波动的稳定，一般要求调压室断面大于托马斯断面，初步分析时可取（1.0~1.1）F Tℎ，作为调压室的设计断面。

这里选取D=7.8m，则系数k为：F k=47.784k=F k/F Tℎ=1.052. 最高涌波水位计算按正常蓄水位时共用同一调压室的三台机组全部满载运行瞬时丢弃全部负荷（即流量由Q max=57减至流量Q=0）作为设计工况。

引水隧洞的糙率取尽可能的最小值（能耗少，涌波高）。

n=0.012引水道损失由表格1和表格2得：ℎw0=ℎw0程+ℎw0局=17.802+0.296=18.098mv0为时段开始时管中流速v0=Qf=3.43m s⁄；f为引水隧洞断面面积。

F为调压井断面面积，145.267m2；引水隧洞长L=12662m，g=9.81m s2⁄得引水道—调压室系统的特性系数。

λ=Lfv022gFℎw0=12662×16.62×3.4322×9.81×47.784×18.098=145.89令X0=ℎw0λ=0.124，X=zλ，则要求最高涌波水位z max，只需要求出X max=z maxλ即可。

X max的符号在静水位以上为负，以下为正。

ln（1+X max）−X max=−X0运用牛顿切线法求解方程的根令φ(x)=ln（1+x）−x+X0牛顿迭代公式为：x k+1=x k−φ(x)φ′(x)取迭代初值x0=−0.5，计算结果见下表表格3迭代计算结果由表格3可以看出精确到0.001，X max=z maxλ=−0.419|z max|=−λX max=61.128m根据《水电站调压室设计规范》调压室最高涌波水位以上的安全超高不宜小于1m。

所以调压室的顶高程：Z=1279+61.128+1.5=1341.63m3. 最低涌波水位计算3.1 增加负荷时的最低涌波水位上游水位取死水位，下游取最低水位。

引水隧洞糙率取最大值，取n=0.016。

水轮发电机组由2台增至3台满负荷运行。

由表1和表2引水隧洞水头损失为：ℎw0=ℎw0程+ℎw0局=31.649+0.296=31.945m|z min |ℎw0=1+(√ε−0.275√m +0.05ε⁄−0.9)(1−m)(1−m ε0.62)⁄ 式中ε=2λℎw0⁄=9.134，m = 23 ，解得：z min =50.44m所以当增加负荷时水位降到静水位下23.74m 3.2 丢弃全负荷时产生的第二波动振幅X 2+ln (1−X 2)=X max +ln (1−X max )解之X 2+ln (1−X 2)=−0.069X 2=0.327z min =λX 2=47.70m所以当丢弃全负荷时，水位降到静水位下47.70m 综上所述调压室的最低涌波水位应为静水位下50.44m 根据《水电站调压室设计规范》上游调压室最低涌波水位与调压室处压力引水道顶部之间的安全高度应不小于2～3m 。

所以压力引水道顶部高程为：Z=1274-50.44-2.5=1221.06m4. 四阶龙格库塔法4.1调压室水位的基本微分方程dz dt =(Q−Q m)F=f1(t,z,Q)dQ dt =(H R−z−KQ s|Q s|−RQ|Q|)gAL=f2(t,z,Q)已知压力管道的流量变化规律，则调压室流量等于隧洞流量减压力管道流量，即：Q s=Q−Q m。

上述方程两个未知量，Q和Z，采用四阶龙格库塔法进行逐步求解。

4.2四阶龙格库塔法计算公式如果已知t时刻的Q t，Z t值，采用以下公式求解t+∆t时刻的Q t+∆t，Z t+∆t。

Z t+∆t=Z t+16(K1+2K2+2K3+K4)K1=∆tf1(t,Z t,Q t)K2=∆tf1(t+∆t2,Z t+K12,Q t+L12)K3=∆tf1(t+∆t2,Z t+K22,Q t+L22)K4=∆tf1(t+∆t,Z t+K3,Q t+L3)Q t+∆t=Q t+16(L1+2L2+2L3+L4)L1=∆tf2(t,Z t,Q t)L2=∆tf2(t+∆t2,Z t+K12,Q t+L12)L3=∆tf2(t+∆t2,Z t+K22,Q t+L22)L4=∆tf2(t+∆t,Z t+K3,Q t+L3)4.3程序框图4.4 VBA程序代码定义函数f1：Public Function f1(t As Double, z As Double, q As Double)Dim a As Doublea = 5 'a为导叶直线关闭时间If t < a Thenf1 = (q - (57 - 57 * t / a)) / 47.784Elsef1 = q / 47.784End IfEnd Function定义函数f2：Public Function f2(t As Double, z As Double, q As Double)Dim r As Doubler = 0.0055 'r为隧洞的沿程损失和局部损失系数f2 = 9.81 * 16.61 * (1279 - z - r * Abs(q) * q) / 12662End Function调压室甩负荷时水位波动程序Sub tyssw()Dim t As DoubleDim z As DoubleDim q As Double '将t,z,q定义为实数Dim k1, k2, k3, k4 As DoubleDim l1, l2, l3, l4 As Double '将k1, k2, k3, k4,l1, l2, l3, l4定义为实数Dim dt, stoptime As Double 'dt为时间步长，stoptime为计算的调压室水位波动的时间，单位（s）stoptime = 200 '计算从导叶开始关闭到t=200s的水位波动过程i = 2 'i用于计算结果计数dt = 1t = 0q = 57z = 1260.9 't,q,z稳态时的初始值Do While t < stoptime '当时间t小于stoptime时执行以下循环过程k1 = dt * f1(t, z, q)l1 = dt * f2(t, z, q)'k2 = dt * f1(t + dt / 2, z + k1 / 2, q + l1 / 2)l2 = dt * f2(t + dt / 2, z + k1 / 2, q + l1 / 2)'k3 = dt * f1(t + dt / 2, z + k2 / 2, q + l2 / 2)l3 = dt * f2(t + dt / 2, z + k2 / 2, q + l2 / 2)'k4 = dt * f1(t + dt, z + k3, q + l3)l4 = dt * f2(t + dt, z + k3, q + l3)'t = t + dtz = z + (k1 + 2 * k2 + 2 * k3 + k4) / 6q = q + (l1 + 2 * l2 + 2 * l3 + l4) / 6 '如上是四阶龙格库塔法计算过程i = i + 1Cells(i, 1) = tCells(i, 2) = zCells(i, 3) = qLoop '循环结束标志End Sub调压室增加负荷时水位波动程序主程序相同，只需要将初始值改为Q=38，Z=1259.8，其中自定义函数改为如下：Public Function f1(t As Double, z As Double, q As Double)Dim a As Doublea = 5 'a为导叶直线关闭时间If t < a Thenf1 = (q - (38 + 19 * t / a)) / 47.784Elsef1 = (q -57)/ 47.784End IfEnd Function定义函数f2：Public Function f2(t As Double, z As Double, q As Double)Dim r As Doubler = 0.001 'r为隧洞的沿程损失和局部损失系数f2 = 9.81 * 16.61 * (1279 - z - r * Abs(q) * q) / 12662End Function4.5计算结果及调压室水位振荡曲线甩负荷增负荷 57减至0 57减至38 57减至19 19增至3819增至5738增至57最高上升水位 相应水位 61.14 1340.14 52.22 1331.22 57.32 1336.32 最低下降水位 相应水位-47.73 1231.27-46.47 1232.53-47.25 1231.75-22.86 1251.14-48.28 1225.72-30.95 1243.051245125012551260126512701275128012850200400600800100012001400增负荷时调压室水位振荡曲线122012401260128013001320134013600200400600800100012001400甩负荷时调压室水位振荡曲线Z min =1231.27T s =400sZ max =1340.14图 1流量由57减至0调压室水位振荡曲线图 2流量由38增至57调压室水位振荡曲线Z min =1243.05T s =200s5. 解析法与四阶龙格库塔法对比表格4计算结果对比单位（m）在计算最低涌波水位时。