damyarou

python, GMT などのプログラム

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

以　上

Python 軸対称応力解析プログラム

はじめに

軸対称モデルとして扱う構造物として、水平トンネルなど水平方向に回転軸を持つ構造物と、調圧水槽など鉛直方向に回転軸を持つ構造物がある。 これまで、四角形要素の２次元応力解析プログラムを軸対称解析用に書き換えたプログラムを使ってきたが、座標系（座標軸の方向）を意識せず回転軸を鉛直にして計算を行ったところ、荷重に対して変位及び応力が本来あるべき方向と逆方向の解が計算されてきた。 このため、回転軸（ｚ軸）が水平の場合と、回転軸が鉛直の場合で、場合分けして解析するようプログラムの見直しを行った。

はじめから下図Case-2の座標系を反時計回りに90度回転させてｚ軸を上向きとした座標を与えればいいのだが、実際のメッシュ切では、２次元応力解析用のメッシャーで作成した要素座標を使うのが便利なため、今回の対応策を導入することにしたものである。

プログラムの改良

以下に座標系と要素のイメージ図を示す。

Case-1およびCase-2では、回転軸をｚ軸、半径方向軸をｒ軸としている。

f:id:damyarou:20200615110019j:plain

Reference

Reference は通常の２次元応力解析用の座標である。ｘ軸を水平右向き、ｙ軸を鉛直上向きに設定している。 この場合、通常、要素番号は、節点番号を反時計回りに指定して要素を定義する。

Case-1

Case-1は、ｚ軸（回転軸）を水平右向き、ｒ軸（半径方向）を鉛直上向きに設定している。 この場合も、Referenceの２次元応力解析と同様に、要素番号は、節点番号を反時計回りに指定して要素を定義する。（そうするようにプログラムが組まれている）

Case-2

Case-2は、ｚ軸（回転軸）を鉛直上向き、ｒ軸（半径方向）を水平右向きに設定している。 この場合、Reference の２次元応力解析と同様に反時計回りの節点番号で要素を定義するすると、Case-1と比べてｚ軸とｒ軸が入れ替わっているため、要素の節点番号が逆周りで定義されたことになり、剛性マトリックスや物体力ベクトル、温度荷重ベクトルなど要素座標を用いて計算する行列やベクトルが通常と逆符号で計算されることとなる。

対応策

このため、今回整備した軸対称応力解析プログラムでは、変数 nzdir を導入し、座標の方向に応じて、剛性マトリックス、物体力ベクトル、温度荷重ベクトルの方向を逆転させることにした。 すなわち、上図Case-1では、nzdir=1、上図Case-2ではnzdir=-1とし、メインプログラム上で、求まった剛性マトリックス、物体力ベクトル、温度荷重ベクトルの符号を調整するようにしている。

        sm=sm*float(nzdir) # 剛性マトリックス
        tfe=tfe*float(nzdir) # 温度荷重ベクトル
        bfe=bfe*float(nzdir) # 物体力（慣性力）ベクトル

入力データ・フォーマット

npoin  nele  nsec  npfix  nlod  nzdir # Basic values for analysis
E  po  alpha  gamma  gkz              # Material properties
    ..... (1 to nsec) .....           #
node1  node2  node3  node4  isec      # Element connectivity, material set number
    ..... (1 to nele) .....           #   (counter-clockwise order of node number)
z  r  deltaT                          # Coordinate, Temperature change of node
    ..... (1 to npoin) .....          #
node  koz  kor  rdisz  rdisr          # Restricted node and restrict conditions
    ..... (1 to npfix) .....          #
node  fz  fr                          # Loaded node and loading conditions
    ..... (1 to nlod) .....           #   (omit data input if nlod=0)
変数説明
npoin, nele, nsec節点数，要素数，材料特性数
npfix, nlod拘束節点数，載荷節点数
nzdirｚ軸方向（水平右向き: １，鉛直上向き: −１）
E, po, alpha弾性係数，ポアソン比，線膨張係数
gamma, gkz単位体積重量，ｚ方向加速度（gの比）
node1, node2, node3, node4, isec節点1，節点2，節点3，節点4，材料特性番号
z, r, deltaT節点z座標，節点r座標，節点温度変化
node, koz, kor拘束節点番号，z及びr方向拘束の有無（拘束: 1，自由: 0）
rdisz, rdisrz及びr方向変位（無拘束でも0を入力）
node, fz, fr載荷重節点番号，z方向荷重，r方向荷重

（注）回転軸ｚの方向が水平（nzdir=1）であっても鉛直（nzdir=-1）であっても、座標入力は (z , r) の順に行う。

プログラム全文

################################
# Axisymmetric Stress Analysis
################################
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import csr_matrix
import sys
import time


