damyarou

python, GMT などのプログラム

LG 27インチ４Kディスプレイ購入

９月７日注文、９月８日到着で、LGの27inch 4K display を購入しました。￥57,980也。ちょっと奮発しすぎたか。

４月より在宅勤務しています。会社に行くのは月に２〜３回ですかね。 しかし７月よりお客さんから設計の仕事をいただき、今はかなり忙しい。まあ色々検討してレポートを書いているのですが、ディスプレイ上の作業領域が「もっと欲しい！」ということで、サブディスプレイを購入したのでした。LGにした理由は、USB-Cポートがついていたからというのが大きいです。このポートがあると、ディスプレイとiMacをUSB-Cケーブル１本でつなぐことができます。 最初は１万円くらいのやつでいいかと思っていたのですが、色々Webを見ながら考えていたら、「目も悪くなってきたし、せっかくなら４Kにするか！」と思い立ったのでした。実際には27インチで4Kだと、100%のscaleでは文字が小さくて読めません。なので200%とかに拡大しています。「拡大したら4Kの意味ないじゃん」という向きもありますが、画面のドットそのものは精細なので、文字はとてもきれいです。 買った当初は縦置きにしようと思い、試しに縦置きで使ってみたのですが、上の方がとても遠く感じます。そこで通常の横置きに。図面などを見るにも視野が広くなり便利です。このディスプレイは横から縦へ回転させることができるので、縦置きにするのは、南北に長い地形図やGoogle earthを見るときだけとし、普段は横置きとすることにしました。

はじめは縦置きにしてみましたが、あまり使いよくない。 f:id:damyarou:20200910154018j:plain

基本は横置きです。図面を見る時など視野が広くて便利。 f:id:damyarou:20200910154049j:plain

Google Earthなどで縦に画像を表示させたい場合は便利。 f:id:damyarou:20200910154122j:plain

とりあえず持っている機械類の集合写真。この写真はiPhoneで撮っているので、かわいそうに、iPhoneは写っていません。 正面が新しいディスプレイ、左側のWin機は会社支給のPC. f:id:damyarou:20200910154154j:plain

つなぎ方のイメージ図。iMac以外はHDMIでつないでいます。 f:id:damyarou:20200910154221p:plain

以　上

Python 逆調整池運用解析修正版（Flood Routine）

逆調整池運用解析プログラムの修正版

  • 関数FLOODの中の収束計算で、前回貯水池水位のみを既知としていたが、本来貯水池水位と放流量はペアなので、前回放流量も既知とした。プログラム作成時になぜ収束計算の中で前回流量を既知として確定していなかったのかは思い出せない。
  • 関連出力は、データフレームにして一気に書き出すことにした。

以下最新版プログラム（2020.08.08）

# ==============================================
# Flood Routine Analysis
# ==============================================
import sys
import numpy as np
import pandas as pd
from scipy import interpolate


def FLOOD(iOFC,fnameR1,fnameR2,fnameR3,ram,qrel):
    nr,rh,rv=INP_HV(fnameR1)                   # reservoir H-V
    nd,ti,q_in,ELini,Qoini=INP_INFLOW(fnameR2) # inflow time history
    no,elof,qof=INP_OUTFLOW(fnameR3)           # outflow capacity

    pEL=np.zeros(nd,dtype=np.float64) # storage of reervoir water level
    pQo=np.zeros(nd,dtype=np.float64) # storage of outflow
    pDE=np.zeros(nd,dtype=np.float64) # storage of error
    pVO=np.zeros(nd,dtype=np.float64) # storage of reservoir volume
    pRA=np.zeros(nd,dtype=np.float64) # storage of ramp up rate
    pUD=np.zeros(nd,dtype=int)        # storage of indicator od water level
    
    Vin=np.sum(q_in)*(ti[2]-ti[1])/2 # total volume of inflow
    Qimax=np.max(q_in)
    # Initial setting
    elv=ELini
#    elv=WLINI(q_in[0],elof,qof,no)     # Initial reservoir water level
    EL=elv
    VOL=RET_V(nr,rh,rv,elv)           # Initial reservoir volume
    q_out=Qoini
    qout0=Qoini
    print('q_in0={0:10.3f}'.format(q_in[0]))
    print('elv_ini={0:10.3f}'.format(EL))
    print('q_ini={0:10.3f}'.format(q_out))
    i=0
    iUD=0
    pUD[i]=iUD
    pEL[i]=EL
    pQo[i]=q_out
    pDE[i]=EL-elv
    pVO[i]=VOL
    pRA[i]=0
    # Iterative calculation
    iUD=1        # Assume reservoir water level increase
    dh=0.0005    # Water level increment for searching equilibrium point
    eps=0.001    # Convergence criteria
    itmax=int(1.0/dh)*100
    Vf=RET_V(nr,rh,rv,60.0)
    Vz=0
    Vin2=0
    Vot2=0
    q1=q_out
    for i in range(0,nd-1):
        dt=ti[i+1]-ti[i]             # Time interval
        Qin=0.5*(q_in[i]+q_in[i+1])  # Average inflow rate during dt
        hh=0.0                       # Zero set of water level increment
        k=0
        j=0
        while 1:
            k=k+1
            j=j+1;
            hh=hh+float(iUD)*dh      # Update waler level increment
            elv1=EL
            elv2=EL+hh
            if iOFC==8: q2=OFC8(q1,qout0,elv2,elof,qof,no,ti[i],ram,qrel)
            Qout=0.5*(q1+q2)         # Average outflow rate during dt
            dS=(Qin-Qout)*dt*3600.0  # Storage accumulated during dt
            R=VOL+dS                 # Total reservoir volume
            elv=RET_H(nr,rh,rv,R)    # Reservoir water level at reservoir volume R
            if iUD==1 and j==10:  # 10 times trial assuming reservoir water level increase
                if EL+hh>elv:     # if reservoir water level is less than assumed water level,
                    iUD=-1        # assume water level decrease
                    hh=0.0
                    j=0
            if iUD==-1 and j==10: # 10 times trial assuming reservoir water level decrease
                if EL+hh<elv:     # if reservoir water level is greater than assumed water level,
                    iUD=1         # assume water level decrease
                    hh=0.0
                    j=0
            if np.abs(EL+hh-elv)<eps: break  # Judgement of convergence
            if itmax<k: break
        VOL=R    # Cumulative volume
        EL=EL+hh # Elevation
        q_out=q2 # Outflow
        q1=q2
        
        pUD[i+1]=iUD
        pEL[i+1]=EL
        pQo[i+1]=q_out
        pDE[i+1]=EL-elv
        pVO[i+1]=VOL
        pRA[i+1]=(pQo[i+1]-pQo[i])/(ti[i+1]-ti[i])

    hmax=max(pEL)-ELini
    sys.stdout.write('Time: %10.3f\n'% np.max(ti))
    sys.stdout.write('h   : %10.3f\n'% hmax)
    sys.stdout.write('Qin : %10.3f %10.3f\n'% (np.min(q_in),np.max(q_in)))
    sys.stdout.write('Qout: %10.3f %10.3f\n'% (np.min(pQo),np.max(pQo)))
    sys.stdout.write('EL  : %10.3f %10.3f\n'% (np.min(pEL),np.max(pEL)))
    sys.stdout.write('\n')

    # TWC of 80m downstream of Re-Regulating Dam
    fnameR='inp_twc.txt'
    data=np.loadtxt(fnameR,skiprows=0,usecols=(0,1))
    qq=data[:,0]
    ee=data[:,1]
    ff=interpolate.interp1d(qq,ee)
    el_ds=ff(pQo)
    
    df = pd.DataFrame({ 'i'  : np.arange(0,nd), # numper
                        'iUD': pUD,   # indicator for water level
                        'ti' : ti,    # time
                        'EL' : pEL,   # water level of reservoir
                        'dEL': pDE,   # error
                        'VOL': pVO,   # reservoir volume
                        'Qi' : q_in,  # inflow to reservoir
                        'Qo' : pQo,   # outflow from reservoir
                        'RA' : pRA,   # ramp up rate
                        'ELD': el_ds  # downstream water level
                      })
    df = df.set_index('i')
    return df 


def WLINI(q,elof,qof,no):
    # To obtain initial water level at oitflow discharge q
    for i in range(0,no-1):
        if qof[i]<=q and q<=qof[i+1]: break
    if qof[no-1]<q: i=no-2
    x1=qof[i]
    y1=elof[i]
    x2=qof[i+1]
    y2=elof[i+1]
    a=(y2-y1)/(x2-x1)
    b=(x2*y1-x1*y2)/(x2-x1)
    elv=a*q+b
    if q<0.0: elv=elof[0]
    return elv


