damyarou

python, GMT などのプログラム

Python 開水路トンネルの不等流解析プログラム

開水路トンネルの不等流解析（常流）を行ったので、そのプログラムをアップしておく。 全区間、常流であるため、下流端で固定水位を与え、上流に向かって逐次水位を求めていく。

成果図

f:id:damyarou:20200126133415p:plain:w300

f:id:damyarou:20200126133433p:plain:w600

f:id:damyarou:20200126133328p:plain

f:id:damyarou:20200126133348p:plain

このプログラムでのTips

与えられた固定点数点に対し計算点を等間隔に指定する

from scipy import interpolate

scipy interpolateで線形補間により計算点を指定。 ｘ座標と同一点でｚ座標と流量が与えられているので、それらも計算点に合わせて補間する。

al=1076.0
dx=1.0
xx=np.arange(0,al+dx,dx)          # distance
f1 = interpolate.interp1d([0,983,996,1017,1041,1076],[888.222,891.300,891.300,891.300,891.300,891.500])
zz = f1(xx)
f1 = interpolate.interp1d([0,983,996,1017,1041,1076],[68.64,68.64,68.64,45.76,22.88,22.88])
qq = f1(xx)

インバート線の下を灰色で塗りつぶす

# invert line
plt.fill_between(xx,zz,zz-0.5,facecolor='#aaaaaa',alpha=0.5)
plt.plot(xx,zz,color='#000000',lw=2,label='Invert')

四角の中に文字列を記述する

 ss=''
    ss=ss+'* Input coordinates & discharge               \n'
    ss=ss+'dL(m)   x(m)  Floor level  Q(m3/s)  Remarks   \n'
    ss=ss+'----------------------------------------------\n'
    ss=ss+'   0   1076   EL.891.500   22.88   PS wall(#1)\n'
    ss=ss+' +35   1041   EL.891.300   22.88              \n'
    ss=ss+' +24   1017   EL.891.300   45.76   #2 flow-in \n'
    ss=ss+' +21    996   EL.891.300   68.64   #3 flow-in \n'
    ss=ss+' +13    983   EL.891.300   68.64              \n'
    ss=ss+'+983      0   EL.888.222   68.64   Outlet     '
    sinp=ss

上記のように表のような文字列を定義し、これを四角の中に書き込む。 表形式なのでプロポーショナルではなく等幅フォント（monospace）を指定する。

bbox_props = dict(boxstyle='square,pad=0.3', fc='#ffffff', ec='#000000', lw=1)
plt.text(-80,899.5, sinp, ha='right', va='top', rotation=0, size=fsz, bbox=bbox_props,fontfamily='monospace')

不等流の基本方程式

\begin{equation*} \left(\cfrac{Q^2}{2 g A_2{}^2}+h_2+z_2\right) - \left(\cfrac{Q^2}{2 g A_1{}^2}+h_1+z_1\right) = \cfrac{1}{2}\left(\cfrac{n^2 Q^2}{R_1{}^{4/3} A_1{}^2} + \cfrac{n^2 Q^2}{R_2{}^{4/3} A_2{}^2}\right)\cdot\Delta x \end{equation*}
  • Q : discharge
  • n : Manning's roughness coefficient
  • A_2, h_2, z_2, R_2 : flow section area, water depth, channel bottom elevation and hydraulic radius at upstream section
  • A_1, h_1, z_1, R_1 : flow section area, water depth, channel bottom elevation and hydraulic radius at downstream section
  • B : channel width

合成粗度係数

\begin{equation} n=\left\{\cfrac{\sum\left(n_i^{3/2}\cdot s_i\right)}{\sum s_i}\right\}^{2/3} \end{equation}
  • n : combined roughness coefficient
  • n_i : roughness coefficient of zone i
  • s_i : perimeter of zone i

Manning-Strickler equation

\begin{equation} n=\cfrac{k^{1/6}}{7.66 \sqrt{g}} \end{equation}
  • n : roughness coefficient
  • k : equivalent roughness in meter
  • g : gravity acceleration (=9.8 m/s^2)

矩形断面の限界水深計算式

\begin{equation*} H_c=\left(\cfrac{Q^2}{g B^2}\right)^{1/3} \end{equation*}