def inpdata_ax4(fnameR,nod,nfree):
    f=open(fnameR,'r')
    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
    nzdir =int(text[5]) # direction of z-axis (=1 or =-1)
    # array declanation
    ae    =np.zeros([5,nsec],dtype=np.float64)      # Section characteristics
    node  =np.zeros([nod+1,nele],dtype=np.int)      # Node-element relationship
    x     =np.zeros([2,npoin],dtype=np.float64)     # Coordinates of nodes
    deltaT=np.zeros(npoin,dtype=np.float64)         # Temperature change of node
    mpfix =np.zeros([nfree,npoin],dtype=np.int)     # Ristrict conditions
    rdis  =np.zeros([nfree,npoin],dtype=np.float64) # Ristricted displacement
    fp    =np.zeros(nfree*npoin,dtype=np.float64)   # External force vector
   # section characteristics
    for i in range(0,nsec):
        text=f.readline()
        text=text.strip()
        text=text.split()
        ae[0,i]=float(text[0]) #E    : Elastic modulus
        ae[1,i]=float(text[1]) #po   : Poisson's ratio
        ae[2,i]=float(text[2]) #alpha: Thermal expansion coefficient
        ae[3,i]=float(text[3]) #gamma: Unit weight of material
        ae[4,i]=float(text[4]) #gkz  : Acceleration in z-direction
    # element-node
    for i in range(0,nele):
        text=f.readline()
        text=text.strip()
        text=text.split()
        node[0,i]=int(text[0]) #node_1
        node[1,i]=int(text[1]) #node_2
        node[2,i]=int(text[2]) #node_3
        node[3,i]=int(text[3]) #node_4
        node[4,i]=int(text[4]) #section characteristic number
    # node coordinates
    for i in range(0,npoin):
        text=f.readline()
        text=text.strip()
        text=text.split()
        x[0,i]   =float(text[0]) # z-coordinate
        x[1,i]   =float(text[1]) # r-coordinate
        deltaT[i]=float(text[2]) # Temperature change of node
    # boundary conditions (0:free, 1:restricted)
    if 0<npfix:
        for i in range(0,npfix):
            text=f.readline()
            text=text.strip()
            text=text.split()
            lp=int(text[0])             #fixed node
            mpfix[0,lp-1]=int(text[1])  #fixed in z-direction
            mpfix[1,lp-1]=int(text[2])  #fixed in r-direction
            rdis[0,lp-1]=float(text[3]) #fixed displacement in z-direction
            rdis[1,lp-1]=float(text[4]) #fixed displacement in r-direction
    # load
    if 0<nlod:
        for i in range(0,nlod):
            text=f.readline()
            text=text.strip()
            text=text.split()
            lp=int(text[0])           #loaded node
            fp[2*lp-2]=float(text[1]) #load in z-direction
            fp[2*lp-1]=float(text[2]) #load in r-direction
    f.close()
    return npoin,nele,nsec,npfix,nlod,nzdir,ae,node,x,deltaT,mpfix,rdis,fp


def prinp_ax4(fnameW,npoin,nele,nsec,npfix,nlod,nzdir,ae,x,fp,deltaT,mpfix,rdis,node):
    fout=open(fnameW,'w')
    # print out of input data
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s}'
    .format('npoin','nele','nsec','npfix','nlod','nzdir'),file=fout)
    print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d}'
    .format(npoin,nele,nsec,npfix,nlod,nzdir),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s}'
    .format('sec','E','po','alpha','gamma','gkz'),file=fout)
    for i in range(0,nsec):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e}'
        .format(i+1,ae[0,i],ae[1,i],ae[2,i],ae[3,i],ae[4,i]),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>5s} {7:>5s}'
    .format('node','z','r','fz','fr','deltaT','koz','kor'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:5d} {7:5d}'
        .format(lp,x[0,i],x[1,i],fp[2*i],fp[2*i+1],deltaT[i],mpfix[0,i],mpfix[1,i]),file=fout)
    if 0<npfix:
        print('{0:>5s} {1:>5s} {2:>5s} {3:>15s} {4:>15s}'
        .format('node','koz','kor','rdis_z','rdis_r'),file=fout)
        for i in range(0,npoin):
            if 0<mpfix[0,i]+mpfix[1,i]:
                lp=i+1
                print('{0:5d} {1:5d} {2:5d} {3:15.7e} {4:15.7e}'
                .format(lp,mpfix[0,i],mpfix[1,i],rdis[0,i],rdis[1,i]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s}'
    .format('elem','i','j','k','l','sec'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d}'
        .format(ne+1,node[0,ne],node[1,ne],node[2,ne],node[3,ne],node[4,ne]),file=fout)
    fout.close()


def prout_ax4(fnameW,npoin,nele,disg,strs):
    fout=open(fnameW,'a')
    # displacement
    print('{0:>5s} {1:>15s} {2:>15s}'.format('node','dis-z','dis-r'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        print('{0:5d} {1:15.7e} {2:15.7e}'.format(lp,disg[2*lp-2],disg[2*lp-1]),file=fout)
    # section force
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>15s} {6:>15s} {7:>15s}'
    .format('elem','sig_z','sig_r','sig_t','tau_zr','p1','p2','ang'),file=fout)
    for ne in range(0,nele):
        sigz =strs[0,ne]
        sigr =strs[1,ne]
        sigt =strs[2,ne]
        tauzr=strs[3,ne]
        ps1,ps2,ang=pst_ax4(sigz,sigr,tauzr)
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:15.7e} {6:15.7e} {7:15.7e}'
        .format(ne+1,sigz,sigr,sigt,tauzr,ps1,ps2,ang),file=fout)
    fout.close()


def dmat_ax4(E,po):
    d=np.zeros([4,4],dtype=np.float64)
    d[0,0]=1.0-po; d[0,1]=po    ; d[0,2]=po    ; d[0,3]=0.0
    d[1,0]=po    ; d[1,1]=1.0-po; d[1,2]=po    ; d[1,3]=0.0
    d[2,0]=po    ; d[2,1]=po    ; d[2,2]=1.0-po; d[2,3]=0.0
    d[3,0]=0.0   ; d[3,1]=0.0   ; d[3,2]=0.0   ; d[3,3]=0.5*(1.0-2.0*po)
    d=d*E/(1.0+po)/(1.0-2.0*po)
    return d


def bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4):
    bm=np.zeros([4,8],dtype=np.float64)
    #[N]
    nm1=0.25*(1.0-a)*(1.0-b)
    nm2=0.25*(1.0+a)*(1.0-b)
    nm3=0.25*(1.0+a)*(1.0+b)
    nm4=0.25*(1.0-a)*(1.0+b)
    rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4
    #dN/da,dN/db
    dn1a=-0.25*(1.0-b)
    dn2a= 0.25*(1.0-b)
    dn3a= 0.25*(1.0+b)
    dn4a=-0.25*(1.0+b)
    dn1b=-0.25*(1.0-a)
    dn2b=-0.25*(1.0+a)
    dn3b= 0.25*(1.0+a)
    dn4b= 0.25*(1.0-a)
    #Jacobi matrix and det(J)
    J11=dn1a*z1 + dn2a*z2 + dn3a*z3 + dn4a*z4
    J12=dn1a*r1 + dn2a*r2 + dn3a*r3 + dn4a*r4
    J21=dn1b*z1 + dn2b*z2 + dn3b*z3 + dn4b*z4
    J22=dn1b*r1 + dn2b*r2 + dn3b*r3 + dn4b*r4
    detJ=J11*J22-J12*J21
    #[B]=[dN/dx][dN/dy]
    bm[0,0]= J22*dn1a-J12*dn1b
    bm[0,2]= J22*dn2a-J12*dn2b
    bm[0,4]= J22*dn3a-J12*dn3b
    bm[0,6]= J22*dn4a-J12*dn4b
    bm[1,1]=-J21*dn1a+J11*dn1b
    bm[1,3]=-J21*dn2a+J11*dn2b
    bm[1,5]=-J21*dn3a+J11*dn3b
    bm[1,7]=-J21*dn4a+J11*dn4b
    bm[2,1]=nm1/rr*detJ
    bm[2,3]=nm2/rr*detJ
    bm[2,5]=nm3/rr*detJ
    bm[2,7]=nm4/rr*detJ
    bm[3,0]=-J21*dn1a+J11*dn1b
    bm[3,1]= J22*dn1a-J12*dn1b
    bm[3,2]=-J21*dn2a+J11*dn2b
    bm[3,3]= J22*dn2a-J12*dn2b
    bm[3,4]=-J21*dn3a+J11*dn3b
    bm[3,5]= J22*dn3a-J12*dn3b
    bm[3,6]=-J21*dn4a+J11*dn4b
    bm[3,7]= J22*dn4a-J12*dn4b
    bm=bm/detJ
    return bm,detJ


def ab_ax4(kk):
    if kk==0:
        a=-0.5773502692
        b=-0.5773502692
    if kk==1:
        a= 0.5773502692
        b=-0.5773502692
    if kk==2:
        a= 0.5773502692
        b= 0.5773502692
    if kk==3:
        a=-0.5773502692
        b= 0.5773502692
    return a,b


def sm_ax4(E,po,z1,r1,z2,r2,z3,r3,z4,r4):
    sm=np.zeros([8,8],dtype=np.float64)
    #Stiffness matrix [sm]=[B]T[D][B]*t*det(J)
    d=dmat_ax4(E,po)
    for kk in range(0,4):
        a,b=ab_ax4(kk)
        bm,detJ=bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4)
        nm1=0.25*(1.0-a)*(1.0-b)
        nm2=0.25*(1.0+a)*(1.0-b)
        nm3=0.25*(1.0+a)*(1.0+b)
        nm4=0.25*(1.0-a)*(1.0+b)
        rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4
        sm=sm+np.dot(bm.T,np.dot(d,bm))*rr*detJ
    return sm


def calst_ax4(E,po,alpha,dt1,dt2,dt3,dt4,wd,z1,r1,z2,r2,z3,r3,z4,r4):
    eps0  =np.zeros(4,dtype=np.float64)
    stress=np.zeros(4,dtype=np.float64)
    #stress vector {stress}=[D][B]{u}
    d=dmat_ax4(E,po)
    for kk in range(0,4):
        a,b=ab_ax4(kk)
        bm,detJ=bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4)
        eps=np.dot(bm,wd)
        # Strain due to temperature change
        nm1=0.25*(1.0-a)*(1.0-b)
        nm2=0.25*(1.0+a)*(1.0-b)
        nm3=0.25*(1.0+a)*(1.0+b)
        nm4=0.25*(1.0-a)*(1.0+b)
        tem=nm1*dt1 + nm2*dt2 + nm3*dt3 + nm4*dt4
        #Thermal strain
        eps0[0]=alpha*tem
        eps0[1]=alpha*tem
        eps0[2]=alpha*tem
        eps0[3]=0.0
        stress=stress+np.dot(d,(eps-eps0))
    stress=0.25*stress
    return stress


def pst_ax4(sigz,sigr,tauzr):
    ps1=0.5*(sigz+sigr)+np.sqrt(0.25*(sigz-sigr)*(sigz-sigr)+tauzr*tauzr)
    ps2=0.5*(sigz+sigr)-np.sqrt(0.25*(sigz-sigr)*(sigz-sigr)+tauzr*tauzr)
    if sigz==sigr:
        if tauzr >0.0: ang= 45.0
        if tauzr <0.0: ang=135.0
        if tauzr==0.0: ang=  0.0
    else:
        ang=0.5*np.arctan(2.0*tauzr/(sigz-sigr))
        ang=180.0/np.pi*ang
        if sigz>sigr and tauzr>=0.0: ang=ang
        if sigz>sigr and tauzr <0.0: ang=ang+180.0
        if sigz<sigr:                ang=ang+90.0
    return ps1,ps2,ang