def OFC8(q1,qout0,elv2,elof,qof,no,tii,ram,qrel):
    Qlim=300.0
    if   5.0<tii and tii<= 29.0 :qout0=qrel
    if  29.0<tii and tii<= 53.0 :qout0=qrel
    if  53.0<tii and tii<= 77.0 :qout0=qrel
    if  77.0<tii and tii<=101.0 :qout0=qrel
    if 101.0<tii: qout0=39.0
    q2=OFC(elv2,elof,qof,no)
    rq2=q2
    if q2<qout0: qout0=q2
    q2=qout0
    st=0.01
    ts1=0.0
    ts2=24.0
    ts3=48.0
    ts4=72.0
    ts5=96.0
    if ts1+st<=tii and tii<=ts1+5.0+st:
        if elv2<68.0:
            q2=q1+ram*st
            if rq2<q2: q2=rq2
            if Qlim<=q2: q2=Qlim
        else:
            q2=rq2
    if ts2+st<=tii and tii<=ts2+5.0+st:
        if elv2<68.0:
            q2=q1+ram*st
            if rq2<q2: q2=rq2
            if Qlim<=q2: q2=Qlim
        else:
            q2=rq2
    if ts3+st<=tii and tii<=ts3+5.0+st:
        if elv2<68.0:
            q2=q1+ram*st
            if rq2<q2: q2=rq2
            if Qlim<=q2: q2=Qlim
        else:
            q2=rq2
    if ts4+st<=tii and tii<=ts4+5.0+st:
        if elv2<68.0:
            q2=q1+ram*st
            if rq2<q2: q2=rq2
            if Qlim<=q2: q2=Qlim
        else:
            q2=rq2
    if ts5+st<=tii and tii<=ts5+5.0+st:
        if elv2<68.0:
            q2=q1+ram*st
            if rq2<q2: q2=rq2
            if Qlim<=q2: q2=Qlim
        else:
            q2=rq2

    return q2


def OFC(elv,elof,qof,no):
    # Free overflow discharge
    for i in range(0,no-1):
        if elof[i]<=elv and elv<=elof[i+1]: break
    if elof[no-1]<elv: i=no-2
    x1=elof[i]
    y1=qof[i]
    x2=elof[i+1]
    y2=qof[i+1]
    a=(y2-y1)/(x2-x1)
    b=(x2*y1-x1*y2)/(x2-x1)
    q=a*elv+b
    if elv<=elof[0]: q=0.0
    return q

def RET_V(nr,rh,rv,elv):
    # To obtain reservoir cumulative volume from the water level
    for i in range(0,nr-1):
        if rh[i]<=elv and elv<=rh[i+1]: break
    if rh[nr-1]<elv: i=nr-2
    x1=rv[i]
    y1=rh[i]
    x2=rv[i+1]
    y2=rh[i+1]
    a=(y2-y1)/(x2-x1)
    b=(x2*y1-x1*y2)/(x2-x1)
    v=(elv-b)/a
    return v

def RET_H(nr,rh,rv,v):
    # To obtain reservoir water level from cumulative volume
    for i in range(0,nr-1):
        if rv[i]<=v and v<=rv[i+1]: break
    if rv[nr-1]<v: i=nr-2
    x1=rv[i]
    y1=rh[i]
    x2=rv[i+1]
    y2=rh[i+1]
    a=(y2-y1)/(x2-x1)
    b=(x2*y1-x1*y2)/(x2-x1)
    elv=a*v+b
    return elv

def INP_HV(fnameR1):
    # Input reservoir H-V data
    fin=open(fnameR1,'r')
    dat=fin.readlines()
    fin.close()
    n=len(dat)
    text=dat[0]
    text=text.strip()
    text=text.split()
    vcoef=float(text[0])
    rh=np.zeros(n-1,dtype=np.float64)
    rv=np.zeros(n-1,dtype=np.float64)
    for i in range(0,n-1):
        text=dat[i+1]
        text=text.strip()
        text=text.split()
        rh[i]=float(text[0])
        rv[i]=float(text[1])*vcoef
    nr=len(rh)
    return nr,rh,rv

def INP_INFLOW(fnameR2):
    # Input time sequence of inflow
    fin=open(fnameR2,'r')
    dat=fin.readlines()
    fin.close()
    n=len(dat)
    text=dat[0]
    text=text.strip()
    text=text.split()
    ELini=float(text[0])
    Qoini=float(text[1])
    ti  =np.zeros(n-1,dtype=np.float64)
    q_in=np.zeros(n-1,dtype=np.float64)
    for i in range(0,n-1):
        text=dat[i+1]
        text=text.strip()
        text=text.split()
        ti[i]  =float(text[0])
        q_in[i]=float(text[1])
    nd=len(ti)
    return nd,ti,q_in,ELini,Qoini

def INP_OUTFLOW(fnameR3):
    # Input outflow capacity (water level-discharge curve)
    fin=open(fnameR3,'r')
    dat=fin.readlines()
    fin.close()
    n=len(dat)
    elof=np.zeros(n,dtype=np.float64)
    qof =np.zeros(n,dtype=np.float64)
    for i in range(0,n):
        text=dat[i]
        text=text.strip()
        text=text.split()
        elof[i]=float(text[0])
        qof[i] =float(text[1])
    no=len(elof)
    return no,elof,qof


def main():
    lout=['out_150.csv','out_175.csv','out_200.csv','out_225.csv','out_250.csv']
    lram=[150,175,200,225,250]
    lqrel=[89.0,87.1,82.9,81.5,80.7]
    for out,ram,qrel in zip(lout,lram,lqrel):
        iOFC=8
        fnameR1='inp_RHV.txt' # H-V data of reservoir
        fnameR2='inp_c00.txt' # Inflow time history (Hydrograph)
        fnameR3='inp_OF8.txt' # Outflow capacity
        fnameW = out # Output file name
        df=FLOOD(iOFC,fnameR1,fnameR2,fnameR3,ram,qrel)
        df.to_csv(fnameW, sep=",")
    

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

以　上

Relationship between critical depth and discharge coefficient for free overflow on dam spillway

In general, the discharge for free overflow on dam spillway can be expressed as follow:

\begin{equation} Q=C\cdot B\cdot H^{3/2} \end{equation}
\begin{align} &Q & & \text{discharge}\\ &C & & \text{discharge coefficient}\\ &B & & \text{length of overflow crest}\\ &H & & \text{design head including approach velocity head} \end{align}

On the other hand, a critical depth of rectangular cross section channel can be expressed as follow:

\begin{equation} h_c=\left(\cfrac{Q^2}{g\cdot B^2}\right)^{1/3} \end{equation}
\begin{align} &h_c & & \text{critical depth}\\ &Q & & \text{discharge}\\ &B & & \text{length of overflow crest}\\ &g & & \text{gravity acceleration} \end{align}

When it is assumed the water depth at the overflow crest of the dam is equal to the critical depth, the design head including approach velocity head can be expressed as follow:

\begin{equation} H=h_c+\cfrac{1}{2g}\left(\cfrac{Q}{B\cdot h_c}\right)^2 \end{equation}

Considering a unit length of overflow crest which means B=1.0 m,

\begin{equation} Q=C\cdot H^{3/2} \qquad H=h_c+\cfrac{1}{2g}\left(\cfrac{Q}{h_c}\right)^2 \qquad h_c=\left(\cfrac{Q^2}{g}\right)^{1/3} \end{equation}

From above, H can be expressed as follow:

\begin{equation} H=\cfrac{3}{2}\cdot\left(\cfrac{Q^2}{g}\right)^{1/3} \end{equation}

Therefore, discharge coefficient C becomes a constant and it can be expressed as follow:

\begin{equation} C=\cfrac{Q}{H^{3/2}}=\left(\cfrac{8}{27}\cdot g\right)^{1/2}=1.704 \qquad (g=9.8 m/s^2) \end{equation}

That is all.

Python：計算結果をエクセル出力し加工する（月平均値の算出とエクセル化）

能書き

報告書を作成しているとき、Pythonで計算した結果をWordに表形式で貼り付けたい場合がある。 これが、TeXやhtmlであれば、プロブラム内でタグを埋め込んだ形で標準出力し、コピー＆ペーストでTeXもしくはhtml文書に貼り付けるのだが、相手がWordとなるとなかなか面倒である。PythonからWordを制御することができることも知っている（使ったことはない）が、そのためのみのコードを書くのも大変だし、後処理計算での活用などを考えるとエクセル出力しておき、それをwordに貼り付けるほうが、実用的である。

私がよく使う数表をwordに貼り付ける方法は以下の３つである。

(1) 計算結果をコンマ区切りで標準出力し、Wordにコピー＆ペースト。その後はWordで[Insert]=>[table]=>[Convert Text to Table]=>[Separate test at Commas]=.[ok]の手順で、Word上で表に加工する方法。Word上で加工できるので仕上がりは結構綺麗にいく。

(2) 計算結果をエクセルに出力しフォーマットを整える。その後はエクセル上で範囲指定しクリップボードコピー、Word上で[Paste Special]=>[As Picture(PNG)]=>[ok]の手順で、画像としてWordに貼り付ける方法。画像としてwordに貼り付けるためWord上で加工できず、結果がぼやける場合もあり、最近は使わない。

(3)計算結果の表をhtml出力してブラウザで表示し、それをコピー＆ペーストでwordに貼り付け、仕上げはword上で行う。ブラウザはsafariで行うのがやりやすいようである、

方法(3)は最近使い出したのだが、面倒なようでなかなか良い。 というのも計算結果は一回出力するのみでなく、次の処理で使う場合も多く、そうなると１ファイルに複数シートを配置できるエクセルで保存しておくのも後工程での使い勝手がよい。

エクセルからhtmlへの変換

エクセルからhtmlへの出力は、ネット上で便利なものが有り、以下のものを使わせてもらっている。

ExcelからHTMLテーブルタグ変換

エクセルから変換したhtmlの表を貼り付けたものが下の表。二重線が消えているがhtml出力が最終成果ではないので気にしない。 必要あれば装飾をかければ良い。