基本方程式の解法

基本方程式において、常流のため、下流水深がわかっており、上流に向かって水深を求めていく。 式中の上流側水深 h2 が求める水深であるが、断面積 A、動水半径 R が h の関数でありhに関する非線形方程式なので自前の二分法で解いている。

def sec(h,x,nn1,nn2):
    r0=2.9
    b0=5.5
    h0=4.0
    a1=b0*h
    s1=b0+2*h
    a2=0
    s2=0
    if 40.0<=x and h0<=h:
        a1=b0*h0
        s1=b0+2*h0
        theta=np.arcsin((h-h0)/r0)
        s2=2*r0*theta
        a2=r0**2*theta+r0**2*np.sin(theta)*np.cos(theta)
    aa=a1+a2
    rr=(a1+a2)/(s1+s2)
    n=((nn1**(3/2)*s1+nn2**(3/2)*s2)/(s1+s2))**(2/3)
    return aa,rr,n


def func(h2,h1,z1,z2,x1,x2,dx,q,nn1,nn2):
    g=9.8
    a1,r1,n1=sec(h1,x1,nn1,nn2)
    a2,r2,n2=sec(h2,x2,nn1,nn2)
    e2=q**2/2/g/a2**2+h2+z2
    e1=q**2/2/g/a1**2+h1+z1
    e3=0.5*(n1**2*q**2/r1**(4/3)/a1**2 + n2**2*q**2/r2**(4/3)/a2**2)*dx
    return e2-e1-e3


def bisection(h1,z1,z2,x1,x2,dx,q,nn1,nn2):
    ha=z1+h1-z2-0.001
    hb=6.9
    for k in range(100):
        hi=0.5*(ha+hb)
        fa=func(ha,h1,z1,z2,x1,x2,dx,q,nn1,nn2)
        fb=func(hb,h1,z1,z2,x1,x2,dx,q,nn1,nn2)
        fi=func(hi,h1,z1,z2,x1,x2,dx,q,nn1,nn2)
        if fa*fi<0: hb=hi
        if fb*fi<0: ha=hi
        #print(fa,fi,fb)
        if np.abs(hb-ha)<1e-10: break
    return hi

......

h2=bisection(h1,z1,z2,x1,x2,dx,q,nn1,nn2)

......

不等流計算プログラム

# Non-Uniform Flow Analysis (Subcritical flow)
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as plt