def tfvec_ax4(E,po,alpha,dt1,dt2,dt3,dt4,z1,r1,z2,r2,z3,r3,z4,r4):
    eps0=np.zeros(4,dtype=np.float64)
    tfe =np.zeros(8,dtype=np.float64)
    # {tfe=[B]T[D]{eps0}
    d=dmat_ax4(E,po)
    for kk in range(0,4):
        a,b=ab_ax4(kk)
        bm,detJ=bmat_ax4(a,b,z1,r1,z2,r2,z3,r3,z4,r4)
        # Strain due to temperature change
        nm1=0.25*(1.0-a)*(1.0-b)
        nm2=0.25*(1.0+a)*(1.0-b)
        nm3=0.25*(1.0+a)*(1.0+b)
        nm4=0.25*(1.0-a)*(1.0+b)
        tem=nm1*dt1 + nm2*dt2 + nm3*dt3 + nm4*dt4
        #Thermal strain
        eps0[0]=alpha*tem
        eps0[1]=alpha*tem
        eps0[2]=alpha*tem
        eps0[3]=0.0
        rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4
        tfe=tfe+np.dot(bm.T,np.dot(d,eps0))*rr*detJ
    return tfe


def bfvec_ax4(gamma,gkz,z1,r1,z2,r2,z3,r3,z4,r4):
    ntn=np.zeros([8,8],dtype=np.float64)
    nm =np.zeros([2,8],dtype=np.float64)
    bfe=np.zeros(8,dtype=np.float64)
    for kk in range(0,4):
        a,b=ab_ax4(kk)
        nm1=0.25*(1.0-a)*(1.0-b)
        nm2=0.25*(1.0+a)*(1.0-b)
        nm3=0.25*(1.0+a)*(1.0+b)
        nm4=0.25*(1.0-a)*(1.0+b)
        rr=nm1*r1 + nm2*r2 + nm3*r3 + nm4*r4
        dn1a=-0.25*(1.0-b)
        dn2a= 0.25*(1.0-b)
        dn3a= 0.25*(1.0+b)
        dn4a=-0.25*(1.0+b)
        dn1b=-0.25*(1.0-a)
        dn2b=-0.25*(1.0+a)
        dn3b= 0.25*(1.0+a)
        dn4b= 0.25*(1.0-a)
        #Jacobi matrix and det(J)
        J11=dn1a*z1 + dn2a*z2 + dn3a*z3 + dn4a*z4
        J12=dn1a*r1 + dn2a*r2 + dn3a*r3 + dn4a*r4
        J21=dn1b*z1 + dn2b*z2 + dn3b*z3 + dn4b*z4
        J22=dn1b*r1 + dn2b*r2 + dn3b*r3 + dn4b*r4
        detJ=J11*J22-J12*J21
        nm[0,0]=nm1
        nm[0,2]=nm2
        nm[0,4]=nm3
        nm[0,6]=nm4
        nm[1,1]=nm[0,0]
        nm[1,3]=nm[0,2]
        nm[1,5]=nm[0,4]
        nm[1,7]=nm[0,6]
        ntn=ntn+np.dot(nm.T,nm)*gamma*rr*detJ
    w=np.array([gkz,0,gkz,0,gkz,0,gkz,0],dtype=np.float64)
    bfe=np.dot(ntn,w)
    return bfe