Year Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Ave
2001 19.6 17.8 16.1 14.5 21.3 36.1 35.9 84.0 36.0 47.7 39.4 26.7 33.1
2002 22.1 20.4 18.5 16.7 24.1 19.9 43.1 36.7 44.8 32.9 29.6 23.8 27.8
2003 21.4 19.7 18.0 16.2 14.6 13.9 20.6 22.9 30.7 29.9 28.1 21.5 21.5
2004 19.2 17.6 15.9 14.3 21.1 37.7 47.0 48.7 110.6 47.1 32.0 24.2 36.2
2005 22.0 20.0 18.1 16.3 14.9 15.0 44.5 49.7 92.6 45.4 32.2 25.6 33.0
2006 22.3 20.5 18.6 16.7 15.2 16.1 32.7 65.6 104.4 62.3 34.7 25.6 36.3
2007 22.8 20.8 18.8 16.9 17.7 20.3 24.5 35.7 29.5 33.4 26.8 21.8 24.1
2008 20.0 18.3 16.6 15.0 13.7 25.0 56.6 83.9 74.4 41.0 36.4 28.9 35.9
2009 24.1 22.2 20.1 18.1 16.4 15.3 18.0 37.2 31.9 29.1 25.4 20.2 23.2
2010 18.6 17.0 15.3 13.8 12.4 11.8 34.5 114.9 89.1 45.6 32.9 24.4 36.0
2011 21.8 19.9 18.0 16.4 17.9 28.9 42.7 129.4 111.2 55.2 36.3 27.8 43.9
2012 24.9 22.7 20.4 18.4 17.8 25.9 41.6 60.0 72.2 42.4 33.8 26.5 33.9
2013 23.6 21.6 19.6 17.6 15.9 16.2 28.0 103.1 82.5 52.8 38.4 28.2 37.4
2014 23.7 21.8 19.8 17.8 16.0 18.3 57.8 56.8 58.2 34.4 28.8 23.1 31.5
2015 21.2 19.3 17.5 15.7 14.6 18.5 20.7 43.4 29.2 28.2 24.8 20.3 22.8
2016 18.8 17.1 15.5 14.0 12.6 13.6 39.1 48.7 58.3 35.8 29.9 23.9 27.3
2017 21.0 19.5 17.7 16.0 14.5 14.5 29.3 36.6 51.0 33.2 29.9 23.8 25.6
2018 21.1 19.6 17.8 16.1 14.8 33.8 40.6 58.9 41.9 64.7 35.2 26.6 32.7
2019 23.2 21.3 19.4 17.4 15.6
2020
Ave 21.7 19.8 18.0 16.2 16.4 21.2 36.5 62.0 63.8 42.3 31.9 24.6 31.2
RGS1 (CA=3062km2) tank model monthly and annual discharge, unit: m3/s

プログラム上の処理

Pythonプログラム上での出力したい表の処理は以下の手順で行う。

  • １．pandasのデータフレーム作成
  • ２．一度データフレームをそのままエクセルに保存
  • ３．保存したエクセルファイルを再度pythonで呼び出し、エクセうの書式指定、罫線などの可能を行い、再保存する。

データフレーム作成

    df = pd.DataFrame({ 'Year' : lyy,
                        'Jan' : qqm[0:n,0],
                        'Feb' : qqm[0:n,1],
                        'Mar' : qqm[0:n,2],
                        'Apr' : qqm[0:n,3],
                        'May' : qqm[0:n,4],
                        'Jun' : qqm[0:n,5],
                        'Jul' : qqm[0:n,6],
                        'Aug' : qqm[0:n,7],
                        'Sep' : qqm[0:n,8],
                        'Oct' : qqm[0:n,9],
                        'Nov' : qqm[0:n,10],
                        'Dec' : qqm[0:n,11],
                        'Ave' : qqm[0:n,12]
                      })
    df = df.set_index('Year')

データフレーム保存

関連する内容を、１ファイルに、複数シートに分割して保存する。

    fnameW='out_dis_month.xlsx'
    with pd.ExcelWriter(fnameW) as writer:
        df1_tank.to_excel(writer, sheet_name='RGS1_tank')
        df2_tank.to_excel(writer, sheet_name='RGS2_tank')
        df1_mlp.to_excel(writer, sheet_name='RGS1_mlp')
        df2_mlp.to_excel(writer, sheet_name='RGS2_mlp')

エクセルファイルの再読み込みと保存

この間に処理を記載する。

    wb =openpyxl.load_workbook(fnameW)
    #
    ....
    #
    wb.save(fnameW)

フォーマット指定

    # format
    for i in range(2, n+1):
        for j in range(2, m+1):
            ws.cell(row=i, column=j).number_format = '0.0'

罫線（前半は１重線、後半は２重線を指定）

from openpyxl.styles.borders import Border, Side

    # border line
    bcc = Border(top=Side(style='thin', color='000000'),
                 bottom=Side(style='thin', color='000000'),
                 left=Side(style='thin', color='000000'),
                 right=Side(style='thin', color='000000')
                 )
    for i in range(1, n+1):
        for j in range(1, m+1):
            ws.cell(row=i, column=j).border = bcc
    bcc = Border(top=Side(style='double', color='000000'),
                 bottom=Side(style='double', color='000000'),
                 left=Side(style='double', color='000000'),
                 right=Side(style='double', color='000000')
                 )
    for j in range(1, m+1):
        ws.cell(row=n, column=j).border = bcc
    for i in range(1, n+1):
        ws.cell(row=i, column=m).border = bcc

セルの結合（長い文字列を入力しセルを結合する）

    ws.cell(row=n+1,column=1).value='RGS2 (CA=7860km2) tank model monthly and annual discharge, unit: m3/s'
    ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss)

セルの塗りつぶし

このプログラムには含まれていないが、塗りつぶしを行う場合は以下の要領で行う。

from openpyxl.styles import PatternFill

# cell color
fill = PatternFill(patternType='solid', fgColor='00ffff')
for i in [6,11,16,21,26,31]:
    for j in range(3,17):
        ws.cell(row=i,column=j).fill = fill

主要計算部分の説明

例外処理を使っているので、その部分を解説。 このプログラムでは、１日でも欠測があれば、その月および年の平均は欠測としている。

mainで数表にする年と月を、リスト lyy と lmm で定義している。

    lyy=['2001','2002','2003','2004','2005','2006','2007','2008','2009','2010',
         '2011','2012','2013','2014','2015','2016','2017','2018','2019','2020']
    lmm=['01','02','03','04','05','06','07','08','09','10','11','12']

以下の関数で各月および各年の平均値を算出している。 引数の意味は以下の通り。

  • rf: 日データを示すseries.
  • lyyy: リストlyyと同じ
  • lmmm: リスト lmm と同じ。
def cal_m(rf,lyyy,lmmm):
    qqm=np.zeros((len(lyyy)+1,len(lmmm)+1),dtype=np.float64)
    kda=np.zeros((len(lyyy)+1,len(lmmm)+1),dtype=np.float64)
    nac=0
    for i,yy in enumerate(lyyy):
        for j,mm in enumerate(lmmm):
            ss=yy+'/'+mm # 処理する年/月を指定
            try:
                # Series: rf[ss]に対する処理（合計およびデータ数カウント）を行う
                qqm[i,j]=rf[ss].sum() # sum in a month
                kda[i,j]=rf[ss].count() # availavle days
                nac=np.count_nonzero(np.isnan(rf[ss])) # unavailable days
            except KeyError:
                # もし[ss]に合致する年/月のデータがなければ配列にnp.nanを代入
                qqm[i,j]=np.nan
                kda[i,j]=np.nan
            if 0<nac:
                # 仮に指定した年/月のデータがあってもnp.nanがあれば、すなわち月の全日のデータが揃っていなければnp.nanを代入
                qqm[i,j]=np.nan
                kda[i,j]=np.nan
        qqm[i,-1]=np.sum(qqm[i,:]) # sum in a year 年間合計値
        kda[i,-1]=np.sum(kda[i,:]) # sum in a year of available days 年間合計値
    for j in range(len(lmmm)+1):
        qqm[-1,j]=np.nansum(qqm[:,j]) # nanを除外したデータの合計（各月合計値）
        kda[-1,j]=np.nansum(kda[:,j]) # nanを除外したデータ数の合計（各月合計値）
    qmean=qqm/kda # 月および年のデータ合計をデータ数で除すことにより平均値を算出
    return qmean

プログラム例（全文）

プログラムは、ファイルから値を読み取り、上に示した表を作るもの。 実際のプログラムの出力はエクセル。 以下にPythonによるプログラム全文を示す。

import pandas as pd
import numpy as np
import datetime
import openpyxl
from openpyxl.styles.borders import Border, Side


def formx(ws,n,m):
    # format
    for i in range(2, n+1):
        for j in range(2, m+1):
            ws.cell(row=i, column=j).number_format = '0.0'
    # border line
    bcc = Border(top=Side(style='thin', color='000000'),
                 bottom=Side(style='thin', color='000000'),
                 left=Side(style='thin', color='000000'),
                 right=Side(style='thin', color='000000')
                 )
    for i in range(1, n+1):
        for j in range(1, m+1):
            ws.cell(row=i, column=j).border = bcc
    bcc = Border(top=Side(style='double', color='000000'),
                 bottom=Side(style='double', color='000000'),
                 left=Side(style='double', color='000000'),
                 right=Side(style='double', color='000000')
                 )
    for j in range(1, m+1):
        ws.cell(row=n, column=j).border = bcc
    for i in range(1, n+1):
        ws.cell(row=i, column=m).border = bcc


            