def drawfig(xx,zz,hh1,hh2,walltop,crown,sinp,kkk):
    if kkk==1: fnameF='fig_tailrace1.png'
    if kkk==2: fnameF='fig_tailrace2.png'

    fsz=12
    xmin=-100
    xmax=1200
    dx=100
    ymin=885
    ymax=903
    dy=1
    fig=plt.figure(figsize=(16,10),facecolor='w')
    plt.rcParams["font.size"] = fsz
    plt.xlabel('Distance x (m)')
    plt.ylabel('Water Level (EL.m)')
    plt.xlim(xmax,xmin)
    plt.ylim(ymin,ymax)
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#777777',linestyle='--')
    # invert (powerhouse)
    xs=xx[-1]; ys=zz[-1]; ss='EL.{0:7.3f}'.format(ys)
    plt.plot([xs+100,xs],[ys,ys],'-',color='#000000',lw=1)
    plt.text(xs+100,ys,ss,va='bottom',ha='left',fontsize=fsz)
    # invert (EL.891.3)
    ix=np.where(zz==891.3)[0][0]
    xs=xx[ix]; ys=zz[ix]; ss='EL.{0:7.3f}'.format(ys)
    plt.plot([xs-100,xs],[ys,ys],'-',color='#000000',lw=1)
    plt.text(xs-100,ys,ss,va='bottom',ha='right',fontsize=fsz)
    # invert (outlet)
    xs=xx[0]; ys=zz[0]; ss='EL.{0:7.3f}'.format(ys)
    plt.plot([xs-100,xs],[ys,ys],'-',color='#000000',lw=1)
    plt.text(xs-100,ys,ss,va='bottom',ha='right',fontsize=fsz)
    # No.1
    els=899
    i1=np.where(xx==1076)[0]
    plt.plot([xx[i1],xx[i1]],[zz[i1],els],'-',color='#000000',lw=1)
    plt.text(xx[i1],els,'PS wall (No.1)',ha='center',va='bottom',fontsize=fsz,rotation=90)
    # No.2    
    i2=np.where(xx==1017)[0]
    plt.plot([xx[i2],xx[i2]],[zz[i2],els],'-',color='#000000',lw=1)
    plt.text(xx[i2],els,'No.2 flow-in',ha='center',va='bottom',fontsize=fsz,rotation=90)
    # No.3
    i3=np.where(xx==996)[0]
    plt.plot([xx[i3],xx[i3]],[zz[i3],els],'-',color='#000000',lw=1)
    plt.text(xx[i3],els,'No.3 flow-in',ha='center',va='bottom',fontsize=fsz,rotation=90)
    # outlet
    plt.plot([xx[0],xx[0]],[zz[0],els-3],'-',color='#000000',lw=1)
    plt.text(xx[0],els-3,'Outlet',ha='center',va='bottom',fontsize=fsz,rotation=90)
    
    # wall top (powerhouse)
    xs=xx[-1]; ys=zz[-1]+4; ss='EL.{0:7.3f}'.format(ys)
    plt.plot([xs+100,xs],[ys,ys],'-',color='#000000',lw=1)
    plt.text(xs+100,ys,ss,va='bottom',ha='left',fontsize=fsz)
    # wall top (outlet)
    xs=xx[0]; ys=walltop[0]; ss='Invert+5.5m\nEL.{0:7.3f}'.format(ys)
    plt.plot([xs-100,xs],[ys,ys],'-',color='#000000',lw=1)
    plt.text(xs-100,ys,ss,va='bottom',ha='right',fontsize=fsz)
    # crest level of MT
    xs1=xx[0]; ys1=zz[0]
    xs2=xx[0]; ys2=887
    plt.plot([xs1,xs2],[ys1,ys2],'-',color='#000000',lw=1)
    plt.plot([xs1-100,xs],[ys2,ys2],'-',color='#000000',lw=1)
    ss='MT Wire crest level: EL.{0:7.3f}'.format(ys2)
    plt.text(xs1-100,ys2,ss,va='bottom',ha='right',fontsize=fsz)

    # water level at powerhouse
    xs=xx[-1]; ys=zz[-1]+hh1[-1]; ss='EL.{0:7.3f}'.format(ys)
    plt.plot([xs+100,xs],[ys,ys],'-',color='#000000',lw=1)
    plt.text(xs+100,ys,ss,va='bottom',ha='left',fontsize=fsz)
    xs=xx[-1]; ys=zz[-1]+hh2[-1]; ss='EL.{0:7.3f}'.format(ys)
    plt.plot([xs+100,xs],[ys,ys],'-',color='#000000',lw=1)
    plt.text(xs+100,ys,ss,va='bottom',ha='left',fontsize=fsz)

    # water level at outlet
    if kkk==1: s1='Critical depth'
    if kkk==2: s1='1:10 flood'
    xs=xx[0]; ys=zz[0]+hh1[0]; ss='EL.{0:7.3f}'.format(ys)+'\n'+s1
    plt.plot([xs-100,xs],[ys,ys],'-',color='#000000',lw=1)
    plt.text(xs-100,ys,ss,va='top',ha='right',fontsize=fsz)
        
    # tunnel crown, walltop
    plt.plot(xx,crown,'--',color='#aaaaaa',lw=2,label='Tunnel crown (invert+6.9m)')
    plt.plot(xx,walltop,'-',color='#aaaaaa',lw=2,label='Wall top (invert+4.0m)')
    # water level
    plt.plot(xx,zz+hh2,'--',color='#0000ff',lw=2,label='water level (n1=0.020,n2=0.032)')
    plt.plot(xx,zz+hh1,'-',color='#0000ff',lw=2,label='water level (n1=0.016,n2=0.028)')
    # invert line
    plt.fill_between(xx,zz,zz-0.5,facecolor='#aaaaaa',alpha=0.5)
    plt.plot(xx,zz,color='#000000',lw=2,label='Invert')
    
    bbox_props = dict(boxstyle='square,pad=0.3', fc='#ffffff', ec='#000000', lw=1)
    plt.text(-80,ymax-0.5, sinp, ha='right', va='top', rotation=0, size=fsz, bbox=bbox_props,fontfamily='monospace')
    tstr='Tailrace Water Surface Profile (channel width: B=5.5m, design discharge: Q=68.64m$^3$/s)'
    plt.title(tstr,loc='left',fontsize=fsz)
    plt.legend(shadow=True, loc='lower left', handlelength=3)
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
    plt.show()


