はじめに

制水口式調圧水槽の設計において用いられるVogt-Forchheimer式は、指定した制水口径および立坑径に対し貯水池水位からの最大上昇水位を求めるものであり、制水口径および立坑径をパラメータとして最大上昇水位を計算し、条件を満足するよう設計諸元を定めるために用いることができる。

Vogt-Forchheimer式は非線形方程式であり、二分法あるいはNewton-Raphson法で解くことができるが、ここでは、フローが単純な二分法を用いたプログラムを紹介する。二分法も自前実装とScipyのBrent法を用いる方法が考えられるが、エラー処理を簡単に行うため、紹介するプログラムでは自前実装としている。

具体的には、ポート径が非常に小さい場合、設定した範囲に解が存在しない。この場合Scipyではエラーとなり計算がストップするため、エラーを回避するための事前処理が必要となる。これに対し、自前関数であれば解が存在しないことがわかった段階で収束計算のループを抜け、解が存在しないことを示す仮の値を出力させることにより、後処理を簡単に行うことができる。

出力

画面出力では、参考データとして、各立坑直径に対する自由サージの最高上昇水位と、水面振動の安定性確保のための要求値を表示している。下の例では、立坑径は8.9mから32.8mの間にする必要があり、立坑径設計値として、9mを選定している。また制水口径（ポート径）は、最適値より少し大きい2.1ｍを選定している。

* Free surge water rise
Ds= 7    z*=40.002
Ds= 8    z*=35.002
Ds= 9    z*=31.113
Ds=10    z*=28.002
Ds=11    z*=25.456
Ds=12    z*=23.335
* Requirements for stability
Hg/6=32.083m  # for static stabity
Dsr1= 5.1m    # for dynamic stability (1)
Dsr2= 8.9m    # for dynamic stability (2)
Dsr3=32.8m    # for critical discharge

f:id:damyarou:20190517102113j:plain

設計理論

Basic Equations and Conditions for Design of Restricted Orifice Surge Tank

Vogt-Forchheimer's Equation for Estimation of the Maximum Water Level Rise

For the estimation of the maximum water level rise from reservoir water level for the total load interception in the restricted orifice surge tank, Vogt-Forchheimer's equations have been used. Vogt-Forchheimer's equations shown below are solutions of the fundamental differential equations of restricted orifice surge tank which are shown in a post of https://damyarou.hatenablog.com/entry/2019/05/18/110312.

$\begin{equation} \begin{cases} & m'\cdot k_0 < 1 & & (1+m'\cdot z_m)-\ln(1+m'\cdot z_m) &=(1+m'\cdot h_0)-\ln(1-m'\cdot k_0) \\ & m'\cdot k_0 > 1 & & (m'\cdot |z_m|-1)+\ln(m'\cdot |z_m|-1)&=\ln(m'\cdot k_0-1)-(m'\cdot h_0+1) \end{cases} \end{equation}$

$\begin{align} & h_0=c\cdot v_0{}^2 \\ & k_0=\cfrac{1}{2g}\left(\cfrac{Q_0}{C_d F_p}\right)^2 \\ & m'=\cfrac{2g F (h_0+k_0)}{L f v_0{}^2} \end{align}$

$z_m$	maximum water level rise from reservoir water level (positive is in the downward direction)
$h_0$	total head loss of the pressure tunnel (including the velocity head at the bottom of the surge tank)
$k_0$	resistance head at the port
$v_0$	flow velocity of the pressure tunnel
$c$	head loss coefficient ( $h_0=c\cdot v_0{}^2$ )
$Q_0$	maximum discharge
$F_p$	section area of the port
$C_d$	discharge coefficient of the port
$L$	length of the pressure tunnel
$f$	section area of the pressure tunnel
$F$	section area of the surge tank shaft
$g$	gravity acceleration

To obtain the maximum water rise of $z_m$ , Bisection method can be applied considering below equations.

$\begin{align} f(z)= \begin{cases} &\{(1+m'\cdot z)-\ln(1+m'\cdot z)\} &-\{(1+m'\cdot h_0)-\ln(1-m'\cdot k_0)\} &=0 & & (m'\cdot k_0 < 1) \\ &\{(m'\cdot |z|-1)+\ln(m'\cdot |z|-1)\}&-\{\ln(m'\cdot k_0-1)-(m'\cdot h_0+1)\} &=0 & & (m'\cdot k_0 > 1) \end{cases} \end{align}$

In above equation, the following conditions have to be satisfied everytime. These are for maintaining the positive values in the logarithm paragraph in the function of $f(z)$ .

$\begin{equation} \begin{cases} &0 < 1+m'\cdot z & & (m'\cdot k_0 < 1) \\ &0 < m'\cdot |z|-1 & & (m'\cdot k_0 > 1) \end{cases} \end{equation}$