def main_ax4():
    start=time.time()
    args = sys.argv
    fnameR=args[1] # input data file
    fnameW=args[2] # output data file
    nod=4          # Number of nodes per element
    nfree=2        # Degree of freedom per node
    # data input
    npoin,nele,nsec,npfix,nlod,nzdir,ae,node,x,deltaT,mpfix,rdis,fp=inpdata_ax4(fnameR,nod,nfree)
    # print out of input data
    prinp_ax4(fnameW,npoin,nele,nsec,npfix,nlod,nzdir,ae,x,fp,deltaT,mpfix,rdis,node)
    # array declanation
    ir    =np.zeros(nod*nfree,dtype=np.int)         # Work vector for matrix assembly
    gk    =np.zeros([nfree*npoin,nfree*npoin],dtype=np.float64)         # Global stiffness matrix
    # assembly of stiffness matrix & load vector
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        k=node[2,ne]-1
        l=node[3,ne]-1
        m=node[4,ne]-1
        z1=x[0,i]; r1=x[1,i]
        z2=x[0,j]; r2=x[1,j]
        z3=x[0,k]; r3=x[1,k]
        z4=x[0,l]; r4=x[1,l]
        E     =ae[0,m] #E    : Elastic modulus
        po    =ae[1,m] #po   : Poisson's ratio
        alpha =ae[2,m] #alpha: Thermal expansion coefficient
        gamma =ae[3,m] #gamma: Unit weight of material
        gkz   =ae[4,m] #gkz  : Acceleration in z-direction
        dt1=deltaT[i]; dt2=deltaT[j]; dt3=deltaT[k]; dt4=deltaT[l] # temperature change
        sm=sm_ax4(E,po,z1,r1,z2,r2,z3,r3,z4,r4) # element stiffness matrix
        tfe=tfvec_ax4(E,po,alpha,dt1,dt2,dt3,dt4,z1,r1,z2,r2,z3,r3,z4,r4) # termal load vector
        bfe=bfvec_ax4(gamma,gkz,z1,r1,z2,r2,z3,r3,z4,r4) # body force vector
        sm=sm*float(nzdir)
        tfe=tfe*float(nzdir)
        bfe=bfe*float(nzdir)
        ir[7]=2*l+1; ir[6]=ir[7]-1
        ir[5]=2*k+1; ir[4]=ir[5]-1
        ir[3]=2*j+1; ir[2]=ir[3]-1
        ir[1]=2*i+1; ir[0]=ir[1]-1
        for i in range(0,nod*nfree):
            it=ir[i]
            fp[it]=fp[it]+tfe[i]+bfe[i]
            for j in range(0,nod*nfree):
                jt=ir[j]
                gk[it,jt]=gk[it,jt]+sm[i,j]
    # treatment of boundary conditions
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp[iz]=0.0
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                fp=fp-rdis[j,i]*gk[:,iz]
                gk[:,iz]=0.0
                gk[iz,iz]=1.0
    # solution of simultaneous linear equations
    #disg = np.linalg.solve(gk, fp)
    gk = csr_matrix(gk)
    disg = spsolve(gk, fp, use_umfpack=True)
    # recovery of restricted displacements
    for i in range(0,npoin):
        for j in range(0,nfree):
            if mpfix[j,i]==1:
                iz=i*nfree+j
                disg[iz]=rdis[j,i]
    # calculation of element stress
    wd    =np.zeros(8,dtype=np.float64)
    strs  =np.zeros([4,nele],dtype=np.float64)      # Section force vector
    for ne in range(0,nele):
        i=node[0,ne]-1
        j=node[1,ne]-1
        k=node[2,ne]-1
        l=node[3,ne]-1
        m=node[4,ne]-1
        z1=x[0,i]; r1=x[1,i]
        z2=x[0,j]; r2=x[1,j]
        z3=x[0,k]; r3=x[1,k]
        z4=x[0,l]; r4=x[1,l]
        wd[0]=disg[2*i]; wd[1]=disg[2*i+1]
        wd[2]=disg[2*j]; wd[3]=disg[2*j+1]
        wd[4]=disg[2*k]; wd[5]=disg[2*k+1]
        wd[6]=disg[2*l]; wd[7]=disg[2*l+1]
        E    =ae[0,m]
        po   =ae[1,m]
        alpha=ae[2,m]
        dt1=deltaT[i]; dt2=deltaT[j]; dt3=deltaT[k]; dt4=deltaT[l]
        strs[:,ne]=calst_ax4(E,po,alpha,dt1,dt2,dt3,dt4,wd,z1,r1,z2,r2,z3,r3,z4,r4)
    # print out of result
    prout_ax4(fnameW,npoin,nele,disg,strs)
    # information
    dtime=time.time()-start
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec')
    fout=open(fnameW,'a')
    print('n={0}  time={1:.3f}'.format(nfree*npoin,dtime)+' sec',file=fout)
    fout.close()


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

以　上

トランジション（矩形断面から円形断面）での損失水頭

圧力水路の設計をしていると、取水口周辺で矩形断面から円形断面に変化する部分が必要となるが、その時の形状変化による損失（摩擦損失を除く）の計算について書いておく。 これまでは、矩形断面の断面積に等しい等価円を考え、漸縮管の損失水頭として扱ってきた。

今回同じ損失水頭計算が必要となり、しらべたら「矩形管の等価円」という考え方があったので、それについても記載しておく。

等価円の計算法は以下のページによる、

http://www.mekatoro.net/digianaecatalog/ebara-fan50/book/ebara-fan50-P0380.pdf

上記ページによれば等価円の算定式は以下の通り。

\begin{equation} D_e=1.3 \left[\cfrac{(a b)^5}{(a+b)^2}\right]^{0.125} \end{equation}
\begin{align} &D_e & & \text{equivalent diameter of circular duct}\\ &a & & \text{width of lectangular duct}\\ &b & & \text{height of rectangular duct} \end{align}

上記式および慣用的に用いてきた矩形断面と同一断面積を有する円の直径を計算してみる。

Upstream edge Length Downstream edge
Rectangular shape
(a x b = 8.4m x 8.4m)		 L=5.5m Circular shape
(diameter D=8.4m)
換算式を用いた等価半径
\begin{equation} D_e=1.3 \left[\cfrac{(8.4 \cdot 8.4)^5}{(8.4+8.4)^2}\right]^{0.125}=9.183 m \end{equation}
矩形断面と同一断面積を有する円管の直径
\begin{align} & \text{Section area of upstream rectangular duct} & A_1=8.4 \times 8.4 = 70.560 \\ & \text{Section area of downstream circular duct} & A_2=\pi \times 8.4^2 / 4 = 55.390 \end{align}

上記より、矩形断面と同一断面積を有する円管の直径は以下の通りとなる。

\begin{equation} D_1=\sqrt{\cfrac{4 \cdot 70.560}{\pi}}=9.481 m \end{equation}
損失水頭計算

この事例では、矩形断面から換算された円形断面の直径は、２つの手法でほぼ同じとなった。 よってここでは、これまで慣用的に用いてきた、矩形断面と同一断面積を有する円の直径を用いて、漸縮管としての損失水頭計算を行う。

漸縮管の損失水頭は、以下に示すGardelの図（出典：発電水力演習、千秋信一著）による。

f:id:damyarou:20200607123959p:plain

計算条件を整理すると以下のとおり。ここで設計流量はQ=290m3/sとする。

Q (m3/s) D1 (m) A1 (m2) v1 (m/s) D2 (m) A2 (m2) v2 (m/s) A2/A1 L (m)
290 9.481 70.560 4.110 8.400 55.390 5.236 0.785 5.5

角度 \theta は、次式で計算する。