def wexcel(df1_tank,df2_tank,df1_mlp,df2_mlp):
    fnameW='out_dis_month.xlsx'
    with pd.ExcelWriter(fnameW) as writer:
        df1_tank.to_excel(writer, sheet_name='RGS1_tank')
        df2_tank.to_excel(writer, sheet_name='RGS2_tank')
        df1_mlp.to_excel(writer, sheet_name='RGS1_mlp')
        df2_mlp.to_excel(writer, sheet_name='RGS2_mlp')


    n=len(df1_tank.index)+1
    m=len(df1_tank.columns)+1
    wb =openpyxl.load_workbook(fnameW)
    #
    ws = wb.get_sheet_by_name('RGS1_tank')
    formx(ws,n,m)
    ws.cell(row=n+1,column=1).value='RGS1 (CA=3062km2) tank model monthly and annual discharge, unit: m3/s'
    ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss)
    #
    ws = wb.get_sheet_by_name('RGS2_tank')
    formx(ws,n,m)
    ws.cell(row=n+1,column=1).value='RGS2 (CA=7860km2) tank model monthly and annual discharge, unit: m3/s'
    ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss)
    #
    ws = wb.get_sheet_by_name('RGS1_mlp')
    formx(ws,n,m)
    ws.cell(row=n+1,column=1).value='RGS1 (CA=3062km2) MLP monthly and annual discharge, unit: m3/s'
    ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss)
    #
    ws = wb.get_sheet_by_name('RGS2_mlp')
    formx(ws,n,m)
    ws.cell(row=n+1,column=1).value='RGS2 (CA=7860km2) MLP monthly and annual discharge, unit: m3/s'
    ss='A{0}:N{1}'.format(n+1,n+1); ws.merge_cells(ss)
    #
    wb.save(fnameW)


def rdata(fnameR):
    df=pd.read_csv(fnameR, header=0, index_col=0) # read excel data
    df.index = pd.to_datetime(df.index, format='%Y/%m/%d')
    return df    
    

def cal_m(rf,lyyy,lmmm):
    qqm=np.zeros((len(lyyy)+1,len(lmmm)+1),dtype=np.float64)
    kda=np.zeros((len(lyyy)+1,len(lmmm)+1),dtype=np.float64)
    nac=0
    for i,yy in enumerate(lyyy):
        for j,mm in enumerate(lmmm):
            ss=yy+'/'+mm
            try:
                qqm[i,j]=rf[ss].sum() # sum in a month
                kda[i,j]=rf[ss].count() # availavle days
                nac=np.count_nonzero(np.isnan(rf[ss])) # unavailable days
            except KeyError:
                qqm[i,j]=np.nan
                kda[i,j]=np.nan
            if 0<nac:
                qqm[i,j]=np.nan
                kda[i,j]=np.nan
        qqm[i,-1]=np.sum(qqm[i,:]) # sum in a year
        kda[i,-1]=np.sum(kda[i,:]) # sum in a year of available days
    for j in range(len(lmmm)+1):
        qqm[-1,j]=np.nansum(qqm[:,j])
        kda[-1,j]=np.nansum(kda[:,j])
    qmean=qqm/kda
    return qmean


def qq_mon(lyy,lmm,qq):
    qqm=cal_m(qq,lyy,lmm)
    lyy=lyy+['Ave']
    n=len(lyy)
    df = pd.DataFrame({ 'Year' : lyy,
                        'Jan' : qqm[0:n,0],
                        'Feb' : qqm[0:n,1],
                        'Mar' : qqm[0:n,2],
                        'Apr' : qqm[0:n,3],
                        'May' : qqm[0:n,4],
                        'Jun' : qqm[0:n,5],
                        'Jul' : qqm[0:n,6],
                        'Aug' : qqm[0:n,7],
                        'Sep' : qqm[0:n,8],
                        'Oct' : qqm[0:n,9],
                        'Nov' : qqm[0:n,10],
                        'Dec' : qqm[0:n,11],
                        'Ave' : qqm[0:n,12]
                      })
    df = df.set_index('Year')
    return df


def main():
    lyy=['2001','2002','2003','2004','2005','2006','2007','2008','2009','2010',
         '2011','2012','2013','2014','2015','2016','2017','2018','2019','2020']
    lmm=['01','02','03','04','05','06','07','08','09','10','11','12']

    fname_list=[
        'df_rgs1_tank_result.csv',
        'df_rgs2_tank_result.csv',
        'df_rgs1_mlp_result.csv',
        'df_rgs2_mlp_result.csv'
    ]
    for iii,fnameR in enumerate(fname_list):
        df0=rdata(fnameR) # daily discharge
        df0=df0['2001/01/01':'2019/05/31']
        qq=pd.Series(df0['Q_pred'], index=df0.index)
        df=qq_mon(lyy,lmm,qq)
        print(fnameR)
        if iii==0: df1_tank=df
        if iii==1: df2_tank=df
        if iii==2: df1_mlp=df
        if iii==3: df2_mlp=df

    wexcel(df1_tank,df2_tank,df1_mlp,df2_mlp)


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

以　上

GraphvizでMLPと通常の重回帰分析のイメージ図を作る

作例

GraphvizでMLPと通常重回帰分析のイメージ図（conceptual diagram あるいは schematic diagram かな）を作ってみた。 作例は以下の通り。

f:id:damyarou:20200816080901p:plain:w600

f:id:damyarou:20200816080923p:plain:w300

プログラミング環境は以下の通り。

Graphvizのインストール

まずはGraphvizのインストール。 brewとpipで２回行う。

brew install graphviz

pip install graphviz

MLPのイメージ図作成プログラム

node の表示順を制御するため、gi.node('1','x[0]' ,pos='0,0!') というように pos を追記している。

from graphviz import Digraph

g = Digraph(format='png')
g.attr(rankdir='LR')
g.attr(splines='false')
g.attr(dpi='300')

# Input subgraph
gi=Digraph(name='cluster_i')
gi.attr(label='Input')
gi.attr(penwidth='0')
gi.node('1','x[0]'  ,pos='0,0!')
gi.node('2','x[1]'  ,pos='0,1!')
gi.node('3','x[..]' ,pos='0,2!')
gi.node('4','x[n]',pos='0,3!')
# Hidden subgraph 1
gh1=Digraph(name='cluster_h1')
gh1.attr(label='Hidden layer 1')
gh1.attr(penwidth='0')
gh1.node('5', 'h1[0]'  ,pos='0,0!')
gh1.node('6', 'h1[1]'  ,pos='0,1!')
gh1.node('7', 'h1[..]' ,pos='0,2!')
gh1.node('8', 'h1[n]',pos='0,3!')
# Hidden subgraph 2
gh2=Digraph(name='cluster_h2')
gh2.attr(label='Hidden layer 2')
gh2.attr(penwidth='0')
gh2.node('9', 'h2[0]'  ,pos='0,0!')
gh2.node('10','h2[1]'  ,pos='0,1!')
gh2.node('11','h2[..]' ,pos='0,2!')
gh2.node('12','h2[n]',pos='0,3!')
# Output subgraph
go=Digraph(name='cluster_o')
go.attr(label='Output')
go.attr(penwidth='0')
go.node('13','y',pos='0,1.5!')

g.subgraph(gi)
g.subgraph(gh1)
g.subgraph(gh2)
g.subgraph(go)

g.edge('1', '5')
g.edge('1', '6')
g.edge('1', '7')
g.edge('1', '8')
g.edge('2', '5')
g.edge('2', '6')
g.edge('2', '7')
g.edge('2', '8')
g.edge('3', '5')
g.edge('3', '6')
g.edge('3', '7')
g.edge('3', '8')
g.edge('4', '5')
g.edge('4', '6')
g.edge('4', '7')
g.edge('4', '8')

g.edge('5', '9')
g.edge('5', '10')
g.edge('5', '11')
g.edge('5', '12')
g.edge('6', '9')
g.edge('6', '10')
g.edge('6', '11')
g.edge('6', '12')
g.edge('7', '9')
g.edge('7', '10')
g.edge('7', '11')
g.edge('7', '12')
g.edge('8', '9')
g.edge('8', '10')
g.edge('8', '11')
g.edge('8', '12')

g.edge('9', '13')
g.edge('10', '13')
g.edge('11', '13')
g.edge('12', '13')


g.render('fig_3_mlp', view=True)

重回帰のイメージ図作成プログラム

edge（矢印線）に回帰係数を示す w[..] を表示するため、g.edge('3', '5', label='w[..]') というように label を追記している。

from graphviz import Digraph

g = Digraph(format='png')
g.attr(rankdir='LR')
g.attr(splines='false')
g.attr(dpi='300')

# Input subgraph
gi=Digraph(name='cluster_i')
gi.attr(label='Input')
gi.attr(penwidth='0')
gi.node('1','x[0]'  ,pos='0,0!')
gi.node('2','x[1]'  ,pos='0,1!')
gi.node('3','x[..]' ,pos='0,2!')
gi.node('4','x[n]',pos='0,3!')
# Output subgraph
go=Digraph(name='cluster_o')
go.attr(label='Output')
go.attr(penwidth='0')
go.node('5','y',pos='0,0!')

g.subgraph(gi)
g.subgraph(go)