Considering above conditions, the minimum limit of $z$ in Bisection method can be set as followings.

$\begin{equation} \begin{cases} &z_1=-\cfrac{1}{m'}+0.0001 & & (m'\cdot k_0 < 1)\\ &|z_1|=\cfrac{1}{m'}+0.0001 & & (m'\cdot k_0 > 1) \end{cases} \end{equation}$

And the maximum water rise in the restricted oriffice surge tank have to be smaller than the maximum water rise of free surging. Therefore, the maximum water rise of free surging shown below can be set as the maximum limit of $z$ in Bisection method.

$\begin{equation} z_2=v_0\cdot \sqrt{\cfrac{L f}{g F}} \end{equation}$

Using $z_1$ and $z_2$ as initial values for Bisection method, Vogt-Forchheimer's Equation can be solved.

Conditions for Reservoir Water Level and Roughness in the calculation

When a headrace surge tank is considered as an example, followings can be said:

(1) for load interception (turbine discharge decreasing)
In this case, lower head loss of headrace tunnel gives higher USWL (up-surge water level), because of higher inertia force with direction to the surge tank from reservoir.

(2) for load increase (turbine discharge increasing)
In this case, a surge tank provides the water to the penstock. Therefore, higher head loss of headrace tunnel gives lower DSWL (down-surge water level), because of higher resistance for pulling the water to the surge tank.

Considering above, following conditions are required for hydraulic design of surge tank in order to achieve the safety design in Japan standard.

Conditions for Reservoir water Level and Roughness
Discharge	water level and roughness	headrace Surge tank	Tailrace Surge Tank
Load interception	Reservoir water level Checked water level in the Tank Variation of roughness	HWL of Upper res. Up surge WL. -0.0015	LWL of Lower res. Down surge WL. -0.0015
Load rapidly increase	Reservoir water level Checked water level in the Tank Variation of roughness	LWL of Upper res. Down surge WL. +0.0015	HWL of Lower res. Up surge WL. +0.0015
Input interception	Reservoir water level Checked water level in the Tank Variation of roughness	LWL of Upper res. Down surge WL. -0.0015	HWL of Lower res. Up surge WL. -0.0015
(note) + Above values of roughness variable are for concrete. + In case of steel lining, roughness variable is 0.001 instead of 0.0015 for concrete. + In case of no lining, roughness variable is 0.003 instead of 0.0015 for concrete.

Thoma-Schuller's Stability Conditions against water Level Vibration

For the stability against the water level vibration, following condition shall be satisfied.

$\begin{align} &\text{Static stability condition} & &h_0 < \cfrac{H_g}{3} \sim \cfrac{H_g}{6} \\ &\text{Dynamic stability condition} & &F > \cfrac{L f}{c (1+\eta) g H_g} \sim \cfrac{L f}{2 c g (H_g-z_m)} \end{align}$

where, $H_g$ means the gross head, and $\eta=k_0\,/\,h_0$ .

Critical Discharge by Calame-Garden

The critical discharge by calame-Garden is defined as follow.

$\begin{equation} Q_c=\cfrac{1}{c}\left(\cfrac{1}{2g}\cdot \cfrac{L f^3}{F \eta}\right)^{1/2} \end{equation}$

If the critical discharge of $Q_c$ is smaller than the maximum discharge of $Q_0$ , the maximum water level rise for the critical discharge of $Q_c$ shall be checked. Usually, the surge tank shaft diameter is defined to satisfy the condition of $Q_0 \lt Q_c$ .

Required Conditions for Section Area of Surge Tank Shaft