def cal_hc(q,b):
    g=9.8
    hc=(q**2/g/b**2)**(1/3)
    return hc


def sec(h,x,nn1,nn2):
    r0=2.9
    b0=5.5
    h0=4.0
    a1=b0*h
    s1=b0+2*h
    a2=0
    s2=0
    if 40.0<=x and h0<=h:
        a1=b0*h0
        s1=b0+2*h0
        theta=np.arcsin((h-h0)/r0)
        s2=2*r0*theta
        a2=r0**2*theta+r0**2*np.sin(theta)*np.cos(theta)
    aa=a1+a2
    rr=(a1+a2)/(s1+s2)
    n=((nn1**(3/2)*s1+nn2**(3/2)*s2)/(s1+s2))**(2/3)
    return aa,rr,n


def func(h2,h1,z1,z2,x1,x2,dx,q,nn1,nn2):
    g=9.8
    a1,r1,n1=sec(h1,x1,nn1,nn2)
    a2,r2,n2=sec(h2,x2,nn1,nn2)
    e2=q**2/2/g/a2**2+h2+z2
    e1=q**2/2/g/a1**2+h1+z1
    e3=0.5*(n1**2*q**2/r1**(4/3)/a1**2 + n2**2*q**2/r2**(4/3)/a2**2)*dx
    return e2-e1-e3


def bisection(h1,z1,z2,x1,x2,dx,q,nn1,nn2):
    ha=z1+h1-z2-0.001
    hb=6.9
    for k in range(100):
        hi=0.5*(ha+hb)
        fa=func(ha,h1,z1,z2,x1,x2,dx,q,nn1,nn2)
        fb=func(hb,h1,z1,z2,x1,x2,dx,q,nn1,nn2)
        fi=func(hi,h1,z1,z2,x1,x2,dx,q,nn1,nn2)
        if fa*fi<0: hb=hi
        if fb*fi<0: ha=hi
        #print(fa,fi,fb)
        if np.abs(hb-ha)<1e-10: break
    return hi