g.edge('1', '5', label='w[0]')
g.edge('2', '5', label='w[1]')
g.edge('3', '5', label='w[..]')
g.edge('4', '5', label='w[n]')

g.render('fig_3_reg', view=True)

以　上

Structural design formulas for penstocks embedded in rock

The structural design formulas for penstocks embedded in rock which are shown in 'Technical standards for gates and penstock' in Japan are described.

Allowable stresss of steel material

\begin{equation} \sigma_a = \min(\sigma_Y\: /\: SF_Y,\; \sigma_B\: /\: SF_B) \end{equation}
\begin{align} &\sigma_a & &\text{allowable stress of steel plate}\\ &\sigma_Y & &\text{yield strength of steel plate}\\ &\sigma_B & &\text{tensile strehgth of steel plate}\\ &SF_Y & &\text{safety factor for yield strength (= 1.80)}\\ &SF_B & &\text{safety factor for tensile strength (= 2.35)} \end{align}

Welded joint efficiency

Location of welding RT or UT carrying out
for more than 5% of
welded joint length 		RT or UT carrying out
for less than 5% of
welded joint length
Factory welding 95% 85%
Site welding 90% 80%

Required safety factor for buckling due to external pressure

Pipe shell

\begin{equation} \cfrac{p_k}{P_o} \geq 1.5 \end{equation}
\begin{align} &p_k & &\text{critical buckling pressure of steel pipe}\\ &P_o & &\text{external pressure}\\ \end{align}

Stiffener

\begin{equation} \cfrac{\sigma_{cr}}{\sigma_c} \geq 1.5 \end{equation}
\begin{align} &\sigma_{cr} & &\text{critical buckling stress of stiffener}\\ &\sigma_c & &\text{average compressive stress of stiffener}\\ \end{align}

Design internal pressure

f:id:damyarou:20200625160211j:plain

Definition of dimensions

f:id:damyarou:20200625160239j:plain

\begin{align} &D_0 & &\text{design internal diameter}\\ &D_0' & &\text{design external diameter}\\ &D & &\text{internal diameter adding 1/2 the corrosion allowance}\\ & & &\text{from the internal surface of the pipe shell ($D=D_0 + \epsilon$})\\ &D' & &\text{external diameter subtracting 1/2 the corrosion allowance}\\ & & &\text{from the external surface of the pipe shell ($D'=D_0' - \epsilon$)}\\ &D_m & &\text{diameter of the center of plate thickness ($D_m=2 r_m$)}\\ &r_m & &\text{radius of the center of plate thickness $(r_m=(D_0 + t_0)\:/\:2$)}\\ &t_0 & &\text{design shell thickness}\\ &t & &\text{shell thickness excluding corrosion allowance ($t=t_0 - \epsilon$)}\\ &\epsilon & &\text{allowance thickness for corrosion and wear}\\ \end{align}

Corrosion allowance

An allowance of more than or equal to 1.5mm shall be provided against corrossion and wear.

Minimum shell thickness

The minimum shell thickness used for the pressure lining part shall be more than or equal to those determined from following formula, if stiffeners are not used.

\begin{equation} t_0=\cfrac{D_0 + 800}{400} \end{equation}

The minimum shell thickness shall not be less than 6mm even if the pipe diameter is small and stiffeners are used.

Stress calculation

Tensile stress in circumferential direction due to internal pressure

General

\begin{equation} \sigma=\cfrac{P\cdot D}{2\cdot t} \quad \text{or} \quad \sigma=\cfrac{P\cdot D}{2\cdot (t_0 - \epsilon)} \end{equation}

Internal pressure shared design by bedrock