Considering the conditions for the stability and the critical discharge, the section area of surge tank shaft of $F$ is defined to satisfy the following conditions.

$\begin{align} &F > F_1=\cfrac{h_0 L f}{c (h_0+k_0) g H_g} \\ &F > F_2=\cfrac{L f}{2 c g (H_g-z_m)} \\ &F < F_3=\cfrac{h_0 L f^3}{2 g c^2 k_0 Q_0{}^2} \end{align}$

Resistance of the Port

The resistance of the port $|z_{rm}'|$ can be calculated as follow.

$\begin{equation} |z_{rm}'|=k_0-h_0 \end{equation}$

Relationship between Port Resistance and Maximum water Level Rise

The relationship between Port Resistance $|z_{rm}'|$ and Maximum water Level Rise $|z_m|$ are shown below.

$\|z_m\| < \|z_{rm}'\|$	the port diameter is small because of large port resistance
$\|z_m\| = \|z_{rm}'\|$	the port diameter is optimal
$\|z_m\| > \|z_{rm}'\|$	the port diameter is large because of small port resistance

where, $|z_m|$ is calculated from Vogt-Forchheimer's equations.

Usually, the port diameter is defined as a value which is little bit larger than the optimal point.

f:id:damyarou:20190518100947p:plain

入力データ作成

下の表は、サージング解析のための入力データ作成・整理のためのものですが、Vogt-Forchheimer式を解くためには、少なくとも負荷遮断（load interception）の $h_0$ までと、流量係数 $C_d$ の欄までを完成させておかなければなりません。

負荷遮断（load interception）の条件で、Vogt-Forchheimer式を解くことにより、下表のサージタンク立坑径（断面積： $F$ ）とポート径（断面積： $F_p$ ）を決定します。

事前準備

Case	Load interception	Load increase
* Conditions for Loss calculation
Discharge	49.0 m3/s	49.0 m3/s
Reservoir water Level	EL.282.5 (FSL)	EL.275.0 (MOL)
Roughness coefficient	n=0.014	n=0.014
Head loss of pressure tunnel	3.409 m	3.801 m
* Conditions for Surging Analysis
Discharge	49.0 m3/s	49.0 m3/s
Reservoir water Level	EL.282.5 (FSL)	EL.275.0 (MOL)
Modified roughness coefficient	n=0.0125 (0.014-0.0015)	n=0.0155 (0.014+0.0015)
Modified head loss of pressure tunnel $h_l$ *	2.718 m	4.659 m

$\begin{align} &\text{For load interception}& &h_l=3.409\times\left(\cfrac{49.0}{49.0}\right)^2\times\left(\cfrac{0.0125}{0.014}\right)^2=2.718 m \\ &\text{For load increase} & &h_l=3.801\times\left(\cfrac{49.0}{49.0}\right)^2\times\left(\cfrac{0.0155}{0.014}\right)^2=4.659 m \end{align}$

入力データの整理

Case	Load interception	Load increase
Gross head $H_g$	192.5 m	192.5 m
Discharge $Q$	49.0 => 0 m3/s	24.5 => 49.0 m3/s
Closing/Opening time $T$	8 sec	8 sec
Reservoir water level $RWL$	EL.282.5 (FSL)	EL.275.0 (MOL)
Length of pressure tunnel $L$	6576 m	6576 m
Section area of pressure tunnel $f$	26.162 m2	26.162 m2
Head loss of pressure tunnel $h_l$	2.718 m	4.659 m
Velocity head at bottom of surge tank $h_b$	0.179 m (v=1.873m/s)	0.179 m (v=1.873m/s)
Total head loss $h_0=h_l+h_b$	2.897 m	4.838 m
Average velocity of pressure tunnel $v_0$	1.873 m/s	1.873 m/s
Loss coefficient of pressure tunnel $c=h_0/v_0{}^2$	0.826	1.379
Section area of shaft $F$	***
Section area of port $F_p$	***
Discharge coefficient of port $C_d$	0.9	0.9

プログラム

ポイント

二分法の主要部分のコードは以下の通り。

