damyarou

python, GMT などのプログラム

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

はじめに

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

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

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

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

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

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鉄管の熱膨張係数(1/度)
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()

以 上