\begin{equation} \sigma=\cfrac{P\cdot D}{2\cdot t}\cdot (1-\lambda) \end{equation}
\begin{equation} \lambda=\cfrac{1-\cfrac{E_s}{P}\cdot \alpha_s\cdot \Delta T\cdot \cfrac{2\cdot t}{D}} {1+(1+\beta_c)\cdot \cfrac{E_s}{E_c}\cdot \cfrac{2\cdot t}{D}\cdot \log_e\cfrac{D_R}{D} +(1+\beta_g)\cdot \cfrac{E_s}{E_g}\cdot \cfrac{m_g+1}{m_g}\cdot \cfrac{2\cdot t}{D}} \end{equation}
\begin{align} &\sigma & & \text{tensile stress in circumferential direction of pipe}\\ &P & & \text{internal pressure}\\ &\lambda & & \text{sharing ratio of internal pressure by bedrock}\\ &E_s & & \text{elastic modulus of steel}\\ &\alpha_s & & \text{coefficient of linear expansion of steel}\\ &\Delta T & & \text{temperature change of steel penstock}\\ &\beta_c & & \text{coefficient of plastic deformation of concrete}\\ &E_c & & \text{elastic modulus of concrete}\\ &D_R & & \text{excavation diameter of tunnel}\\ &\beta_g & & \text{coefficient of plastic deformation of bedrock}\\ &E_g & & \text{elastic modulus of bedrock}\\ &m_g & & \text{Poisson's number of bedrock} \end{align}
\begin{equation} t=\cfrac{P\cdot D}{2\cdot \eta\cdot \sigma_a} -\cfrac{\eta\cdot \sigma_a-E_s\cdot \alpha_s\cdot \Delta T} {\eta\cdot \sigma_a\cdot \left\{(1+\beta_c)\cdot \cfrac{E_s}{E_c}\cdot \cfrac{2}{D}\cdot \log_e\cfrac{D_R}{D} +(1+\beta_g)\cdot \cfrac{E_s}{E_g}\cdot \cfrac{m_g+1}{m_g}\cdot \cfrac{2}{D}\right\}} \end{equation}
\begin{align} &\sigma_a & &\text{allowable stress of steel}\\ &\eta & &\text{welded joint efficiency} \end{align}

Tensile stress in penstock axis direction

Stress due to restraint by stiffener

\begin{equation} \sigma=1.82\cdot \cfrac{t_r\cdot h_r}{A_r}\cdot \cfrac{P\cdot D}{2\cdot t}\cdot (1-\lambda) \end{equation}
\begin{equation} A_r=t_r\cdot h_r+\left(1.56\cdot \sqrt{r_m\cdot t}+t_r\right)\cdot t \end{equation}

f:id:damyarou:20200625160336j:plain

Stress due to temperature change

\begin{equation} \sigma=\alpha_s\cdot E_s\cdot \Delta T \end{equation}
\begin{align} &\sigma & & \text{stress due to temperature change}\\ &\alpha_s & & \text{coefficient of linear expansion of steel}\\ &E_s & & \text{elastic modulus of steel}\\ &\Delta T & & \text{temperature change} \end{align}

Stress due to Poisson's effect

\begin{equation} \sigma=\nu_s\cdot \sigma_r \end{equation}
\begin{align} &\sigma & & \text{stress due to Poisson's effect}\\ &\nu_s & & \text{Poisson's ratio of steel}\\ &\sigma_r & & \text{circumferential stress} \end{align}

Calculation for external presure

Without stiffener (E. Amstutz's formula)

\begin{equation} p_k=\cfrac{\sigma_N}{\cfrac{r_m}{t}\left(1+0.35\cdot \cfrac{r_m}{t}\cdot \cfrac{\sigma_F^*-\sigma_N}{E_s^*}\right)} \qquad 35 < \cfrac{r_m}{t} \end{equation}
\begin{equation} \left(\cfrac{k_0}{r_m}+\cfrac{\sigma_N}{E_s^*}\right) \cdot \left(1+12\cdot \cfrac{r_m^2}{t^2}\cdot \cfrac{\sigma_N}{E_s^*}\right)^{1.5} =3.36\cdot \cfrac{r_m}{t}\cdot \cfrac{\sigma_F^*-\sigma_N}{E_s^*} \cdot \left(1-\cfrac{1}{2}\cdot \cfrac{r_m}{t}\cdot \cfrac{\sigma_F^*-\sigma_N}{E_s^*}\right) \end{equation}
\begin{equation} k_0=\cfrac{\left(\alpha_s\cdot \Delta T+\beta_g\cdot \cfrac{\sigma_a\cdot \eta}{E_s}\right)\cdot r_0'}{1+\beta_g} \quad \text{or} \quad k_0=0.4\times 10^{-3}\cdot r_m \end{equation}
\begin{equation} E_s^*=\cfrac{E_s}{1-\nu_s{}^2} \end{equation}
\begin{equation} \sigma_F^*=\mu\cdot \cfrac{\sigma_F}{\sqrt{1-\nu_s+\nu_s{}^2}} \qquad \mu=1.5-0.5\cdot \cfrac{1}{\left(1+0.002\cdot \cfrac{E_s}{\sigma_F}\right)^2} \end{equation}
\begin{align} &p_k & &\text{critical buckling pressure of pipe shell without stiffener}\\ &k_0 & &\text{gap between concrete and external surface of pipe}\\ &\sigma_a & &\text{allowable stress of steel}\\ &\eta & &\text{welded joint efficiency}\\ &\sigma_N & &\text{circumferential direct stress at deformed pipe shell portion}\\ &\sigma_F & &\text{yield point of steel}\\ &E_s & &\text{elastic modulus of steel}\\ &\nu_s & &\text{Poisson's ratio of steel} \end{align}

With stiffener

Pipe shell proper (S. Timoshenko's formula)
\begin{equation} \cfrac{(1-\nu_s{}^2)\cdot r_0'\cdot p_k}{E_s\cdot t} =\cfrac{1-\nu_s{}^2}{(n^2-1)\cdot \left(1+\cfrac{n^2\cdot l^2}{\pi^2\cdot r_0'{}^2}\right)^2} +\cfrac{l^2}{12\cdot r_0'{}^2}\cdot \left\{(n^2-1)+\cfrac{2\cdot n^2-1-\nu_s}{1+\cfrac{n^2\cdot l^2}{\pi^2\cdot r_0'{}^2}}\right\} \end{equation}
\begin{align} &p_k & &\text{critical buckling pressure of pipe shell with stiffener}\\ &n & &\text{number of wrinkles}\\ &l & &\text{actual interval of stiffeners} \end{align}
\begin{equation} l'=\left(l+1.56\cdot \sqrt{r_m\cdot t}\cdot \cos^{-1}\lambda\right)\cdot \left(1+0.037\cdot \cfrac{\sqrt{r_m\cdot t}}{l+1.56\cdot \sqrt{r_m\cdot t}\cdot \cos^{-1}\lambda} \cdot \cfrac{t^3}{I_s} \right) \end{equation}
\begin{equation} \lambda=1-(1+T)\cdot \cfrac{1+\cfrac{t_r}{1.56\cdot \sqrt{r_m\cdot t}}}{1+\cfrac{S_0}{1.56\cdot \sqrt{r_m\cdot t}}} \end{equation}
\begin{equation} T=\cfrac{2\cdot C}{t_r+1.56\cdot \sqrt{r_m\cdot t}} \end{equation}
\begin{equation} C=\cfrac{\cfrac{r_0'{}^2}{t}-\cfrac{(t_r+1.56\cdot \sqrt{r_m\cdot t})\cdot r_0'{}^2}{S_0+1.56\cdot t\cdot \sqrt{r_m\cdot t}}} {\cfrac{3}{\{3\cdot (1-\nu_s{}^2)\}^{0.75}}\cdot \left(\cfrac{r_0'}{t}\right)^{1.5}\cdot \cfrac{\sinh(\beta\cdot l)+\sin(\beta\cdot l)}{\cosh(\beta\cdot l)-\cos(\beta\cdot l)} +\cfrac{2\cdot r_0'{}^2}{S_0+1.56\cdot t\cdot \sqrt{r_m\cdot t}} } \end{equation}
\begin{equation} \beta=\cfrac{\{3\cdot (1-\nu_s{}^2)\}^{0.25}}{\sqrt{r_m\cdot t}} \end{equation}
\begin{equation} S_0=t_r\cdot (t+h_r) \qquad I_s=\cfrac{t_r}{12}\cdot (t+h_r)^3 \end{equation}
\begin{align} &l' & &\text{modified interval of stiffeners by Nagashima and Kozuki}\\ &S_0 & &\text{section area of stiffener}\\ &I_s & &\text{moment of inertia of stiffener} \end{align}

f:id:damyarou:20200625160458j:plain

Stiffener (E. Amstutz's formula)
\begin{equation} \sigma_{cr}=\sigma_N\cdot \left\{1-\cfrac{r_0'}{e}\cdot \cfrac{\sigma_F-\sigma_N} {\left(1+\cfrac{3}{2}\cdot \pi\right)\cdot E_s}\right\} \end{equation}
\begin{equation} \left(\cfrac{k_0}{r_m}+\cfrac{\sigma_N}{E_s}\right)\cdot \left(1+\cfrac{r_m{}^2}{i^2}\cdot \cfrac{\sigma_N}{E_s}\right)^{1.5} =1.68\cdot \cfrac{r_m}{e}\cdot \cfrac{\sigma_F-\sigma_N}{E_s}\cdot \left(1-\cfrac{1}{4}\cdot \cfrac{r_m}{e}\cdot \cfrac{\sigma_F-\sigma_N}{E_s}\right) \end{equation}
\begin{align} &\sigma_{cr} & &\text{critical buckling stress of stiffener}\\ &i & &\text{radius of gyration of combined section of stiffeners}\\ &e & &\text{distance from the center of gravity of a combined section}\\ & & &\text{of stiffener to internal surface of the pipe} \end{align}
\begin{equation} \sigma_c=\cfrac{p'\cdot r_0'\cdot (t_r+1.56\cdot \sqrt{r_m\cdot t})} {S_0+1.56\cdot t\cdot \sqrt{r_m\cdot t}} \end{equation}
\begin{equation} p'=\cfrac{p}{t_r+1.56\cdot \sqrt{r_m\cdot t}}\cdot \left\{(t_r+1.56\cdot \sqrt{r_m\cdot t})+2\cdot C\right\} \end{equation}
\begin{align} &\sigma_c & &\text{average compressive stress of stiffener}\\ &p & &\text{external pressure}\\ &p' & &\text{modified external pressure summing up shearing force}\\ & & &\text{acting both sides of a stiffener from a pipe shell}\\ &C & &\text{coefficient defined in Nagashima and Kozuki's formula} \end{align}
\begin{equation} S_f=\cfrac{\sigma_{cr}}{\sigma_c} \qquad \text{(safety factor of stiffener)} \end{equation}

f:id:damyarou:20200625160527j:plain

That's all.

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

はじめに

岩盤内埋設式水圧鉄管プログラムは、過去に一度掲載している（以下参照）。今回のプログラミング内容は工学的には過去のものとダブルが、過去投稿も、エクセルの罫線の書き方やセルの色付などの参考になるため、削除せずにおく。

今回は、以前のプログラムに対し、更に汎用性をもたせるよう若干の変更を加えるとともに、エクセル書出部があまりにも冗長なので、プログラムを「設計計算部」と「エクセル書出部」に分離した。

「設計計算部」ではデータを読み込み、必要諸元を定めたらデータフレームを作成し、それをエクセル形式で出力するまでとした。具体的には、１つの出力用temporaryファイル（Excel）に、３つのシートを作成し、入力データ、耐内圧設計結果、耐外圧設計結果を格納したデータフレームを保存している。

with pd.ExcelWriter('_df_ps.xlsx') as writer:                # temporary file of outpit
        df0.to_excel(writer, sheet_name='Load', index=False) # dataframe for input data
        df1.to_excel(writer, sheet_name='Pin', index=False)  # dataframe for design against internal pressure
        df2.to_excel(writer, sheet_name='Pex', index=False)  # dataframe for design against external pressure

この段階では出力数値のフォーマットが統一されておらず、きれいとは言えない。

このため、「エクセル書出部」で、「設計計算部」で吐き出したエクセル形式のデータフレームを読み込み、数値の書式を揃えてエクセルに書き出すようにした。この段階で、等幅フォントの使用、セル幅の調整、各セルを破線で囲む、というところまで行っている。これ以上の装飾はプログラムが煩雑になるとともに汎用性をもたせるのが難しいため、以降はエクセルから装飾するという考え方にしている。

設計計算式は、日本の水門鉄管技術基準に基づいている。なお、本プログラムでは、管軸方向応力の検討は行っていない。

出力データ事例（file: out_penstock.xlsx）

入力データ出力（sheet: Load）

f:id:damyarou:20200623173158p:plain

耐内圧設計結果（sheet: Pin）

f:id:damyarou:20200623173217p:plain

耐外圧設計結果（sheet: Pout）

f:id:damyarou:20200623173235p:plain

プログラミング

入出力データファイル

プログラム内で指定している入出力ファイルは以下の通り。 設計計算部の入力データファイル「inp_load.xlsx」の内容は、エクセル書出部出力の入力データ出力と同一（ただし記号の凡例は含まず）である。

プログラム入力データファイル出力データファイル
設計計算部inp_load.xlsx_df_ps.xlsx
エクセル書出部_df_ps.xlsxout_penstock.xlsx

エクセル書き出し部の出力は３つのシートからなっており、それぞれ以下の内容が出力されている。

Sheet内容
Load入力データ（板割り・荷重計算・岩盤物性）
Pin耐内圧設計結果
Pex耐外圧設計結果

入出力項目

以下で示す入力データおよび出力データは、考えるべき「ある長さを有する鉄管」の「下流端」を示すことに注意する。

入力データ項目
列名説明
No通し番号
L区間鉄管長さ
Sum(L)鉄管累計長さ（サージタンクからの長さ）
EL鉄管中心標高
D0鉄管設計内径
GWL鉄管中心に対する地表標高
Hst静水頭
Hsgサージング水頭（一定値）
Hwh水撃圧水頭（水車中心で最大、サージタンクでゼロ）
Hin設計内水頭（Hst + Hsg + Hwh）
Hex設計外水頭（GWL - EL）
Drトンネル掘削径（岩盤負担設計用）
Eg岩盤弾性係数（岩盤負担設計用）
Remarks部位説明


耐内圧計算結果出力項目
列名説明
No通し番号
Pi設計内水圧（設計内水頭 x 0.01）
L区間鉄管長さ
D0鉄管設計内径
t0鉄管設計板厚
steel鋼種
Eg岩盤弾性係数
lam岩盤負担率
sig鉄管発生応力
siga鉄管許容応力
weight区間鉄管管胴重量
Remarks部位説明


耐外圧計算結果出力項目
列名説明
No通し番号
Pe設計外水圧（設計外水頭 x 0.01）
L区間鉄管長さ
D0鉄管設計内径
t0鉄管設計板厚
steel鋼種
pk_0限界座屈圧力（補剛材無し）
SF_0外水圧に対する安全率（補剛材無し）
pk_s限界座屈圧力（補剛材有り）
SF_s外水圧に対する安全率（補剛材有り）
sc_r補剛材の限界座屈応力
sc_c補剛材の合成圧縮応力
SF_c補剛材の合成圧縮応力に対する安全率
stiffener補剛材ピッチ


プログラム内で指定している変数

設計計算に必要な定数は、入力データの他、以下のものをプログラム内で定義している。もちろん値の変更は可能である。

Global変数として値を定義している変数
変数説明
E_steel = 206000鉄管弾性係数（MPa）
po_steel = 0.3鉄管ポアソン比
t_eps = 2.0鉄管余裕厚（mm）
ef = 0.85 溶接継手効率
t0_rec = 25推奨最小板厚

(note) 推奨最小板厚 t0_rec は、水門鉄管技術基準で定める最小板厚とは別に、施工性・耐外圧設計のために独自に定める最小板厚である。計算される板厚はこの値以上となる。計算された補剛材ピッチが小さすぎる場合など耐外圧設計において管胴本体の耐外圧性能を高めたい場合、溶接施工用台車吊り下げのための必要最小板厚を指定したい場合に用いる。

岩盤負担設計用として値を定義している変数
変数説明
alpha = 1.2e-5鉄管の熱膨張係数（１／度）
deltaT = 20.0鉄管の温度降下量（度）
Ec = 20600.0充填コンクリートの弾性係数（MPa）
Bc = 0.0充填コンクリートの塑性変形係数
Bg = 1.0岩盤の塑性変形係数
mg = 4.0岩盤のポアソン数（岩盤のポアソン比 = 0.25）

(note) 鉄管と重点コンクリート間の空隙は、k_0=0.4\times 10^{-3}\cdot r_m で計算している。

設計用鋼種と許容応力・降伏点
鋼種適用板厚許容応力降伏点
JIS: SM400t0 < 40mm130 MPa235 MPa
JIS: SM490t0 < 40mm175 MPa315 MPa
JIS: SM570t0 < 40mm240 MPa450 MPa
JIS: SM57040mm < t0235 MPa430 MPa

(note) 設計計算では耐内圧設計において上記鋼種を自動的に選定する。SM400、SM490では許容応力が40mmで切り替わることを考慮し板厚は40mm以下としている。このプログラム内では使用鋼種はJIS: SM570までとしているが、JIS: SHY685を加えることにより、更に高落差への対応が可能である。

補剛材諸元
変数説明
hr=75補剛材高さ（mm）
tr=20補剛材板厚（mm）
el=3000, 1500, 1000補剛材ピッチ（mm）

(note) 補剛材設計では、補剛材寸法を、hr x tr = 75 x 20 と設定し、補剛材ピッチを el=3000, 1500, 1000と変化させて、外水圧に対する安全率が、管胴・補剛材とも1.5以上となるピッチを選定している。所要安全率が満たされない場合、補剛材寸法・ピッチ共に変更可能であるが、Global変数として定義している推奨最小板厚（変数名：t0_rec）を増加させて管胴の外水圧に対する安全率を高める方法もある。なお、補剛材ピッチは3m単位管に設置することを意識しており、また補剛材高さはトンネル掘削径との兼ね合いを考慮して設定することに注意する。

プログラムコード

設計計算プログラム

# ================================
# Embedded penstock Design
# ================================
import numpy as np
import pandas as pd
from scipy import optimize
import openpyxl


# declaration of global variables
E_steel = 206000  # elastic mpdulus of steel (MPa)
po_steel = 0.3  # Poisson's ratio of steel
t_eps = 2.0  # margin thickness (mm)
ef = 0.85  # welding joint efficiency
t0_rec = 25  # recommended minimum plate thickness


def stcr(dd0, t0, hr, tr, el, pp, sigF):
    # Calculation of critical stress and compressive stress of stiffener
    # dd0  : design internal diameter of pipe
    # t0   : design plate thickness
    # hr   : height of stiffener
    # tr   : plate thickness of stiffener
    # el   : stiffener interval
    # pp   : external pressure
    # sigF : yield point of steel plate material
    def calrg(dd0, t0, eps, hr, tr):
        # calculation of section properties
        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):
        # function for Brent-method for obtaining sigN
        # 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
    Es = E_steel     # elastic modulus of steel plate
    pos = po_steel    # Poisson's ratio of steel plate
    eps = t_eps       # margin of plate thickness
    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-pos**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-pos**2))**0.75*(r0d/t)**1.5*(np.sinh(beta*el)+np.sin(beta*el)) / \
        (np.cosh(beta*el)-np.cos(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):
    # Calculation of 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
    # dd0 : design internal diameter of steel pipe
    # t0  : design plate thickness of steel pipe
    pi = np.pi
    Es = E_steel      # elastic modulus of steel
    pos = po_steel    # Poisson's ratio of steel
    eps = t_eps       # margin of plate thickness
    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-pos**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-pos**2))**0.75*(r0d/t)**1.5*(np.sinh(beta*el)+np.sin(beta*el)) / \
        (np.cosh(beta*el)-np.cos(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-pos**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-pos) /
                             (1+n**2*ell**2/pi**2/r0d**2))
        _pk = (c1+c2)*Es*t/(1-pos**2)/r0d
        apk = apk+[_pk]
    pk = min(apk)
    return pk


