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()

以　上

2020-01-10

Python コンター図作成

コンター図を作成する必要があったため、そのプログラムを作成した。

作例

f:id:damyarou:20200110160608j:plain

プログラムソース

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns


def drawfig(xx,yy,zh,zc):
    q=68.64
    vp1=7.0; dp1=np.sqrt(4*q/np.pi/vp1)
    vp2=2.5; dp2=np.sqrt(4*q/np.pi/vp2)
    vh1=5.0; dh1=np.sqrt(4*q/np.pi/vh1)
    vh2=2.0; dh2=np.sqrt(4*q/np.pi/vh2)
    fsz=16
    xmin=3; xmax=6
    ymin=4; ymax=7
    plt.figure(figsize=(20,10),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'

    plt.subplot(121)
    tstr='Total Head loss Contour of Headrace and Penstock (unit: m)'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Penstock diameter (m)')
    plt.ylabel('Headrace diameter (m)')
    plt.grid(color='#999999',linestyle='dashed')
    plt.gca().set_aspect('equal',adjustable='box')
    plt.title(tstr,loc='left',fontsize=fsz)
    
    plt.fill([xmin,xmax,xmax,xmin],[ymin,ymin,ymax,ymax],color='#aaaaaa',alpha=0.3)
    plt.fill([dp1,dp2,dp2,dp1],[dh1,dh1,dh2,dh2],color='#ffffff')
    co0='#555555'
    xs=dp1; ys=0.5*(dh1+dh2); ss='$v_P$=7.0m/s'
    plt.text(xs,ys,ss,color=co0,va='center',ha='left',rotation=90,fontsize=fsz)
    xs=dp2; ys=0.5*(dh1+dh2); ss='$v_P$=2.5m/s'
    plt.text(xs,ys,ss,color=co0,va='center',ha='right',rotation=90,fontsize=fsz)
    xs=4.25; ys=dh1; ss='$v_H$=5.0m/s'
    plt.text(xs,ys,ss,color=co0,va='bottom',ha='center',rotation=0,fontsize=fsz)
    xs=4.25; ys=dh2; ss='$v_H$=2.0m/s'
    plt.text(xs,ys,ss,color=co0,va='top',ha='center',rotation=0,fontsize=fsz)
    co1='black'
    co2='red'
    cont=plt.contour(xx,yy,zh,colors=[co1,co1,co2,co1,co1,co1,co1,co1],levels=[10,15,16.58,20,30,40,60,80])
    cont.clabel(fmt='%1.2f', fontsize=fsz-1)

    plt.subplot(122)
    tstr='Total Cost Contour of Headrace and Penstock (unit: mil.NRP)'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Penstock diameter (m)')
    plt.ylabel('Headrace diameter (m)')
    plt.grid(color='#999999',linestyle='dashed')
    plt.gca().set_aspect('equal',adjustable='box')
    plt.title(tstr,loc='left',fontsize=fsz)
    
    plt.fill([xmin,xmax,xmax,xmin],[ymin,ymin,ymax,ymax],color='#aaaaaa',alpha=0.3)
    plt.fill([dp1,dp2,dp2,dp1],[dh1,dh1,dh2,dh2],color='#ffffff')
    xs=dp1; ys=0.5*(dh1+dh2); ss='$v_P$=7.0m/s'
    plt.text(xs,ys,ss,color=co0,va='center',ha='left',rotation=90,fontsize=fsz)
    xs=dp2; ys=0.5*(dh1+dh2); ss='$v_P$=2.5m/s'
    plt.text(xs,ys,ss,color=co0,va='center',ha='right',rotation=90,fontsize=fsz)
    xs=4.75; ys=dh1; ss='$v_H$=5.0m/s'
    plt.text(xs,ys,ss,color=co0,va='bottom',ha='center',rotation=0,fontsize=fsz)
    xs=4.75; ys=dh2; ss='$v_H$=2.0m/s'
    plt.text(xs,ys,ss,color=co0,va='top',ha='center',rotation=0,fontsize=fsz)

    cont=plt.contour(xx,yy,zc,colors=[co1],levels=[6000,6500,7000,7500,8000,8200,8500,9000])
    cont.clabel(fmt='%1.0f', fontsize=fsz-1)
    plt.contour(xx,yy,zh,colors=[co2],levels=[16.58])
    xs=5.6; ys=5.65; ss='$h_L$=16.58m'
    plt.text(xs,ys,ss,color='#ff0000',va='bottom',ha='center',rotation=0,fontsize=fsz)
    
    
    fnameF='fig_contour.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()


def headrace(dd):
    a1=164619 # NRP/m
    b1=450019 # NRP/m
    q=68.64
    al=10149.0
    n=0.016
    v=4*q/np.pi/dd**2
    hl=124.5*n**2*al/dd**(4/3)*v**2/2/9.8
    cc=(a1*dd**2/6.0**2+b1*dd/6.0)*al
    return hl,cc

def penstock(dd):
    a2=143829
    b2=2437167
    cf=177919810
    q=68.64
    al=688.2
    n=0.011
    rho=20.0
    fbr=0.4
    v=4*q/np.pi/dd**2
    h1=124.5*n**2*al/dd**(4/3)
    h2=4*(0.131+0.1632*(dd/rho)**(7/2))
    h3=fbr
    hl=(h1+h2+h3)*v**2/2/9.8+1.991
    cc=(a2*dd**2/4.0**2+b2*dd/4.0)*al+cf
    return hl,cc


def main():
    ds=0.1
    d1=np.arange(4.0,7.0+ds,ds)
    d2=np.arange(3.0,6.0+ds,ds)
    hlt=np.zeros((len(d1),len(d2)),dtype=np.float64)
    cct=np.zeros((len(d1),len(d2)),dtype=np.float64)
    for i,dd1 in enumerate(d1):
        hl1,cc1=headrace(dd1)
        for j,dd2 in enumerate(d2):
            hl2,cc2=penstock(dd2)
            hlt[i,j]=(hl1+hl2)*1.1
            #cct[i,j]=(cc1+cc2)/(450-hlt[i,j])*433.42/1e6
            cct[i,j]=(cc1+cc2)/1e6
            if 16.0<=hlt[i,j]<=16.58: print('{0:6.1f}{1:6.1f}{2:8.3f}{3:8.3f}{4:8.3f}{5:12.3f}'.format(dd1,dd2,hl1*1.1,hl2*1.1,hlt[i,j],cct[i,j]))
    xx, yy = np.meshgrid(d2,d1)
    zh=hlt
    zc=cct
    drawfig(xx,yy,zh,zc)



if __name__ == '__main__':
    main()

以　上

2020-01-04

Python 岩盤内埋設式水圧鉄管設計プログラム（旧版）

概要

岩盤内埋設式水圧鉄管の設計プログラムをフルで書いてみたので残しておく。

プログラム前半の非常に長い関数「def xlswrite」は、エクセルに書式指定して結果を書き出すもの。最終的にワードの報告書に貼り付けなくてはならないので、エクセルにしておくと便利である。このため、面倒ではあったがこの部分も作っておいた。

入力

プログラムでは 'load2.xlsx' を読み込んでいるが、これは出力例「荷重計算」と同じ内容（欄外の変数凡例行は除く）のファイルである。

出力

プログラムでは出力は'penstock2.xlsx' というエクセルファイルに行われ、以下のシートに結果が出力される。

Sheet	Outout
Load	荷重計算（入力ファイルと同じ）
Pin	内圧設計結果
Pex	外圧設計結果

なお設計では、関数「calc1」の中で、水門鉄管技術基準で定める最小板厚とは別に、施工性や耐外厚設計などを考慮し、これより大きな推奨最小板厚（recommended minimum thickness: t0t）というものを導入し、板厚がこれ以下にならないようにしている。

プログラム上での使用鋼材はJIS：SM400, SM490, SM570である。

荷重計算

f:id:damyarou:20200112090904j:plain

内圧設計

f:id:damyarou:20200112090943j:plain

外圧設計

f:id:damyarou:20200112091003j:plain

プログラム

import numpy as np
import pandas as pd
from scipy import optimize
import openpyxl
from openpyxl.styles.borders import Border, Side
from openpyxl.styles import PatternFill
from openpyxl.styles.fonts import Font


def xlline(ws, cmax):
    def wline(ws, row_num, cmax, bl, bc, br):
        col_num=1; ws.cell(row=row_num ,column=col_num).border = bl
        for col_num in range(2,cmax):
            ws.cell(row=row_num ,column=col_num).border = bc
        col_num=cmax; ws.cell(row=row_num ,column=col_num).border = br
    
    # font
    font = Font(name='Courier New',
            size=12,
            bold=False,
            italic=False,
            underline='none',
            strike=False,
            color='000000')
    for row_num in range(1,53):
        for col_num in range(1,cmax+1):
            ws.cell(row=row_num ,column=col_num).font = font
    # cell height and width
    lcol=['A','B','C','D','E','F','G','H','I','J','K','L','M','N']
    for row_num in range(1,53):
        ws.row_dimensions[row_num].height = 20
    for i in range(0,cmax):
        ws.column_dimensions[lcol[i]].width = 12
    ws.column_dimensions[lcol[cmax-1]].width = 30
    
    # set color for cell filling
    fill = PatternFill(patternType='solid', fgColor='ffffaa')
    # write in sheet
    for row_num in range(6,12):
        for col_num in range(1,cmax+1):
            ws.cell(row=row_num ,column=col_num).fill = fill
    for row_num in range(20,26):
        for col_num in range(1,cmax+1):
            ws.cell(row=row_num ,column=col_num).fill = fill
    
    bul = Border(top=Side(style='medium', color='000000'), 
                 bottom=Side(style='dashed', color='000000'), 
                 left=Side(style='medium', color='000000'),
                 right=Side(style='dashed', color='000000')
                )
    buc = Border(top=Side(style='medium', color='000000'), 
                 bottom=Side(style='dashed', color='000000'), 
                 left=Side(style='dashed', color='000000'),
                 right=Side(style='dashed', color='000000')
                )
    bur = Border(top=Side(style='medium', color='000000'), 
                 bottom=Side(style='dashed', color='000000'), 
                 left=Side(style='dashed', color='000000'),
                 right=Side(style='medium', color='000000')
                )
    bcl = Border(top=Side(style='dashed', color='000000'), 
                 bottom=Side(style='dashed', color='000000'), 
                 left=Side(style='medium', color='000000'),
                 right=Side(style='dashed', color='000000')
                )
    bcc = Border(top=Side(style='dashed', color='000000'), 
                 bottom=Side(style='dashed', color='000000'), 
                 left=Side(style='dashed', color='000000'),
                 right=Side(style='dashed', color='000000')
                )
    bcr = Border(top=Side(style='dashed', color='000000'), 
                 bottom=Side(style='dashed', color='000000'), 
                 left=Side(style='dashed', color='000000'),
                 right=Side(style='medium', color='000000')
                )
    bll = Border(top=Side(style='dashed', color='000000'), 
                 bottom=Side(style='medium', color='000000'), 
                 left=Side(style='medium', color='000000'),
                 right=Side(style='dashed', color='000000')
                )
    blc = Border(top=Side(style='dashed', color='000000'), 
                 bottom=Side(style='medium', color='000000'), 
                 left=Side(style='dashed', color='000000'),
                 right=Side(style='dashed', color='000000')
                )
    blr = Border(top=Side(style='dashed', color='000000'), 
                 bottom=Side(style='medium', color='000000'), 
                 left=Side(style='dashed', color='000000'),
                 right=Side(style='medium', color='000000')
                )

    row_num=1; wline(ws, row_num, cmax, bul, buc, bur)
    row_num=2; wline(ws, row_num, cmax, bul, buc, bur)
    for row_num in range(3,29):
        wline(ws, row_num, cmax, bcl, bcc, bcr)
    row_num=29; wline(ws, row_num, cmax, bul, buc, bur)
    row_num=30; wline(ws, row_num, cmax, bcl, bcc, bcr)
    row_num=31; wline(ws, row_num, cmax, bcl, bcc, bcr)
    row_num=32; wline(ws, row_num, cmax, bcl, bcc, bcr)
    row_num=33; wline(ws, row_num, cmax, bul, buc, bur)
    row_num=34; wline(ws, row_num, cmax, bcl, bcc, bcr)
    row_num=35; wline(ws, row_num, cmax, bcl, bcc, bcr)
    row_num=36; wline(ws, row_num, cmax, bul, buc, bur)
    row_num=37; wline(ws, row_num, cmax, bcl, bcc, bcr)
    row_num=38; wline(ws, row_num, cmax, bll, blc, blr)
    if cmax==11:
        row_num=39; wline(ws, row_num, cmax, bll, blc, blr)
    

def xlswrite(df0,df1,df2):
    # Write results to Excel 
    fnameW='penstock2.xlsx'
    wb = openpyxl.Workbook()
    #
    ws=wb.create_sheet(index=1,title='Load')
    ws.cell(row=1,column=1).value='No'
    ws.cell(row=1,column=2).value='L(m)'
    ws.cell(row=1,column=3).value='Sum(L)'
    ws.cell(row=1,column=4).value='EL(m)'
    ws.cell(row=1,column=5).value='D0(m)'
    ws.cell(row=1,column=6).value='GWL(m)'
    ws.cell(row=1,column=7).value='Hst(m)'
    ws.cell(row=1,column=8).value='Hsg(m)'
    ws.cell(row=1,column=9).value='Hwh(m)'
    ws.cell(row=1,column=10).value='Hin(m)'
    ws.cell(row=1,column=11).value='Hex(m)'
    ws.cell(row=1,column=12).value='Remarks'
    #
    j=1; temp=np.array(df0['No'],dtype=np.int)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=2; temp=np.array(df0['L(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=3; temp=np.array(df0['Sum(L)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=4; temp=np.array(df0['EL(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=5; temp=np.array(df0['D0(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=6; temp=np.array(df0['GWL(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=7; temp=np.array(df0['Hst(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=8; temp=np.array(df0['Hsg(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=9; temp=np.array(df0['Hwh(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=10; temp=np.array(df0['Hin(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=11; temp=np.array(df0['Hex(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=12; temp=df0['Remarks']
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).value = temp[i]
    # border line
    cmax=12; xlline(ws, cmax)
    # legend
    i=38
    i=i+1; ws.cell(row=i,column=1).value = 'No: section number'
    i=i+1; ws.cell(row=i,column=1).value = 'L: length of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 'Sum(L): cumulative length'
    i=i+1; ws.cell(row=i,column=1).value = 'EL: elevation of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 'D0: internal diameter of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 'GWL: ground water level'
    i=i+1; ws.cell(row=i,column=1).value = 'Hst: hydrostatic head'
    i=i+1; ws.cell(row=i,column=1).value = 'Hsg: pressure head rise by surging'
    i=i+1; ws.cell(row=i,column=1).value = 'Hwh: pressure head rise by water hammering'
    i=i+1; ws.cell(row=i,column=1).value = 'Hin: design internal pressure head'
    i=i+1; ws.cell(row=i,column=1).value = 'Hex: design external pressure head'
    #
    ws=wb.create_sheet(index=1,title='Pin')
    ws.cell(row=1,column=1).value='No'
    ws.cell(row=1,column=2).value='Pi(MPa)'
    ws.cell(row=1,column=3).value='L(m)'
    ws.cell(row=1,column=4).value='D0(mm)'
    ws.cell(row=1,column=5).value='t0(mm)'
    ws.cell(row=1,column=6).value='steel'
    ws.cell(row=1,column=7).value='lam'
    ws.cell(row=1,column=8).value='sig(MPa)'
    ws.cell(row=1,column=9).value='siga(MPa)'
    ws.cell(row=1,column=10).value='Weight(t)'
    ws.cell(row=1,column=11).value='Remarks'
    #
    j=1; temp=np.array(df1['No'],dtype=np.int)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=2; temp=np.array(df1['Pi(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=3; temp=np.array(df1['L(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=4; temp=np.array(df1['D0(mm)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=5; temp=np.array(df1['t0(mm)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=6; temp=df1['steel']
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).value = temp[i]
    j=7; temp=np.array(df1['lam'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=8; temp=np.array(df1['sig(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=9; temp=np.array(df1['siga(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=10; temp=np.array(df1['weight(t)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    ws.cell(row=len(temp)+2,column=1).value = 'Sum'
    ws.cell(row=len(temp)+2,column=10).number_format = "0.000" 
    ws.cell(row=len(temp)+2,column=10).value = np.sum(temp)      
    j=11; temp=df0['Remarks']
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).value = temp[i]
    # border line
    cmax=11; xlline(ws, cmax)
    # legend
    i=39
    i=i+1; ws.cell(row=i,column=1).value = 'No: section number'
    i=i+1; ws.cell(row=i,column=1).value = 'Pi: design internal pressure'
    i=i+1; ws.cell(row=i,column=1).value = 'L: length of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 'D0: internal diameter of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 't0: plate thickness of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 'steel: kind of steel'
    i=i+1; ws.cell(row=i,column=1).value = 'lam: internal pressure shared ratio be bedrock'
    i=i+1; ws.cell(row=i,column=1).value = 'sig: stress of pipe shell'
    i=i+1; ws.cell(row=i,column=1).value = 'siga: allowable stress'
    i=i+1; ws.cell(row=i,column=1).value = 'Weight: weight of pipe section'
    #
    ws=wb.create_sheet(index=1,title='Pex')
    ws.cell(row=1,column=1).value='No'
    ws.cell(row=1,column=2).value='Pe(MPa)'
    ws.cell(row=1,column=3).value='L(m)'
    ws.cell(row=1,column=4).value='D0(mm)'
    ws.cell(row=1,column=5).value='t0(mm)'
    ws.cell(row=1,column=6).value='steel'
    ws.cell(row=1,column=7).value='pk_0(MPa)'
    ws.cell(row=1,column=8).value='SF_0'
    ws.cell(row=1,column=9).value='pk_s(MPa)'
    ws.cell(row=1,column=10).value='SF_s'
    ws.cell(row=1,column=11).value='sc_r(MPa)'
    ws.cell(row=1,column=12).value='sc_c(MPa)'
    ws.cell(row=1,column=13).value='SF_c'
    ws.cell(row=1,column=14).value='stiffener'
    #
    j=1; temp=np.array(df2['No'],dtype=np.int)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=2; temp=np.array(df2['Pe(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=3; temp=np.array(df2['L(m)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=4; temp=np.array(df2['D0(mm)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=5; temp=np.array(df2['t0(mm)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=6; temp=df2['steel']
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).value = temp[i]
    j=7; temp=np.array(df2['pk0(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=8; temp=np.array(df2['sf0'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=9; temp=np.array(df2['pks(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=10; temp=np.array(df2['sfs'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=11; temp=np.array(df2['scr(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=12; temp=np.array(df2['scc(MPa)'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=13; temp=np.array(df2['sfc'],dtype=np.float64)
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).number_format = "0.000" 
        ws.cell(row=i+2,column=j).value = temp[i]
    j=14; temp=df2['stiffener']
    for i in range(len(temp)):
        ws.cell(row=i+2,column=j).value = temp[i]
    # border line
    cmax=14; xlline(ws, cmax)
    # legend
    i=38
    i=i+1; ws.cell(row=i,column=1).value = 'No: section number'
    i=i+1; ws.cell(row=i,column=1).value = 'Pe: design external pressure'
    i=i+1; ws.cell(row=i,column=1).value = 'L: length of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 'D0: internal diameter of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 't0: plate thickness of pipe section'
    i=i+1; ws.cell(row=i,column=1).value = 'steel: kind of steel'
    i=i+1; ws.cell(row=i,column=1).value = 'pk_0: critical buckling pressure of pipe section without stiffener'
    i=i+1; ws.cell(row=i,column=1).value = 'SF_0: safety factor against external pressure without pressure'
    i=i+1; ws.cell(row=i,column=1).value = 'pk_s: critical buckling pressure of pipe section with stiffener'
    i=i+1; ws.cell(row=i,column=1).value = 'SF_s: safety factor against external pressure with stiffener'
    i=i+1; ws.cell(row=i,column=1).value = 'sc_r: critical buckling stress of stiffener'
    i=i+1; ws.cell(row=i,column=1).value = 'sc_c: compressive stress of stiffener'
    i=i+1; ws.cell(row=i,column=1).value = 'SF_c: safety factor against compressive stress in stiffener'
    i=i+1; ws.cell(row=i,column=1).value = 'size of stiffener plate: 75mm height x 20mm thickness'
    #    
    wb.save(fnameW)


def stcr(dd0,t0,eps,hr,tr,el,pp,sigF):
    # Calculation of critical stress and compressive stress of stiffener
    def calrg(dd0,t0,eps,hr,tr):
        rm=0.5*(dd0+t0)
        tt=t0-eps
        # change symbols
        b=1.56*np.sqrt(rm*tt)+tr
        d=tt+hr
        s=tt
        t=tr
        # calculation of radius of gyration of area
        e2=(d**2*t+s**2*(b-t))/2/(b*s+t*(d-s))
        e1=d-e2
        ii=1/3*(t*e1**3+b*e2**3-(b-t)*(e2-s)**3)
        aa=s*b+t*(d-s)
        rg=np.sqrt(ii/aa) # radius of gyration of area
        ee=np.abs(e2-s)
        return rg,ee
    def func(x,params):
        # x=sigN
        k0   =params[0]
        rm   =params[1]
        t    =params[2]
        Es   =params[3]
        sigF =params[4]
        rg   =params[5]
        ee   =params[6]
        aa = k0/rm+x/Es
        bb = 1+rm**2/rg**2*x/Es
        cc = 1.68*rm/ee*(sigF-x)/Es*(1-1/4*rm/ee*(sigF-x)/Es)
        f = aa*bb**1.5-cc
        return f
    # main routine
    eps=2.0         # margin of plate thickness
    Es=206000.0     # elastic modulus of steel plate
    ns=0.3         # Poisson's ratio of steel plate
    #D0=2100.0  # design internal diameter of penstock
    #t0=30.0    # design plate thickness of penstock
    #sigF=315.0 # yield stress of steel plate (SM490)
    #siga=175.0 # allowable stress of steel plate (SM490)
    t=t0-eps
    rm=(dd0+t0)/2
    r0d=(dd0+2*t0)/2
    k0=0.4e-3*rm
    # calculation of critical stress
    rg,ee=calrg(dd0,t0,eps,hr,tr)
    params=np.array([k0,rm,t,Es,sigF,rg,ee])
    sigN=optimize.brentq(func,0.0,5*sigF,args=(params))
    scr=sigN*(1-r0d/ee*(sigF-sigN)/(1+3/2*np.pi)/Es)
    # calculation of compressive stress
    s0=tr*(t+hr)
    iis=tr/12*(hr+t)**3
    beta=(3*(1-ns**2))**0.25/np.sqrt(rm*t)
    c1=r0d**2/t-(tr+1.56*np.sqrt(rm*t))*r0d**2/(s0+1.56*t*np.sqrt(rm*t)) 
    c2=3/(3*(1-ns**2))**0.75*(r0d/t)**1.5*(np.sinh(beta*el)+np.sinh(beta*el))/(np.cosh(beta*el)-np.cosh(beta*el))+2*r0d**2/(s0+1.56*t*np.sqrt(rm*t))
    pd=pp/(tr+1.56*np.sqrt(rm*t))*((tr+1.56*np.sqrt(rm*t))+2*c1/c2)
    sc=pd*r0d*(tr+1.56*np.sqrt(rm*t))/(s0+1.56*t*np.sqrt(rm*t))
    return scr,sc
    

def timo(el,hr,tr,dd0,t0,eps):
    # Calculate criticsl buckling pressure
    # of steel pipe with stiffener by Timoshenko's formula
    # el: interval of stiffener
    # hr: height of stiffener
    # tr: plate thickness of stiffener
    # t: plate thickness
    pi=np.pi
    Es=206000 # elastic modulus of steel
    ns=0.3    # Poisson's ratio of steel
    rm=0.5*(dd0+t0)
    r0d=0.5*(dd0+2*t0)
    t=t0-eps
    s0=tr*(t+hr)
    iis=tr/12*(hr+t)**3
    beta=(3*(1-ns**2))**0.25/np.sqrt(rm*t)
    c1=r0d**2/t-(tr+1.56*np.sqrt(rm*t))*r0d**2/(s0+1.56*t*np.sqrt(rm*t)) 
    c2=3/(3*(1-ns**2))**0.75*(r0d/t)**1.5*(np.sinh(beta*el)+np.sinh(beta*el))/(np.cosh(beta*el)-np.cosh(beta*el))+2*r0d**2/(s0+1.56*t*np.sqrt(rm*t))
    lt=2/(tr+1.56*np.sqrt(rm*t))*c1/c2
    lam=1-(1+lt)*(1+tr/(1.56*np.sqrt(rm*t)))/(1+s0/(1.56*t*np.sqrt(rm*t)))
    ell=(el+1.56*np.sqrt(rm*t)*np.arccos(lam))*(1+0.037*np.sqrt(rm*t)/(el+1.56*np.sqrt(rm*t)*np.arccos(lam))*t**3/iis)
    apk=[]
    for n in range(2,30):
        c1=(1-ns**2)/(n**2-1)/(1+n**2*ell**2/pi**2/r0d**2)**2
        c2=t**2/12/r0d**2*((n**2-1)+(2*n**2-1-ns)/(1+n**2*ell**2/pi**2/r0d**2))
        _pk=(c1+c2)*Es*t/(1-ns**2)/r0d
        apk=apk+[_pk]
    pk=min(apk)
    return pk


def ams(dd0,t0,sigF,siga):
    # Calculate criticsl buckling pressure
    # of steel pipe without stiffener by Amstutz's formula
    def func(x,params):
        # x=sigN
        k0   =params[0]
        rm   =params[1]
        t    =params[2]
        Ess  =params[3]
        sigFs=params[4]
        aa = k0/rm+x/Ess
        bb = 1+12*rm**2/t**2*x/Ess
        cc = 3.36*rm/t*(sigFs-x)/Ess*(1-1/2*rm/t*(sigFs-x)/Ess)
        f = aa*bb**1.5-cc
        return f
    def cal_pk(sigN,sigFs,Ess,rm,t):
        aa=rm/t*(1+0.35*rm/t*(sigFs-sigN)/Ess)
        pk = sigN/aa
        return pk
    def cal_k0(siga,eta,Es,alphas,deltaT,r0d,betag):
        aa = alphas*deltaT+betag*siga*eta/Es
        k0 = aa*r0d/(1+betag)
        return k0    
    # main routine
    eps=2.0         # margin of plate thickness
    Es=206000.0     # elastic modulus of steel plate
    pos=0.3         # Poisson's ratio of steel plate
    eta=0.85         # welding efficiency
    alphas=1.2e-5   # thermal expansion coefficient
    deltaT=20.0     # temperature decrease
    betag=1.0       # plastic deformation coefficient of bedrock
    #D0=2100.0  # design internal diameter of penstock
    #t0=30.0    # design plate thickness of penstock
    #sigF=315.0 # yield stress of steel plate (SM490)
    #siga=175.0 # allowable stress of steel plate (SM490)
    Ess=Es/(1-pos**2)
    mu=1.5-0.5*1/(1+0.002*Es/sigF)**2
    sigFs=mu*sigF/(1-pos+pos**2)**0.5
    t=t0-eps
    rm=(dd0+t0)/2
    r0d=(dd0+2*t0)/2
    k0=cal_k0(siga,eta,Es,alphas,deltaT,r0d,betag)
    k0=0.4e-3*rm
    params=np.array([k0,rm,t,Ess,sigFs])
    sigN=optimize.brentq(func,0.0,sigFs,args=(params))
    pk=cal_pk(sigN,sigFs,Ess,rm,t)
    return pk


def datainp():
    df0 = pd.read_excel('load2.xlsx')
    return df0


def calc1(num,pp,dd0,sa,eps):
    # Design for internal pressure
    # Calculate required plate thickness
    Es=206000.0
    alpha=1.2e-5
    deltaT=20.0
    Ec=20600.0
    Bc=0.0
    Eg=2000.0
    Bg=1.0
    mg=4.0
    dr=5.2
    
    t0t=25 # recommended minimum thickness
    t0min=np.ceil((dd0+800)/400)
    dd=dd0-eps
    lam=0
    tt=pp*dd/2/sa
    t0=np.ceil(tt+eps)
    if t0<t0min: t0=t0min
    if t0<t0t: t0=t0t
    # Internal pressure shared design by bedrock
    if 5<=num<=10 or 20<=num<=25:
        c1=sa-Es*alpha*deltaT
        c2=(1+Bc)*Es/Ec*2/dd*np.log(dr/dd)+(1+Bg)*Es/Eg*(mg+1)/mg*2/dd
        tt=pp*dd/2/sa-c1/c2/sa
        lam1=1-Es/pp*alpha*deltaT*2*tt/dd
        lam2=1+(1+Bc)*Es/Ec*2*tt/dd*np.log(dr/dd)+(1+Bg)*Es/Eg*(mg+1)/mg*2*tt/dd
        lam=lam1/lam2
        t0=np.ceil(tt+eps)
        if t0<t0min: t0=t0min
        if t0<t0t: t0=t0t
        tt=t0-eps
        lam1=1-Es/pp*alpha*deltaT*2*tt/dd
        lam2=1+(1+Bc)*Es/Ec*2*tt/dd*np.log(dr/dd)+(1+Bg)*Es/Eg*(mg+1)/mg*2*tt/dd
        lam=lam1/lam2    
    sig=pp*dd/2/(t0-eps)*(1-lam)
    return t0,sig,lam


def main():    
    df0=datainp()
    no=np.array(df0['No'],dtype=np.int)
    al=np.array(df0['L(m)'],dtype=np.float64)        # length (m)
    d0=np.array(df0['D0(m)'],dtype=np.float64)*1000  # internal diameter (mm)
    pi=np.array(df0['Hin(m)'],dtype=np.float64)*0.01 # internal pressure (MPa)
    pe=np.array(df0['Hex(m)'],dtype=np.float64)*0.01 # external pressure (MPa)
    eps=2.0
    ef=0.85

    a_t0=np.zeros(len(no),dtype=np.float64)
    a_sig=np.zeros(len(no),dtype=np.float64)
    a_lam=np.zeros(len(no),dtype=np.float64)
    a_sa=np.zeros(len(no),dtype=np.float64)
    a_sname=[]
    a_pk0=np.zeros(len(no),dtype=np.float64)
    a_sf0=np.zeros(len(no),dtype=np.float64)
    a_pks=np.zeros(len(no),dtype=np.float64)
    a_sfs=np.zeros(len(no),dtype=np.float64)
    a_scr=np.zeros(len(no),dtype=np.float64)
    a_scc=np.zeros(len(no),dtype=np.float64)
    a_sfc=np.zeros(len(no),dtype=np.float64)
    a_stiff=[]
    
    for num,pp,dd0,ppe in zip(no,pi,d0,pe):
        sname='SM400'
        sa=130*ef
        t0,sig,lam=calc1(num,pp,dd0,sa,eps)
        if sname=='SM400' and 40<t0:
            sname='SM490'
            sa=175*ef
            t0,sig,lam=calc1(num,pp,dd0,sa,eps)
        if sname=='SM490' and 40<t0:
            sname='SM570'
            sa=240*ef
            t0,sig,lam=calc1(num,pp,dd0,sa,eps)
        if sname=='SM570' and 40<t0:
            sname='SM570'
            sa=235*ef
            t0,sig,lam=calc1(num,pp,dd0,sa,eps)
        a_t0[num-1]=t0
        a_sig[num-1]=sig/ef
        a_lam[num-1]=lam
        a_sa[num-1]=sa/ef
        a_sname=a_sname+[sname]
        if sname=='SM400': sigF=235
        if sname=='SM490': sigF=315
        if sname=='SM570' and t0<=40: sigF=450
        if sname=='SM400' and 40<t0: sigF=430
        # without stiffener
        pk0=ams(dd0,t0,sigF,sa/ef)
        sf0=pk0/ppe
        if sf0<1.5:
            # stiffener
            hr=75; tr=20
            # stiffener interval=2000mm
            el=3000.0
            pk1=timo(el,hr,tr,dd0,t0,eps)
            sf1=pk1/ppe
            scr1,sc1=stcr(dd0,t0,eps,hr,tr,el,ppe,sigF)
            sfc1=scr1/sc1
            # stiffener interval=1500mm
            el=1500.0
            pk2=timo(el,hr,tr,dd0,t0,eps)
            sf2=pk2/ppe
            scr2,sc2=stcr(dd0,t0,eps,hr,tr,el,ppe,sigF)
            sfc2=scr2/sc2
            # stiffener interval=1000mm
            el=1000.0
            pk3=timo(el,hr,tr,dd0,t0,eps)
            sf3=pk3/ppe
            scr3,sc3=stcr(dd0,t0,eps,hr,tr,el,ppe,sigF)
            sfc3=scr3/sc3
        # store results
        a_pk0[num-1]=pk0
        a_sf0[num-1]=sf0        
        if 1.5<=sf0:
            a_pks[num-1]=0
            a_sfs[num-1]=0
            a_stiff=a_stiff+['no-need']
            a_scr[num-1]=0
            a_scc[num-1]=0
            a_sfc[num-1]=0
        elif sf0<1.5 and 1.5<=sf1:
            a_pks[num-1]=pk1
            a_sfs[num-1]=sf1
            a_stiff=a_stiff+['interval=3000mm']
            a_scr[num-1]=scr1
            a_scc[num-1]=sc1
            a_sfc[num-1]=sfc1
        elif sf1<1.5 and 1.5<=sf2:
            a_pks[num-1]=pk2
            a_sfs[num-1]=sf2
            a_stiff=a_stiff+['interval=1500mm']
            a_scr[num-1]=scr2
            a_scc[num-1]=sc2
            a_sfc[num-1]=sfc2
        elif sf2<1.5 and 1.5<=sf3:
            a_pks[num-1]=pk3
            a_sfs[num-1]=sf3
            a_stiff=a_stiff+['interval=1000mm']
            a_scr[num-1]=scr3
            a_scc[num-1]=sc3
            a_sfc[num-1]=sfc3
    
    ww=0.25*np.pi*((d0+2*a_t0)**2-d0**2)/1000/1000*al*7.85 # weight
    print('Weight:{0:10.3f} ton'.format(np.sum(ww)))
    print('Average thickness: {0:6.1f} mm'.format(np.sum(a_t0*al)/np.sum(al)))
    
    df1=pd.DataFrame({
        'No': no,
        'Pi(MPa)': pi,
        'L(m)' : al,
        'D0(mm)': d0,
        't0(mm)': a_t0,
        'steel': a_sname,
        'lam': a_lam,
        'sig(MPa)': a_sig,
        'siga(MPa)': a_sa,
        'weight(t)': ww,
        'Remarks': df0['Remarks']
    })
    df2=pd.DataFrame({
        'No': no,
        'Pe(MPa)': pe,
        'L(m)' : al,
        'D0(mm)': d0,
        't0(mm)': a_t0,
        'steel': a_sname,
        'pk0(MPa)': a_pk0,
        'sf0': a_sf0,
        'pks(MPa)': a_pks,
        'sfs': a_sfs,
        'scr(MPa)': a_scr,
        'scc(MPa)': a_scc,
        'sfc': a_sfc,
        'stiffener': a_stiff
    })
    with pd.ExcelWriter('penstock.xlsx') as writer:
        df0.to_excel(writer,sheet_name='Load',index=False)
        df1.to_excel(writer,sheet_name='Internal Pressure',index=False)
        df2.to_excel(writer,sheet_name='External Pressure',index=False)

        
    xlswrite(df0,df1,df2)
        
        
if __name__ == '__main__':
    main()

以　上