\begin{equation} \theta=2 \tan^{-1}\left(\cfrac{D_1-D_2}{2 L}\right) \end{equation}
\begin{equation} \theta=2 \tan^{-1}\left(\cfrac{9.481-8.400}{2\cdot 5.5}\right)=11.225 degrees \end{equation}

断面積比 A_2 / A_1=0.785、角度 \theta=11.225 deg. より、損失係数は、上記図より f_{gc}=0.003 となる。 よって、損失水頭は以下の通り計算できる。

\begin{equation} h_{gc}=f_{gc}\cdot\cfrac{v_2{}^{2}}{2 g}=0.003\cdot \cfrac{5.236^2}{2 g}=0.004 m \end{equation}

以　上

パソコン支給

2020年5月21日、会社からノートパソコンというかSurfaceが支給されました、 マシン仕様は以下の通り。

  • プロセッサ：Intel(R) Core(TM) i5-1035G4 CPU @ 1.10GHz 1.50
  • 実装RAM：8GB
  • OS：Windows 10 Pro
  • ストレージ：SSD128GB

支給日の木曜日と金曜日は忙しかったので5月22日土曜日にソフトウエアのインストールを行いました。

ソフト 用途
Firefox ブラウザ
gPad テキストエディタ
VS code プログラムコード用エディタ
Google Earth Pro バーチャル地球儀
Python 3.8.3 プログラミング言語
MinGW-w64 (gcc version 8.1.0) gcc & gfortran コンパイラ
GMT (version 6.0.0 64-bit) Generic Mapping Tools、地図およびグラフ作成

Pythonの関連ライブラリは以下のものを pip でインストール。

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

Jupyterでは、デフォルトのブラウザをFirefoxに指定し、行番号のデフォルト表示にセット、最後にテーマ設定などをコマンドで設定。

デフォルトブラウザの設定

https://qiita.com/acknpop/items/4e5b57e38780068a9155

行番号のデフォルト表示

https://qiita.com/pollenjp/items/88600df83448a8ff5ea6

テーマ・フォントサイズ・表示幅・行間表示設定
jt -t oceans16 -fs 12 -ofs 12 -cellw 1200 -lineh 120 -N -T
集合写真

手持ちのiMac, Macbook pro と会社支給 Surface（一番左）の集合写真が下。

f:id:damyarou:20200524100408p:plain

以　上

Python 表付きグラフ（効率カーブ）

簡単な表と一体化した画像を作った。

作例

f:id:damyarou:20200516115000j:plain

このグラフにおけるTips

作図領域を上下に分ける