def func(z,param):
    # h0 is global variable
    md=param[0]
    k0=param[1]
    if md*k0<1: f=((1+md*z)-np.log(1+md*z))-((1+md*h0)-np.log(1-md*k0))
    if md*k0>1: f=((md*z-1)+np.log(md*z-1))-(np.log(md*k0-1)-(md*h0+1))
    return f


def bisection(x01,x02,param):
    imax=100    # maximum number of iterations
    err=1.0e-10 # allowable error of the zero value
    x1=x01
    x2=x02
    for i in range(0,imax):
        xi=0.5*(x1+x2)
        f1=func(x1,param)
        f2=func(x2,param)
        fi=func(xi,param)
        if abs(x2-x1) < err: break # converged
        if 0<f1*f2: break          # no solution
        if f1*fi < 0.0: x2=xi
        if fi*f2 < 0.0: x1=xi
    return xi

凡例の「Ds: shaft diameter」の部分は、タイトル（title）として表示している。以下を参照。

plt.legend(loc='lower right',fontsize=fsz-2,shadow=True,title='Ds: shaft diameter').get_title().set_fontsize(fsz-2)

プログラム全文

プログラム全文を以下に示す。作図まで行うようにしている。入力データでは、設計値も入力するようにしているが、最初は適当な値を入れておき、解析結果を見て最終的な設計値を確定する。

#==============================
# Vogt-Forchheimer's eqyation
# Bisection Method
#==============================
import numpy as np
import matplotlib.pyplot as plt


# global variables
g=9.8       # gravity acceleration
q0=49       # discharge of pressure tunnel
cd=0.9      # discharge coefficient of port
h0=2.897    # head loss of pressure tunnel
tl=6559+17  # length of pressure tunnel including untake tunnel
ft=26.162   # section area of pressure tunnel
v0=q0/ft    # flow velocity of pressure tunnel
hg=192.5    # gross head (284.5-92)
dds=9.0     # design shaft diameter
ddp=2.1     # design port diameter
shaft=np.arange(7,13,1).astype(np.float64)      # shaft diameter: parameter
port=np.arange(0.1,6.1,0.05).astype(np.float64) # port diameter: parameter
tstr='Surge Tank (h0=2.897m)'
xmin=0; xmax=6; dx=1
ymin=0; ymax=40; dy=10
xsize=8; ysize=6


