設計 仮排水路設計における Flood Routine の活用（3）
ここでは、「設計 仮排水路設計における Flood Routine の活用（２）」で述べた事項に関連する計算・作図プログラムのうち、以下の4つを掲載する。
- トンネル標準断面作図プログラム
- トンネル通水量計算・作図プログラム
- 貯水池水位容量曲線作図プログラム
- 洪水波形作図プログラム
トンネル標準断面作図プログラム
このような説明図を描くプログラムが意外に難しいというか面倒くさい。
断面形は、幅と高さが等しい複合円である。
この断面の場合、トンネル内幅を とすれば、上半円の半径は
、下半円の半径は
となり、インバート幅は
となる。
道路トンネルの場合、インバート幅をきっちりした数字とする場合が多く、下半円の半径が中途半端な数字となることが多いが、水路トンネルの場合のインバートの制約は、施工用重機が往来できればいい。このため、水路トンネルとしてのこの断面形状は、幾何学的にきれいで結構気に入っている。
トンネル形状は、中心を原点として、水平右の点から角度を指定して (x, y) 座標を求めて置き、一気にプロットしている。
寸法線と寸法は、両矢印に記載するものと、片矢印に記載するものがあるため、以下のように３つの関数を作成して描画している。
両矢印を描く関数
def barrow(x1,y1,x2,y2): # (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): # (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 atext(x1,y1,x2,y2,ss,dd,ii): # (x1,y1),(x2,y2)：矢印始点と終点 # 矢印始点と終点座標は文字列の傾きを考慮して決定する # ss：描画文字列 # dd：文字列中央の寸法線からの距離 # ii：１なら寸法線の上側，２なら寸法線の下側にテキストを描画 al=np.sqrt((x2-x1)**2+(y2-y1)**2) cs=(x2-x1)/al sn=(y2-y1)/al rt=np.degrees(np.arccos(cs)) if y2-y1<0: rt=-rt if ii==1: xs=0.5*(x1+x2)-dd*sn ys=0.5*(y1+y2)+dd*cs plt.text(xs,ys,ss,ha='center',va='center',rotation=rt) else: xs=0.5*(x1+x2)+dd*sn ys=0.5*(y1+y2)-dd*cs plt.text(xs,ys,ss,ha='center',va='center',rotation=rt)
プログラム全文
import numpy as np import matplotlib.pyplot as plt def barrow(x1,y1,x2,y2): # (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): # (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 atext(x1,y1,x2,y2,ss,dd,ii): # (x1,y1),(x2,y2)：矢印始点と終点 # 矢印始点と終点座標は文字列の傾きを考慮して決定する # ss：描画文字列 # dd：文字列中央の寸法線からの距離 # ii：１なら寸法線の上側，２なら寸法線の下側にテキストを描画 al=np.sqrt((x2-x1)**2+(y2-y1)**2) cs=(x2-x1)/al sn=(y2-y1)/al rt=np.degrees(np.arccos(cs)) if y2-y1<0: rt=-rt if ii==1: xs=0.5*(x1+x2)-dd*sn ys=0.5*(y1+y2)+dd*cs plt.text(xs,ys,ss,ha='center',va='center',rotation=rt) else: xs=0.5*(x1+x2)+dd*sn ys=0.5*(y1+y2)-dd*cs plt.text(xs,ys,ss,ha='center',va='center',rotation=rt) def tunnel(d0): rr=d0/2 ang=np.linspace(0,np.pi,180,endpoint=False) xx1=rr*np.cos(ang) yy1=rr*np.sin(ang) rr=d0 ang=np.linspace(np.pi,np.pi+np.radians(30),30,endpoint=True) xx2=d0/2+rr*np.cos(ang) yy2=rr*np.sin(ang) ang=np.linspace(np.radians(-30),0,30,endpoint=True) xx3=-d0/2+rr*np.cos(ang) yy3=rr*np.sin(ang) xx=np.hstack([xx1,xx2,xx3]) yy=np.hstack([yy1,yy2,yy3]) return xx,yy def main(): d0=8.0 fsz=16 xmin=-6; xmax=6 ymin=-6; ymax=6 plt.figure(figsize=(6,6),facecolor='w') plt.rcParams['font.size']=fsz plt.rcParams['font.family']='sans-serif' plt.xlim([xmin,xmax]) plt.ylim([ymin,ymax]) plt.axis('off') plt.gca().set_aspect('equal',adjustable='box') # plt.title(tstr,loc='center',fontsize=fsz) xx,yy=tunnel(d0) plt.plot(xx,yy,'-',color='#000000',lw=3) xtc=0; ytc= d0/2 xbc=0; ybc=-d0/2 xcl=-d0/2; ycl=0 xcr= d0/2; ycr=0 xbl=d0/2-d0*np.cos(np.radians(30)); ybl=-d0/2 xbr=d0*np.cos(np.radians(30))-d0/2; ybr=-d0/2 ds=1.0 plt.plot([0,0],[ybc-0.5*ds,ytc+0.5*ds],'-.',color='#000000',lw=1) plt.plot([xcl-0.5*ds,xcr+0.5*ds],[0,0],'-.',color='#000000',lw=1) ss='$(\sqrt{3}-1) D$'; dd=0.5*ds; ii=-1 plt.plot([xbl,xbl],[ybl,ybl-1.5*ds],'-',color='#000000',lw=1) plt.plot([xbr,xbr],[ybr,ybr-1.5*ds],'-',color='#000000',lw=1) x1=xbl; y1=ybl-1.0*ds; x2=xbr; y2=y1 barrow(x1,y1,x2,y2) atext(x1,y1,x2,y2,ss,dd,ii) ss='$D$'; dd=0.3*ds; ii=1 plt.plot([xcl,xcl],[ycl,ytc+1.5*ds],'-',color='#000000',lw=1) plt.plot([xcr,xcr],[ycr,ytc+1.5*ds],'-',color='#000000',lw=1) x1=xcl; y1=ytc+1.0*ds; x2=xcr; y2=y1 barrow(x1,y1,x2,y2) atext(x1,y1,x2,y2,ss,dd,ii) ss='$D$'; dd=0.3*ds; ii=1 plt.plot([xtc,xcl-1.5*ds],[ytc,ytc],'-',color='#000000',lw=1) plt.plot([xbl,xcl-1.5*ds],[ybl,ybl],'-',color='#000000',lw=1) x1=xcl-1.0*ds; y1=ybl; x2=x1; y2=ytc barrow(x1,y1,x2,y2) atext(x1,y1,x2,y2,ss,dd,ii) ss='$D \; / \; 2$'; dd=0.3*ds; ii=1 x1=d0/2*np.cos(np.radians(30)); y1=x2=d0/2*np.sin(np.radians(30)); x2=0; y2=0 sarrow(x1,y1,x2,y2) atext(x2,y2,x1,y1,ss,dd,ii) ss='$D$'; dd=0.3*ds; ii=1 x1=xbr; y1=ybr; x2=xcl; y2=ycl; sarrow(x1,y1,x2,y2) atext(0.5*(x1+x2),0.5*(y1+y2),x1,y1,ss,dd,ii) x1=xbl; y1=ybl; x2=xcr; y2=ycr; sarrow(x1,y1,x2,y2) atext(x1,y1,0.5*(x1+x2),0.5*(y1+y2),ss,dd,ii) plt.tight_layout() fnameF='fig_model.jpg' plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1) plt.show() #--------------- # Execute #--------------- if __name__ == '__main__': main()
トンネル通水量計算・作図プログラム
等流水深部（自由水面）での計算は、通水断面積・潤辺は厳密な計算をせず、通水部を小さな台形に近似して計算している。
import numpy as np import matplotlib.pyplot as plt # global variables dh=0.1 # increment of reservoir water level g=9.8 # gravity acceleration n=0.015 # Manning's roughness coefficient elvs=64.0 # invert elevation of entrance elve=62.5 # invert elevation of exit def drawfig(): fsz=16 xmin=0 ; xmax=5000; dx=500 ymin=60; ymax=110; dy=10 tstr='Discharge Capacity of Diversion Tunnel' plt.figure(figsize=(10,6),facecolor='w') plt.rcParams['font.size']=fsz plt.rcParams['font.family']='sans-serif' plt.xlim([xmin,xmax]) plt.ylim([ymin,ymax]) plt.xlabel('Discharge Capacity (m$^3$/s)') plt.ylabel('Reservoir water Level (EL.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) a_d0=np.array([5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0]) l_col=[ '#dddddd', '#cccccc', '#bbbbbb', '#aaaaaa', '#999999', '#666666', '#333333', '#000000', ] for d0,col in zip(a_d0,l_col): fnameR='inp_sec_{0:0>2d}.txt'.format(int(d0)) data = np.loadtxt(fnameR, usecols=(0,1) , skiprows=0) el=data[:,0] qq=data[:,1] ss='D={0:.0f}m'.format(d0) plt.plot(qq,el,'-',color=col,lw=2,label=ss) plt.plot([xmin,xmax],[elvs,elvs],'-',color='#0000ff',lw=1) xs=10; ys=elvs; ss='Tunnel Entrance Invert EL.{0:.1f}m'.format(elvs) plt.text(xs,ys,ss,color='#0000ff',ha='left',va='top',fontsize=fsz) plt.legend(fontsize=fsz-2) fnameF='fig_sec.jpg' plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.2) plt.show() def loss_calc(ell,rho1,theta1,rho2,theta2,d0): a0=-0.5*d0 r1=0.5*d0 hd=r1 nn=int((hd+r1)/dh)+1 sel=[] sqq=[] saa=[] w1=d0*(np.sqrt(3)-1) i=(elvs-elve)/ell el=elvs qq=0.0 aa=0.0 ss=w1 rr=aa/ss sel=sel+[el] sqq=sqq+[qq] saa=saa+[aa] # free flow for j in range(1,nn): h=dh*float(j) y=h-hd if y<0.0: r=d0; a=a0 if 0.0<=y and y<=r1: r=r1; a=0.0 x=a+np.sqrt(r**2-y**2) w2=2*x da=0.5*(w1+w2)*dh ds=2*np.sqrt(dh**2+np.abs(0.5*(w1-w2))**2) aa=aa+da ss=ss+ds rr=aa/ss vv=1/n*rr**(2/3)*np.sqrt(i) qq=aa*vv el=elvs+h w1=w2 sel=sel+[el] sqq=sqq+[qq] saa=saa+[aa] # pressure flow fe=0.25 fo=1.0 fb1=(0.131+0.1632*(d0/rho1))*np.sqrt(theta1/90.0) fb2=(0.131+0.1632*(d0/rho2))*np.sqrt(theta2/90.0) ff=2*g*n*n/rr**(1/3) while 1: if 110.0<el: break h=h+dh el=elvs+h head=el-elve-hd vv=np.sqrt(2*g*head/(fe+fo+ff*ell/rr+fb1+fb2)) qq=aa*vv sel=sel+[el] sqq=sqq+[qq] saa=saa+[aa] a_el=np.array(sel) a_qq=np.array(sqq) a_aa=np.array(saa) return a_el,a_qq,a_aa def main(): for d0 in [5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0]: fnameW='inp_sec_{0:0>2d}.txt'.format(int(d0)) # inside tunnel ell=400.0 rho1=230; theta1=39.5 rho2=200; theta2=48.2 sel1,sqq1,saa1=loss_calc(ell,rho1,theta1,rho2,theta2,d0) # outside tunnel ell=450.0 rho1=260; theta1=39.5 rho2=230; theta2=48.2 sel2,sqq2,saa2=loss_calc(ell,rho1,theta1,rho2,theta2,d0) # total flow capacity sel=0.5*(sel1+sel2) sqq=sqq1+sqq2 saa=saa1+saa2 fout=open(fnameW,'w') for j in range(0,len(sel)): print('{0:10.3f}{1:10.3f}{2:10.3f}'.format(sel[j],sqq[j],saa[j]),file=fout) fout.close() drawfig() #============== # Execution #============== if __name__ == '__main__': main()
貯水池水位容量曲線作図プログラム
水位と容量の関係は、オリジナルデータでは標高10ｍピッチで与えられているため、３次スプラインで内挿し、標高1mピッチのデータとして出力・作図している。
import numpy as np from scipy.interpolate import interp1d import matplotlib.pyplot as plt def drawfig(el1,vv1,el2,vv2): fsz=16 xmin=0 ; xmax=500; dx=50 ymin=60; ymax=110; dy=5 tstr='Storage Capacity Curve' plt.figure(figsize=(10,6),facecolor='w') plt.rcParams['font.size']=fsz plt.rcParams['font.family']='sans-serif' plt.xlim([xmin,xmax]) plt.ylim([ymin,ymax]) plt.xlabel('Storage ($\\times 10^6 \; m^3$)') plt.ylabel('Reservoir water Level (EL.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) plt.plot(vv1,el1,'-',color='#0000ff',lw=2,clip_on=False) plt.plot(vv2,el2,'o',color='#0000ff',clip_on=False) fnameF='fig_hv.jpg' plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.2) plt.show() def main(): # vv x 1e6 (m3) data=np.array([[60, 0], [70, 2.485], [80, 21.670], [90, 88.555], [100, 223.959], [110, 446.817], [120, 778.284], [130, 1238.519], [140, 1865.721], [150, 2375.546], [160, 3416.319], [170, 4773.768]]) _el=data[:,0] _vv=data[:,1] f = interp1d(_el, _vv, kind='cubic') xmin=60 xmax=170 dx =1 ndx=int((xmax-xmin)/dx)+1 el = np.linspace(xmin,xmax,ndx) vv =f(el) fnameW='inp_hv.txt' fout=open(fnameW,'w') for i in range(0,ndx): print('{0:5.1f} {1:12.3f}'.format(el[i],vv[i]),file=fout) fout.close() ii1=np.where(el<=110) ii2=np.where(_el<=110) el1=el[ii1]; vv1=vv[ii1] el2=_el[ii2]; vv2=_vv[ii2] drawfig(el1,vv1,el2,vv2) #============== # Execution #============== if __name__ == '__main__': main()
洪水波形作図プログラム
1時間ピッチで与えられた流量をそのままプロットしている。 Flood Routine における水位追跡の時間ピッチは、ここで指定されている洪水波形の時間ピッチとしている。
import numpy as np import matplotlib.pyplot as plt def drawfig(tt,qq): fsz=16 xmin=0 ; xmax=132; dx=12 ymin=0; ymax=3500; dy=500 tstr='Hydrograph (20 years return period flood)' plt.figure(figsize=(10,6),facecolor='w') plt.rcParams['font.size']=fsz plt.rcParams['font.family']='sans-serif' plt.xlim([xmin,xmax]) plt.ylim([ymin,ymax]) plt.xlabel('Time (hours)') plt.ylabel('Discharge (m$^3$/s)') 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) plt.plot(tt,qq,'-',color='#0000ff',lw=2) i=np.argmax(qq) plt.text(tt[i],qq[i],'Qmax={0:.0f}'.format(qq[i]),va='bottom',ha='center',color='#000000',fontsize=fsz) fnameF='fig_hydro.jpg' plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.2) plt.show() def main(): fnameR='inp_hydro.txt' data = np.loadtxt(fnameR, usecols=(0,1) , skiprows=0) tt=data[:,0] qq=data[:,1] drawfig(tt,qq) #============== # Execution #============== if __name__ == '__main__': main()