def main():
    #Main routine
    #Solution of non-linear equation
    ss=''
    ss=ss+'* Input coordinates & discharge               \n'
    ss=ss+'dL(m)   x(m)  Floor level  Q(m3/s)  Remarks   \n'
    ss=ss+'----------------------------------------------\n'
    ss=ss+'   0   1076   EL.891.500   22.88   PS wall(#1)\n'
    ss=ss+' +35   1041   EL.891.300   22.88              \n'
    ss=ss+' +24   1017   EL.891.300   45.76   #2 flow-in \n'
    ss=ss+' +21    996   EL.891.300   68.64   #3 flow-in \n'
    ss=ss+' +13    983   EL.891.300   68.64              \n'
    ss=ss+'+983      0   EL.888.222   68.64   Outlet     '
    sinp=ss
    g=9.8
    al=1076.0
    dx=1.0
    xx=np.arange(0,al+dx,dx)          # distance
    f1 = interpolate.interp1d([0,983,996,1017,1041,1076],[888.222,891.300,891.300,891.300,891.300,891.500])
    zz = f1(xx)
    f1 = interpolate.interp1d([0,983,996,1017,1041,1076],[68.64,68.64,68.64,45.76,22.88,22.88])
    qq = f1(xx)
    walltop=zz+4.0
    crown=zz+6.9
    ii=np.where(xx<40)
    walltop[ii]=zz[ii]+5.5
    crown[ii]=zz[ii]+5.5
    hh1=np.zeros(len(xx),dtype=np.float64) # normal water level, normal roughness
    hh2=np.zeros(len(xx),dtype=np.float64) # normal water level, maximum roughness
    hh3=np.zeros(len(xx),dtype=np.float64) # 10 years return period flood, normal roughness
    hh4=np.zeros(len(xx),dtype=np.float64) # 10 years return period flood, maximum roughness

    nnn1=[0.016,0.020] # Manning's roughness coefficient of concrete
    nnn2=[0.028,0.032] # Manning's roughness coefficient of shotcrete
    
    #print('{0:5.0f}{1:8.3f}{2:8.3f}{3:8.3f}{4:8.3f}{5:8.3f}{6:8.3f}'.format(xx[0],q,b1,z1,h1,v1,z1+h1))
    for iii in [1,2,3,4]:
        _hh=np.zeros(len(xx),dtype=np.float64)
        q=qq[0]
        if iii==1: # normal water level of outlet (normal roughness)
            h1=cal_hc(q,5.5)
            nn1=nnn1[0]
            nn2=nnn2[0]
        if iii==2: # normal water level of outlet (max roughness)
            h1=cal_hc(q,5.5)
            nn1=nnn1[1]
            nn2=nnn2[1]
        if iii==3: # 10 years return period flood (normal roughness)
            h1=893.45-888.222
            nn1=nnn1[0]
            nn2=nnn2[0]
        if iii==4: # 10 yers return period flood (max.roughness)
            h1=893.45-888.222
            nn1=nnn1[1]
            nn2=nnn2[1]
        z1=zz[0]
        _hh[0]=h1
        for i in range(1,len(xx)):
            q=qq[i]
            z1=zz[i-1]
            z2=zz[i]
            x1=xx[i-1]
            x2=xx[i]
            h2=bisection(h1,z1,z2,x1,x2,dx,q,nn1,nn2)
            _hh[i]=h2
            h1=h2
            #print('{0:5.0f}{1:8.3f}{2:8.3f}{3:8.3f}{4:8.3f}{5:8.3f}{6:8.3f}'.format(xx[i],q,b2,z2,h2,v2,z2+h2))
        if iii==1: hh1=_hh
        if iii==2: hh2=_hh
        if iii==3: hh3=_hh
        if iii==4: hh4=_hh

            
    drawfig(xx,zz,hh1,hh2,walltop,crown,sinp,1)
    drawfig(xx,zz,hh3,hh4,walltop,crown,sinp,2)


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

トンネル断面作図プログラム

# Tunnel cross section
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as tick


def barrow(x1,y1,x2,y2):
    col='#333333'
    arst='<->,head_width=3,head_length=10'
    aprop=dict(arrowstyle=arst,connectionstyle='arc3',facecolor=col,edgecolor=col,shrinkA=0,shrinkB=0,lw=1)
    plt.annotate('',
        xy=(x1,y1), xycoords='data',
        xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)


def sarrow(x1,y1,x2,y2):
    col='#333333'
    sv=0
    aprop=dict(shrink=sv,width=1,headwidth=5,headlength=8,
               connectionstyle='arc3',facecolor=col,edgecolor=col)
    plt.annotate('',
        xy=(x1,y1), xycoords='data',
        xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)

    