drawfig1() は下側のグラフ（作図領域の下側70%）、```drawfig2()''' は上川の表を描画する。

    plt.axes((0.0, 0.0, 1.0, 0.7)); drawfig1()
    plt.axes((0.0, 0.7, 1.0, 0.3)); drawfig2()

白抜き記号描画

  • ms=8 ：見号の大きさ
  • color='#000080' ：記号の色
  • markerfacecolor='#ffffff' ：白抜き
  • markeredgewidth=1 ：記号輪郭の太さ
    plt.plot(q_hmax,e_hmax,'s',ms=8,color='#000080',markerfacecolor='#ffffff',markeredgewidth=1,label='H=192.5m')

プログラム全文

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq


def reg(x,y):
    def func(param,x,y):
        a=param[0]
        b=param[1]
        c=param[2]
        return y - (a*(x-b)**2 + c)
    para0 =np.array([-1.0, 20.96, 89.65])
    res = leastsq(func, para0, args=(x,y))
    a=res[0][0]
    b=res[0][1]
    c=res[0][2]
    return a,b,c


def drawfig1():
    q_hmax=np.array([22.37,21.67,20.26,19.50,17.77,17.34,16.32,15.17,15.01,13.67,11.27])
    e_hmax=np.array([89.48,89.45,89.23,88.96,88.18,87.87,87.14,86.23,86.07,84.98,81.82])

    q_hnor=np.array([23.94,23.05,20.96,19.28,18.68,17.27,16.77,15.83,14.67,13.46,10.90])
    e_hnor=np.array([89.53,89.64,89.65,89.17,88.94,88.10,87.73,87.04,85.97,84.86,81.53])
    
    q_hmin=np.array([24.24,22.50,20.46,18.90,18.41,17.18,16.36,15.89,14.79,13.45,10.64])
    e_hmin=np.array([89.27,89.62,89.60,89.12,88.80,88.05,87.39,86.99,85.91,84.78,81.07])

    fsz=14
    xmin=5 ; xmax=25; dx=5
    ymin=75; ymax=95; dy=5
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Turbine discharge per unit Qu (m$^3$/s)')
    plt.ylabel('Combined efficiency ef (%)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#999999',linestyle='solid')
    #plt.gca().set_aspect('equal',adjustable='box')
    #plt.title(tstr,loc='left',fontsize=fsz)

    x=q_hnor
    y=e_hnor
    a,b,c=reg(x,y)
    qmax=24.5
    qmin=qmax*0.35
    xx=np.arange(qmin,qmax,0.005)
    yy=a*(xx-b)**2+c   
    ss='ef={0:6.3f}*(Qu-{1:6.3f})**2+{2:6.3f}'.format(a,b,c)
    print(np.min(yy))
    print(np.max(yy))
    
    plt.plot(xx,yy,'-',color='#000080',lw=2,label=ss)
    plt.plot(q_hmax,e_hmax,'s',ms=8,color='#000080',markerfacecolor='#ffffff',markeredgewidth=1,label='H=192.5m')
    plt.plot(q_hnor,e_hnor,'o',ms=8,color='#000080',markerfacecolor='#ffffff',markeredgewidth=1,label='H=180.0m')
    plt.plot(q_hmin,e_hmin,'^',ms=8,color='#000080',markerfacecolor='#ffffff',markeredgewidth=1,label='H=171.5m')
    
    plt.plot([qmax,qmax],[ymin,ymax],'--',color='#000080',lw=2)
    plt.plot([qmin,qmin],[ymin,ymax],'--',color='#000080',lw=2)
    xs=qmax; ys=ymin+0.02*(ymax-ymin); ss='Qu(max)={0:.1f}m$^3$/s'.format(qmax)
    plt.text(xs,ys,ss,ha='right',va='bottom',fontsize=fsz,rotation=90)
    xs=qmin; ys=ymin+0.02*(ymax-ymin); ss='Qu(min)=Qu(max) x 0.35'
    plt.text(xs,ys,ss,ha='right',va='bottom',fontsize=fsz,rotation=90)
    plt.legend(loc='upper left',shadow=True,fontsize=fsz-2)

    
def drawfig2():
    tstr='Design head setting (Q=49.0m$^3$/s)'
    row4=['Item','RWL','TWL','Margin','head loss','Value of head']
    row3=['Maximum head','EL.284.5','EL.92.0','0.0m','0.0m', 'H=192.5m']
    row2=['Nominal head','EL.284.5','EL.92.0','1.0m','11.5m','H=180.0m']
    row1=['Minimum head','EL.275.0','EL.92.0','0.0m','11.5m','H=171.5m']
    x=np.array([0.5,1.5,2.5,3.5,4.5,5.5])
    y=np.array([0.5,1.5,2.5,3.5])    
    fsz=12
    xmin=0 ; xmax=6
    ymin=-0.5 ; ymax=5
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.axis('off')
    for i in range(len(x)):
        for j in range(len(y)):
            if j==0: _row=row1
            if j==1: _row=row2
            if j==2: _row=row3
            if j==3: _row=row4
            ss=_row[i]
            plt.text(x[i],y[j],ss,va='center',ha='center',fontsize=fsz)
    plt.text(0.5*(xmin+xmax),4.5,tstr,va='center',ha='center',fontsize=fsz)
    for yy in [0,1,2,3,4]:
        plt.plot([xmin,xmax],[yy,yy],'-',color='#000000',lw=1)
    for xx in [0,1,2,3,4,5,6]:
        plt.plot([xx,xx],[0,4],'-',color='#000000',lw=1,clip_on=False)


def main():
    fnameF='fig_efficiency.jpg'
    plt.figure(figsize=(8,6),facecolor='w')
    plt.rcParams['font.family'] ='sans-serif'
    plt.rcParams['font.size']=14
    plt.axes((0.0, 0.0, 1.0, 0.7)); drawfig1()
    plt.axes((0.0, 0.7, 1.0, 0.3)); drawfig2()
    plt.tight_layout()
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()
    
    
#---------------
# Execute
#---------------
if __name__ == '__main__': main()

以　上

Python 月次・年次データ集計（雨量）

日雨量を月・年毎に合計したもの。 各月の雨量の単位は mm/momth （月累計雨量）、 各年の雨量の単位は mm/year（年累計雨量）である。

#==============================
# Rainfall
# Dayly => Monthly, Yearly data
#==============================
import pandas as pd
import numpy as np
import datetime


def rdata():
    fnameR='xls_rdata_20190630.xlsx' # input file name
    sname='daily'
    df=pd.read_excel(fnameR, sheet_name=sname,header=0, index_col=0) # read excel data
    df.index = pd.to_datetime(df.index, format='%Y/%m/%d')
    return df    
    

def cal_y(rf,lyyy):
    rfy=np.zeros(len(lyyy),dtype=np.float)
    for i,ss in enumerate(lyyy):
        rfy[i]=rf[ss].sum() # yearly rainfall
    return rfy


def cal_m(rf,lyyy,lmmm):
    rfm=np.zeros((len(lyyy),len(lmmm)+1),dtype=np.float)
    for i,yy in enumerate(lyyy):
        for j,mm in enumerate(lmmm):
            ss=yy+'/'+mm
            try:
                rfm[i,j]=rf[ss].sum() # monthly rainfall
            except KeyError:
                rfm[i,j]=np.nan
        rfm[i,-1]=np.sum(rfm[i,:])
    return rfm


def rf_year(lyy,la,df0):
    rfy=np.zeros((len(lyy),len(la)),dtype=np.float64)
    for i,ss in enumerate(la):
        rf=pd.Series(df0[ss], index=df0.index)
        rfy[:,i]=cal_y(rf,lyy)

    ss='{0:4s}'.format('Year')
    for k in [1,2,3,4,5,6,7]:
        ss=ss+',{0:>8s}'.format('Area-'+str(k))
    ss=ss+',{0:>8s}'.format('Dam')
    ss=ss+',{0:>8s}'.format('Teromu')
    print(ss)
    for j,yy in enumerate(lyy):
        ss='{0:4s}'.format(yy)
        for i in range(len(la)):
            ss=ss+',{0:8.0f}'.format(rfy[j,i])
        print(ss)               
    ss='{0:4s}'.format('Ave.')
    for i in range(len(la)):
        ss=ss+',{0:8.0f}'.format(np.mean(rfy[:,i]))
    print(ss)


def rf_mon(lyy,lmm,rf):
    rfm=cal_m(rf,lyy,lmm)
    ss='{0:4s}'.format('Year')
    for m in ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aig','Sep','Oct','Nov','Dec','Sum']:
        ss=ss+',{0:>6s}'.format(m)
    print(ss)
    for i,yy in enumerate(lyy):
        ss='{0:4s}'.format(yy)
        for j in range(len(lmm)+1):
            ss=ss+',{0:6.0f}'.format(rfm[i,j])
        print(ss)               
    ss='{0:4s}'.format('Ave.')
    for j in range(len(lmm)+1):
        ss=ss+',{0:6.0f}'.format(np.nanmean(rfm[:,j]))
    print(ss)


def main():
    df0=rdata() # daily rainfall at each area (2000.03-2019.06)
    lyy=['2001','2002','2003','2004','2005','2006','2007','2008','2009',
         '2010','2011','2012','2013','2014','2015','2016','2017','2018']
    la=['area1','area2','area3','area4','area5','area6','area7','Dam','Teromu']
    rf_year(lyy,la,df0)
    
    lyy=['2000','2001','2002','2003','2004','2005','2006','2007','2008','2009',
         '2010','2011','2012','2013','2014','2015','2016','2017','2018','2019']
    lmm=['01','02','03','04','05','06','07','08','09','10','11','12']

    print('* Monthly rainfall at Dam')
    rf=pd.Series(df0['Dam'], index=df0.index)
    rf_mon(lyy,lmm,rf)

    print('* Monthly rainfall at Teromu')
    rf=pd.Series(df0['Teromu'], index=df0.index)
    rf_mon(lyy,lmm,rf)

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

以　上

Python　月次及び年次データ集計（流量）

概要

日平均流量の月次および年次データ集計のプログラム。

ここでは流量の集計を行い、月平均・年平均を算出している。 年平均は、「月平均の平均」ではなく、「１年間全データの平均」としている。 これは365日の流況図から出した年平均と数値を合わせるため。 「月平均の平均」と「年間全データの平均」は結構違う場合があるので注意。

計算のルール

入力

データ取得開始日の月初めから、データ取得完了日の月末まで、データの有無に関わらず、全ての日付が連続的に入力されていること。

出力

  • 各月の中で１日でも欠測があれば、その月の平均は nan とする
  • 各年の平均は、各年の全データ（365日あるいは366日）の平均とする。その年に１日でも欠測があればその年の平均は nan とする。
  • 各月の集計平均（最終行）は、欠測を含まない月の全データの平均値とする。欠測がある月のデータは月平均の計算に用いない。

出力例

  • 年別に各月の平均日流量を表示している。
  • 最終列はその年の平均日流量。
  • 最終行は各月の年をまたいだ平均日流量。
  • Wordに貼り付けて数表を作るため、画面にCSVイメージで出力する。
  • 各月の中で、１日でも欠測があれば月平均は、nanとしている。
  • nan の月を含む年の平均はもちろん nan である。
Year,   Jan,   Feb,   Mar,   Apr,   May,   Jun,   Jul,   Aig,   Sep,   Oct,   Nov,   Dec,  Ave.
2013,  70.6,  83.3, 118.3, 124.3, 117.4, 115.6,  93.4,  53.2,  41.1,  28.4,  76.7, 113.9,  86.3
2014,  46.9,  36.5,  59.2,  71.9,  42.1,  86.1,  79.0,  50.6,  28.1,   nan,  23.0,   nan,   nan
2015,  47.5,  95.7, 101.0,  90.9,  74.9,  37.8,  32.0,  27.8,  15.4,  14.0,   9.1,  37.0,  48.3
2016,  38.2,  72.3, 125.8, 145.9,  72.2,  80.1,  49.4,  54.0,  49.7,  92.2, 112.1,  54.1,  78.7
2017,  50.7,  63.0,  99.2, 108.1, 128.7, 105.6,  86.3,   nan,   nan,   nan,   nan,   nan,   nan
2018,   nan,   nan,   nan,   nan,  84.1,   nan,   nan,   nan,  31.9,  32.6,  37.4,  64.7,   nan
2019,  55.2,  91.4,  68.4, 105.9,  88.5,   nan,   nan,   nan,   nan,   nan,   nan,   nan,   nan
Ave.,  51.5,  73.7,  95.3, 107.8,  86.8,  85.0,  68.0,  46.4,  33.2,  41.8,  51.7,  67.4,  71.1

プログラム

#====================================
# Data aggregation (monthly & yearly)
#====================================
import pandas as pd
import numpy as np
import datetime


def rdata():
    fnameR='df_qq2.csv' # input file name
    df=pd.read_csv(fnameR, header=0, index_col=0) # read excel data
    df.index = pd.to_datetime(df.index, format='%Y/%m/%d')
    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)
    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)
    ss='{0:4s}'.format('Year')
    for m in ['Jan','Feb','Mar','Apr','May','Jun','Jul','Aig','Sep','Oct','Nov','Dec','Ave.']:
        ss=ss+',{0:>6s}'.format(m)
    print(ss)
    lyy=lyy+['Ave.']
    for i,yy in enumerate(lyy):
        ss='{0:4s}'.format(yy)
        for j in range(len(lmm)+1):
            ss=ss+',{0:6.1f}'.format(qqm[i,j])
        print(ss)               


def main():
    df0=rdata() # daily rainfall at each area (2000.03-2019.06)
    lyy=['2013','2014','2015','2016','2017','2018','2019']
    lmm=['01','02','03','04','05','06','07','08','09','10','11','12']

    print('* Monthly discharge at Teromu GS')
    qq=pd.Series(df0['q_tot'], index=df0.index)
    qq_mon(lyy,lmm,qq)

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

以　上