def ams(dd0, t0, sigF):
    # Calculate criticsl buckling pressure
    # of steel pipe without stiffener by Amstutz's formula
    # dd0  : design internal diameter of steel pipe
    # t0   : design plate thickness of steel pipe
    # sigF : yield point of steel material
    def func(x, params):
        # function for Brent-method for obtaining sigN
        # 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):
        # calculation of critical buckling pressure
        aa = rm/t*(1+0.35*rm/t*(sigFs-sigN)/Ess)
        pk = sigN/aa
        return pk
    # main routine
    eps = t_eps         # margin of plate thickness
    Es = E_steel        # elastic modulus of steel plate
    pos = po_steel      # Poisson's ratio of steel plate
    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
    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 calc1(pp, dd0, sa, Eg, dr):
    # Design for internal pressure
    # (Calculate required plate thickness)
    # pp  : internal pressure
    # dd0 : design internal diameter
    # sa  : allowable stress of steel pipe
    # Eg  : elastic modulus of bedrock
    # dr  : excavation diameter of tunnel
    Es = E_steel  # elastic modulus of steel
    eps = t_eps  # margin thickness
    alpha = 1.2e-5  # thermal expansion coefficient of steel
    deltaT = 20.0  # temperatire change of steel
    Ec = 20600.0  # elastic modulus of backfill concrete
    Bc = 0.0  # plastic deformation modulus of cbackfill concrete
    Bg = 1.0  # plastoc deformation modulus of bedrock
    mg = 4.0  # Poisson number of bedrock
    t0t = t0_rec  # 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 1 <= Eg:
        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 datainp():
    df0 = pd.read_excel('inp_load.xlsx')
    return df0


def main():
    df0 = datainp()
    no = np.array(df0['No'], dtype=np.int)              # pipe number
    al = np.array(df0['L(m)'], dtype=np.float64)        # pipe 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)
    dte = np.array(df0['Dr(m)'], dtype=np.float64) * \
        1000  # tunnel excavation diameter (mm)
    # elastic modulus of bedrock (MPa)
    egg = np.array(df0['Eg(MPa)'], dtype=np.float64)

    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, dd0, ppi, ppe, dr, Eg in zip(no, d0, pi, pe, dte, egg):
        sname = 'SM400'
        sa = 130*ef
        t0, sig, lam = calc1(ppi, dd0, sa, Eg, dr)
        if sname == 'SM400' and 40 < t0:
            sname = 'SM490'
            sa = 175*ef
            t0, sig, lam = calc1(ppi, dd0, sa, Eg, dr)
        if sname == 'SM490' and 40 < t0:
            sname = 'SM570'
            sa = 240*ef
            t0, sig, lam = calc1(ppi, dd0, sa, Eg, dr)
        if sname == 'SM570' and 40 < t0:
            sname = 'SM570'
            sa = 235*ef
            t0, sig, lam = calc1(ppi, dd0, sa, Eg, dr)
        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)
        sf0 = pk0/ppe
        if sf0 < 1.5:
            # stiffener size
            hr = 75  # height of stiffener plate
            tr = 20  # thickness of stiffener plate
            # stiffener interval=2000mm
            el = 3000.0
            pk1 = timo(el, hr, tr, dd0, t0)
            sf1 = pk1/ppe
            scr1, sc1 = stcr(dd0, t0, hr, tr, el, ppe, sigF)
            sfc1 = scr1/sc1
            # stiffener interval=1500mm
            el = 1500.0
            pk2 = timo(el, hr, tr, dd0, t0)
            sf2 = pk2/ppe
            scr2, sc2 = stcr(dd0, t0, hr, tr, el, ppe, sigF)
            sfc2 = scr2/sc2
            # stiffener interval=1000mm
            el = 1000.0
            pk3 = timo(el, hr, tr, dd0, t0)
            sf3 = pk3/ppe
            scr3, sc3 = stcr(dd0, t0, 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,
        'Eg(MPa)': egg,
        '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,
        'pk_0(MPa)': a_pk0,
        'SF_0': a_sf0,
        'pk_s(MPa)': a_pks,
        'SF_s': a_sfs,
        'sc_r(MPa)': a_scr,
        'sc_c(MPa)': a_scc,
        'SF_c': a_sfc,
        'stiffener': a_stiff
    })
    with pd.ExcelWriter('_df_ps.xlsx') as writer:
        df0.to_excel(writer, sheet_name='Load', index=False)
        df1.to_excel(writer, sheet_name='Pin', index=False)
        df2.to_excel(writer, sheet_name='Pex', index=False)