def drawfig(port,zz,zd,dzz):
    fsz=16
    fig=plt.figure(figsize=(xsize,ysize),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Port diameter (m)')
    plt.ylabel('Maximum water level rise (m)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#999999',linestyle='solid')
    plt.title(tstr,loc='left',fontsize=fsz)
    for i,sh in enumerate(shaft):
        # Since zm is monotone increasing function for port diameter, 
        # non-converged values of zm are omitted from plot data
        n=np.argmin(zz[i,:])
        plt.plot(port[n:],zz[i,n:],color='#000000',lw=1)
        xs=xmax
        ys=zz[i,-1]
        ss='Ds={0:.0f}m'.format(sh)
        plt.text(xs,ys,ss,color='#000000',ha='right',va='bottom',fontsize=fsz-2)
    plt.plot([xmin],[ymin],'-',color='#000000',lw=1,label="Water level rise")
    plt.plot(port,zd,'--',color='#000000',lw=1,label="$|z_m'|=k_0-h_0$")
    plt.plot([ddp],[dzz],'o',color='#000000',ms=8,label='Design point')   
    plt.plot([xmin,ddp],[dzz,dzz],'-',color='#000000',lw=0.5)
    plt.plot([ddp,ddp],[ymin,dzz],'-',color='#000000',lw=0.5)
    ss1='Dp={0:.1f}m'.format(ddp)
    ss2='Zmax={0:.3f}m'.format(dzz)
    plt.text(ddp,ymin,ss1,color='#000000',ha='left',va='bottom',rotation=90,fontsize=fsz-2)
    plt.text(xmin,dzz,ss2,color='#000000',ha='left',va='bottom',rotation=0,fontsize=fsz-2)
#    plt.legend(loc='lower right',fontsize=fsz-2,shadow=True)
    plt.legend(loc='lower right',fontsize=fsz-2,shadow=True,title='Ds: shaft diameter').get_title().set_fontsize(fsz-2)
    plt.tight_layout()
    fnameF='fig_surge_zmax.jpg'
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
    plt.show()

    
def func(z,param):
    # h0 is global variable
    md=param[0]
    k0=param[1]
    if md*k0<1: f=((1+md*z)-np.log(1+md*z))-((1+md*h0)-np.log(1-md*k0))
    if md*k0>1: f=((md*z-1)+np.log(md*z-1))-(np.log(md*k0-1)-(md*h0+1))
    return f


def bisection(x01,x02,param):
    imax=100    # maximum number of iterations
    err=1.0e-10 # allowable error of the zero value
    x1=x01
    x2=x02
    for i in range(0,imax):
        xi=0.5*(x1+x2)
        f1=func(x1,param)
        f2=func(x2,param)
        fi=func(xi,param)
        if abs(x2-x1) < err: break # converged
        if 0<f1*f2: break          # no solution
        if f1*fi < 0.0: x2=xi
        if fi*f2 < 0.0: x1=xi
    return xi

           
def cal_z(param,zs):
    md=param[0]
    k0=param[1]
    eps1=1.0e-10
    zmax=0
    z2=zs # maximum limit of assumed solution
    if md*k0<1: z1=-1/md+eps1 # minimum limit of assumed solution
    if md*k0>1: z1= 1/md+eps1 # minimum limit of assumed solution
    zmax=bisection(z1,z2,param)
    return zmax


def dreq(dzz):
    # required diameter
    fp=np.pi*ddp**2/4
    cc=h0/v0**2
    k0=1/2/g*(q0/cd/fp)**2
    f1d=h0*tl*ft/cc/(h0+k0)/g/hg;       dr1=np.sqrt(4*f1d/np.pi) # Thoma-Schuller's dynamic stability (1)
    f2d=tl*ft/2/cc/g/(hg-dzz);          dr2=np.sqrt(4*f2d/np.pi) # Thoma-Schuller's dynamic stability (2)
    f3d=h0*tl*ft**3/2/g/cc**2/k0/q0**2; dr3=np.sqrt(4*f3d/np.pi) # Calame-Garden's critical discharge
    hh=hg/6
    print('* Requirements for stability')
    print('Hg/6={0:6.3f}m  # for static stabity'.format(hh))
    print('Dsr1={0:4.1f}m    # for dynamic stability (1)'.format(dr1))
    print('Dsr2={0:4.1f}m    # for dynamic stability (2)'.format(dr2))
    print('Dsr3={0:4.1f}m    # for critical discharge'.format(dr3))


def free_surge(fs):
    zs=v0*np.sqrt(tl*ft/g/fs) # free surge water rise
    print('* Free surge water rise')
    for i in range(len(fs)):
        print('Ds={0:2.0f}    z*={1:6.3f}'.format(shaft[i],zs[i]))
    

def main():
    # h0 is global variable
    fs=np.pi*shaft**2/4 # shaft diameter
    fp=np.pi*port**2/4  # port diameter
    zz=np.zeros((len(shaft),len(port)),dtype=np.float64) # array for |zm|
    zd=np.zeros(len(port),dtype=np.float64)              # array for |zrm'|
    eps=0.0001
    dzz=0
    for i in range(len(fs)):
        for j in range(len(port)):
            k0=1/2/g*(q0/cd/fp[j])**2        # k0
            md=2*g*fs[i]*(h0+k0)/tl/ft/v0**2 # m'
            zs=v0*np.sqrt(tl*ft/g/fs[i])     # free surge water rise
            param=np.array([md,k0])
            zz[i,j]=np.abs(cal_z(param,zs))     # calculation of |zm|
            # set design water level rize at design point
            #     dzz : design water level rise
            #     dds : design shaft diameter 
            #     ddp : design port diameter
            if dds-eps<shaft[i] and shaft[i]<dds+eps:
                if ddp-eps<port[j] and port[j]<ddp+eps:
                    dzz=zz[i,j]
    zd=1/2/g*(q0/cd/fp)**2-h0 # |zrm'|=k0-h0
    free_surge(fs) # free surge water rize (reference)
    dreq(dzz)      # requirement for stability
    drawfig(port,zz,zd,dzz) # drawing


#==============
# Execution
#==============
if __name__ == '__main__': main()