2020-03-30

流況曲線と運転パターン図

はじめに

同じような図は何回か作成しているのだが、また必要になったため、新しく作り直した。最新版のプログラムをアップしておく。

成果図

f:id:damyarou:20200330084330p:plain

プログラム

import numpy as np
import pandas as pd
import datetime
import matplotlib.pyplot as plt


def rdata1():
    fnameR='df_tank_nk.csv' # input file name
    df=pd.read_csv(fnameR, header=0, index_col=0) # read excel data
    df.index = pd.to_datetime(df.index, format='%Y/%m/%d')
    df=df['1985/01/01':'2017/12/31']
    return df


def fig_duration(ppp, qqq):
    qe=5.84
    qt0=7.0
    fsz=14
    xmin=0; xmax=100; dx=10
    ymin=0; ymax=40; dy=5
    plt.figure(figsize=(8,5),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    tstr='Typical operation pattern of units'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Non-exceedance Probability (%)')
    plt.ylabel('River 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)
    qt1=np.zeros(len(qqq),dtype=np.float64)
    qt2=np.zeros(len(qqq),dtype=np.float64)
    qt3=np.zeros(len(qqq),dtype=np.float64)
    for i in range(len(qqq)):
        if qe+qt0*3<=qqq[i]: qt1[i]=qe+qt0*3
        if qe+qt0*2<=qqq[i]<qe+qt0*3: qt1[i]=qqq[i]
        if qqq[i]<qe+qt0*2: qt1[i]=0
    for i in range(len(qqq)):
        if qe+qt0*3<=qqq[i]: qt2[i]=qe+qt0*2
        if qe+qt0*2<=qqq[i]<qt0*3+qe: qt2[i]=qe+qt0+(qqq[i]-qt0-qe)/2.0
        if qe+qt0*1<=qqq[i]<qt0*2+qe: qt2[i]=qqq[i]
        if qqq[i]<qe+qt0*1: qt2[i]=0            
    for i in range(len(qqq)):
        if qe+qt0*3<=qqq[i]: qt3[i]=qt0+qe
        if qe+qt0*2<=qqq[i]<qe+qt0*3: qt3[i]=qe+qt0
        if qe+qt0*1<=qqq[i]<qe+qt0*2: qt3[i]=qe+(qqq[i]-qe)/2
        if qe+0.2*qt0<=qqq[i]<qe+qt0*1: qt3[i]=qqq[i]            
    
    plt.fill_between(ppp,qt1,0,facecolor='#00ffff',alpha=1.0)
    plt.fill_between(ppp,qt2,0,facecolor='#ffffcc',alpha=1.0)
    plt.fill_between(ppp,qt3,0,facecolor='#ffd700',alpha=1.0)
    plt.fill_between(ppp,qqq,0,where=qqq<qe,facecolor='#ff00ff',alpha=1.0)
    plt.fill_between(ppp,qe,0,where=qe<qqq,facecolor='#ff00ff',alpha=1.0)
    plt.plot(ppp,qqq,'-', lw=2,color='#000080')
    plt.plot([xmin,xmax],[qe+qt0*3,qe+qt0*3],'-',color='#000000',lw=1)
    plt.plot([xmin,xmax],[qe+qt0*2,qe+qt0*2],'-',color='#000000',lw=1)
    plt.plot([xmin,xmax],[qe+qt0*1,qe+qt0*1],'-',color='#000000',lw=1)
    plt.plot([xmin,xmax],[qe+qt0*0,qe+qt0*0],'-',color='#000000',lw=1)
    ss1='Qe + 3 x Qu'
    ss2='Qe + 2 x Qu'
    ss3='Qe + 1 x Qu'
    ss4='Qe'
    plt.text(xmax,qe+qt0*3,ss1,va='bottom',ha='right',fontsize=fsz)
    plt.text(xmax,qe+qt0*2,ss2,va='bottom',ha='right',fontsize=fsz)
    plt.text(xmax,qe+qt0*1,ss3,va='bottom',ha='right',fontsize=fsz)
    plt.text(xmax,qe+qt0*0,ss4,va='bottom',ha='right',fontsize=fsz)
    plt.text(10,qe+qt0*2.5,'No.3 unit',va='center',ha='left',fontsize=fsz)
    plt.text(10,qe+qt0*1.5,'No.2 unit',va='center',ha='left',fontsize=fsz)
    plt.text(10,qe+qt0*0.5,'No.1 unit',va='center',ha='left',fontsize=fsz)
    plt.text(10,qe*0.5,'Environmental flow',va='center',ha='left',fontsize=fsz)
    xs=xmin+0.97*(xmax-xmin)
    ys=ymin+0.95*(ymax-ymin)
    ss='Qu : turbine discharge per unit\nQe : environmental flow'
    props = dict(boxstyle='square', facecolor='#ffffff', alpha=1.0)
    plt.text(xs, ys, ss, fontsize=fsz-2,va='top',ha='right',bbox=props)
    fnameF='fig_operation.png'
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
    plt.show()


def main():
    df=rdata1()
    q0 = np.array(df['Qtank'])
    qqq=np.sort(q0)[::-1] # target discharge sort
    ppp=np.linspace(0,100,len(q0)) # non-exceedance probability
    fig_duration(ppp,qqq)
   

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

以上

2020-03-15

サージング解析プログラムの挙動解明（２）

解析条件

解析パラメータを下表に示す。変化させているパラメータは、制水口径および遮断時間でる。

項目	採用値
解析条件	負荷遮断
シャフト内径（断面積）	D=10.0 m (F=78.5 m2)
ポート内径（断面積）	Dp=2.0, 3.0, 4.0, 8.0m (Fp=3.14, 7.065, 12.56, 50.24 m2)
ポート流量係数	Cd=0.9
圧力トンネル延長	L=2553.370 m
圧力トンネル内径（断面積）	d0=8.2 m (f=52.783 m2)
圧力トンネル損失係数	c=0.179
初期流量=>遮断時流量	Q=340 m3/s => 0 m3/s
遮断時間（線形遮断）	t=0.1 〜 1000 sec
貯水池水位	EL.1100
サージタンク頂部標高	EL.1200
サージタンク底部標高	EL.1000

解析結果

下図中、Tは以下で表される自由サージの周期であり、参考のため、自由サージ周期の1/4 (T/4) を一点鎖線で記載してある。

$\begin{equation} T=2\pi \sqrt{\cfrac{F L}{f g}} \end{equation}$

最高上昇水位および制水口抵抗最大値

f:id:damyarou:20200315072327j:plain

上図から以下のことがわかる。

遮断時間を増加させていくと、最高上昇水位は増加していき、ある遮断時間を超えると、最高上昇水位は減少する。
この傾向は、制水口径が大きくなってもあるが、制水口径が大きいほど、遮断時間が小さい場合（ここでは0.1秒）に対する水位上昇量は小さくなる。
制水口の最大抵抗は、遮断時間が小さい間は遮断完了時に発生するが、制水口最大抵抗が発生する時刻は、社団時間が大きくなってくと、自由サージ周期の1/4時間周辺で頭打ちとなる。

水面変動時刻歴

制水口径Dp=2.0m	制水口径Dp=3.0m

制水口径Dp=4.0m	制水口径Dp=8.0m

上図からわかるように、社団時間が150病を超えてくると、水位が調整池水位より高いところで２山を持つようになり、社団時間が200秒となると、２山目の水位が最高水位となる。

計算結果によれば、制水口型サージタンクの場合、遮断時間が長くなるにつれて最高上昇水位が高くなる傾向にあり、ある遮断時間を超えると最高上昇水位は減少していくことが確認できた。
これが実現象を正しく表しているかについて、明言はできないが、解析結果のプロットは不自然とはいえず、あり得るものと考えたほうがよいと考える。
この現象を確認するには、実験によることが好ましいと考える。
この現象の原因については不明であるが、制水口抵抗は遮断時間が大きくなるにつれて単調に減少していくため、制水口抵抗と水面を押し上げる力のバランスによりこのような現象が出ている可能性がある。

以　上

2020-03-12

サージング解析プログラムの挙動解明

きっかけ

昨日、私のwebページをご覧になった方から、サージング解析プログラムに関するメールを頂いた。

下が、その疑惑のページのコピーである。図を見てわかるように。

サージング解析では、「制水口（ポート）径が小さい場合、遮断時間を大きくしていくと、最高上昇水位が大きくなる」とされているが、これは本当なのだろうか？という疑問である。確かに計算結果はそうなっているし、そうコメントしている。これに関してどのように解釈したら良いのか解明することになった。

f:id:damyarou:20200312055605p:plain

試し計算

現象を確認するため、以下の条件でサージング計算を行ってみた。変化させるパラメータは遮断時間であり、ポート径は比較的小さい2.5mとしシャフト径は現象が顕著に出ている10mとした。貯水池水位は区切りよくEL.1100mとし、サージタンク頂部標高EL.1150、サージタンク底部標高EL.1050（シャフト高さ100m）とした。

解析パラメータを下表に示す。

項目	採用値
解析条件	負荷遮断
シャフト内径（断面積）	D=10.0 m (F=78.5 m2)
ポート内径（断面積）	Dp=2.5m (Fp=4.906 m2)
ポート流量係数	Cd=0.9
圧力トンネル延長	L=2553.370 m
圧力トンネル内径（断面積）	d0=8.2 m (f=52.783 m2)
圧力トンネル損失係数	c=0.179
初期流量=>遮断時流量	Q=340 m3/s => 0 m3/s
遮断時間（線形遮断）	t=0.1, 1, 10, 100, 1000 sec
貯水池水位	EL.1100
サージタンク頂部標高	EL.1150
サージタンク底部標高	EL.1050

計算結果は以下の通り。

f:id:damyarou:20200312060929j:plain

f:id:damyarou:20200312060947j:plain

f:id:damyarou:20200312061005j:plain

f:id:damyarou:20200312061031j:plain

f:id:damyarou:20200312061051j:plain

上図から以下のことがわかる。

遮断時間t=0.1secからt=10secまでは遮断時間が増加しているにも関わらず、最高上昇水位が増加している
遮断時間t=100secでの最高上昇水位は遮断時間t=0.1secでのものより小さくなっている
遮断時間t=1000secでは顕著なサージング現象は発生せず、定性的には期待される挙動となっている

よって、計算結果と想定される実現象に致命的な差異はないと予想されるが、遮断時間が比較的小さい領域では、遮断時間を大きくするに連れ最高上昇水位が上昇するという計算結果の解釈が課題として残る。

以　上

2020-01-21

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)

以　上

2020-01-18

macOS Catalinaのクリーンインストールと作業環境更新

2020年1月18日（土）、iMacとMacbook proにmacOS Catalinaをクリーンインストールし、作業環境を更新した。更新内容を記録しておく。

マシン

iMac (Retina 4K, 21.5-inch, 2017)
MacBook Pro (Retina, 13-inch, Mid 2014)

不便な点

この更新を行ってこれまでできていたのにできなくなったことがある。それは、

Nortonでチェックをかけると途中で止まるので、手動で勧めてやる必要があること。
Python matplotlibでpillowがインストールしてあるにも関わらずjpg出力ができないこと。どうしてもjpgが欲しい場合は、pillowでpng=>jpg変換すれば良い。

macOS catalina

ここ：https://st-over.com/pc-environment/macos-catalina/に従ってクリーンインストール。

手動インストール

Norton
Microsoft Office 2016 for Mac
Firefox           # ブラウザ
Google Chrome     # ブラウザ
Ricty Diminished  # フォント (プログラミング用等幅)
IPAex             # フォント (TeX 用日本語)
CotEditor         # テキストエディタ
Atom              # テキストエディタ
Google Earth Pro  # バーチャル地球儀システム
Google 日本語入力 # 日本語入力システム

Atomを起動しようとすると「開けないよ」と言われたら、リンゴマークから「system preferences => Seculity&Privacy => General」を確認。Atomへのアクセスを許すかどうかの表示があったらこれを許可。（鍵マークを開いておく）

Command Line Tools

'gcc -v'すると「command line toolsを入れてね」というようなメッセージが出てcommand line toolsインストール用ダイアログボックスが立ち上がるので、これに従ってインストール。もしくはこれ。

 xcode-select --install

Homebrew

ここ：https://brew.sh/index_jaから以下のインストール用スクリプトをコピーしてターミナルに貼り付けて実行。

/usr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"

Homebrewによるインストール

brew install gcc                    # gfortran含む
brew install gawk                   # GMT使用時の補助プログラミング
brew install ghostscript            # eps取扱用
brew install gmt                    # 作図：Generic Mapping Tools

Python関係

CatelinaにもPython3.7が入っているが、最新版のPythonを使いたいので、pyenvを使うことにする。

pyenvを使う

brew install pyenv
pyenv install -l
pyenv install 3.8.1
pyenv global 3.8.1

.zshrc にパスを書き込む

１行目はプロンプト表示を簡略化するおまじないです。

PROMPT="%# "

export PYENV_ROOT=${HOME}/.pyenv
export PATH=${PYENV_ROOT}/bin:$PATH
eval "$(pyenv init -)"

pipで欲しいライブラリをインストール

pip3 install --upgrade setuptools
pip3 install --upgrade pip
pip3 install numpy       # 数値計算用ライブラリ
pip3 install scipy       # 数値解析用ライブラリ
pip3 install matplotlib  # グラフ作成用ライブラリ
pip3 install pillow      # 画像処理用ライブラリ
pip3 install pandas      # データ加工支援用ライブラリ
pip3 install xlrd        # エクセルデータ読込用ライブラリ
pip3 install xlwt        # エクセルデータ書込用ライブラリ
pip3 install openpyxl    # エクセルデータ読み書き
pip3 install sympy       # 記号計算用ライブラリ
pip3 install scikit-learn #sklearn
pip3 install seaborn      #seaborn

Jupyter Notebook

Jupyterもpipでインストールします。

pip3 install jupyter
pip3 install jupyterthemes

下記によりテーマとフォントサイズ，セル幅，行間を変更する

jt -l
Available Themes:
   chesterish
   grade3
   gruvboxd
   gruvboxl
   monokai
   oceans16
   onedork
   solarizedd
   solarizedl

jt -t oceans16 -fs 12 -ofs 12 -cellw 1200 -lineh 120 -N -T

これ：https://qiita.com/pollenjp/items/88600df83448a8ff5ea6に従って行番号をデフォルト表示にする。

BasicTex

ここ：https://qiita.com/yaplus/items/55fa6957c1b15dbcf387に従ってインストール。パス「/Library/TeX/texbin」は自動的に追加されるようです。「mktexlsr」も自動で実行されるようです。

以　上

2020-01-18

Python 画像処理関係

画像変換

from PIL import Image
import os

files = os.listdir()

for file in files:
    base,ext=os.path.splitext(file)
    if ext=='.png':
        input_im = Image.open(base + ".png")
        rgb_im = input_im.convert('RGB')
        rgb_im.save(base + ".jpg",quality=30)
        print("transaction finished" + base)
        os.remove(base+'.png')

余白削除

from PIL import Image, ImageChops
import glob, os

def trim(im, border):
  bg = Image.new(im.mode, im.size, border)
  diff = ImageChops.difference(im, bg)
  bbox = diff.getbbox()
  if bbox:
    return im.crop(bbox)


def main():
    lfig=[os.path.basename(r) for r in glob.glob('*.jpg')]
    for fig in lfig:
        img_org=Image.open(fig,'r')
        img_new=trim(img_org,'#ffffff')
        img_new.save(fig, 'JPEG', quality=100, optimize=True)
        img_new.show()


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

画像結合

# 画像結合
# https://note.nkmk.me/python-pillow-concat-images/

from PIL import Image
import os
import glob


def comb_h(im1, im2):
    dst = Image.new('RGB', (im1.width + im2.width, im1.height))
    dst.paste(im1, (0, 0))
    dst.paste(im2, (im1.width, 0))
    return dst

def comb_v(im1, im2):
    dst = Image.new('RGB', (im1.width, im1.height + im2.height))
    dst.paste(im1, (0, 0))
    dst.paste(im2, (0, im1.height))
    return dst

def multi_comb_h(im_list):
    _im=im_list.pop(0)
    for im in im_list:
        _im=comb_h(_im, im)
    return _im

def multi_comb_v(im_list):
    _im=im_list.pop(0)
    for im in im_list:
        _im=comb_v(_im, im)
    return _im



fig_01='fig_cont1.jpg'; im_01 = Image.open(fig_01)
fig_02='fig_cont2.jpg'; im_02 = Image.open(fig_02)
fig_03='fig_cont3.jpg'; im_03 = Image.open(fig_03)
fig_04='fig_vect1.jpg'; im_04 = Image.open(fig_04)
fig_05='fig_vect2.jpg'; im_05 = Image.open(fig_05)
fig_06='fig_disp0.jpg'; im_06 = Image.open(fig_06)

fnameF='fig_ax4.jpg'
im_list_1=[im_01, im_02, im_03]
im_list_2=[im_04, im_05, im_06]
_im1=multi_comb_v(im_list_1)
_im2=multi_comb_v(im_list_2)
multi_comb_h([_im1, _im2]).save(fnameF)

file_list = glob.glob('_*.jpg')
for file in file_list:
    print("remove：{0}".format(file))
    os.remove(file)

以　上

2020-01-13

Python RC部材の設計

RC部材の設計をしたので流れを記録しておく。

リンク

平面骨組解析プログラムはこちら（Qiita）。
配筋計算プログラムの自作関数はこちら（はてな）。

出力

モデル図

f:id:damyarou:20200113182801j:plain

断面力図

f:id:damyarou:20200113182824j:plain

配鉄設計

f:id:damyarou:20200125102959p:plain:w300

プログラム

モデル図

import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate


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 drawfig():
    tstr='Model of Tailrace Outlet Channel'
    fsz=12
    xmin=-11
    xmax=11
    ymin=-8
    ymax=2
    iw=12
    ih=(ymax-ymin)/(xmax-xmin)*iw
    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_direction (m)')
    plt.ylabel('y_direction (m)')
    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.title(tstr,loc='left',fontsize=fsz)
    
    x=np.array([-3.75,3.75,3.75,2.75,2.75,-2.75,-2.75,-3.75,-3.75])
    y=np.array([-5.0,-5.0,0.0,0.0,-4.0,-4.0,0.0,0.0,-5.0]) 
    plt.fill(x,y,color='#eeeeee')
    plt.plot(x,y,'-',color='#000000',lw=1) 
    x1=-3.75; y1=-1.0; x2=-2.75; y2=y1; barrow(x1,y1,x2,y2)
    xs=0.5*(x1+x2); ys=y1; plt.text(xs,ys,'{0:.3f}'.format(x2-x1),va='bottom',ha='center',fontsize=fsz)
    x1=-2.75; y1=-1.0; x2= 2.75; y2=y1; barrow(x1,y1,x2,y2)
    xs=0.5*(x1+x2); ys=y1; plt.text(xs,ys,'{0:.3f}'.format(x2-x1),va='bottom',ha='center',fontsize=fsz)
    x1= 2.75; y1=-1.0; x2= 3.75; y2=y1; barrow(x1,y1,x2,y2)
    xs=0.5*(x1+x2); ys=y1; plt.text(xs,ys,'{0:.3f}'.format(x2-x1),va='bottom',ha='center',fontsize=fsz)
    x1=-1.5; y1=-5.0; x2=x1; y2=-4.0; barrow(x1,y1,x2,y2)
    xs=x1; ys=0.5*(y1+y2); plt.text(xs,ys,'{0:.3f}'.format(y2-y1),va='center',ha='right',rotation=90,fontsize=fsz)
    x1=-1.5; y1=-4.0; x2=x1; y2=0.0; barrow(x1,y1,x2,y2)
    xs=x1; ys=0.5*(y1+y2); plt.text(xs,ys,'{0:.3f}'.format(y2-y1),va='center',ha='right',rotation=90,fontsize=fsz)
    x1=-4.5; y1=-1.3; x2=x1; y2=0.0; barrow(x1,y1,x2,y2)
    xs=x1; ys=0.5*(y1+y2); plt.text(xs,ys,'{0:.3f}'.format(y2-y1),va='center',ha='right',rotation=90,fontsize=fsz)
    
    # lateral water pressure
    x1=-3.75; p1=3.7*0.5
    x=np.array([x1,x1,x1-p1,x1])
    y=np.array([-5.0,-1.3,-5.0,-5.0]) 
    plt.fill(x,y,color='#00ffff',alpha=0.4)
    plt.fill(-x,y,color='#00ffff',alpha=0.4)
    yy=np.array([-5.0,-4.5,-4.0,-3.5,-3.0,-2.5,-2.0])
    f1 = interpolate.interp1d([-5.0,-1.3], [3.7,0])
    uu = f1(yy)*0.5
    xx=x1-uu
    vv=np.zeros(len(yy),dtype=np.float64)
    plt.quiver(xx,yy,uu,vv,angles='xy',scale_units='xy',scale=1,width=0.003)
    plt.quiver(-xx,yy,-uu,vv,angles='xy',scale_units='xy',scale=1,width=0.003)
    xs=-3.75; ys=-5.0; ss='{0:.1f}kN/m$^2$'.format(37)
    plt.text(xs,ys,ss,va='top',ha='right',fontsize=fsz)
    xs=3.75; ys=-5.0; ss='{0:.1f}kN/m$^2$'.format(37)
    plt.text(xs,ys,ss,va='top',ha='left',fontsize=fsz)
    xs=-3.75; ys=0; ss='Lateral\nwater\npressure'
    plt.text(xs,ys,ss,va='bottom',ha='right',fontsize=fsz)
    xs=-x1; ys=0
    plt.text(xs,ys,ss,va='bottom',ha='left',fontsize=fsz)
    
    # uplift
    x=np.array([-3.75,3.75,3.75,-3.75,-3.75])
    y=np.array([-5.0-p1,-5.0-p1,-5.0,-5.0,-5.0-p1]) 
    plt.fill(x,y,color='#00ffff',alpha=0.4)
    xx=np.linspace(-3.75,3.75,16)
    yy=np.zeros(len(xx),dtype=np.float64)-5.0-3.7/2
    uu=np.zeros(len(xx),dtype=np.float64)
    vv=np.zeros(len(xx),dtype=np.float64)+3.7/2
    plt.quiver(xx,yy,uu,vv,angles='xy',scale_units='xy',scale=1,width=0.003)
    xs=0; ys=-5.0-p1; ss='Uplift {0:.1f}kN/m$^2$'.format(37)
    plt.text(xs,ys,ss,va='top',ha='center',fontsize=fsz)

    # lateral earth pressure
    k0=0.5; qq=1.0
    db=5.0; pb=k0*(1.3*db+qq)*0.5
    di=1.3; pi=k0*(2.3*di+qq)*0.5
    dt=0.0; pt=k0*(2.3*dt+qq)*0.5
    x1=-6.0
    x=np.array([x1,x1,x1-pt,x1-pi,x1-pb,x1])
    y=np.array([-5.0,0.0,0.0,-di,-5.0,-5.0])    
    plt.fill(x,y,color='#00ff00',alpha=0.4)
    plt.fill(-x,y,color='#00ff00',alpha=0.4)
    yy=np.array([-5.0,-4.5,-4.0,-3.5,-3.0,-2.5,-2.0,-1.5,-1.0,-0.5])
    f1 = interpolate.interp1d([-5.0,-1.3,0], [pb,pi,pt])
    uu = f1(yy)
    xx=x1-uu
    vv=np.zeros(len(yy),dtype=np.float64)
    plt.quiver(xx,yy,uu,vv,angles='xy',scale_units='xy',scale=1,width=0.003)
    plt.quiver(-xx,yy,-uu,vv,angles='xy',scale_units='xy',scale=1,width=0.003)
    xs=x1-pt; ys=dt; ss='{0:.1f}kN/m$^2$'.format(pt*2*10)
    plt.text(xs,ys,ss,va='center',ha='right',fontsize=fsz)
    xs=x1-pi; ys=-di; ss='{0:.1f}kN/m$^2$'.format(pi*2*10)
    plt.text(xs,ys,ss,va='center',ha='right',fontsize=fsz)
    xs=x1-pb; ys=-db; ss='{0:.1f}kN/m$^2$'.format(pb*2*10)
    plt.text(xs,ys,ss,va='center',ha='right',fontsize=fsz)
    xs=-(x1-pt); ys=dt; ss='{0:.1f}kN/m$^2$'.format(pt*2*10)
    plt.text(xs,ys,ss,va='center',ha='left',fontsize=fsz)
    xs=-(x1-pi); ys=-di; ss='{0:.1f}kN/m$^2$'.format(pi*2*10)
    plt.text(xs,ys,ss,va='center',ha='left',fontsize=fsz)
    xs=-(x1-pb); ys=-db; ss='{0:.1f}kN/m$^2$'.format(pb*2*10)
    plt.text(xs,ys,ss,va='center',ha='left',fontsize=fsz)
    xs=x1; ys=-5.5; ss='Lateral\nearth\npressure'
    plt.text(xs,ys,ss,va='top',ha='right',fontsize=fsz)
    xs=-x1; ys=-5.5
    plt.text(xs,ys,ss,va='top',ha='left',fontsize=fsz)
    
    
    
    plt.tight_layout()
    fnameF='fig_model.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()



def main():
    drawfig()
    
        
if __name__ == '__main__':
    main()

断面力図

import matplotlib.pyplot as plt
import numpy as np
import sys



def drawfig(disx,ax1,sh1,mo1,disy,ax2,sh2,mo2,x,y,node,nele):
    nmax=np.max([np.max(np.abs(ax1)),np.max(np.abs(ax2))])
    smax=np.max([np.max(np.abs(sh1)),np.max(np.abs(sh2))])
    mmax=np.max([np.max(np.abs(mo1)),np.max(np.abs(mo2))])
    dmax=np.max([np.max(np.abs(disx)),np.max(np.abs(disy))])
    fsz=10
    xmin=-7
    xmax=7
    ymin=-7
    ymax=4
    scl_dis=1.0
    scl_axi=1.0
    scl_she=1.0
    scl_mom=2.0
    iw=10
    ih=(ymax-ymin)/(xmax-xmin)*iw
    plt.figure(figsize=(iw,ih),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'

    for nnn in [1,2,3,4]:
        nplot=220+nnn
        plt.subplot(nplot)
        plt.xlim([xmin,xmax])
        plt.ylim([ymin,ymax])
        plt.xlabel('x_direction (m)')
        plt.ylabel('y_direction (m)')
        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')

        if nnn==1: # displacement
            tstr='Displacement mode'
            ls1='disx_max={0:8.3f}mm'.format(np.max(disx))
            ls2='disx_min={0:8.3f}mm'.format(np.min(disx))
            ls3='disy_max={0:8.3f}mm'.format(np.max(disy))
            ls4='disy_min={0:8.3f}mm'.format(np.min(disy))
            dx=x+disx/dmax*scl_dis
            dy=y+disy/dmax*scl_dis
            for ne in range(0,nele):
                n1=node[0,ne]-1
                n2=node[1,ne]-1
                plt.plot([x[n1],x[n2]],[y[n1],y[n2]],color='gray',linewidth=0.5)
                plt.plot([dx[n1],dx[n2]],[dy[n1],dy[n2]],color='black',linewidth=1)
        if nnn==2: # axial force diagram
            tstr='Axial force diagram'
            ls1='N_max={0:10.3f}kN'.format(np.max([np.max(ax1),np.max(ax2)]))
            ls2='N_min={0:10.3f}kN'.format(np.min([np.min(ax1),np.min(ax2)]))
            ls3=''
            ls4=''
            d1=ax1/nmax*scl_axi
            d2=ax2/nmax*scl_axi
            for ne in range(0,nele):
                x1,x2,x3,x4,y1,y2,y3,y4=calc(ne,node,x,y,d1,d2)
                if d1[ne]<=0.0: # compression
                    plt.fill([x1,x2,x3,x4],[y1,y2,y3,y4],color='black',alpha=0.1)
                else: # tension
                    plt.fill([x1,x2,x3,x4],[y1,y2,y3,y4],color='black',alpha=0.2)
            for ne in range(0,nele):
                n1=node[0,ne]-1
                n2=node[1,ne]-1
                plt.plot([x[n1],x[n2]],[y[n1],y[n2]],color='black',linewidth=0.5)
        if nnn==3: # shearing force
            tstr='Shear force diagram'
            ls1='S_max={0:10.3f}kN'.format(np.max([np.max(-sh1),np.max(-sh2)]))
            ls2='S_min={0:10.3f}kN'.format(np.min([np.min(-sh1),np.min(-sh2)]))
            ls3=''
            ls4=''
            d1=sh1/smax*scl_she
            d2=sh2/smax*scl_she
            for ne in range(0,nele):
                x1,x2,x3,x4,y1,y2,y3,y4=calc(ne,node,x,y,d1,d2)
                plt.fill([x1,x2,x3,x4],[y1,y2,y3,y4],color='black',alpha=0.1)
            for ne in range(0,nele):
                n1=node[0,ne]-1
                n2=node[1,ne]-1
                plt.plot([x[n1],x[n2]],[y[n1],y[n2]],color='black',linewidth=0.5)
        if nnn==4:
            # moment
            tstr='Bending moment diagram'
            ls1='M_max={0:10.3f}kN-m'.format(np.max([np.max(-mo1),np.max(-mo2)]))
            ls2='M_min={0:10.3f}kN-m'.format(np.min([np.min(-mo1),np.min(-mo2)]))
            ls3=''
            ls4=''
            d1=mo1/mmax*scl_mom
            d2=mo2/mmax*scl_mom
            for ne in range(0,nele):
                x1,x2,x3,x4,y1,y2,y3,y4=calc(ne,node,x,y,d1,d2)
                plt.fill([x1,x2,x3,x4],[y1,y2,y3,y4],color='black',alpha=0.1)
            for ne in range(0,nele):
                n1=node[0,ne]-1
                n2=node[1,ne]-1
                plt.plot([x[n1],x[n2]],[y[n1],y[n2]],color='black',linewidth=0.5)

        plt.plot(xmin,ymin,'.',label=ls1)
        plt.plot(xmin,ymin,'.',label=ls2)
        plt.plot(xmin,ymin,'.',label=ls3)
        plt.plot(xmin,ymin,'.',label=ls4)
        plt.legend(loc='upper right',numpoints=1,markerscale=0, frameon=False,prop={'family':'monospace','size':12})
        plt.title(tstr,loc='left',fontsize=fsz)
    
    plt.tight_layout()
    fnameF='fig_force.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()



def calc(ne,node,x,y,d1,d2):
    i=node[0,ne]-1
    j=node[1,ne]-1
    x1=x[i]
    y1=y[i]
    x2=x[j]
    y2=y[j]
    al=np.sqrt((x2-x1)**2+(y2-y1)**2)
    theta=np.arccos((x2-x1)/al)
    x4=x1-d1[ne]*np.sin(theta)
    y4=y1+d1[ne]*np.cos(theta)
    x3=x2-d2[ne]*np.sin(theta)
    y3=y2+d2[ne]*np.cos(theta)
    return x1,x2,x3,x4,y1,y2,y3,y4


def main():
    # data input
    #args = sys.argv
    #fnameR=args[1] # input data file
    fnameR='out_fem_outlet.txt'
    
    f=open(fnameR,'r')
    text=f.readline()
    text=f.readline()
    text=text.strip()
    text=text.split()
    npoin=int(text[0]) # Number of nodes
    nele =int(text[1]) # Number of elements
    nsec =int(text[2]) # Number of sections
    npfix=int(text[3]) # Number of restricted nodes
    nlod =int(text[4]) # Number of loaded nodes
    x   =np.zeros(npoin,dtype=np.float64) # Coordinates of nodes
    y   =np.zeros(npoin,dtype=np.float64) # Coordinates of nodes
    node=np.zeros([3,nele],dtype=np.int)  # Node-element relationship
    disx=np.zeros(npoin,dtype=np.float64) # Coordinates of nodes
    disy=np.zeros(npoin,dtype=np.float64) # Coordinates of nodes
    ax1  =np.zeros(nele,dtype=np.float64)  # Section force vector
    sh1  =np.zeros(nele,dtype=np.float64)  # Section force vector
    mo1  =np.zeros(nele,dtype=np.float64)  # Section force vector
    ax2  =np.zeros(nele,dtype=np.float64)  # Section force vector
    sh2  =np.zeros(nele,dtype=np.float64)  # Section force vector
    mo2  =np.zeros(nele,dtype=np.float64)  # Section force vector
    text=f.readline()
    for i in range(0,nsec):
        text=f.readline()
    text=f.readline()
    for i in range(0,npoin):
        text=f.readline()
        text=text.strip()
        text=text.split()
        x[i]=float(text[1]) # x-coordinate
        y[i]=float(text[2]) # y-coordinate
    text=f.readline()
    for i in range(0,npfix):
        text=f.readline()
    text=f.readline()
    for i in range(nele):
        text=f.readline()
        text=text.strip()
        text=text.split()
        node[0,i]=int(text[1]) #node_1
        node[1,i]=int(text[2]) #node_2
        node[2,i]=int(text[2]) #mat-no
    text=f.readline()
    for i in range(0,npoin):
        text=f.readline()
        text=text.strip()
        text=text.split()
        disx[i]=float(text[1]) # displacement in x-direction
        disy[i]=float(text[2]) # displacement in y-direction
    text=f.readline()
    for i in range(0,nele):
        text=f.readline()
        text=text.strip()
        text=text.split()
        ax1[i]=-float(text[1]) # axial force at node-1
        sh1[i]= float(text[2]) # shear force at node-1
        mo1[i]= float(text[3]) # moment at node-1
        ax2[i]= float(text[4]) # axial force at node-2
        sh2[i]=-float(text[5]) # shear force at node-2
        mo2[i]=-float(text[6]) # moment at node-2
    f.close()

    disx=disx*1000
    disy=disy*1000
    ax1=ax1*10
    sh1=sh1*10
    mo1=mo1*10
    ax2=ax2*10
    sh2=sh2*10
    mo2=mo2*10
    drawfig(disx,ax1,sh1,mo1,disy,ax2,sh2,mo2,x,y,node,nele)
    
        
if __name__ == '__main__':
    main()

配筋設計

プログラム中で用いている、自作関数 im_rebar.py については、ここ：https://damyarou.hatenablog.com/entry/2019/06/02/172405 に紹介しています。

#=============================================
# Ultimate strehgth of Rectangular RC section
# (use module: im_rebar)
#=============================================
import numpy as np
import matplotlib.pyplot as plt
import time
import sys; sys.path.append('/Users/kk/data_kk/pypro')
import im_rebar


def drawfig(fnp,fmp,fm0,nreq,mreq,strt,fnameF):
    fsz=16
    fig=plt.figure(figsize=(6,6),facecolor='w')
    plt.rcParams["font.size"] = fsz
    plt.rcParams['font.family'] ='sans-serif'
    plt.xlabel('Bending moment $M\: / \:bh^2$ (MPa)')
    plt.ylabel('Axial force $N\: / \:bh$ (MPa)')
    plt.grid(which='major',color='#777777',linestyle='solid')
    plt.plot(fmp,fnp,'-',color='#000080',label='Ultimate strength')
    plt.plot(mreq,nreq,'o',color='#ff0000',label='required')
    xs=mreq[0]; ys=nreq[0]+0.5; ss='{0:.3f}'.format(xs)
    plt.text(xs,ys,ss,va='bottom',ha='center',fontsize=fsz,rotation=90)    
    plt.legend(loc='upper right')
    plt.title(strt,loc='left',fontsize=fsz-2)
    plt.tight_layout()
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
    plt.show()


def main():
    start=time.time()
    a32=32*32*3.14/4
    a25=25*25*3.14/4
    a20=20*20*3.14/4
    fcu=35.0   # design strength of concrete
    for iii in [1]:
        if iii==1:
            bb=1000
            hh=1000
            As=np.array([a20*5,a20*5])
            ys=np.array([100,hh-100])
            strt='Outlet channel B1000xH1000,C35-T20-200(D)'
            fnameF='fig_rc_channel.png'
        mh=int(hh/10)
        nreq=np.array([0])        # kN
        mreq=np.array([257.561*1.4]) # kN-m
        nreq=nreq*1000/bb/hh
        mreq=mreq*1000*1000/bb/hh/hh
        fnp,fmp,fm0=im_rebar.calmain(mh,fcu,bb,hh,As,ys)
        drawfig(fnp,fmp,fm0,nreq,mreq,strt,fnameF)
    print(time.time()-start)
    
    
#==============
# Execution
#==============
if __name__ == '__main__': main()

以　上