設計 仮排水路設計における Flood Routine の活用（２）
ここでは、「設計 仮排水路設計における Flood Routine の活用（２）」で述べた事項に関連する計算・作図プログラムのうち、以下の２つを掲載する。
- Flood Routine 解析プログラム
- Flood Routine 作図プログラム
Flood Routine 解析プログラム
import numpy as np def flood(elini,outlet,rh,rv,nr,ti,q_in,nd,elof,qof,no,fout): p_el=np.zeros(nd-1,dtype=np.float64) p_qo=np.zeros(nd-1,dtype=np.float64) # Initial outflow elv=elini el=elv vol=ret_v(nr,rh,rv,elv) q_out=ofc(elv,elof,qof,no)+outlet i=0 iud=0 print('{0:>5s} {1:>5s} {2:>10s} {3:>10s} {4:>15s} {5:>15s} {6:>15s} {7:>15s}'.format('i','iud','time','el','el-elv','vol','q_in','q_out'),file=fout) print('{0:5d} {1:5d} {2:10.3f} {3:10.3f} {4:15e} {5:15e} {6:15e} {7:15e}'.format(i,iud,ti[i],el,el-elv,vol,q_in[i],q_out),file=fout) # Iterative calculation iud=1 dh=0.0005 eps=0.001 itmax=int(1.0/dh)*100 for i in range(0,nd-1): qqin=0.5*(q_in[i]+q_in[i+1]) tim=0.5*(ti[i+1]-ti[i]) qin=0.5*(q_in[i]+q_in[i+1])*(ti[i+1]-ti[i])*3600.0 hh=0.0 k=0 j=0 while 1: k=k+1 j=j+1; hh=hh+float(iud)*dh elv=el; q1=ofc(elv,elof,qof,no)+outlet elv=el+hh; q2=ofc(elv,elof,qof,no)+outlet qout=0.5*(q1+q2)*(ti[i+1]-ti[i])*3600.0 R=vol+qin-qout elv=ret_h(nr,rh,rv,R) if iud==1 and j==10: if el+hh>elv: iud=-1 hh=0.0 j=0 if iud==-1 and j==10: if el+hh<elv: iud=1 hh=0.0 j=0 if abs(el+hh-elv)<eps: break if itmax<k: break vol=R # Cumulative volume el=el+hh # Elevation q_out=q2 # Outflow print('{0:5d} {1:5d} {2:10.3f} {3:10.3f} {4:15e} {5:15e} {6:15e} {7:15e}'.format(i+1,iud,ti[i+1],el,el-elv,vol,q_in[i+1],q_out),file=fout) p_el[i]=el p_qo[i]=q_out hmax=np.max(p_el)-elini print('Time: {0:10.3f}'.format(np.max(ti))) print('hmax: {0:10.3f}'.format(hmax)) print('Qin : {0:10.3f} {1:10.3f}'.format(np.min(q_in),np.max(q_in))) print('Wout: {0:10.3f} {1:10.3f}'.format(np.min(p_qo),np.max(p_qo))) print('EL : {0:10.3f} {1:10.3f}'.format(np.min(p_el),np.max(p_el))) del p_el,p_qo def ofc(elv,elof,qof,no): for i in range(0,no-1): if elof[i]<=elv and elv<=elof[i+1]: break if elof[no-1]<elv: i=no-2 x1=elof[i] y1=qof[i] x2=elof[i+1] y2=qof[i+1] a=(y2-y1)/(x2-x1) b=(x2*y1-x1*y2)/(x2-x1) q=a*elv+b return q def ret_v(nr,rh,rv,elv): # To obtain the cumulative volume from the water level for i in range(0,nr-1): if rh[i]<=elv and elv<=rh[i+1]: break if rh[nr-1]<elv: i=nr-2 x1=rv[i] y1=rh[i] x2=rv[i+1] y2=rh[i+1] a=(y2-y1)/(x2-x1) b=(x2*y1-x1*y2)/(x2-x1) v=(elv-b)/a return v def ret_h(nr,rh,rv,v): # To obtain the water level from cumulative volume for i in range(0,nr-1): if rv[i]<=v and v<=rv[i+1]: break if rv[nr-1]<v: i=nr-2 x1=rv[i] y1=rh[i] x2=rv[i+1] y2=rh[i+1] a=(y2-y1)/(x2-x1) b=(x2*y1-x1*y2)/(x2-x1) elv=a*v+b return elv def inp_hv(fnameR1): # Input H-V data vcoef=1e6 data = np.loadtxt(fnameR1, usecols=(0,1) , skiprows=0) rh=data[:,0] rv=data[:,1]*vcoef nr=len(rh) return nr,rh,rv def inp_inflow(fnameR2): # Input time sequence of inflow data = np.loadtxt(fnameR2, usecols=(0,1) , skiprows=0) ti =data[:,0] q_in=data[:,1] nd=len(ti) elini=64.0 outlet=0.0 return nd,ti,q_in,elini,outlet def inp_outflow(fnameR3): # Input outflow capacity data = np.loadtxt(fnameR3, usecols=(0,1) , skiprows=0) elof=data[:,0] qof =data[:,1] no=len(elof) return no,elof,qof def main(): #Main routine fnameR1='inp_hv.txt' fnameR2='inp_hydro.txt' nr,rh,rv=inp_hv(fnameR1) nd,ti,q_in,elini,outlet=inp_inflow(fnameR2) for ss in ['05','06','07','08','09','10','11','12']: fnameR3='inp_sec_'+ss+'.txt' # Outflow capacity fnameW ='out_'+ss+'.txt' # Output file name no,elof,qof=inp_outflow(fnameR3) fout=open(fnameW,'w') flood(elini,outlet,rh,rv,nr,ti,q_in,nd,elof,qof,no,fout) fout.close() #============== # Execution #============== if __name__ == '__main__': main()
Flood Routine 作図プログラム
import numpy as np import matplotlib.pyplot as plt fsz=16 xmin =0.0 ; xmax =132.0 ; dx =12 # time ymin1=0.0 ; ymax1=5000.0; dy1=1000 # discharge ymin2=60.0; ymax2=110.0 ; dy2=10 # water level def drawfig(nnn,d0): ss='{0:0>2}'.format(d0) fnameR='out_'+ss+'.txt' # calculation result data=np.loadtxt(fnameR,skiprows=1,usecols=(2,3,6,7)) ti=data[:,0] el=data[:,1] qi=data[:,2] qo=data[:,3] n1=np.argmax(qi) n2=np.argmax(qo) n3=np.argmax(el) plt.subplot(nnn) plt.xlim([xmin,xmax]) plt.ylim([ymin1,ymax1]) plt.xlabel('Time (hour)') plt.ylabel('Discharge (m$^3$/s)') plt.xticks(np.arange(xmin,xmax+dx,dx)) plt.yticks(np.arange(ymin1,ymax1+dy1,dy1)) plt.grid(color='#999999',linestyle='solid') plt.plot(ti,qi,'--',color='#000000',lw=2.0,label='Q (inflow)') plt.plot(ti,qo,'-.',color='#000000',lw=2.0,label='Q (outflow)') plt.text(ti[n1],qi[n1],'max:{0:.0f}'.format(qi[n1]),fontsize=fsz,color='#000000',ha='center',va='bottom') plt.text(ti[n2],qo[n2],'max:{0:.0f}'.format(qo[n2]),fontsize=fsz,color='#000000',ha='center',va='bottom') plt.twinx() plt.ylim([ymin2,ymax2]) plt.ylabel('Water Level (EL.m)') plt.plot(ti,el,'-',color='#0000ff',lw=2.0,label='Water Level') plt.text(ti[n3],el[n3],'ELmax:{0:.3f}'.format(el[n3]),fontsize=fsz,color='#000000',ha='left',va='bottom') plt.text(xmin+5,ymax2-2,'D{0:.1f}m x 2'.format(d0),fontsize=fsz+4,color='#000000',ha='left',va='top') elti=64.0 eltc=elti+d0 plt.plot([xmin,xmax],[elti,elti],'-',color='#0000ff',lw=1) plt.plot([xmin,xmax],[eltc,eltc],'-',color='#0000ff',lw=1) plt.text(xmax-5,elti,'Tunnel Invert:EL.{0:.1f}'.format(elti),fontsize=fsz,color='#0000ff',ha='right',va='top') plt.text(xmax-5,eltc,'Tunnel Crown:EL.{0:.1f}'.format(eltc) ,fontsize=fsz,color='#0000ff',ha='right',va='bottom') plt.plot([0],[0],'--',color='#000000',lw=2.0,label='Q (inflow)') plt.plot([0],[0],'-.',color='#000000',lw=2.0,label='Q (outflow)') plt.legend(shadow=True,loc='upper right',handlelength=2,fontsize=fsz-2) return qo[n2],el[n3] def main(): # Main routine fnameF='fig_diversion.jpg' # image file ww=8; hh=5 plt.figure(figsize=(ww*2,hh*4),facecolor='w') plt.rcParams['font.size']=fsz plt.rcParams['font.family']='sans-serif' l_d0=[5,6,7,8,9,10,11,12] qomax=[] elmax=[] nnn=420 for d0 in l_d0: nnn=nnn+1 qo,el=drawfig(nnn,d0) qomax=qomax+[qo] elmax=elmax+[el] print(qomax) print(elmax) plt.tight_layout() plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.2) plt.show() #============== # Execution #============== if __name__ == '__main__': main()
