damyarou

python, GMT などのプログラム

設計 岩盤内埋設式水圧鉄管の限界座屈圧力

設計

記事の最後に行く

設計理論

Amstutzの式として知られる方法を紹介する（水門鉄管技術基準　水圧鉄管・鉄鋼構造物編　第29条解説による）

限界座屈圧力 (p_k) 算定式

\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)} \tag{1} \end{equation}

座屈発生時の周方向一様圧縮応力 (\sigma_N) 算定式

\begin{equation} \left(\cfrac{k_0}{r_m}+\cfrac{\sigma_N}{E_s^*}\right)\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^*}\left(1-\cfrac{1}{2}\cdot\cfrac{r_m}{t}\cdot\cfrac{\sigma_F^*-\sigma_N}{E_s^*}\right) \tag{2} \end{equation}

間隙量 (k_0) 算定式

\begin{equation} k_0=\cfrac{\left(\alpha_s\cdot\Delta T+\beta_g\cfrac{\sigma_a\cdot\eta}{E_s}\right)r_0'}{1+\beta_g} \tag{3} \end{equation}

鋼材の力学特性

\begin{equation} E_s^*=\cfrac{E_s}{1-\nu_s{}^2} \qquad \sigma_F^*=\mu\cfrac{\sigma_F}{\sqrt{1-\nu_s+\nu_s{}^2}} \qquad \mu=1.5-0.5\cfrac{1}{\left(1+0.002\cfrac{E_s}{\sigma_F}\right)^2} \tag{4} \end{equation}

記号の定義

記号説明１説明２変数
p_k 鉄管の限界座屈圧力 最終解 pk
\sigma_N 座屈発生時の周方向一様圧縮応力 中間解 sigN
k_0 鉄管外面とコンクリート間の間隙量 パラメータ k0
D_0 設計鉄管内径 入力値 D0
t_0 設計鉄管板厚 入力値 t0
\sigma_F 鉄管材料の降伏点 入力値 sigF
\sigma_a 鉄管材料の許容応力 入力値 siga
\epsilon 板厚の余裕代 標準値（=1.5 mm） eps
E_s 鉄管材料の弾性係数 標準値（=206,000 MPa） Es
\nu_s 手間材料のポアソン比 標準値（=0.3） pos
\eta 溶接継手効率 標準値（=1.0） eta
\alpha_s 鉄管材料の線膨張係数 標準値（=1.2\times 10^{-5} /^{\circ}C） alphas
\Delta T 鉄管の温度降下量 標準値（=20 ^{\circ}C） deltaT
\beta_g 周辺岩盤の塑性変形係数 標準値（=1.0） betag
t 余裕代を差し引いた板厚（=t_0-\epsilon） 途中計算値 t
r_m 鉄管の板厚中心半径（=(D_0+t_0)/2） 途中計算値 rm
r_0' 鉄管の外半径（=(D_0+2 t_0)/2） 途中計算値 r0d
E_s^* 鉄管材料のポアソン効果を考慮した弾性係数 途中計算値 Ess
\sigma_F^* 鉄管材料のポアソン効果を考慮した降伏点 途中計算値 sigFs
\mu 支持効果を表す係数 途中計算値 mu 
 Material properties
 ----------------------------------------------------------
 Kind of steel           HT100  HT-80  SM570  SM490  SM400 
 ----------------------------------------------------------
 Yield strength   (MPa)   885    685    450    315    235
 Allowable stress (MPa)   400    330    240    175    130
 ----------------------------------------------------------

出力画像

f:id:damyarou:20190430155257p:plain

プログラム

設計理論中の\sigma_Nは非線形方程式を解いて求める必要がある。 \sigma_Nの値は0と\sigma_F^*の間にあるため、二分法により求めることができる。

二分法プログラムとして、自作関数によるものと scipy brentq を用いたものを作成した。 計算時間は以下の通りであり、scipy brentq のほうが圧倒的に早い。収束判定基準はいずれも 1e-12 である。

  • 二分法自作関数：0.262秒
  • scipy brentq：0.029秒

以下にプログラムを示す。

二分法自作関数

import するのはnumpyと計算時間計測用のtimeのみである。 二分法の計算は、関数：funcを呼び出しながら、関数：bisectionで行っている。

プログラム全文

import numpy as np
import time


def func(x,params):
    # x=sigN
    k0   =params[0]
    rm   =params[1]
    t    =params[2]
    Ess  =params[3]
    sigFs=params[4]
    aa = k0/rm+x/Ess
    bb = 1+12*rm**2/t**2*x/Ess
    cc = 3.36*rm/t*(sigFs-x)/Ess*(1-1/2*rm/t*(sigFs-x)/Ess)
    f = aa*bb**1.5-cc
    return f


def cal_pk(sigN,sigFs,Ess,rm,t):
    aa=rm/t*(1+0.35*rm/t*(sigFs-sigN)/Ess)
    pk = sigN/aa
    return pk


def cal_k0(siga,eta,Es,alphas,deltaT,r0d,betag):
    aa = alphas*deltaT+betag*siga*eta/Es
    k0 = aa*r0d/(1+betag)
    return k0


def bisection(x01,x02,params):
    imax=100
    err=2.0e-12
    x1=x01
    x2=x02
    for i in range(0,imax):
        xi=0.5*(x1+x2)
        f1=func(x1,params)
        f2=func(x2,params)
        fi=func(xi,params)
        if abs(x2-x1) < err: break
        if f1*fi < 0.0: x2=xi
        if fi*f2 < 0.0: x1=xi
    return xi


def calc():
    # main calculation
    eps=1.5         # margin of plate thickness
    Es=206000.0     # elastic modulus of steel plate
    pos=0.3         # Poisson's ratio of steel plate
    eta=1.0         # welding efficiency
    alphas=1.2e-5   # thermal expansion coefficient
    deltaT=20.0     # temperature decrease
    betag=1.0       # plastic deformation coefficient of bedrock

    t0=30.0    # design plate thickness of penstock

    # Material properties
    # ----------------------------------------------------------
    # Kind of steel           HT100  HT-80  SM570  SM490  SM400 
    # ----------------------------------------------------------
    # Yield strength   (MPa)   885    685    450    315    235
    # Allowable stress (MPa)   400    330    240    175    130
    # ----------------------------------------------------------

    lnm=[1,2,3,4,5]           # number
    lys=[885,685,450,315,235] # yield strength
    las=[400,330,240,175,130] # allowable stress
    asl=np.arange(35,141,1)   # range of slenderness ratio
    data=np.zeros([len(asl),len(lnm)+1],dtype=np.float64) # array for data store
    data[:,0]=asl[:]

    for i,sigF,siga in zip(lnm,lys,las):
        for k,slr in enumerate(asl):
            D0=2*t0*slr  # design internal diameter of penstock
            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=(D0+t0)/2
            r0d=(D0+2*t0)/2
            k0=cal_k0(siga,eta,Es,alphas,deltaT,r0d,betag)
            params=np.array([k0,rm,t,Ess,sigFs])
            sigN=bisection(0.0,sigFs,params) # self made
            pk=cal_pk(sigN,sigFs,Ess,rm,t)
            data[k,i]=pk
    return data
        
        
def main():
    start=time.time()
    data=calc()
    fnameW='out.txt'
    np.savetxt(fnameW,data)
    dtime=time.time()-start
    print(dtime)


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

Scipy Brent法

import するのは numpy と計算時間計測用の time に加え、brentq を使うための scipy.optimize である。 二分法の計算を自作関数(bisection)を用いて行うプログラムとの違いは、インポートするモジュールと、\sigma_Nを計算するための以下の１文のみである。

# 自作関数bisectionの呼び出し
sigN=bisection(0.0,sigFs,params) # self made

# scipy brenyqの呼び出し
sigN=optimize.brentq(func,0.0,sigFs,args=(params)) # scipy

プログラム全文

import numpy as np
from scipy import optimize
import time


def func(x,params):
    # x=sigN
    k0   =params[0]
    rm   =params[1]
    t    =params[2]
    Ess  =params[3]
    sigFs=params[4]
    aa = k0/rm+x/Ess
    bb = 1+12*rm**2/t**2*x/Ess
    cc = 3.36*rm/t*(sigFs-x)/Ess*(1-1/2*rm/t*(sigFs-x)/Ess)
    f = aa*bb**1.5-cc
    return f


def cal_pk(sigN,sigFs,Ess,rm,t):
    aa=rm/t*(1+0.35*rm/t*(sigFs-sigN)/Ess)
    pk = sigN/aa
    return pk


def cal_k0(siga,eta,Es,alphas,deltaT,r0d,betag):
    aa = alphas*deltaT+betag*siga*eta/Es
    k0 = aa*r0d/(1+betag)
    return k0


def calc():
    # main calculation
    eps=1.5         # margin of plate thickness
    Es=206000.0     # elastic modulus of steel plate
    pos=0.3         # Poisson's ratio of steel plate
    eta=1.0         # welding efficiency
    alphas=1.2e-5   # thermal expansion coefficient
    deltaT=20.0     # temperature decrease
    betag=1.0       # plastic deformation coefficient of bedrock

    t0=30.0    # design plate thickness of penstock

    # Material properties
    # ----------------------------------------------------------
    # Kind of steel           HT100  HT-80  SM570  SM490  SM400 
    # ----------------------------------------------------------
    # Yield strength   (MPa)   885    685    450    315    235
    # Allowable stress (MPa)   400    330    240    175    130
    # ----------------------------------------------------------

    lnm=[1,2,3,4,5]           # number
    lys=[885,685,450,315,235] # yield strength
    las=[400,330,240,175,130] # allowable stress
    asl=np.arange(35,141,1)   # range of slenderness ratio
    data=np.zeros([len(asl),len(lnm)+1],dtype=np.float64) # array for data store
    data[:,0]=asl[:]

    for i,sigF,siga in zip(lnm,lys,las):
        for k,slr in enumerate(asl):
            D0=2*t0*slr  # design internal diameter of penstock
            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=(D0+t0)/2
            r0d=(D0+2*t0)/2
            k0=cal_k0(siga,eta,Es,alphas,deltaT,r0d,betag)
            params=np.array([k0,rm,t,Ess,sigFs])
            sigN=optimize.brentq(func,0.0,sigFs,args=(params)) # scipy
            pk=cal_pk(sigN,sigFs,Ess,rm,t)
            data[k,i]=pk
    return data


def main():
    start=time.time()
    data=calc()
    fnameW='out.txt'
    np.savetxt(fnameW,data)
    dtime=time.time()-start
    print(dtime)


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

作図プログラム

計算プログラムでは計算結果を numpy 配列に格納・ファイル書き出しし、作図プログラムでは書き出した numpy 配列をそのまま読み込んで配列データを作成している。

# numpy配列書き出し・保存（計算プログラム）
fnameW='out.txt'
np.savetxt(fnameW,data)

# numpy配列読み込み（作図プログラム）
fnameR='out.txt'
data = np.loadtxt(fnameR, usecols=(0,1,2,3,4,5),skiprows=0)

プログラム全文

import numpy as np
import matplotlib.pyplot as plt


def datainp():
    fnameR='out.txt'
    data = np.loadtxt(fnameR, usecols=(0,1,2,3,4,5),skiprows=0)
    return data


def sarrow(x1,y1,x2,y2,ss,fsz):
    sv=0
    plt.annotate(
        ss,
        xy=(x1,y1),
        xycoords='data',
        xytext=(x2,y2),
        textcoords='data',
        fontsize=fsz,
        ha='left', va='center',
        arrowprops=dict(shrink=sv,width=0.1,headwidth=5,headlength=10,
                        connectionstyle='arc3',facecolor='#777777',edgecolor='#777777'))

    
def drawfig(data):
    # Drawing
    fsz=16
    xmin=30; xmax=140; dx=10
    ymin=0; ymax=12; dy=1
    plt.figure(figsize=(10,8),facecolor='w')
    plt.rcParams["font.size"] = fsz
    plt.rcParams['font.family'] ='sans-serif'
    plt.xlim(xmin,xmax)
    plt.ylim(ymin,ymax)
    plt.xlabel('$D_0 / 2 t_0$')
    plt.ylabel('Critical buckling pressure (MPa)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(which='both',color='#999999',linestyle='solid')
    plt.title('Critical Buckling Pressure of Embeded Penstock',loc='left',fontsize=fsz)
    for i in range(1,6):
        plt.plot(data[:,0],data[:,i],'-',color='#000000',lw=1.5)
    plt.plot([35,35],[ymin,ymax],'--',color='#000000',lw=1.0)
    plt.text(35,1,'Application limit ($D_0/2t_0=35$)',ha='right',va='bottom',rotation=90,fontsize=fsz)

    _i=np.where(data[:,0]==40); i1=_i[0][0]
    _i=np.where(data[:,0]==45); i2=_i[0][0]
    _i=np.where(data[:,0]==50); i3=_i[0][0]
    _i=np.where(data[:,0]==55); i4=_i[0][0]
    _i=np.where(data[:,0]==60); i5=_i[0][0]
    kk=[1,2,3,4,5]
    ii=[i1,i2,i3,i4,i5]
    st=[
        '$\sigma_F$=885 MPa (HT100)',
        '$\sigma_F$=685 MPa (HT-80)',
        '$\sigma_F$=450 MPa (SM590)',
        '$\sigma_F$=315 MPa (SM490)',
        '$\sigma_F$=235 MPa (SM400)'
    ]
    for k,i,ss in zip(kk,ii,st):
        x1=data[i,0]; y1=data[i,k]
        x2=x1+10    ; y2=y1+2
        sarrow(x1,y1,x2,y2,ss,fsz)

    fnameF='fig_ams.png'
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
    plt.show()

    
def main():
    data=datainp()
    drawfig(data)
    
    
#==============
# Execution
#==============
if __name__ == '__main__': main()

Thank you.

記事の先頭に行く