if __name__ == '__main__':
    main()

エクセル書出プログラム

# ================================
# Makeing formated excel file
# ================================
import numpy as np
import pandas as pd
import openpyxl
from openpyxl.styles.borders import Border, Side
from openpyxl.styles import PatternFill
from openpyxl.styles.fonts import Font


def xlline(ws, row_m, col_m):
    # border line
    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')
                 )
    for i in range(0, row_m):
        for j in range(0, col_m):
            ws.cell(row=i+1, column=j+1).border = bcc


def xlswrite(fnameW, df0, df1, df2):
    # Write results to Excel
    wb = openpyxl.Workbook()
    row_height = 20
    font = Font(name='Courier New',
                size=12,
                bold=False,
                italic=False,
                underline='none',
                strike=False,
                color='000000')
    #
    #
    ws = wb.create_sheet(index=1, title='Load')
    df = df0
    pcol = [
        ['No', '0', 'A', 12],
        ['L(m)', '0.000', 'B', 12],
        ['Sum(L)', '0.000', 'C', 12],
        ['EL(m)', '0.000', 'D', 12],
        ['D0(m)', '0.000', 'E', 12],
        ['GWL(m)', '0.000', 'F', 12],
        ['Hst(m)', '0.000', 'G', 12],
        ['Hsg(m)', '0.000', 'H', 12],
        ['Hwh(m)', '0.000', 'I', 12],
        ['Hin(m)', '0.000', 'J', 12],
        ['Hex(m)', '0.000', 'K', 12],
        ['Dr(m)', '0.000', 'L', 12],
        ['Eg(MPa)', '0', 'M', 12],
        ['Remarks', '', 'N', 30]
    ]
    for j in range(0, len(pcol)):
        ws.cell(row=1, column=j+1).value = pcol[j][0]
        ws.cell(row=1, column=j+1).font = font
    for j in range(0, len(pcol)):
        temp = np.array(df[pcol[j][0]])
        for i in range(len(temp)):
            ws.cell(row=i+2, column=j+1).number_format = pcol[j][1]
            ws.cell(row=i+2, column=j+1).value = temp[i]
            ws.cell(row=i+2, column=j+1).font = font
    # legend
    snote = [
        'No: section number',
        'L: length of pipe section',
        'Sum(L): cumulative length',
        'EL: elevation of pipe section',
        'D0: internal diameter of pipe section',
        'GWL: ground water level',
        'Hst: hydrostatic head',
        'Hsg: pressure head rise by surging',
        'Hwh: pressure head rise by water hammering',
        'Hin: design internal pressure head',
        'Hex: design external pressure head',
        'Dr: excavation diameter of tunnel',
        'Eg: elastic modulus of bedrock'
    ]
    k = 1+len(temp)+1
    for i, ss in enumerate(snote):
        ws.cell(row=k+i, column=1).value = ss
        ws.cell(row=k+i, column=1).font = font
    # cell height and cell width
    row_max = k+len(snote)
    col_max = len(pcol)
    for i in range(0, row_max):
        ws.row_dimensions[i+1].height = row_height
    for j in range(0, col_max):
        ws.column_dimensions[pcol[j][2]].width = pcol[j][3]
    # line and font
    xlline(ws, row_max-len(snote)-1, col_max)
    #
    #
    ws = wb.create_sheet(index=1, title='Pin')
    df = df1
    pcol = [
        ['No', '0', 'A', 12],
        ['Pi(MPa)', '0.000', 'B', 12],
        ['L(m)', '0.000', 'C', 12],
        ['D0(mm)', '0', 'D', 12],
        ['t0(mm)', '0', 'E', 12],
        ['steel', '', 'F', 12],
        ['Eg(MPa)', '0', 'G', 12],
        ['lam', '0.000', 'H', 12],
        ['sig(MPa)', '0.000', 'I', 12],
        ['siga(MPa)', '0',  'J', 12],
        ['weight(t)', '0.000', 'K', 12],
        ['Remarks', '',  'L', 30]
    ]
    for j in range(0, len(pcol)):
        ws.cell(row=1, column=j+1).value = pcol[j][0]
        ws.cell(row=1, column=j+1).font = font
    for j in range(0, len(pcol)):
        temp = np.array(df[pcol[j][0]])
        for i in range(len(temp)):
            ws.cell(row=i+2, column=j+1).number_format = pcol[j][1]
            ws.cell(row=i+2, column=j+1).value = temp[i]
            ws.cell(row=i+2, column=j+1).font = font
    k = 1+len(temp)+1
    totalw = np.sum(np.array(df['weight(t)']))
    ws.cell(row=k, column=1).number_format = ''
    ws.cell(row=k, column=1).value = 'Sum'
    ws.cell(row=k, column=1).font = font
    ws.cell(row=k, column=11).number_format = '0.000'
    ws.cell(row=k, column=11).value = totalw
    ws.cell(row=k, column=11).font = font
    # legend
    snote = [
        'No: section number',
        'Pi: design internal pressure',
        'L: length of pipe section',
        'D0: internal diameter of pipe section',
        't0: plate thickness of pipe section',
        'steel: kind of steel',
        'Eg: elastic modulus of bedrock',
        'lam: internal pressure shared ratio be bedrock',
        'sig: stress of pipe shell',
        'siga: allowable stress',
        'Weight: weight of pipe section'
    ]
    k = k+1
    for i, ss in enumerate(snote):
        ws.cell(row=k+i, column=1).value = ss
        ws.cell(row=k+i, column=1).font = font
    # cell height and cell width
    row_max = k+len(snote)
    col_max = len(pcol)
    for i in range(0, row_max):
        ws.row_dimensions[i+1].height = row_height
    for j in range(0, col_max):
        ws.column_dimensions[pcol[j][2]].width = pcol[j][3]
    # line
    xlline(ws, row_max-len(snote)-1, col_max)
    #
    #
    ws = wb.create_sheet(index=1, title='Pex')
    df = df2
    pcol = [
        ['No', '0', 'A', 12],
        ['Pe(MPa)', '0.000', 'AB', 12],
        ['L(m)', '0.000', 'C', 12],
        ['D0(mm)', '0', 'D', 12],
        ['t0(mm)', '0', 'E', 12],
        ['steel', '', 'F', 12],
        ['pk_0(MPa)', '0.000', 'G', 12],
        ['SF_0', '0.000', 'H', 12],
        ['pk_s(MPa)', '0.000', 'I', 12],
        ['SF_s', '0.000', 'J', 12],
        ['sc_r(MPa)', '0.000', 'K', 12],
        ['sc_c(MPa)', '0.000', 'L', 12],
        ['SF_c', '0.000', 'M', 12],
        ['stiffener', '', 'N', 30]
    ]
    for j in range(0, len(pcol)):
        ws.cell(row=1, column=j+1).value = pcol[j][0]
        ws.cell(row=1, column=j+1).font = font
    for j in range(0, len(pcol)):
        temp = np.array(df[pcol[j][0]])
        for i in range(len(temp)):
            ws.cell(row=i+2, column=j+1).number_format = pcol[j][1]
            ws.cell(row=i+2, column=j+1).value = temp[i]
            ws.cell(row=i+2, column=j+1).font = font
    # legend
    snote = [
        'No: section number',
        'Pe: design external pressure',
        'L: length of pipe section',
        'D0: internal diameter of pipe section',
        't0: plate thickness of pipe section',
        'steel: kind of steel',
        'pk_0: critical buckling pressure of pipe section without stiffener',
        'SF_0: safety factor against external pressure without pressure',
        'pk_s: critical buckling pressure of pipe section with stiffener',
        'SF_s: safety factor against external pressure with stiffener',
        'sc_r: critical buckling stress of stiffener',
        'sc_c: compressive stress of stiffener',
        'SF_c: safety factor against compressive stress in stiffener',
        'size of stiffener plate: 75mm height x 20mm thickness'
    ]
    k = 1+len(temp)+1
    for i, ss in enumerate(snote):
        ws.cell(row=k+i, column=1).value = ss
        ws.cell(row=k+i, column=1).font = font
    # cell height and cell width
    row_max = k+len(snote)
    col_max = len(pcol)
    for i in range(0, row_max):
        ws.row_dimensions[i+1].height = row_height
    for j in range(0, col_max):
        ws.column_dimensions[pcol[j][2]].width = pcol[j][3]
    # line
    xlline(ws, row_max-len(snote)-1, col_max)
    #
    #
    wb.save(fnameW)


def main():
    fnameR = '_df_ps.xlsx'
    fnameW = 'out_penstock.xlsx'
    sn0 = 'Load'
    sn1 = 'Pin'
    sn2 = 'Pex'
    df0 = pd.read_excel(fnameR, sheet_name=sn0, header=0)
    df1 = pd.read_excel(fnameR, sheet_name=sn1, header=0)
    df2 = pd.read_excel(fnameR, sheet_name=sn2, header=0)
    xlswrite(fnameW, df0, df1, df2)


if __name__ == '__main__':
    main()

以　上