def drawfig(x1,y1,x2,y2,x3,y3):
    fsz=14
    xmin=-4; xmax=5; dx=1
    ymin=-6; ymax=5; dy=1
    iw=6
    ih=iw*(ymax-ymin)/(xmax-xmin)
    plt.figure(figsize=(iw,ih),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('x_distance (m)')
    plt.ylabel('y_distance (m)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.gca().spines['right'].set_visible(False)
    plt.gca().spines['top'].set_visible(False)
    plt.gca().yaxis.set_ticks_position('left')
    plt.gca().xaxis.set_ticks_position('bottom')
    plt.gca().set_aspect('equal',adjustable='box')
    
    plt.fill(x1,y1,facecolor='#aaaaaa',edgecolor='#000000',lw=1)
    plt.fill(x2,y2,facecolor='#ffffff',edgecolor='#000000',lw=1)
    plt.fill(x3,y3,color='#ffffff')
    plt.fill(x3,y3,fill=False, hatch='..',lw=2)
    
    plt.plot([-3.5,4.1],[0,0],'-.',color='#000000',lw=1)
    plt.plot([0,0],[-4.7,3.5],'-.',color='#000000',lw=1)
    
    plt.plot([0,4.1],[3.1,3.1],'-',color='#000000',lw=1)
    plt.plot([0,4.1],[2.9,2.9],'-',color='#000000',lw=1)
    plt.plot([0,4.1],[-4.0,-4.0],'-',color='#000000',lw=1)
    plt.plot([0,4.1],[-4.5,-4.5],'-',color='#000000',lw=1)
    x1=4.0; y1=-4.5; x2=x1; y2=y1-0.5; sarrow(x1,y1,x2,y2)
    x1=4.0; y1=-4.0; x2=x1; y2=0     ; barrow(x1,y1,x2,y2)
    x1=4.0; y1=0   ; x2=x1; y2=2.9   ; barrow(x1,y1,x2,y2)
    x1=4.0; y1=3.1 ; x2=x1; y2=y1+0.5; sarrow(x1,y1,x2,y2)
    plt.text(4.1,-4.0-0.25,'500',ha='left',va='center',fontsize=fsz,rotation=90)
    plt.text(4.1,-4.0/2,'4000',ha='left',va='center',fontsize=fsz,rotation=90)
    plt.text(4.1,2.9/2,'2900',ha='left',va='center',fontsize=fsz,rotation=90)
    plt.text(4.1,3.0,'200',ha='left',va='center',fontsize=fsz,rotation=90)

    plt.plot([-3.05,-3.05],[-4.5,-5.1],'-',color='#000000',lw=1)
    plt.plot([-2.75,-2.75],[-4.5,-5.1],'-',color='#000000',lw=1)
    plt.plot([ 2.75, 2.75],[-4.5,-5.1],'-',color='#000000',lw=1)
    plt.plot([ 3.05, 3.05],[-4.5,-5.1],'-',color='#000000',lw=1)
    x1=-3.05; y1=-5.0; x2=x1-0.5; y2=y1; sarrow(x1,y1,x2,y2)
    x1=-2.75; y1=-5.0; x2=2.75  ; y2=y1; barrow(x1,y1,x2,y2)
    x1= 3.05; y1=-5.0; x2=x1+0.5; y2=y1; sarrow(x1,y1,x2,y2)
    plt.text(-3.05+0.3/2,-5.2,'300',ha='center',va='top',fontsize=fsz,rotation=0)
    plt.text( 3.05-0.3/2,-5.2,'300',ha='center',va='top',fontsize=fsz,rotation=0)
    plt.text( 0,-5.2,'5500',ha='center',va='top',fontsize=fsz,rotation=0)

    plt.plot([-3.1,-3.1],[0,4.1],'-',color='#000000',lw=1)
    plt.plot([ 3.1, 3.1],[0,4.1],'-',color='#000000',lw=1)
    x1=-3.1; y1=4.0; x2=3.1  ; y2=y1; barrow(x1,y1,x2,y2)
    plt.text(0,4.0,'6200',ha='center',va='bottom',fontsize=fsz,rotation=0)
    
    bbox_args = dict(boxstyle="round", fc='#eeeeee')
    arrow_args = dict(arrowstyle="->")
    xs=-2.9*np.cos(np.radians(45))
    ys= 2.9*np.sin(np.radians(45))
    plt.annotate('Shotcrete', xy=(xs,ys), xycoords='data',
             xytext=(xs+0.5, ys-0.5), textcoords='data',
             ha='left', va='top',
             bbox=bbox_args,
             arrowprops=arrow_args)
    xs=-2.75
    ys=-2.00
    plt.annotate('Concrete', xy=(xs,ys), xycoords='data',
             xytext=(xs+0.5, ys-0.5), textcoords='data',
             ha='left', va='top',
             bbox=bbox_args,
             arrowprops=arrow_args)
    
    plt.tight_layout()
    fnameF='fig_section.png'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()

    
def main():
    r=3.1
    b=2*r
    h=4.5
    x=[]
    y=[]
    x=x+[b/2]; y=y+[-h]
    for i in range(0,181):
        theta=np.radians(float(i))
        x=x+[r*np.cos(theta)] 
        y=y+[r*np.sin(theta)]
    x=x+[-b/2]; y=y+[-h]
    x1=np.array(x)
    y1=np.array(y)
        
    r=2.9
    b=2*r
    h=4.5
    x=[]
    y=[]
    x=x+[b/2]; y=y+[-h]
    for i in range(0,181):
        theta=np.radians(float(i))
        x=x+[r*np.cos(theta)] 
        y=y+[r*np.sin(theta)]
    x=x+[-b/2]; y=y+[-h]
    x2=np.array(x)
    y2=np.array(y)

    x=[]
    y=[]
    x=x+[3.05]; y=y+[-4.5]
    x=x+[3.05]; y=y+[0]
    x=x+[2.75]; y=y+[0]
    x=x+[2.75]; y=y+[-4.0]
    x=x+[-2.75]; y=y+[-4.0]
    x=x+[-2.75]; y=y+[0]
    x=x+[-3.05]; y=y+[0]
    x=x+[-3.05]; y=y+[-4.5]
    x3=np.array(x)
    y3=np.array(y)

    drawfig(x1,y1,x2,y2,x3,y3)


#---------------
# Execute
#---------------
if __name__ == '__main__': main()

断面水理特性作図プログラム

import numpy as np
import matplotlib.pyplot as plt


def drawfig(hh,aa,rr,n1,n2):
    fsz=12
    ymin=0; ymax=8; dy=1
    fig=plt.figure(figsize=(12,6),facecolor='w')
    plt.rcParams["font.size"] = fsz
    tc=6.9
    wt=4.0

    plt.subplot(131)
    xmin=0; xmax=40; dx=10
    tstr='Flow area'
    plt.xlabel('Flow area $A$ (m$^2$)')
    plt.ylabel('Water depth $h$ (m)')
    plt.xlim(xmin,xmax)
    plt.ylim(ymin,ymax)
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#777777',linestyle='--')
    plt.title(tstr,loc='left',fontsize=fsz)

    plt.fill([xmin,xmax,xmax,xmin],[tc,tc,tc+0.3,tc+0.3],fill=False, hatch='///',lw=0)
    plt.plot([xmin,xmax],[tc,tc],color='#000000',lw=2)
    plt.text(0.05*(xmax-xmin),tc-0.05,'Inner height=6.9m',ha='left',va='top',fontsize=fsz)
    plt.plot([xmin,xmax],[wt,wt],'--',color='#000000',lw=2)
    plt.text(0.05*(xmax-xmin),wt,'Shotcrete\nConcrete',ha='left',va='center',fontsize=fsz)
    plt.plot(aa,hh,'-',color='#0000ff',lw=2)

    plt.subplot(132)
    xmin=0; xmax=2; dx=0.5
    tstr='Hydraulic radius'
    plt.xlabel('Hydraulic radius $R$ (m)')
    plt.ylabel('Water depth $h$ (m)')
    plt.xlim(xmin,xmax)
    plt.ylim(ymin,ymax)
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#777777',linestyle='--')
    plt.title(tstr,loc='left',fontsize=fsz)

    plt.fill([xmin,xmax,xmax,xmin],[tc,tc,tc+0.3,tc+0.3],fill=False, hatch='///',lw=0)
    plt.plot([xmin,xmax],[tc,tc],color='#000000',lw=2)
    plt.text(0.05*(xmax-xmin),tc-0.05,'Inner height=6.9m',ha='left',va='top',fontsize=fsz)
    plt.plot([xmin,xmax],[wt,wt],'--',color='#000000',lw=2)
    plt.text(0.05*(xmax-xmin),wt,'Shotcrete\nConcrete',ha='left',va='center',fontsize=fsz)
    plt.plot(rr,hh,'-',color='#0000ff',lw=2)

    plt.subplot(133)
    xmin=0; xmax=0.03; dx=0.01
    tstr='Manning\'s roughness coefficient'
    plt.xlabel('Combined roughness coefficient $n$')
    plt.ylabel('Water depth $h$ (m)')
    plt.xlim(xmin,xmax)
    plt.ylim(ymin,ymax)
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#777777',linestyle='--')
    plt.title(tstr,loc='left',fontsize=fsz)

    plt.text(n1[0],0.5,'Average roughness',ha='left',va='bottom',fontsize=fsz,rotation=90)
    plt.text(n2[0],0.5,'Maximum roughness',ha='left',va='bottom',fontsize=fsz,rotation=90)
    plt.fill([xmin,xmax,xmax,xmin],[tc,tc,tc+0.3,tc+0.3],fill=False, hatch='///',lw=0)
    plt.plot([xmin,xmax],[tc,tc],color='#000000',lw=2)
    plt.text(0.05*(xmax-xmin),tc-0.05,'Inner height=6.9m',ha='left',va='top',fontsize=fsz)
    plt.plot([xmin,xmax],[wt,wt],'--',color='#000000',lw=2)
    s1=''
    s1=s1+'$n_{ave}=0.028$\n'
    s1=s1+'$n_{max}=0.032$\n'
    s1=s1+'Shotcrete'
    s2=''
    s2=s2+'Concrete\n'
    s2=s2+'$n_{ave}=0.016$\n'
    s2=s2+'$n_{max}=0.020$'
    plt.text(0.05*(xmax-xmin),wt,s1+'\n'+s2,ha='left',va='center',fontsize=fsz)
    plt.plot(n1,hh,'-',color='#0000ff',lw=2,label='Average roughness')
    plt.plot(n2,hh,'--',color='#0000ff',lw=2,label='Maximum roughness')

    fnameF='fig_sec_property.png'
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
    plt.show()




def sec(h,x,nn1,nn2):
    r0=2.9
    b0=5.5
    h0=4.0
    a1=b0*h
    s1=b0+2*h
    a2=0
    s2=0
    if 40.0<=x and h0<=h:
        a1=b0*h0
        s1=b0+2*h0
        theta=np.arcsin((h-h0)/r0)
        s2=2*r0*theta
        a2=r0**2*theta+r0**2*np.sin(theta)*np.cos(theta)
    aa=a1+a2
    rr=(a1+a2)/(s1+s2)
    n=((nn1**(3/2)*s1+nn2**(3/2)*s2)/(s1+s2))**(2/3)
    return aa,rr,n



def main():
    x=50
    hh=np.arange(0,6.9,0.1)
    aa=np.zeros(len(hh),dtype=np.float64)
    rr=np.zeros(len(hh),dtype=np.float64)
    n1=np.zeros(len(hh),dtype=np.float64)
    n2=np.zeros(len(hh),dtype=np.float64)


    for i,h in enumerate(hh):
        nn1=0.016
        nn2=0.028
        a,r,n=sec(h,x,nn1,nn2)
        #print('{0:8.3f}{1:8.3f}{2:8.3f}{3:8.4f}'.format(h,a,r,n))
        aa[i]=a
        rr[i]=r
        n1[i]=n
        nn1=0.020
        nn2=0.032
        a,r,n=sec(h,x,nn1,nn2)
        n2[i]=n

    drawfig(hh,aa,rr,n1,n2)
        

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

おまけ（二分法による等流水深計算）

二分法を自前で実装。常流を対象とするため、限界水深を計算し、そこから上で水位の検索をかける。

import numpy as np

def func(h,q,b,n,i):
    a=b*h
    r=b*h/(b+2*h)
    v=1/n*r**(2/3)*i**(1/2)
    return q-a*v


g=9.8                    # gravity acceleration
n=0.016                  # Manning's roughness coefficient
b=5.5                    # channel width
i=(891.300-888.222)/983  # gradient of invert
q=68.64                  # discharge
hc=(q**2/g/b**2)**(1/3)  # critical depth

# Bisection method
h1=hc # initial value of water depth-1
h2=5  # initial value of water depth-2
for k in range(100):
    hi=0.5*(h1+h2)
    f1=func(h1,q,b,n,i)
    f2=func(h2,q,b,n,i)
    fi=func(hi,q,b,n,i)
    if f1*fi<0: h2=hi
    if f2*fi<0: h1=hi
    print(h1,h2)
    if np.abs(h2-h1)<1e-6: break
print(k,hc,hi,i)

以　上