2019-05-29

設計仮排水路設計における Flood Routine の活用（3）

設計

ここでは、「設計仮排水路設計における Flood Routine の活用（２）」で述べた事項に関連する計算・作図プログラムのうち、以下の4つを掲載する。

トンネル標準断面作図プログラム
トンネル通水量計算・作図プログラム
貯水池水位容量曲線作図プログラム
洪水波形作図プログラム

トンネル標準断面作図プログラム

このような説明図を描くプログラムが意外に難しいというか面倒くさい。

断面形は、幅と高さが等しい複合円である。この断面の場合、トンネル内幅を $D$ とすれば、上半円の半径は $D/2$ 、下半円の半径は $D$ となり、インバート幅は $(\sqrt{3}-1) D=0.732 D$ となる。道路トンネルの場合、インバート幅をきっちりした数字とする場合が多く、下半円の半径が中途半端な数字となることが多いが、水路トンネルの場合のインバートの制約は、施工用重機が往来できればいい。このため、水路トンネルとしてのこの断面形状は、幾何学的にきれいで結構気に入っている。

トンネル形状は、中心を原点として、水平右の点から角度を指定して (x, y) 座標を求めて置き、一気にプロットしている。

寸法線と寸法は、両矢印に記載するものと、片矢印に記載するものがあるため、以下のように３つの関数を作成して描画している。

両矢印を描く関数

def barrow(x1,y1,x2,y2):
    # (x1,y1),(x2,y2)：両矢印始点と終点
    col='#333333'
    arst='<->,head_width=3,head_length=10'
    aprop=dict(arrowstyle=arst,connectionstyle='arc3',facecolor=col,edgecolor=col,shrinkA=0,shrinkB=0,lw=1)
    plt.annotate('',
        xy=(x1,y1), xycoords='data',
        xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)

片矢印を描く関数

def sarrow(x1,y1,x2,y2):
    # (x1,y1)：矢印先端座標
    # (x2,y2)：矢印始点座標
    col='#333333'
    arst='->,head_width=3,head_length=10'
    aprop=dict(arrowstyle=arst,connectionstyle='arc3',facecolor=col,edgecolor=col,shrinkA=0,shrinkB=0,lw=1)
    plt.annotate('',
        xy=(x1,y1), xycoords='data',
        xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)

指定した直線に対し文字列を描画する関数

def atext(x1,y1,x2,y2,ss,dd,ii):
    # (x1,y1),(x2,y2)：矢印始点と終点
    # 矢印始点と終点座標は文字列の傾きを考慮して決定する
    # ss：描画文字列
    # dd：文字列中央の寸法線からの距離
    # ii：１なら寸法線の上側，２なら寸法線の下側にテキストを描画
    al=np.sqrt((x2-x1)**2+(y2-y1)**2)
    cs=(x2-x1)/al
    sn=(y2-y1)/al
    rt=np.degrees(np.arccos(cs))
    if y2-y1<0: rt=-rt
    if ii==1:
        xs=0.5*(x1+x2)-dd*sn
        ys=0.5*(y1+y2)+dd*cs
        plt.text(xs,ys,ss,ha='center',va='center',rotation=rt)
    else:
        xs=0.5*(x1+x2)+dd*sn
        ys=0.5*(y1+y2)-dd*cs
        plt.text(xs,ys,ss,ha='center',va='center',rotation=rt)

プログラム全文

import numpy as np
import matplotlib.pyplot as plt


def barrow(x1,y1,x2,y2):
    # (x1,y1),(x2,y2)：両矢印始点と終点
    col='#333333'
    arst='<->,head_width=3,head_length=10'
    aprop=dict(arrowstyle=arst,connectionstyle='arc3',facecolor=col,edgecolor=col,shrinkA=0,shrinkB=0,lw=1)
    plt.annotate('',
        xy=(x1,y1), xycoords='data',
        xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)

    
def sarrow(x1,y1,x2,y2):
    # (x1,y1)：矢印先端座標
    # (x2,y2)：矢印始点座標
    col='#333333'
    arst='->,head_width=3,head_length=10'
    aprop=dict(arrowstyle=arst,connectionstyle='arc3',facecolor=col,edgecolor=col,shrinkA=0,shrinkB=0,lw=1)
    plt.annotate('',
        xy=(x1,y1), xycoords='data',
        xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)

    
def atext(x1,y1,x2,y2,ss,dd,ii):
    # (x1,y1),(x2,y2)：矢印始点と終点
    # 矢印始点と終点座標は文字列の傾きを考慮して決定する
    # ss：描画文字列
    # dd：文字列中央の寸法線からの距離
    # ii：１なら寸法線の上側，２なら寸法線の下側にテキストを描画
    al=np.sqrt((x2-x1)**2+(y2-y1)**2)
    cs=(x2-x1)/al
    sn=(y2-y1)/al
    rt=np.degrees(np.arccos(cs))
    if y2-y1<0: rt=-rt
    if ii==1:
        xs=0.5*(x1+x2)-dd*sn
        ys=0.5*(y1+y2)+dd*cs
        plt.text(xs,ys,ss,ha='center',va='center',rotation=rt)
    else:
        xs=0.5*(x1+x2)+dd*sn
        ys=0.5*(y1+y2)-dd*cs
        plt.text(xs,ys,ss,ha='center',va='center',rotation=rt)

        
def tunnel(d0):
    rr=d0/2
    ang=np.linspace(0,np.pi,180,endpoint=False)
    xx1=rr*np.cos(ang)
    yy1=rr*np.sin(ang)
    rr=d0
    ang=np.linspace(np.pi,np.pi+np.radians(30),30,endpoint=True)
    xx2=d0/2+rr*np.cos(ang)
    yy2=rr*np.sin(ang)
    ang=np.linspace(np.radians(-30),0,30,endpoint=True)
    xx3=-d0/2+rr*np.cos(ang)
    yy3=rr*np.sin(ang)
    xx=np.hstack([xx1,xx2,xx3])
    yy=np.hstack([yy1,yy2,yy3])
    return xx,yy
    

def main():
    d0=8.0
    fsz=16
    xmin=-6; xmax=6
    ymin=-6; ymax=6
    plt.figure(figsize=(6,6),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.axis('off')
    plt.gca().set_aspect('equal',adjustable='box')
#    plt.title(tstr,loc='center',fontsize=fsz)
    xx,yy=tunnel(d0)
    plt.plot(xx,yy,'-',color='#000000',lw=3)

    xtc=0; ytc= d0/2
    xbc=0; ybc=-d0/2
    xcl=-d0/2; ycl=0
    xcr= d0/2; ycr=0
    xbl=d0/2-d0*np.cos(np.radians(30)); ybl=-d0/2
    xbr=d0*np.cos(np.radians(30))-d0/2; ybr=-d0/2
    ds=1.0
    plt.plot([0,0],[ybc-0.5*ds,ytc+0.5*ds],'-.',color='#000000',lw=1)
    plt.plot([xcl-0.5*ds,xcr+0.5*ds],[0,0],'-.',color='#000000',lw=1)

    ss='$(\sqrt{3}-1) D$'; dd=0.5*ds; ii=-1
    plt.plot([xbl,xbl],[ybl,ybl-1.5*ds],'-',color='#000000',lw=1)
    plt.plot([xbr,xbr],[ybr,ybr-1.5*ds],'-',color='#000000',lw=1)
    x1=xbl; y1=ybl-1.0*ds; x2=xbr; y2=y1
    barrow(x1,y1,x2,y2)
    atext(x1,y1,x2,y2,ss,dd,ii)

    ss='$D$'; dd=0.3*ds; ii=1
    plt.plot([xcl,xcl],[ycl,ytc+1.5*ds],'-',color='#000000',lw=1)
    plt.plot([xcr,xcr],[ycr,ytc+1.5*ds],'-',color='#000000',lw=1)
    x1=xcl; y1=ytc+1.0*ds; x2=xcr; y2=y1
    barrow(x1,y1,x2,y2)
    atext(x1,y1,x2,y2,ss,dd,ii)

    ss='$D$'; dd=0.3*ds; ii=1
    plt.plot([xtc,xcl-1.5*ds],[ytc,ytc],'-',color='#000000',lw=1)
    plt.plot([xbl,xcl-1.5*ds],[ybl,ybl],'-',color='#000000',lw=1)
    x1=xcl-1.0*ds; y1=ybl; x2=x1; y2=ytc
    barrow(x1,y1,x2,y2)
    atext(x1,y1,x2,y2,ss,dd,ii)
 
    ss='$D \; / \; 2$'; dd=0.3*ds; ii=1
    x1=d0/2*np.cos(np.radians(30)); y1=x2=d0/2*np.sin(np.radians(30)); x2=0; y2=0
    sarrow(x1,y1,x2,y2)
    atext(x2,y2,x1,y1,ss,dd,ii)
    
    ss='$D$'; dd=0.3*ds; ii=1
    x1=xbr; y1=ybr; x2=xcl; y2=ycl;
    sarrow(x1,y1,x2,y2)
    atext(0.5*(x1+x2),0.5*(y1+y2),x1,y1,ss,dd,ii)
    x1=xbl; y1=ybl; x2=xcr; y2=ycr;
    sarrow(x1,y1,x2,y2)
    atext(x1,y1,0.5*(x1+x2),0.5*(y1+y2),ss,dd,ii)
    
    plt.tight_layout()
    fnameF='fig_model.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()


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

トンネル通水量計算・作図プログラム

等流水深部（自由水面）での計算は、通水断面積・潤辺は厳密な計算をせず、通水部を小さな台形に近似して計算している。

import numpy as np
import matplotlib.pyplot as plt


# global variables
dh=0.1 # increment of reservoir water level
g=9.8       # gravity acceleration
n=0.015     # Manning's roughness coefficient
elvs=64.0   # invert elevation of entrance     
elve=62.5   # invert elevation of exit


def drawfig():
    fsz=16
    xmin=0 ; xmax=5000; dx=500
    ymin=60; ymax=110; dy=10
    tstr='Discharge Capacity of Diversion Tunnel'
    plt.figure(figsize=(10,6),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Discharge Capacity (m$^3$/s)')
    plt.ylabel('Reservoir water Level (EL.m)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#999999',linestyle='solid')
    plt.title(tstr,loc='left',fontsize=fsz)
    a_d0=np.array([5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0])
    l_col=[
        '#dddddd',
        '#cccccc',
        '#bbbbbb',
        '#aaaaaa',
        '#999999',
        '#666666',
        '#333333',
        '#000000',
    ]
    for d0,col in zip(a_d0,l_col):
        fnameR='inp_sec_{0:0>2d}.txt'.format(int(d0))
        data = np.loadtxt(fnameR, usecols=(0,1) , skiprows=0)
        el=data[:,0]
        qq=data[:,1]
        ss='D={0:.0f}m'.format(d0)
        plt.plot(qq,el,'-',color=col,lw=2,label=ss)
    plt.plot([xmin,xmax],[elvs,elvs],'-',color='#0000ff',lw=1)
    xs=10; ys=elvs; ss='Tunnel Entrance Invert EL.{0:.1f}m'.format(elvs)
    plt.text(xs,ys,ss,color='#0000ff',ha='left',va='top',fontsize=fsz)
    plt.legend(fontsize=fsz-2)
    fnameF='fig_sec.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.2)
    plt.show()

    
def loss_calc(ell,rho1,theta1,rho2,theta2,d0):
    a0=-0.5*d0
    r1=0.5*d0
    hd=r1
    nn=int((hd+r1)/dh)+1

    sel=[]
    sqq=[]
    saa=[]
    w1=d0*(np.sqrt(3)-1)
    i=(elvs-elve)/ell
    el=elvs
    qq=0.0
    aa=0.0
    ss=w1
    rr=aa/ss
    sel=sel+[el]
    sqq=sqq+[qq]
    saa=saa+[aa]
    # free flow
    for j in range(1,nn):
        h=dh*float(j)
        y=h-hd
        if y<0.0:
            r=d0; a=a0
        if 0.0<=y and y<=r1:
            r=r1; a=0.0
        x=a+np.sqrt(r**2-y**2)
        w2=2*x
        da=0.5*(w1+w2)*dh
        ds=2*np.sqrt(dh**2+np.abs(0.5*(w1-w2))**2)
        aa=aa+da
        ss=ss+ds
        rr=aa/ss
        vv=1/n*rr**(2/3)*np.sqrt(i)
        qq=aa*vv
        el=elvs+h
        w1=w2
        sel=sel+[el]
        sqq=sqq+[qq]
        saa=saa+[aa]
    # pressure flow
    fe=0.25
    fo=1.0
    fb1=(0.131+0.1632*(d0/rho1))*np.sqrt(theta1/90.0)
    fb2=(0.131+0.1632*(d0/rho2))*np.sqrt(theta2/90.0)
    ff=2*g*n*n/rr**(1/3)
    while 1:
        if 110.0<el: break
        h=h+dh
        el=elvs+h
        head=el-elve-hd
        vv=np.sqrt(2*g*head/(fe+fo+ff*ell/rr+fb1+fb2))
        qq=aa*vv
        sel=sel+[el]
        sqq=sqq+[qq]
        saa=saa+[aa]
        a_el=np.array(sel)
        a_qq=np.array(sqq)
        a_aa=np.array(saa)
    return a_el,a_qq,a_aa


def main():
    for d0 in [5.0,6.0,7.0,8.0,9.0,10.0,11.0,12.0]:
        fnameW='inp_sec_{0:0>2d}.txt'.format(int(d0))
        # inside tunnel
        ell=400.0
        rho1=230; theta1=39.5
        rho2=200; theta2=48.2
        sel1,sqq1,saa1=loss_calc(ell,rho1,theta1,rho2,theta2,d0)
        # outside tunnel
        ell=450.0
        rho1=260; theta1=39.5
        rho2=230; theta2=48.2
        sel2,sqq2,saa2=loss_calc(ell,rho1,theta1,rho2,theta2,d0)
        # total flow capacity
        sel=0.5*(sel1+sel2)
        sqq=sqq1+sqq2
        saa=saa1+saa2
        fout=open(fnameW,'w')
        for j in range(0,len(sel)):
            print('{0:10.3f}{1:10.3f}{2:10.3f}'.format(sel[j],sqq[j],saa[j]),file=fout)
        fout.close()
    drawfig()

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

貯水池水位容量曲線作図プログラム

水位と容量の関係は、オリジナルデータでは標高10ｍピッチで与えられているため、３次スプラインで内挿し、標高1mピッチのデータとして出力・作図している。

import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt


def drawfig(el1,vv1,el2,vv2):
    fsz=16
    xmin=0 ; xmax=500; dx=50
    ymin=60; ymax=110; dy=5
    tstr='Storage Capacity Curve'
    plt.figure(figsize=(10,6),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Storage ($\\times 10^6 \; m^3$)')
    plt.ylabel('Reservoir water Level (EL.m)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#999999',linestyle='solid')
    plt.title(tstr,loc='left',fontsize=fsz)
    plt.plot(vv1,el1,'-',color='#0000ff',lw=2,clip_on=False)
    plt.plot(vv2,el2,'o',color='#0000ff',clip_on=False)
    fnameF='fig_hv.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.2)
    plt.show()


def main():
    # vv x 1e6 (m3)
    data=np.array([[60, 0],
               [70, 2.485],
               [80, 21.670],
               [90, 88.555],
               [100, 223.959],
               [110, 446.817],
               [120, 778.284],
               [130, 1238.519],
               [140, 1865.721],
               [150, 2375.546],
               [160, 3416.319],
               [170, 4773.768]])
    _el=data[:,0]
    _vv=data[:,1]
    f = interp1d(_el, _vv, kind='cubic')
    xmin=60
    xmax=170
    dx  =1
    ndx=int((xmax-xmin)/dx)+1
    el = np.linspace(xmin,xmax,ndx)
    vv =f(el)
    fnameW='inp_hv.txt'
    fout=open(fnameW,'w')
    for i in range(0,ndx):
        print('{0:5.1f}  {1:12.3f}'.format(el[i],vv[i]),file=fout)
    fout.close()

    ii1=np.where(el<=110)
    ii2=np.where(_el<=110)
    el1=el[ii1]; vv1=vv[ii1]
    el2=_el[ii2]; vv2=_vv[ii2]
    drawfig(el1,vv1,el2,vv2)

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

洪水波形作図プログラム

1時間ピッチで与えられた流量をそのままプロットしている。 Flood Routine における水位追跡の時間ピッチは、ここで指定されている洪水波形の時間ピッチとしている。

import numpy as np
import matplotlib.pyplot as plt


def drawfig(tt,qq):
    fsz=16
    xmin=0 ; xmax=132; dx=12
    ymin=0; ymax=3500; dy=500
    tstr='Hydrograph (20 years return period flood)'
    plt.figure(figsize=(10,6),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Time (hours)')
    plt.ylabel('Discharge (m$^3$/s)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#999999',linestyle='solid')
    plt.title(tstr,loc='left',fontsize=fsz)
    plt.plot(tt,qq,'-',color='#0000ff',lw=2)
    i=np.argmax(qq)
    plt.text(tt[i],qq[i],'Qmax={0:.0f}'.format(qq[i]),va='bottom',ha='center',color='#000000',fontsize=fsz)
    fnameF='fig_hydro.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.2)
    plt.show()


def main():
    fnameR='inp_hydro.txt'
    data = np.loadtxt(fnameR, usecols=(0,1) , skiprows=0)
    tt=data[:,0]
    qq=data[:,1]
    drawfig(tt,qq)

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

Thank you.

記事の先頭に行く

2019-05-29

設計仮排水路設計における Flood Routine の活用（２）

設計

記事の最後に行く

ここでは、「設計仮排水路設計における Flood Routine の活用（２）」で述べた事項に関連する計算・作図プログラムのうち、以下の２つを掲載する。

Flood Routine 解析プログラム
Flood Routine 作図プログラム

Flood Routine 解析プログラム

import numpy as np


def flood(elini,outlet,rh,rv,nr,ti,q_in,nd,elof,qof,no,fout):
    p_el=np.zeros(nd-1,dtype=np.float64)
    p_qo=np.zeros(nd-1,dtype=np.float64)
    # Initial outflow
    elv=elini
    el=elv
    vol=ret_v(nr,rh,rv,elv)
    q_out=ofc(elv,elof,qof,no)+outlet
    i=0
    iud=0
    print('{0:>5s} {1:>5s} {2:>10s} {3:>10s} {4:>15s} {5:>15s} {6:>15s} {7:>15s}'.format('i','iud','time','el','el-elv','vol','q_in','q_out'),file=fout)
    print('{0:5d} {1:5d} {2:10.3f} {3:10.3f} {4:15e} {5:15e} {6:15e} {7:15e}'.format(i,iud,ti[i],el,el-elv,vol,q_in[i],q_out),file=fout)
    # Iterative calculation
    iud=1
    dh=0.0005
    eps=0.001
    itmax=int(1.0/dh)*100
    for i in range(0,nd-1):
        qqin=0.5*(q_in[i]+q_in[i+1])
        tim=0.5*(ti[i+1]-ti[i])
        qin=0.5*(q_in[i]+q_in[i+1])*(ti[i+1]-ti[i])*3600.0
        hh=0.0
        k=0
        j=0
        while 1:
            k=k+1
            j=j+1;
            hh=hh+float(iud)*dh
            elv=el;    q1=ofc(elv,elof,qof,no)+outlet
            elv=el+hh; q2=ofc(elv,elof,qof,no)+outlet
            qout=0.5*(q1+q2)*(ti[i+1]-ti[i])*3600.0
            R=vol+qin-qout
            elv=ret_h(nr,rh,rv,R)
            if iud==1 and j==10:
                if el+hh>elv:
                    iud=-1
                    hh=0.0
                    j=0
            if iud==-1 and j==10:
                if el+hh<elv:
                    iud=1
                    hh=0.0
                    j=0
            if abs(el+hh-elv)<eps: break
            if itmax<k: break
        vol=R    # Cumulative volume
        el=el+hh # Elevation
        q_out=q2 # Outflow
        print('{0:5d} {1:5d} {2:10.3f} {3:10.3f} {4:15e} {5:15e} {6:15e} {7:15e}'.format(i+1,iud,ti[i+1],el,el-elv,vol,q_in[i+1],q_out),file=fout)
        p_el[i]=el
        p_qo[i]=q_out
    hmax=np.max(p_el)-elini
    print('Time: {0:10.3f}'.format(np.max(ti)))
    print('hmax: {0:10.3f}'.format(hmax))
    print('Qin : {0:10.3f} {1:10.3f}'.format(np.min(q_in),np.max(q_in)))
    print('Wout: {0:10.3f} {1:10.3f}'.format(np.min(p_qo),np.max(p_qo)))
    print('EL  : {0:10.3f} {1:10.3f}'.format(np.min(p_el),np.max(p_el)))
    del p_el,p_qo

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


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


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


def inp_hv(fnameR1):
    # Input H-V data
    vcoef=1e6
    data = np.loadtxt(fnameR1, usecols=(0,1) , skiprows=0)
    rh=data[:,0]
    rv=data[:,1]*vcoef
    nr=len(rh)
    return nr,rh,rv


def inp_inflow(fnameR2):
    # Input time sequence of inflow
    data = np.loadtxt(fnameR2, usecols=(0,1) , skiprows=0)
    ti  =data[:,0]
    q_in=data[:,1]
    nd=len(ti)
    elini=64.0
    outlet=0.0
    return nd,ti,q_in,elini,outlet


def inp_outflow(fnameR3):
    # Input outflow capacity
    data = np.loadtxt(fnameR3, usecols=(0,1) , skiprows=0)
    elof=data[:,0]
    qof =data[:,1]
    no=len(elof)
    return no,elof,qof


def main():
    #Main routine
    fnameR1='inp_hv.txt'
    fnameR2='inp_hydro.txt'
    nr,rh,rv=inp_hv(fnameR1)
    nd,ti,q_in,elini,outlet=inp_inflow(fnameR2)
    for ss in ['05','06','07','08','09','10','11','12']:
        fnameR3='inp_sec_'+ss+'.txt' # Outflow capacity
        fnameW ='out_'+ss+'.txt'     # Output file name
        no,elof,qof=inp_outflow(fnameR3)
        fout=open(fnameW,'w')
        flood(elini,outlet,rh,rv,nr,ti,q_in,nd,elof,qof,no,fout)
        fout.close()


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

Flood Routine 作図プログラム

import numpy as np
import matplotlib.pyplot as plt


fsz=16
xmin =0.0 ; xmax =132.0 ; dx =12    # time
ymin1=0.0 ; ymax1=5000.0; dy1=1000  # discharge
ymin2=60.0; ymax2=110.0 ; dy2=10    # water level


def drawfig(nnn,d0):
    ss='{0:0>2}'.format(d0)
    fnameR='out_'+ss+'.txt' # calculation result
    data=np.loadtxt(fnameR,skiprows=1,usecols=(2,3,6,7))
    ti=data[:,0]
    el=data[:,1]
    qi=data[:,2]
    qo=data[:,3]
    n1=np.argmax(qi)
    n2=np.argmax(qo)
    n3=np.argmax(el)
    plt.subplot(nnn)
    plt.xlim([xmin,xmax])
    plt.ylim([ymin1,ymax1])
    plt.xlabel('Time (hour)')
    plt.ylabel('Discharge (m$^3$/s)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin1,ymax1+dy1,dy1))
    plt.grid(color='#999999',linestyle='solid')
    plt.plot(ti,qi,'--',color='#000000',lw=2.0,label='Q (inflow)')
    plt.plot(ti,qo,'-.',color='#000000',lw=2.0,label='Q (outflow)')
    plt.text(ti[n1],qi[n1],'max:{0:.0f}'.format(qi[n1]),fontsize=fsz,color='#000000',ha='center',va='bottom')
    plt.text(ti[n2],qo[n2],'max:{0:.0f}'.format(qo[n2]),fontsize=fsz,color='#000000',ha='center',va='bottom')
    plt.twinx()
    plt.ylim([ymin2,ymax2])
    plt.ylabel('Water Level (EL.m)')
    plt.plot(ti,el,'-',color='#0000ff',lw=2.0,label='Water Level')
    plt.text(ti[n3],el[n3],'ELmax:{0:.3f}'.format(el[n3]),fontsize=fsz,color='#000000',ha='left',va='bottom')
    plt.text(xmin+5,ymax2-2,'D{0:.1f}m x 2'.format(d0),fontsize=fsz+4,color='#000000',ha='left',va='top')
    elti=64.0
    eltc=elti+d0
    plt.plot([xmin,xmax],[elti,elti],'-',color='#0000ff',lw=1)
    plt.plot([xmin,xmax],[eltc,eltc],'-',color='#0000ff',lw=1)
    plt.text(xmax-5,elti,'Tunnel Invert:EL.{0:.1f}'.format(elti),fontsize=fsz,color='#0000ff',ha='right',va='top')
    plt.text(xmax-5,eltc,'Tunnel Crown:EL.{0:.1f}'.format(eltc) ,fontsize=fsz,color='#0000ff',ha='right',va='bottom')
    plt.plot([0],[0],'--',color='#000000',lw=2.0,label='Q (inflow)')
    plt.plot([0],[0],'-.',color='#000000',lw=2.0,label='Q (outflow)')
    plt.legend(shadow=True,loc='upper right',handlelength=2,fontsize=fsz-2)
    return qo[n2],el[n3]
    

def main():
    # Main routine
    fnameF='fig_diversion.jpg' # image file
    ww=8; hh=5
    plt.figure(figsize=(ww*2,hh*4),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    l_d0=[5,6,7,8,9,10,11,12]
    qomax=[]
    elmax=[]
    nnn=420
    for d0 in l_d0:
        nnn=nnn+1      
        qo,el=drawfig(nnn,d0)
        qomax=qomax+[qo]
        elmax=elmax+[el]

    print(qomax)
    print(elmax)
    plt.tight_layout()        
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.2)
    plt.show()

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

Thank you.

記事の先頭に行く

2019-05-29

設計仮排水路設計における Flood Routine の活用（１）

設計

記事の最後に行く

概要

仮排水路トンネル断面寸法と仮締切堤高さを決定するための設計検討を行う。具体的には、トンネル径を変化させ、Flood Routine により、洪水流入波形、流出流量、貯水池水位の関係を求める。

Conditions for Flood Routine Analysis

For flood routine analyssis, following conditions are required.

Hydrograph as unflow characteristic.
Reservoir storage capacity curve as storage characteristics
Discharge capacity curve of diversion tunnel as outfow characteristics

Hydrograph

Following hydrograph for 20 years return period flood is used as a given condition.

f:id:damyarou:20190529110557j:plain

Reservoir Storage Capacity Curve

Following reservoir capacity curve is used as a given condition.

f:id:damyarou:20190529112722j:plain

Discharge Capacity Curve of Diversion Tunnel

In this project, two lanes diversion tunnels are planed, and the characteristics of alignments and standard tunnel shape are shown below.

Information of Diversion Tunnel Alignment
Items	Inside	Outside
Length	L=400m	L=450m
Invert level of entrance	EL.64.0m	EL.64.0m
Invert level of exit	EL.62.5m	EL.62.5m
Gradient	i=0.00375	i=0.00333
Curve-1	Bending radius: $\rho$ =230m Bending angle: $\theta$ =39.5deg.	Bending radius: $\rho$ =260m Bending angle: $\theta$ =39.5deg.
Curve-2	Bending radius: $\rho$ =200m Bending angle: $\theta$ =48.2deg.	Bending radius: $\rho$ =230m Bending angle: $\theta$ =48.2deg.
The distance between tunnels are set to 30m (center to center)

f:id:damyarou:20190529110509j:plain

The assumptions for the calculation of flow capacity curve of diversion tunnels are shown below.

In case that the reservoir water level is lower than or equal to the crown level of the diversion tunnels, the water flows down with uniform flow state.
In case that the reservoir water levei is higher than the crown level of diversion tunnels, the water flows down with pressured tunnel flow state.
The discharge capacity of diversion is defined as the sum of the discharges of two tunnels.

Calculation of Discharge Capacity for Uniform Flow

$\begin{equation} Q=A\cdot v \qquad v=\cfrac{1}{n}\cdot R^{2/3}\cdot I^{1/2} \end{equation}$

$Q$	discharge
$A$	flow section area
$v$	average velocity
$n$	manning's roughness coefficient (=0.015)
$R$	hydraulic radius
$I$	invert gradient

Calculation of Discharge Capacity for Pressued Tunnel Flow

$\begin{gather} Q=A\cdot v \qquad v=\sqrt{\cfrac{2 g \Delta H}{f_e+f_o+f_{b1}+f_{b2}+f\cdot L\;/\;R}} \\ f_b=\left\{0.131+0.1632\cdot\left(\cfrac{D}{\rho}\right)^{7/2}\right\}\cdot\sqrt{\cfrac{\theta}{90}} \\ f=\cfrac{2 g n^2}{R^{1/3}} \end{gather}$

$Q$	discharge
$A$	section area of diversion tunnel
$v$	average velocity
$\Delta H$	difference of water head between upstream and downstream of diversion tunnel
$f_e$	head loss coefficient of entrance (=0.25)
$f_o$	head loss coefficient of exit (=1.00)
$f_{b1}$	head loss coefficient of bending for curve-1
$f_{b2}$	head loss coefficient of bending for curve-2
$f$	friction head loss coefficient
$n$	Manning's roughness coefficient (=0.015)
$L$	length of tunnel
$R$	hydraulic radius of tunnel
$D$	internal diameter of tunnel
$\rho$	radius of curvature of bending
$\theta$	inter angle of bending
$g$	gravity acceleration (=9.8 m/s $^2$ )

The calculated discharge capacity of diversion tunnels are shown below.

f:id:damyarou:20190529112757j:plain

Result of Flood Routine Analysis

The result of flood routine analysis with parameters of tunnel diameter from 5m to 12m is shown below.

Summary of Flood Routine Analysis Result (Max. Inflow: 3130 m3/s)
Tunnel	D5.0m x 2	D6.0m x 2	D7.0m x 2	D8.0m x 2	D9.0m x 2	D10.0m x 2	D11.0m x 2	D12.0m x 2
Max. Outflow (m3/s)	660	959	1287	1628	1964	2278	2548	2756
Max. Water Level (EL.m)	99.641	97.061	94.262	91.364	88.503	85.735	83.079	80.618

f:id:damyarou:20190529110714j:plain

Thank you.

記事の先頭に行く

2019-05-27

matplotlib RC円形圧力トンネルモデル図

matplotlib

記事の最後に行く

はじめに

以下に示す、「設計 RC円形圧力トンネルの配筋設計（１）」で示した図を作成するプログラムである。

f:id:damyarou:20190527132711j:plain

ポイント

アノーテーション用矢印とテキスト描画

説明用のボックス内テキストと矢印を描画する。
matplotlib では、annotate により、矢印とテキストを同時に描画できるが、矢印とテキストの関係が思うようにいかないため、ここでは、矢印とボックス内テキストを別々に描画し、連結させている。
矢印の色は濃い灰色とし、線の太さは細目に設定。
bbox_props を定義し、plt.text の中で bbox=bbox_props とすることにより定義したボックス内に文字を描画する。
この事例では、描画する文字列が２行にわたっているが、行間を詰めるため、plt.text の中で linespacing=1 を指定している。
引数の説明は以下の通り。

x1, y1	矢印先端座標
x2, y2	矢印線始点座標（文字列ボックスとの接点座標）
ss	描画文字列
lstr	ボックスから出る矢印の位置（文字列 l, t, r, b, n のいずれかで指定）
fszz	描画文字列のフォントサイズ

lstr は、１文字で指定し、その意味は以下の通り。

'l'	矢印線はボックスの左からスタート
't'	矢印線はボックスの上からスタート
'r'	矢印線はボックスの右からスタート
'b'	矢印線はボックスの下からスタート
'n'	矢印は描かず、(x1, y1)、(x2, y2)で指定した２点の中央にボックス入り文字列を描画

def sarrow_a(x1,y1,x2,y2,ss,lstr,fszz):
    # arrow for annotation
    if lstr=='n':
        x2=0.5*(x1+x2)
        y2=0.5*(y1+y2)
    else:
        col='#777777'
        sv=0
        aprop=dict(shrink=sv,width=0.3,headwidth=4,headlength=8,
                   connectionstyle='arc3',facecolor=col,edgecolor=col)
        plt.annotate('',
            xy=(x1,y1), xycoords='data',
            xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)
    # text drawing
    bbox_props = dict(boxstyle='square,pad=0.1', fc='#ffffff', ec='#777777', lw=1)
    if lstr=='l':
        hstr='right'; vstr='center'
    if lstr=='t':
        hstr='center'; vstr='top'
    if lstr=='r':
        hstr='left'; vstr='center'
    if lstr=='b':
        hstr='center'; vstr='bottom'
    if lstr=='n':
        hstr='center'; vstr='center'
    plt.text(x2,y2,ss,ha=hstr,va=vstr,fontsize=fszz,bbox=bbox_props,linespacing=1)

荷重用矢印

矢印の色は黒とし、線は太めに設定。

def sarrow_p(x1,y1,x2,y2):
    # arrow for load
    col='#000000'
    sv=0
    aprop=dict(shrink=sv,width=1,headwidth=5,headlength=8,
               connectionstyle='arc3',facecolor=col,edgecolor=col)
    plt.annotate('',
        xy=(x1,y1), xycoords='data',
        xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)

ひび割れ描画

コンクリート内のひび割れは、サインカーブで定義。まず、座標 $(r_a, 0)$ - $(r_b, 0)$ 間でサインカーブを定義し、これを座標変換しながら anga で定めた角度に描画していく。

def crack(ra,rb):
    ds=0.2
    m=4
    ell=(rb-ra)/m
    x0=np.linspace(ra,rb,101)
    y0=ds*np.sin(2*np.pi/ell*x0)
    anga=np.linspace(0,2*np.pi,13)
    for ang in anga:
        x=x0*np.cos(ang)-y0*np.sin(ang)
        y=x0*np.sin(ang)+y0*np.cos(ang)
        plt.plot(x,y,'-',color='#999999')

塗りつぶし

ハッチングの場合

color でハッチングの色を指定できる。ここでは濃い灰色を指定。ハッチングはクロス線で行う。

plt.fill(xr,yr,fill=False, color='#999999',hatch='xx',lw=0)

単色で塗りつぶす場合

facecolor で塗りつぶす領域の色を、edgecolor で境界の色を指定する。

plt.fill(xb,yb,facecolor='#eeeeee',edgecolor='#000000',lw=1)

円を描く

np.linspace で、 $0$ から $2\pi$ を分割して角度を指定して $(x, y)$ 座標を求め、これを折れ線で結んでいく。

angc=np.linspace(0,2*np.pi,181) # for concrete outline
xfa=rfa*np.cos(angc)
yfa=rfa*np.sin(angc)
plt.plot(xfa,yfa,color='#000000',lw=3)

太字の描画

この場合はタイトル描画のために使っているが、fontweight='bold' で太字を描画できる。

https://matplotlib.org/examples/pylab_examples/fonts_demo.html

plt.title(tstr,loc='center',fontsize=fsz,fontweight='bold')

プログラム全文

import numpy as np
import matplotlib.pyplot as plt


# Global variables
fsz=12
xmin=-10 ; xmax=10; dx=5
ymin=-10; ymax=10; dy=5


def sarrow_a(x1,y1,x2,y2,ss,lstr,fszz):
    # arrow for annotation
    if lstr=='n':
        x2=0.5*(x1+x2)
        y2=0.5*(y1+y2)
    else:
        col='#777777'
        sv=0
        aprop=dict(shrink=sv,width=0.3,headwidth=4,headlength=8,
                   connectionstyle='arc3',facecolor=col,edgecolor=col)
        plt.annotate('',
            xy=(x1,y1), xycoords='data',
            xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)
    # text drawing
    bbox_props = dict(boxstyle='square,pad=0.1', fc='#ffffff', ec='#777777', lw=1)
    if lstr=='l':
        hstr='right'; vstr='center'
    if lstr=='t':
        hstr='center'; vstr='top'
    if lstr=='r':
        hstr='left'; vstr='center'
    if lstr=='b':
        hstr='center'; vstr='bottom'
    if lstr=='n':
        hstr='center'; vstr='center'
    plt.text(x2,y2,ss,ha=hstr,va=vstr,fontsize=fszz,bbox=bbox_props,linespacing=1)

    
def sarrow_p(x1,y1,x2,y2):
    # arrow for load
    col='#000000'
    sv=0
    aprop=dict(shrink=sv,width=1,headwidth=5,headlength=8,
               connectionstyle='arc3',facecolor=col,edgecolor=col)
    plt.annotate('',
        xy=(x1,y1), xycoords='data',
        xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)

    
def crack(ra,rb):
    ds=0.2
    m=4
    ell=(rb-ra)/m
    x0=np.linspace(ra,rb,101)
    y0=ds*np.sin(2*np.pi/ell*x0)
    anga=np.linspace(0,2*np.pi,13)
    for ang in anga:
        x=x0*np.cos(ang)-y0*np.sin(ang)
        y=x0*np.sin(ang)+y0*np.cos(ang)
        plt.plot(x,y,'-',color='#999999')    

        
def drawfig(nnn):
    ra=4 # inner surface
    rb=7 # outer surface
    pl=1 # length of load arrow
    br=rb+1.5 # bedrock area
    rfa=ra+0.6 # inner reinforcement
    rfb=rb-0.6 # outer reinforcement
    angc=np.linspace(0,2*np.pi,181) # for concrete outline
    angp=np.linspace(0,2*np.pi,19) # for load arrow
    
    if nnn==121: tstr='RC Pressure Tunnel under Internal Pressure\n(Double reinforcement model)'
    if nnn==122: tstr='RC Pressure Tunnel under External Pressure\n(Double reinforcement model)'
    plt.subplot(nnn)
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.axis('off')
    plt.gca().set_aspect('equal',adjustable='box')
    plt.title(tstr,loc='center',fontsize=fsz,fontweight='bold')
    if nnn==121:
        # hatching of bedrock area
        xr=br*np.cos(angc)
        yr=br*np.sin(angc)
        plt.fill(xr,yr,fill=False, color='#999999',hatch='xx',lw=0)
    # outer surface line of concrete
    xb=rb*np.cos(angc)
    yb=rb*np.sin(angc)
    plt.fill(xb,yb,facecolor='#eeeeee',edgecolor='#000000',lw=1)
    # inner surface line of concrete
    xa=ra*np.cos(angc)
    ya=ra*np.sin(angc)
    plt.fill(xa,ya,facecolor='#ffffff',edgecolor='#000000',lw=1)
    # drawing cracks
    if nnn==121: crack(ra,rb)
    # inner reinforcement
    xfa=rfa*np.cos(angc)
    yfa=rfa*np.sin(angc)
    plt.plot(xfa,yfa,color='#000000',lw=3)
    # outer reinforcement
    xfb=rfb*np.cos(angc)
    yfb=rfb*np.sin(angc)
    plt.plot(xfb,yfb,color='#000000',lw=3)
    # pressure load drawing
    if nnn==121:
        r1=ra; r2=ra-pl
    if nnn==122:
        r1=rb; r2=rb+pl
    xp1=r1*np.cos(angp) # x-coordinate of end of arrow line
    yp1=r1*np.sin(angp) # y-coordinate of end of arrow line
    xp2=r2*np.cos(angp) # x-coordinate of start of arrow line
    yp2=r2*np.sin(angp) # y-coordinate of start of arrow line
    for x1,y1,x2,y2 in zip(xp1,yp1,xp2,yp2):
        sarrow_p(x1,y1,x2,y2)
    if nnn==121: plt.text(ra-pl-0.2,0,'$P_a$',va='center',ha='right',fontsize=fsz+4)
    if nnn==122: plt.text(rb+pl+0.2,0,'$P_b$',va='center',ha='left',fontsize=fsz+4)
    # explanation
    if nnn==121: ss='Concrete\nwith crack'
    if nnn==122: ss='Concrete\nwithout crack'
    x1=0.5*(ra+rb)*np.cos(np.radians(135)); x2=x1-3.5
    y1=0.5*(ra+rb)*np.sin(np.radians(135)); y2=y1+3.5
    sarrow_a(x1,y1,x2,y2,ss,'b',fsz)
    ss='Outer\nReinforcement'
    x1=rfb*np.cos(np.radians(225)); x2=x1-3
    y1=rfb*np.sin(np.radians(225)); y2=y1-3
    sarrow_a(x1,y1,x2,y2,ss,'t',fsz)
    ss='Inner\nReinforcement'
    x1=rfa*np.cos(np.radians(315)); x2=x1+4
    y1=rfa*np.sin(np.radians(315)); y2=y1-4
    sarrow_a(x1,y1,x2,y2,ss,'t',fsz)
    if nnn==121: 
        ss='Bedrock'
        x1=0 ; x2=x1
        y1=br; y2=y1
        sarrow_a(x1,y1,x2,y2,ss,'n',fsz)
    
        
def main():
    plt.figure(figsize=(10,5),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    nnn=121; drawfig(nnn)
    nnn=122; drawfig(nnn)
    plt.tight_layout()
    fnameF='fig_model.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()


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

Thank you.

記事の先頭に行く

2019-05-26

設計 RC円形圧力トンネルの配筋設計（２）

設計

記事の最後に行く

内水圧を受けるRC圧力トンネル

出力

ra	覆工内半径 (mm)
rb	覆工外半径 (mm)
r0	岩盤外縁半径位 (mm)
ta	内側鉄筋等価板厚 (mm)
tb	外側鉄筋等価板厚 (mm)
da	内側鉄筋かぶり (mm)
db	外側鉄筋かぶり (mm)
Pa	内水圧 (MPa)
T	温度変化量（マイナスは温度低下）
Loc.	場所
u_r	半径方向変位 (mm)
sig_r	半径方向直応力 (MPa)
sig_t	円周方向直応力 (MPa)
sig_z	トンネル軸方向直応力 (MPa)

* Double Reinforcement
    ra     rb     r0     ta     tb     da     db     Pa      T
  4000   4800 100000  4.021  4.021    100    100  1.000 -10.000
Loc.           u_r      sig_r      sig_t      sig_z
r7_rock      0.225     -0.000      0.002      0.001
r6_rock      3.123     -0.519      0.521      0.001
r6_conc      3.123     -0.519      0.000     -0.104
r5_conc      3.135     -0.530      0.000     -0.106
r5_rbar      3.135     -0.530    174.965     52.331
r4_rbar      3.137     -0.680    175.116     52.331
r4_conc      3.137     -0.680      0.000     -0.136
r3_conc      3.214     -0.778      0.000     -0.156
r3_rbar      3.214     -0.778    200.345     59.870
r2_rbar      3.216     -0.976    200.542     59.870
r2_conc      3.216     -0.976      0.000     -0.195
r1_conc      3.230     -1.000      0.000     -0.200
* Single Reinforcement
    ra     rb     r0     ta     tb     da     db     Pa      T
  4000   4600 100000  4.021      0    100      0  1.000 -10.000
Loc.           u_r      sig_r      sig_t      sig_z
r5_rock      0.264      0.000      0.003      0.001
r4_rock      3.823     -0.663      0.666      0.001
r4_conc      3.823     -0.663      0.000     -0.133
r3_conc      3.887     -0.743      0.000     -0.149
r3_rbar      3.887     -0.743    236.399     70.697
r2_rbar      3.889     -0.976    236.632     70.697
r2_conc      3.889     -0.976      0.000     -0.195
r1_conc      3.903     -1.000      0.000     -0.200

プログラム

# Displacements and Stresses
# of RC Circular Tunnel under Internal Pressure
import numpy as np


def bdr(Eg,ng,r,cg1,cg2):
    # disp. and stresses of bedrock
    u=cg1*r+cg2/r
    sr=Eg/(1+ng)/(1-2*ng)*cg1-Eg/(1+ng)*cg2/r**2
    st=Eg/(1+ng)/(1-2*ng)*cg1+Eg/(1+ng)*cg2/r**2
    sz=ng*(sr+st)
    return u,sr,st,sz


def con(Ec,nc,alpc,r,r0,tt,cc1,cc2):
    # disp. and stresses of concrete (no-tension)
    u=cc1+cc2*np.log(r)+alpc*tt*(r-r0)
    sr=Ec*cc2/r
    st=0
    sz=nc*(sr+st)
    return u,sr,st,sz


def reb(Es,ns,alps,r,r0,tt,cs1,cs2):
    # disp. and stresses of reinforcement
    u=(1+ns)/(1-ns)*alps*tt*(r**2-r0**2)/2/r+cs1*r+cs2/r
    sr=-Es*alps*tt/(1-ns)*(r**2-r0**2)/2/r**2+Es/(1+ns)/(1-2*ns)*cs1-Es/(1+ns)*cs2/r**2
    st=-Es*alps*tt/(1-ns)*(r**2+r0**2)/2/r**2+Es/(1+ns)/(1-2*ns)*cs1+Es/(1+ns)*cs2/r**2
    sz=ns*(sr+st)
    return u,sr,st,sz


def pi_calc_d(rr,Eg,ng,Es,ns,alps,Ec,nc,alpc,tt,pa):
    # double reinforcement
    r7=rr[6] # outer boundary of bedrock
    r6=rr[5] # external surface of concrete lining
    r5=rr[4] # boundary of outer cover concrete and outer reinforcement
    r4=rr[3] # boundary of outer reinforcement and middle concrete
    r3=rr[2] # boundary of middle concrete and inner reinforcement
    r2=rr[1] # boundary of inner reinforcement and inner cover concrete
    r1=rr[0] # internal surface of concrete lining
    
    n=12
    aa=np.zeros((n,n),dtype=np.float64)
    bb=np.zeros(n,dtype=np.float64)
    aa[ 0, 0]=Eg/(1+ng)/(1-2*ng); aa[ 0,1]=-Eg/(1+ng)/r7**2
    aa[ 1, 0]=r6                ; aa[ 1,1]=1/r6            ; aa[ 1, 2]=-1                 ; aa[ 1, 3]=-np.log(r6)
    aa[ 2, 0]=Eg/(1+ng)/(1-2*ng); aa[ 2,1]=-Eg/(1+ng)/r6**2; aa[ 2, 2]=0                  ; aa[ 2, 3]=-Ec/r6
    aa[ 3, 2]=1                 ; aa[ 3,3]=np.log(r5)      ; aa[ 3, 4]=-r5                ; aa[ 3, 5]=-1/r5
    aa[ 4, 2]=0                 ; aa[ 4,3]=Ec/r5           ; aa[ 4, 4]=-Es/(1+ns)/(1-2*ns); aa[ 4, 5]=Es/(1+ns)/r5**2
    aa[ 5, 4]=r4                ; aa[ 5,5]=1/r4            ; aa[ 5, 6]=-1                 ; aa[ 5, 7]=-np.log(r4)
    aa[ 6, 4]=Es/(1+ns)/(1-2*ns); aa[ 6,5]=-Es/(1+ns)/r4**2; aa[ 6, 6]=0                  ; aa[ 6, 7]=-Ec/r4
    aa[ 7, 6]=1                 ; aa[ 7,7]=np.log(r3)      ; aa[ 7, 8]=-r3                ; aa[ 7, 9]=-1/r3
    aa[ 8, 6]=0                 ; aa[ 8,7]=Ec/r3           ; aa[ 8, 8]=-Es/(1+ns)/(1-2*ns); aa[ 8, 9]=Es/(1+ns)/r3**2
    aa[ 9, 8]=r2                ; aa[ 9,9]=1/r2            ; aa[ 9,10]=-1                 ; aa[ 9,11]=-np.log(r2)
    aa[10, 8]=Es/(1+ns)/(1-2*ns); aa[10,9]=-Es/(1+ns)/r2**2; aa[10,10]=0                  ; aa[10,11]=-Ec/r2
    aa[11,11]=Ec/r1
    
    bb[1]=alpc*tt*(r6-r5)
    bb[3]=(1+ns)/(1-ns)*alps*tt*(r5**2-r4**2)/2/r5
    bb[4]=-Es*alps*tt/(1-ns)*(r5**2-r4**2)/2/r5**2
    bb[5]=alpc*tt*(r4-r3)
    bb[7]=(1+ns)/(1-ns)*alps*tt*(r3**2-r2**2)/2/r3
    bb[8]=-Es*alps*tt/(1-ns)*(r3**2-r2**2)/2/r3**2
    bb[9]=alpc*tt*(r2-r1)
    bb[11]=-pa

    cc = np.linalg.solve(aa, bb)

    uu=np.zeros(n,dtype=np.float64)
    sr=np.zeros(n,dtype=np.float64)
    st=np.zeros(n,dtype=np.float64)
    sz=np.zeros(n,dtype=np.float64)
    ss=[]
    ss=ss+['r7_rock']
    ss=ss+['r6_rock']
    ss=ss+['r6_conc']
    ss=ss+['r5_conc']
    ss=ss+['r5_rbar']
    ss=ss+['r4_rbar']
    ss=ss+['r4_conc']
    ss=ss+['r3_conc']
    ss=ss+['r3_rbar']
    ss=ss+['r2_rbar']
    ss=ss+['r2_conc']
    ss=ss+['r1_conc']
    i=0;  uu[i],sr[i],st[i],sz[i]=bdr(Eg,ng,rr[6],cc[0],cc[1])
    i=1;  uu[i],sr[i],st[i],sz[i]=bdr(Eg,ng,rr[5],cc[0],cc[1])
    i=2;  uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[5],rr[4],tt,cc[2],cc[3])
    i=3;  uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[4],rr[4],tt,cc[2],cc[3])
    i=4;  uu[i],sr[i],st[i],sz[i]=reb(Es,ns,alps,rr[4],rr[3],tt,cc[4],cc[5])
    i=5;  uu[i],sr[i],st[i],sz[i]=reb(Es,ns,alps,rr[3],rr[3],tt,cc[4],cc[5])
    i=6;  uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[3],rr[2],tt,cc[6],cc[7])
    i=7;  uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[2],rr[2],tt,cc[6],cc[7])
    i=8;  uu[i],sr[i],st[i],sz[i]=reb(Es,ns,alps,rr[2],rr[1],tt,cc[8],cc[9])
    i=9;  uu[i],sr[i],st[i],sz[i]=reb(Es,ns,alps,rr[1],rr[1],tt,cc[8],cc[9])
    i=10; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[1],rr[0],tt,cc[10],cc[11])
    i=11; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[0],rr[0],tt,cc[10],cc[11])
    print('* Double Reinforcement')
    print('{0:>6s} {1:>6s} {2:>6s} {3:>6} {4:>6s} {5:>6s} {6:>6s} {7:>6s} {8:>6s}'
          .format('ra','rb','r0','ta','tb','da','db','Pa','T'))
    print('{0:6.0f} {1:6.0f} {2:6.0f} {3:6.3f} {4:6.3f} {5:6.0f} {6:6.0f} {7:6.3f} {8:6.3f}'
          .format(r1,r6,r7,r3-r2,r5-r4,r2-r1,r6-r5,pa,tt))
    print('{0:7s} {1:>10s} {2:>10s} {3:>10s} {4:>10s}'.format('Loc.','u_r','sig_r','sig_t','sig_z'))
    for i in range(n):
        print('{0:7s} {1:10.3f} {2:10.3f} {3:10.3f} {4:10.3f}'.format(ss[i],uu[i],sr[i],st[i],sz[i]))

        
def pi_calc_s(rr,Eg,ng,Es,ns,alps,Ec,nc,alpc,tt,pa):
    # single reinforcement
    r5=rr[4] # outer boundary of bedrock
    r4=rr[3] # external surface of concrete lining
    r3=rr[2] # boundary of outer cover concrete and reinforcement
    r2=rr[1] # boundary of reinforcement and inner cover concrete
    r1=rr[0] # inner surface of concrete lining
    
    n=8
    aa=np.zeros((n,n),dtype=np.float64)
    bb=np.zeros(n,dtype=np.float64)
    aa[0,0]=Eg/(1+ng)/(1-2*ng); aa[0,1]=-Eg/(1+ng)/r5**2
    aa[1,0]=r4                ; aa[1,1]=1/r4            ; aa[1,2]=-1                 ; aa[1,3]=-np.log(r4)
    aa[2,0]=Eg/(1+ng)/(1-2*ng); aa[2,1]=-Eg/(1+ng)/r4**2; aa[2,2]=0                  ; aa[2,3]=-Ec/r4
    aa[3,2]=1                 ; aa[3,3]=np.log(r3)      ; aa[3,4]=-r3                ; aa[3,5]=-1/r3
    aa[4,2]=0                 ; aa[4,3]=Ec/r3           ; aa[4,4]=-Es/(1+ns)/(1-2*ns); aa[4,5]=Es/(1+ns)/r3**2
    aa[5,4]=r2                ; aa[5,5]=1/r2            ; aa[5,6]=-1                 ; aa[5,7]=-np.log(r2)
    aa[6,4]=Es/(1+ns)/(1-2*ns); aa[6,5]=-Es/(1+ns)/r2**2; aa[6,6]=0                  ; aa[6,7]=-Ec/r2
    aa[7,7]=Ec/r1
    
    bb[1]=alpc*tt*(r4-r3)
    bb[3]=(1+ns)/(1-ns)*alps*tt*(r3**2-r2**2)/2/r3
    bb[4]=-Es*alps*tt/(1-ns)*(r3**2-r2**2)/2/r3**2
    bb[5]=alpc*tt*(r2-r1)
    bb[7]=-pa

    cc = np.linalg.solve(aa, bb)

    uu=np.zeros(n,dtype=np.float64)
    sr=np.zeros(n,dtype=np.float64)
    st=np.zeros(n,dtype=np.float64)
    sz=np.zeros(n,dtype=np.float64)
    ss=[]
    ss=ss+['r5_rock']
    ss=ss+['r4_rock']
    ss=ss+['r4_conc']
    ss=ss+['r3_conc']
    ss=ss+['r3_rbar']
    ss=ss+['r2_rbar']
    ss=ss+['r2_conc']
    ss=ss+['r1_conc']
    i=0; uu[i],sr[i],st[i],sz[i]=bdr(Eg,ng,rr[4],cc[0],cc[1])
    i=1; uu[i],sr[i],st[i],sz[i]=bdr(Eg,ng,rr[3],cc[0],cc[1])
    i=2; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[3],rr[2],tt,cc[2],cc[3])
    i=3; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[2],rr[2],tt,cc[2],cc[3])
    i=4; uu[i],sr[i],st[i],sz[i]=reb(Es,ns,alps,rr[2],rr[1],tt,cc[4],cc[5])
    i=5; uu[i],sr[i],st[i],sz[i]=reb(Es,ns,alps,rr[1],rr[1],tt,cc[4],cc[5])
    i=6; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[1],rr[0],tt,cc[6],cc[7])
    i=7; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[0],rr[0],tt,cc[6],cc[7])
    print('* Single Reinforcement')
    print('{0:>6s} {1:>6s} {2:>6s} {3:>6} {4:>6s} {5:>6s} {6:>6s} {7:>6s} {8:>6s}'
          .format('ra','rb','r0','ta','tb','da','db','Pa','T'))
    print('{0:6.0f} {1:6.0f} {2:6.0f} {3:6.3f} {4:6.0f} {5:6.0f} {6:6.0f} {7:6.3f} {8:6.3f}'
          .format(r1,r4,r5,r3-r2,0,r2-r1,0,pa,tt))
    print('{0:7s} {1:>10s} {2:>10s} {3:>10s} {4:>10s}'.format('Loc.','u_r','sig_r','sig_t','sig_z'))
    for i in range(n):
        print('{0:7s} {1:10.3f} {2:10.3f} {3:10.3f} {4:10.3f}'.format(ss[i],uu[i],sr[i],st[i],sz[i]))
        

def main():
    # properties of materials
    Es=200000   # elastic modulus of steel
    ns=0.3      # Poisson's ratio of steel
    alps=1.0e-5 # thermal expansion coefficient of steel
    Ec=25000    # elastic modulus of concrete
    nc=0.2      # Poisson's ratio of concrete
    alpc=1.0e-5 # thermal expansion coefficient of concrete
    Eg=1000    # elastic modulus of bedrock
    ng=0.25     # Poisson's ratio of bedrock
    T32=32**2*np.pi/4
    T25=25**2*np.pi/4
    T20=20**2*np.pi/4

    # double reinforcement
    ro=100000.0 # Outer boundary of bedrock (mm)
    ra=4000.0 # internal radius of pressure tunnel (mm)
    rb=4800.0 # external radius of pressure tunnel (mm)
    ta=T32*5/1000   # equivalent steel thickness (mm)
    tb=T32*5/1000   # equivalent steel thickness (mm)
    da=100    # cover fro internal surface of concrete (mm)
    db=100    # cover from external surface of concrete (mm)
    tt=-10    # temperature change (negative: temperature drop)
    pa=1.0    # internal pressure (MPa)
    rr=np.array([ra,ra+da,ra+da+ta,rb-db-ta,rb-db,rb,ro])
    pi_calc_d(rr,Eg,ng,Es,ns,alps,Ec,nc,alpc,tt,pa)

    # single reinforcement
    ro=100000.0 # Outer boundary of bedrock (mm)
    ra=4000.0 # internal radius of pressure tunnel (mm)
    rb=4600.0 # external radius of pressure tunnel (mm)
    ta=T32*5/1000   # equivalent steel thickness (mm)
    da=100    # cover fro internal surface of concrete (mm)
    tt=-10    # temperature change (negative: temperature drop)
    pa=1.0    # internal pressure (MPa)
    rr=np.array([ra,ra+da,ra+da+ta,rb,ro])
    pi_calc_s(rr,Eg,ng,Es,ns,alps,Ec,nc,alpc,tt,pa)
    
    
#==============
# Execution
#==============
if __name__ == '__main__': main()

外水圧を受けるRC圧力トンネル

出力

ra	覆工内半径 (mm)
rb	覆工外半径 (mm)
ta	内側鉄筋等価板厚 (mm)
tb	外側鉄筋等価板厚 (mm)
da	内側鉄筋かぶり (mm)
db	外側鉄筋かぶり (mm)
Pb	外水圧 (MPa)
T	温度変化量（マイナスは温度低下）
Loc.	場所
u_r	半径方向変位 (mm)
sig_r	半径方向直応力 (MPa)
sig_t	円周方向直応力 (MPa)
sig_z	トンネル軸方向直応力 (MPa)

* Double Reinforcement
    ra     rb     r0     ta     tb     da     db     Pb      T
  4000   4800      0  4.021  4.021    100    100  1.000  0.000
Loc.           u_r      sig_r      sig_t      sig_z
r6_conc     -0.912     -1.000     -5.199     -1.240
r5_conc     -0.914     -0.910     -5.289     -1.240
r5_rbar     -0.914     -0.910    -40.722     -8.326
r4_rbar     -0.914     -0.876    -40.756     -8.326
r4_conc     -0.914     -0.876     -5.286     -1.232
r3_conc     -0.933     -0.194     -5.968     -1.232
r3_rbar     -0.933     -0.194    -47.406     -9.520
r2_rbar     -0.933     -0.147    -47.453     -9.520
r2_conc     -0.933     -0.147     -5.964     -1.222
r1_conc     -0.939      0.000     -6.111     -1.222
* Single Reinforcement
    ra     rb     r0     ta     tb     da     db     Pb      T
  4000   4800      0  4.021      0    100      0  1.000  0.000
Loc.           u_r      sig_r      sig_t      sig_z
r4_conc     -0.940     -1.000     -5.351     -1.270
r3_conc     -0.962     -0.200     -6.152     -1.270
r3_rbar     -0.962     -0.200    -48.865     -9.813
r2_rbar     -0.962     -0.152    -48.913     -9.813
r2_conc     -0.962     -0.152     -6.147     -1.260
r1_conc     -0.968      0.000     -6.299     -1.260

プログラム

# Displacements and Stresses
# of RC Circular Tunnel under External Pressure
import numpy as np


def con(Ec,nc,alpc,r,r0,tt,cc1,cc2):
    # disp. and stresses of concrete
    u=(1+nc)/(1-nc)*alpc*tt*(r**2-r0**2)/2/r+cc1*r+cc2/r
    sr=-Ec*alpc*tt/(1-nc)*(r**2-r0**2)/2/r**2+Ec/(1+nc)/(1-2*nc)*cc1-Ec/(1+nc)*cc2/r**2
    st=-Ec*alpc*tt/(1-nc)*(r**2+r0**2)/2/r**2+Ec/(1+nc)/(1-2*nc)*cc1+Ec/(1+nc)*cc2/r**2
    sz=nc*(sr+st)
    return u,sr,st,sz


def reb(Es,ns,alps,r,r0,tt,cs1,cs2):
    # disp. and stresses of reinforcement
    u=(1+ns)/(1-ns)*alps*tt*(r**2-r0**2)/2/r+cs1*r+cs2/r
    sr=-Es*alps*tt/(1-ns)*(r**2-r0**2)/2/r**2+Es/(1+ns)/(1-2*ns)*cs1-Es/(1+ns)*cs2/r**2
    st=-Es*alps*tt/(1-ns)*(r**2+r0**2)/2/r**2+Es/(1+ns)/(1-2*ns)*cs1+Es/(1+ns)*cs2/r**2
    sz=ns*(sr+st)
    return u,sr,st,sz


def pe_calc_d(rr,Es,ns,alps,Ec,nc,alpc,tt,pb):
    # double reinforcement
    r6=rr[5] # external surface of concrete lining
    r5=rr[4] # boundary of outer cover concrete and outer reinforcement
    r4=rr[3] # boundary of outer reinforcement and middle concrete
    r3=rr[2] # boundary of middle concrete and inner reinforcement
    r2=rr[1] # boundary of inner reinforcement and inner cover concrete
    r1=rr[0] # inner surface of concrete lining
    
    n=10
    aa=np.zeros((n,n),dtype=np.float64)
    bb=np.zeros(n,dtype=np.float64)
    aa[0,0]=Ec/(1+nc)/(1-2*nc); aa[0,1]=-Ec/(1+nc)/r6**2
    aa[1,0]=r5                ; aa[1,1]=1/r5            ; aa[1,2]=-r5                ; aa[1,3]=-1/r5
    aa[2,0]=Ec/(1+nc)/(1-2*nc); aa[2,1]=-Ec/(1+nc)/r5**2; aa[2,2]=-Es/(1+ns)/(1-2*ns); aa[2,3]=Es/(1+ns)/r5**2
    aa[3,2]=r4                ; aa[3,3]=1/r4            ; aa[3,4]=-r4                ; aa[3,5]=-1/r4
    aa[4,2]=Es/(1+ns)/(1-2*ns); aa[4,3]=-Es/(1+ns)/r4**2; aa[4,4]=-Ec/(1+nc)/(1-2*nc); aa[4,5]=Ec/(1+nc)/r4**2
    aa[5,4]=r3                ; aa[5,5]=1/r3            ; aa[5,6]=-r3                ; aa[5,7]=-1/r3
    aa[6,4]=Ec/(1+nc)/(1-2*nc); aa[6,5]=-Ec/(1+nc)/r3**2; aa[6,6]=-Es/(1+ns)/(1-2*ns); aa[6,7]=Es/(1+ns)/r3**2
    aa[7,6]=r2                ; aa[7,7]=1/r2            ; aa[7,8]=-r2                ; aa[7,9]=-1/r2
    aa[8,6]=Es/(1+ns)/(1-2*ns); aa[8,7]=-Es/(1+ns)/r2**2; aa[8,8]=-Ec/(1+nc)/(1-2*nc); aa[8,9]=Ec/(1+nc)/r2**2
    aa[9,8]=Ec/(1+nc)/(1-2*nc); aa[9,9]=-Ec/(1+nc)/r1**2
    
    bb[0]=-pb+Ec*alpc*tt/(1-nc)*(r6**2-r5**2)/2/r6**2
    bb[1]=(1+ns)/(1-ns)*alps*tt*(r5**2-r4**2)/2/r5
    bb[2]=-Es*alps*tt/(1-ns)*(r5**2-r4**2)/2/r5**2
    bb[3]=(1+nc)/(1-nc)*alpc*tt*(r4**2-r3**2)/2/r4
    bb[4]=-Ec*alpc*tt/(1-nc)*(r4**2-r3**2)/2/r4**2
    bb[5]=(1+ns)/(1-ns)*alps*tt*(r3**2-r2**2)/2/r3
    bb[6]=-Es*alps*tt/(1-ns)*(r3**2-r2**2)/2/r3**2
    bb[7]=(1+nc)/(1-nc)*alpc*tt*(r2**2-r1**2)/2/r2
    bb[8]=-Ec*alpc*tt/(1-nc)*(r2**2-r1**2)/2/r2**2
    bb[9]=0

    cc = np.linalg.solve(aa, bb)

    uu=np.zeros(n,dtype=np.float64)
    sr=np.zeros(n,dtype=np.float64)
    st=np.zeros(n,dtype=np.float64)
    sz=np.zeros(n,dtype=np.float64)
    ss=[]
    ss=ss+['r6_conc']
    ss=ss+['r5_conc']
    ss=ss+['r5_rbar']
    ss=ss+['r4_rbar']
    ss=ss+['r4_conc']
    ss=ss+['r3_conc']
    ss=ss+['r3_rbar']
    ss=ss+['r2_rbar']
    ss=ss+['r2_conc']
    ss=ss+['r1_conc']
    i=0; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[5],rr[4],tt,cc[0],cc[1])
    i=1; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[4],rr[4],tt,cc[0],cc[1])
    i=2; uu[i],sr[i],st[i],sz[i]=reb(Es,ns,alps,rr[4],rr[3],tt,cc[2],cc[3])
    i=3; uu[i],sr[i],st[i],sz[i]=reb(Es,ns,alps,rr[3],rr[3],tt,cc[2],cc[3])
    i=4; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[3],rr[2],tt,cc[4],cc[5])
    i=5; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[2],rr[2],tt,cc[4],cc[5])
    i=6; uu[i],sr[i],st[i],sz[i]=reb(Es,ns,alps,rr[2],rr[1],tt,cc[6],cc[7])
    i=7; uu[i],sr[i],st[i],sz[i]=reb(Es,ns,alps,rr[1],rr[1],tt,cc[6],cc[7])
    i=8; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[1],rr[0],tt,cc[8],cc[9])
    i=9; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[0],rr[0],tt,cc[8],cc[9])
    print('* Double Reinforcement')
    print('{0:>6s} {1:>6s} {2:>6s} {3:>6} {4:>6s} {5:>6s} {6:>6s} {7:>6s} {8:>6s}'
          .format('ra','rb','r0','ta','tb','da','db','Pb','T'))
    print('{0:6.0f} {1:6.0f} {2:6.0f} {3:6.3f} {4:6.3f} {5:6.0f} {6:6.0f} {7:6.3f} {8:6.3f}'
          .format(r1,r6,0,r3-r2,r5-r4,r2-r1,r6-r5,pb,tt))
    print('{0:7s} {1:>10s} {2:>10s} {3:>10s} {4:>10s}'.format('Loc.','u_r','sig_r','sig_t','sig_z'))
    for i in range(n):
        print('{0:7s} {1:10.3f} {2:10.3f} {3:10.3f} {4:10.3f}'.format(ss[i],uu[i],sr[i],st[i],sz[i]))

        
def pe_calc_s(rr,Es,ns,alps,Ec,nc,alpc,tt,pb):
    # single reinforcement
    r4=rr[3] # external surface of concrete lining
    r3=rr[2] # boundary of reinforcement and outer cover concrete
    r2=rr[1] # boundary of reinforcement and inner cover concrete
    r1=rr[0] # inner surface of concrete lining
    
    n=6
    aa=np.zeros((n,n),dtype=np.float64)
    bb=np.zeros(n,dtype=np.float64)
    aa[0,0]=Ec/(1+nc)/(1-2*nc); aa[0,1]=-Ec/(1+nc)/r4**2
    aa[1,0]=r3                ; aa[1,1]=1/r3            ; aa[1,2]=-r3                ; aa[1,3]=-1/r3
    aa[2,0]=Ec/(1+nc)/(1-2*nc); aa[2,1]=-Ec/(1+nc)/r3**2; aa[2,2]=-Es/(1+ns)/(1-2*ns); aa[2,3]=Es/(1+ns)/r3**2
    aa[3,2]=r2                ; aa[3,3]=1/r2            ; aa[3,4]=-r2                ; aa[3,5]=-1/r2
    aa[4,2]=Es/(1+ns)/(1-2*ns); aa[4,3]=-Es/(1+ns)/r2**2; aa[4,4]=-Ec/(1+nc)/(1-2*nc); aa[4,5]=Ec/(1+nc)/r2**2
    aa[5,4]=Ec/(1+nc)/(1-2*nc); aa[5,5]=-Ec/(1+nc)/r1**2
    
    bb[0]=-pb+Ec*alpc*tt/(1-nc)*(r4**2-r3**2)/2/r4**2
    bb[1]=(1+ns)/(1-ns)*alps*tt*(r3**2-r2**2)/2/r3
    bb[2]=-Es*alps*tt/(1-ns)*(r3**2-r2**2)/2/r3**2
    bb[3]=(1+nc)/(1-nc)*alpc*tt*(r2**2-r1**2)/2/r2
    bb[4]=-Ec*alpc*tt/(1-nc)*(r2**2-r1**2)/2/r2**2
    bb[5]=0

    cc = np.linalg.solve(aa, bb)

    uu=np.zeros(n,dtype=np.float64)
    sr=np.zeros(n,dtype=np.float64)
    st=np.zeros(n,dtype=np.float64)
    sz=np.zeros(n,dtype=np.float64)
    ss=[]
    ss=ss+['r4_conc']
    ss=ss+['r3_conc']
    ss=ss+['r3_rbar']
    ss=ss+['r2_rbar']
    ss=ss+['r2_conc']
    ss=ss+['r1_conc']
    i=0; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[3],rr[2],tt,cc[0],cc[1])
    i=1; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[2],rr[2],tt,cc[0],cc[1])
    i=2; uu[i],sr[i],st[i],sz[i]=reb(Es,ns,alps,rr[2],rr[1],tt,cc[2],cc[3])
    i=3; uu[i],sr[i],st[i],sz[i]=reb(Es,ns,alps,rr[1],rr[1],tt,cc[2],cc[3])
    i=4; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[1],rr[0],tt,cc[4],cc[5])
    i=5; uu[i],sr[i],st[i],sz[i]=con(Ec,nc,alpc,rr[0],rr[0],tt,cc[4],cc[5])
    print('* Single Reinforcement')
    print('{0:>6s} {1:>6s} {2:>6s} {3:>6} {4:>6s} {5:>6s} {6:>6s} {7:>6s} {8:>6s}'
          .format('ra','rb','r0','ta','tb','da','db','Pb','T'))
    print('{0:6.0f} {1:6.0f} {2:6.0f} {3:6.3f} {4:6.0f} {5:6.0f} {6:6.0f} {7:6.3f} {8:6.3f}'
          .format(r1,r4,0,r3-r2,0,r2-r1,0,pb,tt))
    print('{0:7s} {1:>10s} {2:>10s} {3:>10s} {4:>10s}'.format('Loc.','u_r','sig_r','sig_t','sig_z'))
    for i in range(n):
        print('{0:7s} {1:10.3f} {2:10.3f} {3:10.3f} {4:10.3f}'.format(ss[i],uu[i],sr[i],st[i],sz[i]))
        

def main():
    # properties of materials
    Es=200000   # elastic modulus of steel
    ns=0.2      # Poisson's ratio of steel
    alps=1.0e-5 # thermal expansion coefficient of steel
    Ec=25000    # elastic modulus of concrete
    nc=0.2      # Poisson's ratio of concrete
    alpc=1.0e-5 # thermal expansion coefficient of concrete
    T32=32**2*np.pi/4
    T25=25**2*np.pi/4
    T20=20**2*np.pi/4

    # double reinforcement
    ra=4000.0 # internal radius of pressure tunnel (mm)
    rb=4800.0 # external radius of pressure tunnel (mm)
    ta=T32*5/1000   # equivalent steel thickness (mm)
    tb=T32*5/1000   # equivalent steel thickness (mm)
    da=100    # cover fro internal surface of concrete (mm)
    db=100    # cover from external surface of concrete (mm)
    tt=0    # temperature change (negative: temperature drop)
    pb=3.0    # external pressure (MPa)
    rr=np.array([ra,ra+da,ra+da+ta,rb-db-ta,rb-db,rb])
    pe_calc_d(rr,Es,ns,alps,Ec,nc,alpc,tt,pb)

    # single reinforcement
    ra=4000.0 # internal radius of pressure tunnel (mm)
    rb=4800.0 # external radius of pressure tunnel (mm)
    ta=T32*5/1000   # equivalent steel thickness (mm)
    da=100    # cover fro internal surface of concrete (mm)
    tt=0    # temperature change (negative: temperature drop)
    pb=3.0    # external pressure (MPa)
    rr=np.array([ra,ra+da,ra+da+ta,rb])
    pe_calc_s(rr,Es,ns,alps,Ec,nc,alpc,tt,pb)
    
    
#==============
# Execution
#==============
if __name__ == '__main__': main()

Thank you.

記事の先頭に行く

2019-05-23

設計 RC円形圧力トンネルの配筋設計（１）

設計

記事の最後に行く

概要

均等な内水圧および外水圧を受けるRC円形圧力トンネルを、トンネル軸方向に均一な厚肉円筒として、平面ひずみ状態でモデル化する。モデル化においては、以下の考え方を採用した。

圧力水路の構造は鉄筋コンクリート構造とする。
内水圧を受ける場合のモデル化範囲は、圧力水路の覆工（鉄筋含む）および周辺岩盤とする。
内水圧を受ける覆工コンクリートは、ひび割れが発生しているものとし、周方向の応力は分担しない。鉄筋は完全弾性体として扱う。
外水圧は覆工コンクリートの外面に作用するものとし、岩盤はモデル化しない。
外水圧を受ける覆工には半径方向・周方向ともに圧縮応力が作用するため、コンクリート・鉄筋とも完全弾性体として扱う。
温度変化量は、簡略化のため、覆工内で均一とし、岩盤内での温度変化は考慮しない。
水路における鉄筋のかぶりは、断面寸法に対し比較的大きいので、これを考慮する。

f:id:damyarou:20190527132711j:plain

Design of RC Circular Pressure Tunnel

Basic Equations for Elastic Cylinder Model in Plane Strain State

In this discussion, a long circular cylinder including surrounding area is considered as a pressure tunnel model using polar coordinates. Symbols used are shown below in this discussion.

$\sigma_r$	normal stress in the radial direction
$\epsilon_r$	normal strain in the radial direction
$\sigma_{\theta}$	normal stress in the circumferential direction
$\epsilon_{\theta}$	normal strain in the circumferential direction
$\sigma_z$	normal stress in the axial direction
$\epsilon_z$	normal strain in the axial direction
$u$	displacement in the radial direction
$r$	distance to a point in a cylinder from origin of polar coordinates in the radial direction
$a$	distance to inner boundary of a cylinder from origin of polar coordinates in the radial direction
$E$	elastic modulus of material
$\nu$	Poisson's ratio of material
$\alpha$	thermal expansion coefficient of material
$T$	temperature change of material

Basic Equations for Isotropic Elastic Material in Plane Strain State

Stress - Strain Relationship

Stress - strain relationships for an elastic material in polar coordinates are shown below.

$\begin{equation} \begin{cases} &\epsilon_{r} -\alpha T = \cfrac{1}{E}\{\sigma_{r} -\nu (\sigma_{\theta}+\sigma_{z})\} \\ &\epsilon_{\theta}-\alpha T = \cfrac{1}{E}\{\sigma_{\theta}-\nu (\sigma_{z} +\sigma_{r})\} \\ &\epsilon_{z} -\alpha T = \cfrac{1}{E}\{\sigma_{z} -\nu (\sigma_{r} +\sigma_{\theta})\} \end{cases} \end{equation}$

Considering a condition os $\epsilon_z=0$ , followings can be obtained from above.

$\begin{equation} \begin{cases} &\epsilon_{r} -\alpha T = \cfrac{1}{E}\{(1-\nu^2) \sigma_{r} -\nu (1+\nu) \sigma_{\theta}\} \\ &\epsilon_{\theta}-\alpha T = \cfrac{1}{E}\{(1-\nu^2) \sigma_{\theta}-\nu (1+\nu) \sigma_{r}\} \\ &\sigma_z=\nu (\sigma_r+\sigma_{\theta}) - E \alpha T \end{cases} \end{equation}$

$\begin{equation} \begin{cases} \sigma_r =\cfrac{E}{(1+\nu)(1-2\nu)}\{(1-\nu) \epsilon_r + \nu\epsilon_{\theta} - (1+\nu)\alpha T\} \\ \sigma_{\theta}=\cfrac{E}{(1+\nu)(1-2\nu)}\{\nu \epsilon_r + (1-\nu) \epsilon_{\theta} - (1+\nu)\alpha T\} \end{cases} \end{equation}$

Strain - Displacement Relationship

Relationships between strain and displacement can be expressed as follows.

$\begin{equation} \epsilon_r=\cfrac{du}{dr} \qquad \epsilon_{\theta}=\cfrac{u}{r} \end{equation}$

Equilibrium Equation

An equilibrium equation of stress can be expressed as follow.

$\begin{equation} \cfrac{d\sigma_r}{dr}+\cfrac{\sigma_r-\sigma_{\theta}}{r}=0 \end{equation}$

General Expressions of Displacements and Stresses

$\begin{equation} \cfrac{d\sigma_r}{dr}+\cfrac{\sigma_r-\sigma_{\theta}}{r}=\cfrac{E(1-\nu)}{(1+\nu)(1-2\nu)} \left\{\cfrac{d^2u}{dr^2}+\cfrac{1}{r}\cfrac{du}{dr}-\cfrac{u}{r^2}-\cfrac{1+\nu}{1-\nu} \alpha \cfrac{dT}{dr}\right\}=0 \end{equation}$

$\begin{equation} \cfrac{d^2u}{dr^2}+\cfrac{1}{r}\cfrac{du}{dr}-\cfrac{u}{r^2}=\cfrac{1+\nu}{1-\nu} \alpha \cfrac{dT}{dr} \quad \rightarrow \quad \cfrac{d}{dr}\left[\cfrac{1}{r}\cfrac{d(ru)}{dr}\right]=\cfrac{1+\nu}{1-\nu} \alpha \cfrac{dT}{dr} \end{equation}$

$\begin{align} &\cfrac{d}{dr}\left[\cfrac{1}{r}\cfrac{d(ru)}{dr}\right]=0 & \rightarrow & & &u=C_1\cdot r + \cfrac{C_2}{r} & &\text{(general solution)}\\ &\cfrac{d}{dr}\left[\cfrac{1}{r}\cfrac{d(ru)}{dr}\right]=\cfrac{1+\nu}{1-\nu} \alpha \cfrac{dT}{dr} & \rightarrow & & &u=\cfrac{1+\nu}{1-\nu} \alpha \cfrac{1}{r} \int_a^r T r dr & &\text{(paticular solution)} \end{align}$

From above, following general expressions of displacement and stresses can be obtained. (These equations are shown in 'Theory of Elasticity, S. Timoshenko and J. N. Goodier, Chapter 14. Thermal Stress, 135. The Long Circular Cylinder.')

$\begin{align} &u =\ \ \cfrac{1+\nu}{1-\nu} \alpha \cfrac{1}{r} \int_a^r T r dr & &+ C_1\cdot r + \cfrac{C_2}{r} \\ &\sigma_r = -\cfrac{\alpha E}{1-\nu} \cfrac{1}{r^2} \int_a^r T r dr & &+ \cfrac{E}{(1+\nu)(1-2\nu)}\cdot C_1-\cfrac{E}{(1+\nu)}\cfrac{C_2}{r^2} \\ &\sigma_{\theta}=\ \ \cfrac{\alpha E}{1-\nu} \cfrac{1}{r^2} \int_a^r T r dr - \cfrac{\alpha E T}{1-\nu} & &+ \cfrac{E}{(1+\nu)(1-2\nu)}\cdot C_1+\cfrac{E}{(1+\nu)}\cfrac{C_2}{r^2} \\ &\sigma_z = -\cfrac{\alpha E T}{1-\nu} & &+ \cfrac{2 \nu E}{(1+\nu)(1-2\nu)}\cdot C_1 \end{align}$

In case of uniform distribution of temperature in a material, thermal items can be expressed as results of integrations.

$\begin{align} &\cfrac{1+\nu}{1-\nu} \alpha \cfrac{1}{r} \int_a^r T r dr & &=\ \ \cfrac{1+\nu}{1-\nu} \alpha T \cfrac{r^2-a^2}{2 r} \\ &\cfrac{\alpha E}{1-\nu} \cfrac{1}{r^2} \int_a^r T r dr & &=-\cfrac{E \alpha T}{1-\nu} \cfrac{r^2-a^2}{2 r^2} \\ &\cfrac{\alpha E}{1-\nu} \cfrac{1}{r^2} \int_a^r T r dr - \cfrac{\alpha E T}{1-\nu} & &=-\cfrac{E \alpha T}{1-\nu} \cfrac{r^2+a^2}{2 r^2} \end{align}$

Therefore, displacement and stresses can be expressed as follows.

$\begin{align} &u =\ \ \cfrac{1+\nu}{1-\nu} \alpha T \cfrac{r^2-a^2}{2 r} & &+ C_1\cdot r + \cfrac{C_2}{r} \\ &\sigma_r =-\cfrac{E \alpha T}{1-\nu} \cfrac{r^2-a^2}{2 r^2} & &+ \cfrac{E}{(1+\nu)(1-2\nu)}\cdot C_1-\cfrac{E}{(1+\nu)}\cfrac{C_2}{r^2} \\ &\sigma_{\theta}=-\cfrac{E \alpha T}{1-\nu} \cfrac{r^2+a^2}{2 r^2} & &+ \cfrac{E}{(1+\nu)(1-2\nu)}\cdot C_1+\cfrac{E}{(1+\nu)}\cfrac{C_2}{r^2} \\ &\sigma_z = -\cfrac{\alpha E T}{1-\nu} & &+ \cfrac{2 \nu E}{(1+\nu)(1-2\nu)}\cdot C_1 \end{align}$

Basic Equations for No-tension Material in Circumferential Direction

For no-tension material in circumferential direction such as concrete under the internal pressure, the equilibrium equation becomes shown below considering a condition of $\sigma_{\theta}=0$ .

$\begin{equation} \cfrac{d \sigma_r}{dr} + \cfrac{\sigma_r}{r} = 0 \end{equation}$

Using an assumption of Poisson's ratio of zero ( $\nu=0$ ), the stress in the radial direction can be expressed as follow.

$\begin{equation} \sigma_r=E\epsilon_r-E \alpha T \quad \rightarrow \quad \sigma_r=E\cfrac{du}{dr}-E \alpha T \quad \left(\epsilon_r=\cfrac{du}{dr}\right) \end{equation}$

$\begin{equation} \cfrac{d \sigma_r}{dr} + \cfrac{\sigma_r}{r} = E \left\{ \cfrac{d^2 u}{dr^2} + \cfrac{1}{r} \cfrac{du}{dr} - \alpha \left( \cfrac{dT}{dr} + \cfrac{T}{r} \right) \right\}=0 \end{equation}$

$\begin{align} &\cfrac{d^2 u}{dr^2} + \cfrac{1}{r}\cfrac{du}{dr}=0 & \rightarrow & & &u=C_1 + C_2\cdot \ln(r) & &\text{(general solution)}\\ &\cfrac{d^2 u}{dr^2} + \cfrac{1}{r}\cfrac{du}{dr}=\alpha \left(\cfrac{dT}{dr}+\cfrac{T}{r}\right) & \rightarrow & & &u=\alpha \int_a^r T dr & &\text{(paticular solution)} \end{align}$

From above, following basic equations for no-tension material in circumferential direction can be obtained.

$\begin{align} &u=C_1+C_2\cdot \ln(r) + \alpha \int_a^r T dr \\ &\sigma_r=E \cfrac{C_2}{r} \end{align}$

Using an assumption of uniform distribution of temperature, a integration of thermal item becoms $\alpha T (r-a)$ . Therefore, displacement and stresses can be expressed as follows.

$\begin{align} &u=C_1+C_2\cdot \ln(r) + \alpha T (r-a) \\ &\sigma_r=E \cfrac{C_2}{r} \\ &\sigma_{\theta}=0 \end{align}$

Model of RC Circular Tunnel under Internal Pressure

Double Reinforcement Section

Components of Model

Components of Model
Bedrock	Elastic material. Thermal stress is ignored.
Concrete (outer cover)	No-tension material. Thermal stress is considered.
Outer Reinforcement	Elastic material. Thermal stress is considered.
Concrete (middle location)	No-tension material. Thermal stress is considered.
Inner Reinforcement	Elastic material. Thermal stress is considered.
Concrete (Inner cover)	No-tension material. Thermal stress is considered.

Coordinate of Boundary in Radial Direction
$r_7$	Outer boundary of bedrock
$r_6$	External surface of concrete lining
$r_5$	Boundary of outer cover concrete and outer reinforcement
$r_4$	Boundary of outer reinforcement and middle concrete
$r_3$	Boundary of middle concrete and inner reinforcement
$r_2$	Boundary of inner reinforcement and inner cover concrete
$r_1$	Inner surface of concrete lining

Basic Equations for Each Material

Bedrock ( $r_6 \leqq r \leqq r_7$ )

$\begin{align} % displacement of bedrock &u=C_{g1}\cdot r + \cfrac{C_{g2}}{r} \\ % stress of bedrock &\sigma_r=\cfrac{E_g}{(1+\nu_g)(1-2\nu_g)}\cdot C_{g1}-\cfrac{E_g}{(1+\nu_g)}\cfrac{C_{g2}}{r^2} \\ &\sigma_{\theta}=\cfrac{E_g}{(1+\nu_g)(1-2\nu_g)}\cdot C_{g1}+\cfrac{E_g}{(1+\nu_g)}\cfrac{C_{g2}}{r^2} \\ \end{align}$

Outer Cover Concrete ( $r_5 \leqq r \leqq r_6$ )

$\begin{align} % displacement of outer cover &u=C_{co1}+C_{co2}\cdot \ln(r) + \alpha_c T (r-r_5) \\ % stress of outer cover &\sigma_r=E_c \cfrac{C_{co2}}{r} \\ &\sigma_{\theta}=0 \end{align}$

Outer Reinforcement ( $r_4 \leqq r \leqq r_5$ )

$\begin{align} % displacement of outer rebar &u =\ \ \cfrac{1+\nu_s}{1-\nu_s} \alpha_s T \cfrac{r^2-r_4{}^2}{2 r} & &+ C_{so1}\cdot r + \cfrac{C_{so2}}{r} \\ % stress of outer rebar &\sigma_r =-\cfrac{E_s \alpha_s T}{1-\nu_s} \cfrac{r^2-r_4{}^2}{2 r^2} & &+ \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{so1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{so2}}{r^2} \\ &\sigma_{\theta}=-\cfrac{E_s \alpha_s T}{1-\nu_s} \cfrac{r^2+r_4{}^2}{2 r^2} & &+ \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{so1}+\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{so2}}{r^2} \end{align}$

Middle Concrete ( $r_3 \leqq r \leqq r_4$ )

$\begin{align} % displacement of middle concrete &u=C_{cm1}+C_{cm2}\cdot \ln(r) + \alpha_c T (r-r_3) \\ % stress of middle rebar &\sigma_r=E_c \cfrac{C_{cm2}}{r} \\ &\sigma_{\theta}=0 \end{align}$

Inner Reinforcement ( $r_2 \leqq r \leqq r_3$ )

$\begin{align} % displacement of inner rebar &u =\ \ \cfrac{1+\nu_s}{1-\nu_s} \alpha_s T \cfrac{r^2-r_2{}^2}{2 r} & &+ C_{si1}\cdot r + \cfrac{C_{si2}}{r} \\ % stress of inner rebar &\sigma_r =-\cfrac{E_s \alpha_s T}{1-\nu_s} \cfrac{r^2-r_2{}^2}{2 r^2} & &+ \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{si1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{si2}}{r^2} \\ &\sigma_{\theta}=-\cfrac{E_s \alpha_s T}{1-\nu_s} \cfrac{r^2+r_2{}^2}{2 r^2} & &+ \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{si1}+\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{si2}}{r^2} \end{align}$

Inner Cover Concrete ( $r_1 \leqq r \leqq r_2$ )

$\begin{align} % displacement of inner cover &u=C_{ci1}+C_{ci2}\cdot \ln(r) + \alpha_c T (r-r_1) \\ % stress of inner cover &\sigma_r= E_c \cfrac{C_{ci2}}{r} \\ &\sigma_{\theta}=0 \end{align}$

Simultaneous linear equations

To fix un-known parameters, simultaneous linear equations will be created considering boundary conditions shown below.

Stress in the radial direction at outer boundary of bedrock is equal to zero.
Displacement and stress in the radial direction have continuity at the boundaries of each material.
Stress in the radial direction at inner surface of concrete lining is equal to the internal pressure.
Positive sign of the internal pressure $P_a$ is toward outer direction.

$\begin{align} &1.&r=r_7 \;\text{(stress)}\quad &\left\{ % stress of bedrock \cfrac{E_g}{(1+\nu_g)(1-2\nu_g)}\cdot C_{g1}-\cfrac{E_g}{(1+\nu_g)}\cfrac{C_{g2}}{r_7{}^2} \right\} & &=0 \\ % &2.&r=r_6 \;\text{(disp.)}\quad &\left\{ % displacement of bedrock C_{g1}\cdot r_6 + \cfrac{C_{g2}}{r_6} \right\}-\left\{ % displacement of outer cover C_{co1}+C_{co2}\cdot \ln(r_6) \right\} & &=\alpha_c T (r_6-r_5) \\ % &3.&r=r_6 \;\text{(stress)}\quad &\left\{ % stress of bedrock \cfrac{E_g}{(1+\nu_g)(1-2\nu_g)}\cdot C_{g1}-\cfrac{E_g}{(1+\nu_g)}\cfrac{C_{g2}}{r_6{}^2} \right\}-\left\{ % stress of outer cover E_c \cfrac{C_{co2}}{r_6} \right\} & &=0 \\ % &4.&r=r_5 \;\text{(disp.)}\quad &\left\{ % displacement of outer cover C_{co1}+C_{co2}\cdot \ln(r_5) \right\}-\left\{ % displacement of outer rebar C_{so1}\cdot r_5 + \cfrac{C_{so2}}{r_5} \right\} & &=\cfrac{1+\nu_s}{1-\nu_s} \alpha_s T \cfrac{r_5{}^2-r_4{}^2}{2 r_5} \\ % &5.&r=r_5 \;\text{(stress)}\quad &\left\{ % stress of outer cover E_c \cfrac{C_{co2}}{r_5} \right\}-\left\{ % stress of outer rebar \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{so1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{so2}}{r_5{}^2} \right\} & &=-\cfrac{E_s \alpha_s T}{1-\nu_s} \cfrac{r_5{}^2-r_4{}^2}{2 r_5{}^2} \\ % &6.&r=r_4 \;\text{(disp.)}\quad &\left\{ % displacement of outer rebar C_{so1}\cdot r_4 + \cfrac{C_{so2}}{r_4} \right\}-\left\{ % displacement of middle concrete C_{cm1}+C_{cm2}\cdot \ln(r_4) \right\} & &=\alpha_c T (r_4-r_3) \\ &7.&r=r_4 \;\text{(stress)}\quad &\left\{ % stress of outer rebar \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{so1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{so2}}{r_4{}^2} \right\}-\left\{ % stress of middle rebar E_c \cfrac{C_{cm2}}{r_4} \right\} & &=0 \\ % &8.&r=r_3 \;\text{(disp.)}\quad &\left\{ % displacement of middle concrete C_{cm1}+C_{cm2}\cdot \ln(r_3) \right\}-\left\{ % displacement of inner rebar C_{si1}\cdot r_3 + \cfrac{C_{si2}}{r_3} \right\} & &=\cfrac{1+\nu_s}{1-\nu_s} \alpha_s T \cfrac{r_3{}^2-r_2{}^2}{2 r_3} \\ &9.&r=r_3 \;\text{(stress)}\quad &\left\{ % stress of middle rebar E_c \cfrac{C_{cm2}}{r_3} \right\}-\left\{ % stress of inner rebar \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{si1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{si2}}{r_3{}^2} \right\} & &=-\cfrac{E_s \alpha_s T}{1-\nu_s} \cfrac{r_3{}^2-r_2{}^2}{2 r_3{}^2} \\ % &10.&r=r_2 \;\text{(disp.)}\quad &\left\{ % displacement of inner rebar C_{si1}\cdot r_2 + \cfrac{C_{si2}}{r_2} \right\}-\left\{ % displacement of inner cover C_{ci1}+C_{ci2}\cdot \ln(r_2) \right\} & &=\alpha_c T (r_2-r_1) \\ &11.&r=r_2 \;\text{(stress)}\quad &\left\{ % stress of inner rebar \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{si1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{si2}}{r_2{}^2} \right\}-\left\{ % stress of inner cover E_c \cfrac{C_{ci2}}{r_2} \right\} & &=0 \\ % &12.&r=r_1 \;\text{(stress)}\quad &\left\{ E_c \cfrac{C_{ci2}}{r_1} \right\} & &=-P_a \end{align}$

Matrix Expression of Simultaneous linear equations for Computer Programing

$\begin{equation} \begin{bmatrix} a_{1,1} & a_{1,2} & 0 & 0 & \dots & 0 & 0 & 0 & 0 \\ a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} & \dots & 0 & 0 & 0 & 0 \\ a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & \dots & 0 & 0 & 0 & 0 \\ 0 & 0 & a_{4,3} & a_{4,4} & \dots & 0 & 0 & 0 & 0 \\ 0 & 0 & a_{5,3} & a_{5,4} & \dots & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \dots & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \dots & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \dots &a_{8,9} & a_{8,10} & 0 & 0 \\ 0 & 0 & 0 & 0 & \dots &a_{9,9} & a_{9,10} & 0 & 0 \\ 0 & 0 & 0 & 0 & \dots &a_{10,9} & a_{10,10} & a_{10,11} & a_{10,12} \\ 0 & 0 & 0 & 0 & \dots &a_{11,9} & a_{11,10} & a_{11,11} & a_{11,12} \\ 0 & 0 & 0 & 0 & \dots & 0 & 0 & a_{12,11} & a_{12,12} \end{bmatrix} \begin{Bmatrix} C_{g1} \\ C_{g2} \\ C_{co1} \\ C_{co2} \\ C_{so1} \\ C_{so2} \\ C_{cm1} \\ C_{cm2} \\ C_{si1} \\ C_{si2} \\ C_{ci1} \\ C_{ci2} \\ \end{Bmatrix} = \begin{Bmatrix} 0 \\ \alpha_c T (r_6-r_5) \\ 0 \\ \frac{1+\nu_s}{1-\nu_s} \alpha_s T \frac{r_5{}^2-r_4{}^2}{2 r_5} \\ -\frac{E_s \alpha_s T}{1-\nu_s} \frac{r_5{}^2-r_4{}^2}{2 r_5{}^2} \\ \alpha_c T (r_4-r_3) \\ 0 \\ \frac{1+\nu_s}{1-\nu_s} \alpha_s T \frac{r_3{}^2-r_2{}^2}{2 r_3} \\ -\frac{E_s \alpha_s T}{1-\nu_s} \frac{r_3{}^2-r_2{}^2}{2 r_3{}^2} \\ \alpha_c T (r_2-r_1) \\ 0 \\ -P_a \end{Bmatrix} \end{equation}$

$\begin{align} &a_{1,1}=\cfrac{E_g}{(1+\nu_g)(1-2\nu_g)} & &a_{1,2}=-\cfrac{E_g}{(1+\nu_g)r_7{}^2} & &a_{1,3}=0 & &a_{1,4}=0 \\ % &a_{2,1}=r_6 & &a_{2,2}=\cfrac{1}{r_6} & &a_{2,3}=-1 & &a_{2,4}=-\ln(r_6) \\ &a_{3,1}=\cfrac{E_g}{(1+\nu_g)(1-2\nu_g)} & &a_{3,2}=-\cfrac{E_g}{(1+\nu_g)r_6{}^2} & &a_{3,3}=0 & &a_{3,4}=-\cfrac{E_c}{r_6} \\ % &a_{4,3}=1 & &a_{4,4}=\ln(r_5) & &a_{4,5}=-r_5 & &a_{4,6}=-\cfrac{1}{r_5} \\ &a_{5,3}=0 & &a_{5,4}=\cfrac{E_c}{r_5} & &a_{5,5}=-\cfrac{E_s}{(1+\nu_s)(1-2\nu_s)} & &a_{5,6}=\cfrac{E_s}{(1+\nu_s)r_5{}^2} \\ % &a_{6,5}=r_4 & &a_{6,6}=\cfrac{1}{r_4} & &a_{6,7}=-1 & &a_{6,8}=-\ln(r_4) \\ &a_{7,5}=\cfrac{E_s}{(1+\nu_s)(1-2\nu_s)} & &a_{7,6}=-\cfrac{E_s}{(1+\nu_s)r_4{}^2} & &a_{7,7}=0 & &a_{7,8}=-\cfrac{E_c}{r_4} \\ % &a_{8,7}=1 & &a_{8,8}=\ln(r_3) & &a_{8,9}=-r_3 & &a_{8,10}=-\cfrac{1}{r_3} \\ &a_{9,7}=0 & &a_{9,8}=\cfrac{E_c}{r_3} & &a_{9,9}=-\cfrac{E_s}{(1+\nu_s)(1-2\nu_s)} & &a_{9,10}=\cfrac{E_s}{(1+\nu_s)r_3{}^2} \\ % &a_{10,9}=r_2 & &a_{10,10}=\cfrac{1}{r_2} & &a_{10,11}=-1 & &a_{10,12}=-\ln(r_2) \\ &a_{11,9}=\cfrac{E_s}{(1+\nu_s)(1-2\nu_s)} & &a_{11,10}=-\cfrac{E_s}{(1+\nu_s)r_2{}^2} & &a_{11,11}=0 & &a_{11,12}=-\cfrac{E_c}{r_2} \\ % &a_{12,9}=0 & &a_{12,10}=0 & &a_{12,11}=0 & &a_{12,12}=\cfrac{E_c}{r_1} \end{align}$

Single Reinforcement Section

Components of Model

Components of Model
Bedrock	Elastic material. Thermal stress is ignored.
Concrete (outer cover)	No-tension material. Thermal stress is considered.
Inner Reinforcement	Elastic material. Thermal stress is considered.
Concrete (Inner cover)	No-tension material. Thermal stress is considered.

Coordinate of Boundary in Radial Direction
$r_5$	Outer boundary of bedrock
$r_4$	External surface of concrete lining
$r_3$	Boundary of outer cover concrete and reinforcement
$r_2$	Boundary of reinforcement and inner cover concrete
$r_1$	Inner surface of concrete lining

Basic Equations for Each Material

Bedrock ( $r_4 \leqq r \leqq r_5$ )

Outer Cover Concrete ( $r_3 \leqq r \leqq r_4$ )

$\begin{align} % displacement of outer cover &u=C_{co1}+C_{co2}\cdot \ln(r) + \alpha_c T (r-r_3) \\ % stress of outer cover &\sigma_r=E_c \cfrac{C_{co2}}{r} \\ &\sigma_{\theta}=0 \end{align}$

Reinforcement ( $r_2 \leqq r \leqq r_3$ )

$\begin{align} % displacement of rebar &u =\ \ \cfrac{1+\nu_s}{1-\nu_s} \alpha_s T \cfrac{r^2-r_2{}^2}{2 r} & &+ C_{si1}\cdot r + \cfrac{C_{si2}}{r} \\ % stress of rebar &\sigma_r =-\cfrac{E_s \alpha_s T}{1-\nu_s} \cfrac{r^2-r_2{}^2}{2 r^2} & &+ \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{si1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{si2}}{r^2} \\ &\sigma_{\theta}=-\cfrac{E_s \alpha_s T}{1-\nu_s} \cfrac{r^2+r_2{}^2}{2 r^2} & &+ \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{si1}+\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{si2}}{r^2} \end{align}$

Inner Cover Concrete ( $r_1 \leqq r \leqq r_2$ )

$\begin{align} % displacement of inner cover &u=C_{ci1}+C_{ci2}\cdot \ln(r) + \alpha_c T (r-r_1) \\ % stress of inner cover &\sigma_r= E_c \cfrac{C_{ci2}}{r} \\ &\sigma_{\theta}=0 \end{align}$

Simultaneous linear equations

To fix un-known parameters, simultaneous linear equations will be created considering boundary conditions shown below.

Stress in the radial direction at outer boundary of bedrock is equal to zero.
Displacement and stress in the radial direction have continuity at the boundaries of each material.
Stress in the radial direction at inner surface of concrete lining is equal to the internal pressure.
Positive sign of the internal pressure $P_a$ is toward outer direction.

$\begin{align} &1.&r=r_5 \;\text{(stress)}\quad &\left\{ % stress of bedrock \cfrac{E_g}{(1+\nu_g)(1-2\nu_g)}\cdot C_{g1}-\cfrac{E_g}{(1+\nu_g)}\cfrac{C_{g2}}{r_5{}^2} \right\} & &=0 \\ % &2.&r=r_4 \;\text{(disp.)}\quad &\left\{ % displacement of bedrock C_{g1}\cdot r_4 + \cfrac{C_{g2}}{r_4} \right\}-\left\{ % displacement of outer cover C_{co1}+C_{co2}\cdot \ln(r_4) \right\} & &=\alpha_c T (r_4-r_3) \\ % &3.&r=r_4 \;\text{(stress)}\quad &\left\{ % stress of bedrock \cfrac{E_g}{(1+\nu_g)(1-2\nu_g)}\cdot C_{g1}-\cfrac{E_g}{(1+\nu_g)}\cfrac{C_{g2}}{r_4{}^2} \right\}-\left\{ % stress of outer cover E_c \cfrac{C_{co2}}{r_4} \right\} & &=0 \\ % &4.&r=r_3 \;\text{(disp.)}\quad &\left\{ % displacement of middle concrete C_{co1}+C_{co2}\cdot \ln(r_3) \right\}-\left\{ % displacement of inner rebar C_{si1}\cdot r_3 + \cfrac{C_{si2}}{r_3} \right\} & &=\cfrac{1+\nu_s}{1-\nu_s} \alpha_s T \cfrac{r_3{}^2-r_2{}^2}{2 r_3} \\ % &5.&r=r_3 \;\text{(stress)}\quad &\left\{ % stress of middle rebar E_c \cfrac{C_{co2}}{r_3} \right\}-\left\{ % stress of inner rebar \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{si1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{si2}}{r_3{}^2} \right\} & &=-\cfrac{E_s \alpha_s T}{1-\nu_s} \cfrac{r_3{}^2-r_2{}^2}{2 r_3{}^2} \\ % &6.&r=r_2 \;\text{(disp.)}\quad &\left\{ % displacement of inner rebar C_{si1}\cdot r_2 + \cfrac{C_{si2}}{r_2} \right\}-\left\{ % displacement of inner cover C_{ci1}+C_{ci2}\cdot \ln(r_2) \right\} & &=\alpha_c T (r_2-r_1) \\ % &7.&r=r_2 \;\text{(stress)}\quad &\left\{ % stress of inner rebar \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{si1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{si2}}{r_2{}^2} \right\}-\left\{ % stress of inner cover E_c \cfrac{C_{ci2}}{r_2} \right\} & &=0 \\ % &8.&r=r_1 \;\text{(stress)}\quad &\left\{ E_c \cfrac{C_{ci2}}{r_1} \right\} & &=-P_a \end{align}$

Matrix Expression of Simultaneous linear equations for Computer Programing

$\begin{equation} \begin{bmatrix} a_{1,1} & a_{1,2} & 0 & 0 & 0 & 0 & 0 & 0 \\ a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} & 0 & 0 & 0 & 0 \\ a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & 0 & 0 & 0 & 0 \\ 0 & 0 & a_{4,3} & a_{4,4} &a_{4,5} & a_{4,6} & 0 & 0 \\ 0 & 0 & a_{5,3} & a_{5,4} &a_{5,5} & a_{5,6} & 0 & 0 \\ 0 & 0 & 0 & 0 &a_{6,5} & a_{6,6} & a_{6,7} & a_{6,8} \\ 0 & 0 & 0 & 0 &a_{7,5} & a_{7,6} & a_{7,7} & a_{7,8} \\ 0 & 0 & 0 & 0 & 0 & 0 & a_{8,7} & a_{8,8} \end{bmatrix} \begin{Bmatrix} C_{g1} \\ C_{g2} \\ C_{co1} \\ C_{co2} \\ C_{si1} \\ C_{si2} \\ C_{ci1} \\ C_{ci2} \\ \end{Bmatrix} = \begin{Bmatrix} 0 \\ \alpha_c T (r_4-r_3) \\ 0 \\ \frac{1+\nu_s}{1-\nu_s} \alpha_s T \frac{r_3{}^2-r_2{}^2}{2 r_3} \\ -\frac{E_s \alpha_s T}{1-\nu_s} \frac{r_3{}^2-r_2{}^2}{2 r_3{}^2} \\ \alpha_c T (r_2-r_1) \\ 0 \\ -P_a \end{Bmatrix} \end{equation}$

$\begin{align} &a_{1,1}=\cfrac{E_g}{(1+\nu_g)(1-2\nu_g)} & &a_{1,2}=-\cfrac{E_g}{(1+\nu_g)r_5{}^2} & &a_{1,3}=0 & &a_{1,4}=0 \\ % &a_{2,1}=r_4 & &a_{2,2}=\cfrac{1}{r_4} & &a_{2,3}=-1 & &a_{2,4}=-\ln(r_4) \\ &a_{3,1}=\cfrac{E_g}{(1+\nu_g)(1-2\nu_g)} & &a_{3,2}=-\cfrac{E_g}{(1+\nu_g)r_4{}^2} & &a_{3,3}=0 & &a_{3,4}=-\cfrac{E_c}{r_4} \\ % &a_{4,3}=1 & &a_{4,4}=\ln(r_3) & &a_{4,5}=-r_3 & &a_{4,6}=-\cfrac{1}{r_3} \\ &a_{5,3}=0 & &a_{5,4}=\cfrac{E_c}{r_3} & &a_{5,5}=-\cfrac{E_s}{(1+\nu_s)(1-2\nu_s)} & &a_{5,6}=\cfrac{E_s}{(1+\nu_s)r_3{}^2} \\ % &a_{6,5}=r_2 & &a_{6,6}=\cfrac{1}{r_2} & &a_{6,7}=-1 & &a_{6,8}=-\ln(r_2) \\ &a_{7,5}=\cfrac{E_s}{(1+\nu_s)(1-2\nu_s)} & &a_{7,6}=-\cfrac{E_s}{(1+\nu_s)r_2{}^2} & &a_{7,7}=0 & &a_{7,8}=-\cfrac{E_c}{r_2} \\ % &a_{8,5}=0 & &a_{8,6}=0 & &a_{8,7}=0 & &a_{8,8}=\cfrac{E_c}{r_1} \end{align}$

Model of RC Circular Tunnel under External Pressure

Double Reinforcement Section

Components of Model

Components of Model
Concrete (outer cover)	Elastic material. Thermal stress is considered.
Outer Reinforcement	Elastic material. Thermal stress is considered.
Concrete (middle location)	Elastic material. Thermal stress is considered.
Inner Reinforcement	Elastic material. Thermal stress is considered.
Concrete (Inner cover)	Elastic material. Thermal stress is considered.

Coordinate of Boundary in Radial Direction
$r_6$	External surface of concrete lining
$r_5$	Boundary of outer cover concrete and outer reinforcement
$r_4$	Boundary of outer reinforcement and middle concrete
$r_3$	Boundary of middle concrete and inner reinforcement
$r_2$	Boundary of inner reinforcement and inner cover concrete
$r_1$	Inner surface of concrete lining

Basic Equations for Each Material

Outer Cover Concrete ( $r_5 \leqq r \leqq r_6$ )

$\begin{align} % displacement of outer concrete &u =\ \ \cfrac{1+\nu_c}{1-\nu_c} \alpha_c T \cfrac{r^2-r_5{}^2}{2 r} & &+ C_{co1}\cdot r + \cfrac{C_{co2}}{r} \\ % stress of outer concrete &\sigma_r =-\cfrac{E_c \alpha_c T}{1-\nu_c} \cfrac{r^2-r_5{}^2}{2 r^2} & &+ \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{co1}-\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{co2}}{r^2} \\ &\sigma_{\theta}=-\cfrac{E_c \alpha_c T}{1-\nu_c} \cfrac{r^2+r_5{}^2}{2 r^2} & &+ \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{co1}+\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{co2}}{r^2} \end{align}$

Outer Reinforcement ( $r_4 \leqq r \leqq r_5$ )

Middle Concrete ( $r_3 \leqq r \leqq r_4$ )

$\begin{align} % displacement of middle concrete &u =\ \ \cfrac{1+\nu_c}{1-\nu_c} \alpha_c T \cfrac{r^2-r_3{}^2}{2 r} & &+ C_{cm1}\cdot r + \cfrac{C_{cm2}}{r} \\ % stress of middle concrete &\sigma_r =-\cfrac{E_c \alpha_c T}{1-\nu_c} \cfrac{r^2-r_3{}^2}{2 r^2} & &+ \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{cm1}-\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{cm2}}{r^2} \\ &\sigma_{\theta}=-\cfrac{E_c \alpha_c T}{1-\nu_c} \cfrac{r^2+r_3{}^2}{2 r^2} & &+ \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{cm1}+\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{cm2}}{r^2} \end{align}$

Inner Reinforcement ( $r_2 \leqq r \leqq r_3$ )

Inner Cover Concrete ( $r_1 \leqq r \leqq r_2$ )

$\begin{align} % displacement of inner concrete &u =\ \ \cfrac{1+\nu_c}{1-\nu_c} \alpha_c T \cfrac{r^2-r_1{}^2}{2 r} & &+ C_{ci1}\cdot r + \cfrac{C_{ci2}}{r} \\ % stress of inner concrete &\sigma_r =-\cfrac{E_c \alpha_c T}{1-\nu_c} \cfrac{r^2-r_1{}^2}{2 r^2} & &+ \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{ci1}-\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{ci2}}{r^2} \\ &\sigma_{\theta}=-\cfrac{E_c \alpha_c T}{1-\nu_c} \cfrac{r^2+r_1{}^2}{2 r^2} & &+ \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{ci1}+\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{ci2}}{r^2} \end{align}$

Simultaneous linear equations

To fix un-known parameters, simultaneous linear equations will be created considering boundary conditions shown below.

Stress in the radial direction at outer surface of concrete lining is equal to the external pressure.
Displacement and stress in the radial direction have continuity at the boundaries of each material.
Stress in the radial direction at inner surface of concrete lining is equal to zero.
Positive sign of the external pressure $P_b$ is toward inner direction.

$\begin{align} &1.\quad r=r_6 \;\text{(stress)}& &\quad \\ &\left\{ % stress of outer concrete \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{co1}-\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{co2}}{r_6{}^2} \right\} & &=-P_b+\cfrac{E_c \alpha_c T}{1-\nu_c} \cfrac{r_6{}^2-r_5{}^2}{2 r_6{}^2} \\ % &2.\quad r=r_5 \;\text{(disp.)}& &\quad \\ &\left\{ % displacement of outer concrete C_{co1}\cdot r_5 + \cfrac{C_{co2}}{r_5} \right\}-\left\{ % displacement of outer rebar C_{so1}\cdot r_5 + \cfrac{C_{so2}}{r_5} \right\} & &=\cfrac{1+\nu_s}{1-\nu_s} \alpha_s T \cfrac{r_5{}^2-r_4{}^2}{2 r_5} \\ % &3.\quad r=r_5 \;\text{(stress)}& &\quad \\ &\left\{ % stress of outer concrete \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{co1}-\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{co2}}{r_5{}^2} \right\}-\left\{ % stress of outer rebar \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{so1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{so2}}{r_5{}^2} \right\} & &=-\cfrac{E_s \alpha_s T}{1-\nu_s} \cfrac{r_5{}^2-r_4{}^2}{2 r_5{}^2} \\ % &4.\quad r=r_4 \;\text{(disp.)}& &\quad \\ &\left\{ % displacement of outer rebar C_{so1}\cdot r_4 + \cfrac{C_{so2}}{r_4} \right\}-\left\{ % displacement of middle concrete C_{cm1}\cdot r_4 + \cfrac{C_{cm2}}{r_4} \right\} & &=\cfrac{1+\nu_c}{1-\nu_c} \alpha_c T \cfrac{r_4{}^2-r_3{}^2}{2 r_4} \\ % &5.\quad r=r_4 \;\text{(stress)}& &\quad \\ &\left\{ % stress of outer rebar \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{so1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{so2}}{r_4{}^2} \right\}-\left\{ % stress of middle concrete \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{cm1}-\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{cm2}}{r_4{}^2} \right\} & &=-\cfrac{E_c \alpha_c T}{1-\nu_c} \cfrac{r_4{}^2-r_3{}^2}{2 r_4{}^2} \\ % &6.\quad r=r_3 \;\text{(disp.)}& &\quad \\ &\left\{ % displacement of middle concrete C_{cm1}\cdot r_3 + \cfrac{C_{cm2}}{r_3} \right\}-\left\{ % displacement of inner rebar C_{si1}\cdot r_3 + \cfrac{C_{si2}}{r_3} \right\} & &=\cfrac{1+\nu_s}{1-\nu_s} \alpha_s T \cfrac{r_3{}^2-r_2{}^2}{2 r_3} \\ % &7.\quad r=r_3 \;\text{(stress)}& &\quad \\ &\left\{ % stress of middle concrete \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{cm1}-\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{cm2}}{r_3{}^2} \right\}-\left\{ % stress of inner rebar \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{si1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{si2}}{r_3{}^2} \right\} & &=-\cfrac{E_s \alpha_s T}{1-\nu_s} \cfrac{r_3{}^2-r_2{}^2}{2 r_3{}^2} \\ % &8.\quad r=r_2 \;\text{(disp.)}& &\quad \\ &\left\{ % displacement of inner rebar C_{si1}\cdot r_2 + \cfrac{C_{si2}}{r_2} \right\}-\left\{ % displacement of inner concrete C_{ci1}\cdot r_2 + \cfrac{C_{ci2}}{r_2} \right\} & &=\cfrac{1+\nu_c}{1-\nu_c} \alpha_c T \cfrac{r_2{}^2-r_1{}^2}{2 r_2} \\ % &9.\quad r=r_2 \;\text{(stress)}& &\quad \\ &\left\{ % stress of inner rebar \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{si1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{si2}}{r_2{}^2} \right\}-\left\{ % stress of inner concrete \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{ci1}-\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{ci2}}{r_2{}^2} \right\} & &=-\cfrac{E_c \alpha_c T}{1-\nu_c} \cfrac{r_2{}^2-r_1{}^2}{2 r_2{}^2} \\ % &10.\quad r=r_1 \;\text{(stress)}& &\quad \\ &\left\{ % stress of inner concrete \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{ci1}-\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{ci2}}{r_1^2} \right\} & &=0 \end{align}$

Matrix Expression of Simultaneous linear equations for Computer Programing

$\begin{equation} \begin{bmatrix} a_{1,1} & a_{1,2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} & 0 & 0 & 0 & 0 & 0 & 0 \\ a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & a_{4,3} & a_{4,4} & a_{4,5} & a_{4,6} & 0 & 0 & 0 & 0 \\ 0 & 0 & a_{5,3} & a_{5,4} & a_{5,5} & a_{5,6} & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & a_{6,5} & a_{6,6} & a_{6,7} & a_{6,8} & 0 & 0 \\ 0 & 0 & 0 & 0 & a_{7,5} & a_{7,6} & a_{7,7} & a_{7,8} & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & a_{8,7} & a_{8,8} & a_{ 8,9} & a_{ 8,10} \\ 0 & 0 & 0 & 0 & 0 & 0 & a_{9,7} & a_{9,8} & a_{ 9,9} & a_{ 9,10} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & a_{10,9} & a_{10,10} \end{bmatrix} \begin{Bmatrix} C_{co1} \\ C_{co2} \\ C_{so1} \\ C_{so2} \\ C_{cm1} \\ C_{cm2} \\ C_{si1} \\ C_{si2} \\ C_{ci1} \\ C_{ci2} \\ \end{Bmatrix} = \begin{Bmatrix} -P_b+\frac{E_c \alpha_c T}{1-\nu_c} \frac{r_6{}^2-r_5{}^2}{2 r_6{}^2} \\ \frac{1+\nu_s}{1-\nu_s} \alpha_s T \frac{r_5{}^2-r_4{}^2}{2 r_5} \\ -\frac{E_s \alpha_s T}{1-\nu_s} \frac{r_5{}^2-r_4{}^2}{2 r_5{}^2} \\ \frac{1+\nu_c}{1-\nu_c} \alpha_c T \frac{r_4{}^2-r_3{}^2}{2 r_4} \\ -\frac{E_c \alpha_c T}{1-\nu_c} \frac{r_4{}^2-r_3{}^2}{2 r_4{}^2} \\ \frac{1+\nu_s}{1-\nu_s} \alpha_s T \frac{r_3{}^2-r_2{}^2}{2 r_3} \\ -\frac{E_s \alpha_s T}{1-\nu_s} \frac{r_3{}^2-r_2{}^2}{2 r_3{}^2} \\ \frac{1+\nu_c}{1-\nu_c} \alpha_c T \frac{r_2{}^2-r_1{}^2}{2 r_2} \\ -\frac{E_c \alpha_c T}{1-\nu_c} \frac{r_2{}^2-r_1{}^2}{2 r_2{}^2} \\ 0 \end{Bmatrix} \end{equation}$

$\begin{align} &a_{1,1}=\cfrac{E_c}{(1+\nu_c)(1-2\nu_c)} & &a_{1,2}=-\cfrac{E_c}{(1+\nu_c)r_6{}^2} & &a_{1,3}=0 & &a_{1,4}=0 \\ % &a_{2,1}=r_5 & &a_{2,2}=\cfrac{1}{r_5} & &a_{2,3}=-r_5 & &a_{2,4}=-\cfrac{1}{r_5} \\ &a_{3,1}=\cfrac{E_c}{(1+\nu_c)(1-2\nu_c)} & &a_{3,2}=-\cfrac{E_c}{(1+\nu_c)r_5{}^2} & &a_{3,3}=-\cfrac{E_s}{(1+\nu_s)(1-2\nu_s)} & &a_{3,4}=\cfrac{E_s}{(1+\nu_s)r_5{}^2} \\ % &a_{4,3}=r_4 & &a_{4,4}=\cfrac{1}{r_4} & &a_{4,5}=-r_4 & &a_{4,6}=-\cfrac{1}{r_4} \\ &a_{5,3}=\cfrac{E_s}{(1+\nu_s)(1-2\nu_s)} & &a_{5,4}=-\cfrac{E_s}{(1+\nu_s)r_4{}^2} & &a_{5,5}=-\cfrac{E_c}{(1+\nu_c)(1-2\nu_c)} & &a_{5,6}=\cfrac{E_c}{(1+\nu_c)r_4{}^2} \\ % &a_{6,5}=r_3 & &a_{6,6}=\cfrac{1}{r_3} & &a_{6,7}=-r_3 & &a_{6,8}=-\cfrac{1}{r_3} \\ &a_{7,5}=\cfrac{E_c}{(1+\nu_c)(1-2\nu_c)} & &a_{7,6}=-\cfrac{E_c}{(1+\nu_c)r_3{}^2} & &a_{7,7}=-\cfrac{E_s}{(1+\nu_s)(1-2\nu_s)} & &a_{7,8}=\cfrac{E_s}{(1+\nu_c)r_3{}^2} \\ % &a_{8,7}=r_2 & &a_{8,8}=\cfrac{1}{r_2} & &a_{8,9}=-r_2 & &a_{8,10}=-\cfrac{1}{r_2} \\ &a_{9,7}=\cfrac{E_s}{(1+\nu_s)(1-2\nu_s)} & &a_{9,8}=-\cfrac{E_s}{(1+\nu_s)r_2{}^2} & &a_{9,9}=-\cfrac{E_c}{(1+\nu_c)(1-2\nu_c)} & &a_{9,10}=\cfrac{E_c}{(1+\nu_c)r_2{}^2} \\ % &a_{10,7}=0 & &a_{10,8}=0 & &a_{10,9}=\cfrac{E_c}{(1+\nu_c)(1-2\nu_c)} & &a_{10,10}=-\cfrac{E_c}{(1+\nu_c)r_1{}^2} \end{align}$

Single Reinforcement Section

Components of Model

Components of Model
Concrete (outer cover)	Elastic material. Thermal stress is considered.
Inner Reinforcement	Elastic material. Thermal stress is considered.
Concrete (Inner cover)	Elastic material. Thermal stress is considered.

Coordinate of Boundary in Radial Direction
$r_4$	External surface of concrete lining
$r_3$	Boundary of reinforcement and outer cover concrete
$r_2$	Boundary of reinforcement and inner cover concrete
$r_1$	Inner surface of concrete lining

Basic Equations for Each Material

Outer Cover Concrete ( $r_3 \leqq r \leqq r_4$ )

$\begin{align} % displacement of outer concrete &u =\ \ \cfrac{1+\nu_c}{1-\nu_c} \alpha_c T \cfrac{r^2-r_3{}^2}{2 r} & &+ C_{co1}\cdot r + \cfrac{C_{co2}}{r} \\ % stress of outer concrete &\sigma_r =-\cfrac{E_c \alpha_c T}{1-\nu_c} \cfrac{r^2-r_3{}^2}{2 r^2} & &+ \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{co1}-\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{co2}}{r^2} \\ &\sigma_{\theta}=-\cfrac{E_c \alpha_c T}{1-\nu_c} \cfrac{r^2+r_3{}^2}{2 r^2} & &+ \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{co1}+\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{co2}}{r^2} \end{align}$

Reinforcement ( $r_2 \leqq r \leqq r_3$ )

Inner Cover Concrete ( $r_1 \leqq r \leqq r_2$ )

Simultaneous linear equations

To fix un-known parameters, simultaneous linear equations will be created considering boundary conditions shown below.

Stress in the radial direction at outer surface of concrete lining is equal to the external pressure.
Displacement and stress in the radial direction have continuity at the boundaries of each material.
Stress in the radial direction at inner surface of concrete lining is equal to zero.
Positive sign of the external pressure $P_b$ is toward inner direction.

$\begin{align} &1.\quad r=r_4 \;\text{(stress)}& &\quad \\ &\left\{ % stress of outer concrete \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{co1}-\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{co2}}{r_4{}^2} \right\} & &=-P_b+\cfrac{E_c \alpha_c T}{1-\nu_c} \cfrac{r_4{}^2-r_3{}^2}{2 r_4{}^2} \\ % &2.\quad r=r_3 \;\text{(disp.)}& &\quad\\ &\left\{ % displacement of outer concrete C_{co1}\cdot r_3 + \cfrac{C_{co2}}{r_3} \right\}-\left\{ % displacement of inner rebar C_{si1}\cdot r_3 + \cfrac{C_{si2}}{r_3} \right\} & &=\cfrac{1+\nu_s}{1-\nu_s} \alpha_s T \cfrac{r_3{}^2-r_2{}^2}{2 r_3} \\ % &3.\quad r=r_3 \;\text{(stress)}& &\quad \\ &\left\{ % stress of outer concrete \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{co1}-\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{co2}}{r_3{}^2} \right\}-\left\{ % stress of inner rebar \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{si1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{si2}}{r_3{}^2} \right\} & &=-\cfrac{E_s \alpha_s T}{1-\nu_s} \cfrac{r_3{}^2-r_2{}^2}{2 r_3{}^2} \\ % &4.\quad r=r_2 \;\text{(disp.)}& &\quad \\ &\left\{ % displacement of inner rebar C_{si1}\cdot r_2 + \cfrac{C_{si2}}{r_2} \right\}-\left\{ % displacement of inner concrete C_{ci1}\cdot r_2 + \cfrac{C_{ci2}}{r_2} \right\} & &=\cfrac{1+\nu_c}{1-\nu_c} \alpha_c T \cfrac{r_2{}^2-r_1{}^2}{2 r_2} \\ % &5.\quad r=r_2 \;\text{(stress)}& &\quad \\ &\left\{ % stress of inner rebar \cfrac{E_s}{(1+\nu_s)(1-2\nu_s)}\cdot C_{si1}-\cfrac{E_s}{(1+\nu_s)}\cfrac{C_{si2}}{r_2{}^2} \right\}-\left\{ % stress of inner concrete \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{ci1}-\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{ci2}}{r_2{}^2} \right\} & &=-\cfrac{E_c \alpha_c T}{1-\nu_c} \cfrac{r_2{}^2-r_1{}^2}{2 r_2{}^2} \\ % &6.\quad r=r_1 \;\text{(stress)}& &\quad \\ &\left\{ % stress of inner concrete \cfrac{E_c}{(1+\nu_c)(1-2\nu_c)}\cdot C_{ci1}-\cfrac{E_c}{(1+\nu_c)}\cfrac{C_{ci2}}{r_1^2} \right\} & &=0 \end{align}$

Matrix Expression of Simultaneous linear equations for Computer Programing

$\begin{equation} \begin{bmatrix} a_{1,1} & a_{1,2} & 0 & 0 & 0 & 0 \\ a_{2,1} & a_{2,2} & a_{2,3} & a_{2,4} & 0 & 0 \\ a_{3,1} & a_{3,2} & a_{3,3} & a_{3,4} & 0 & 0 \\ 0 & 0 & a_{4,3} & a_{4,4} & a_{4,5} & a_{4,6} \\ 0 & 0 & a_{5,3} & a_{5,4} & a_{5,5} & a_{5,6} \\ 0 & 0 & 0 & 0 & a_{6,5} & a_{6,6} \end{bmatrix} \begin{Bmatrix} C_{co1} \\ C_{co2} \\ C_{si1} \\ C_{si2} \\ C_{ci1} \\ C_{ci2} \\ \end{Bmatrix} = \begin{Bmatrix} -P_b+\frac{E_c \alpha_c T}{1-\nu_c} \frac{r_4{}^2-r_3{}^2}{2 r_4{}^2} \\ \frac{1+\nu_s}{1-\nu_s} \alpha_s T \frac{r_3{}^2-r_2{}^2}{2 r_3} \\ -\frac{E_s \alpha_s T}{1-\nu_s} \frac{r_3{}^2-r_2{}^2}{2 r_3{}^2} \\ \frac{1+\nu_c}{1-\nu_c} \alpha_c T \frac{r_2{}^2-r_1{}^2}{2 r_2} \\ -\frac{E_c \alpha_c T}{1-\nu_c} \frac{r_2{}^2-r_1{}^2}{2 r_2{}^2} \\ 0 \end{Bmatrix} \end{equation}$

$\begin{align} &a_{1,1}=\cfrac{E_c}{(1+\nu_c)(1-2\nu_c)} & &a_{1,2}=-\cfrac{E_c}{(1+\nu_c)r_4{}^2} & &a_{1,3}=0 & &a_{1,4}=0 \\ % &a_{2,1}=r_3 & &a_{2,2}=\cfrac{1}{r_3} & &a_{2,3}=-r_3 & &a_{2,4}=-\cfrac{1}{r_3} \\ &a_{3,1}=\cfrac{E_c}{(1+\nu_c)(1-2\nu_c)} & &a_{3,2}=-\cfrac{E_c}{(1+\nu_c)r_3{}^2} & &a_{3,3}=-\cfrac{E_s}{(1+\nu_s)(1-2\nu_s)} & &a_{3,4}=\cfrac{E_s}{(1+\nu_c)r_3{}^2} \\ % &a_{4,3}=r_2 & &a_{4,4}=\cfrac{1}{r_2} & &a_{4,5}=-r_2 & &a_{4,6}=-\cfrac{1}{r_2} \\ &a_{5,3}=\cfrac{E_s}{(1+\nu_s)(1-2\nu_s)} & &a_{5,4}=-\cfrac{E_s}{(1+\nu_s)r_2{}^2} & &a_{5,5}=-\cfrac{E_c}{(1+\nu_c)(1-2\nu_c)} & &a_{5,6}=\cfrac{E_c}{(1+\nu_c)r_2{}^2} \\ % &a_{6,3}=0 & &a_{6,4}=0 & &a_{6,5}=\cfrac{E_c}{(1+\nu_c)(1-2\nu_c)} & &a_{6,6}=-\cfrac{E_c}{(1+\nu_c)r_1{}^2} \end{align}$

Thank you.

記事の先頭に行く

2019-05-21

設計圧縮材のせん断と引張に対する安全率

設計

記事の最後に行く

Safety Factor for Shear Strength and Tensile Strength

The definitions of symbols are shown below in this discussion considering compression members such as concrete or rock material.

$c$	cohesion of material
$\phi$	internal friction angle of material
$\sigma_t$	uniaxial tensile strength of material (positive value)
$\sigma_1$	first principal stress
$\sigma_2$	second principal stress

f:id:damyarou:20190521131705j:plain

Referring above figure, following distances on $\sigma$ - $\tau$ plane can be defined.

$\begin{align} & D_1=\left(c-\cfrac{\sigma_1+\sigma_2}{2}\cdot \tan\phi\right)\cdot \cos\phi & \quad & d_1=D_1-\cfrac{\sigma_1-\sigma_2}{2} \\ & D_2=\sigma_t-\cfrac{\sigma_1+\sigma_2}{2} & \quad & d_2=\sigma_t-\sigma_1 \end{align}$

From above, the safety factors for shear strength and for tensile strength can be defined as follows, generally.

$\begin{align} &\text{For shear strength} & &SF_1=\cfrac{D_1}{D_1-d_1}\\ &\text{For tensile strength} & &SF_2=\cfrac{D_2}{D_2-d_2} \end{align}$

Next, a special case is considered. Since there is the following relationship which can be understood from a figure shown above, the safety factors should be defined for a case of $\sigma_1=\sigma_2$ .

$\begin{equation} D_1-d_1=D_2-d_2=\cfrac{\sigma_1-\sigma_2}{2} \end{equation}$

In this case, the safety factors for each strength will be defined as follows.

$\begin{equation} \text{For shear strength} \quad SF_1= \begin{cases} \infty & & (\sigma_1 \lt \sigma_t) \\ 0 & & (\sigma_t \leqq \sigma_1) \end{cases} \end{equation}$

$\begin{equation} \text{For tensile strength} \quad SF_2= \begin{cases} \infty & & (\sigma_1 \leqq 0) \\ \sigma_t / \sigma_1 & & (0 \lt \sigma_1) \end{cases} \end{equation}$

Finally, the minimum safety factor can be obtained as follow.

$\begin{equation} SF_{min}=min(SF_1, SF_2) \end{equation}$

作図プログラム

import numpy as np
import matplotlib.pyplot as plt


# global variables
sigt=3
sigc=-8
phi=np.arcsin((np.abs(sigc)-sigt)/(np.abs(sigc)+sigt))
cc=np.abs(sigc)*sigt/(np.abs(sigc)+sigt)/np.cos(phi)
aa=cc/np.tan(phi)
ps1=0.2*sigc
ps2=0.9*sigc


def barrow(x1,y1,x2,y2,ss,dd,ii):
    # (x1,y1),(x2,y2)：両矢印始点と終点
    # 両矢印始点と終点座標は文字列の傾きを考慮して決定する
    # ss：描画文字列
    # dd：文字列中央の寸法線からの距離
    # ii：１なら寸法線の上側，２なら寸法線の下側にテキストを描画
    col='#333333'
    arst='<->,head_width=3,head_length=10'
    aprop=dict(arrowstyle=arst,connectionstyle='arc3',facecolor=col,edgecolor=col,shrinkA=0,shrinkB=0,lw=1)
    plt.annotate('',
        xy=(x1,y1), xycoords='data',
        xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)
    al=np.sqrt((x2-x1)**2+(y2-y1)**2)
    cs=(x2-x1)/al
    sn=(y2-y1)/al
    rt=np.degrees(np.arccos(cs))
    if y2-y1<0: rt=-rt
    if ii==1:
        xs=0.5*(x1+x2)-dd*sn
        ys=0.5*(y1+y2)+dd*cs
        plt.text(xs,ys,ss,ha='center',va='center',rotation=rt)
    else:
        xs=0.5*(x1+x2)+dd*sn
        ys=0.5*(y1+y2)-dd*cs
        plt.text(xs,ys,ss,ha='center',va='center',rotation=rt)
        
    
def draw_axes(xmin,xmax,ymin,ymax):
    # 座標軸描画
    col='#000000'
    arst='<|-,head_width=0.3,head_length=0.5'
    aprop=dict(arrowstyle=arst,connectionstyle='arc3',facecolor=col,edgecolor=col,lw=1)
    plt.annotate(r'$\sigma$',
        xy=(xmin,0), xycoords='data',
        xytext=(xmax,0), textcoords='data', fontsize=20, ha='left',va='center', arrowprops=aprop)
    plt.annotate(r'$\tau$',
        xy=(0,0), xycoords='data',
        xytext=(0,ymax), textcoords='data', fontsize=20, ha='center',va='bottom', arrowprops=aprop)
    

def drawfig(xmin,xmax,ymin,ymax,fsz):
    # x-y axses
    draw_axes(xmin,xmax,0,ymax)
    # failure criterion
    x1=sigt   ; y1=0
    x2=sigt ; y2=cc-x2*np.tan(phi)
    x3=sigc ; y3=cc-x3*np.tan(phi)
    plt.plot([x1,x2,x3],[y1,y2,y3],'-',color='#ff0000',lw=3)
    xs=0.5*sigc; ys=cc-xs*np.tan(phi)
    strs="Mohr-Coulomb's Failure Criteria\n$\\tau=c-\sigma\cdot\\tan{\phi}$"
    plt.text(xs,ys+0.08*(ymax),strs,ha='center',va='center',rotation=-np.degrees(phi),fontsize=fsz-4)    
    theta=np.linspace(np.pi-phi,np.pi,20,endpoint=True)
    r1=0.07*(xmax-xmin)
    x=r1*np.cos(theta)+x2
    y=r1*np.sin(theta)+y2
    plt.plot(x,y,'-',color='#000000',lw=1)
    r2=r1*0.9
    x=r2*np.cos(theta)+x2
    y=r2*np.sin(theta)+y2
    plt.plot(x,y,'-',color='#000000',lw=1)
    plt.text(sigt-r1,y2,r'$\phi$',ha='right',va='bottom',color='#000000')
    plt.plot([0.5*sigt,x2],[y2,y2],'-',color='#000000',lw=1)
    # Mohr circle
    theta=np.linspace(0,np.pi,180,endpoint=True)
    rr=(ps1-ps2)/2
    psm=(ps1+ps2)/2
    x=rr*np.cos(theta)+psm
    y=rr*np.sin(theta)
    plt.plot(x,y,'-',color='#0000ff',lw=2)
    # points
    x=np.array([ps2,psm,ps1,sigt])
    y=np.array([0,0,0,0])
    plt.plot(x,y,'o',color='#000000',clip_on=False)

    ds=0.01*(ymax)
    plt.text(ps2-ds,0,r'$\sigma_2$',va='bottom',ha='right',fontsize=fsz)
    plt.text(psm,-ds,r'$\frac{\sigma_1+\sigma_2}{2}$',va='top',ha='center',fontsize=fsz)
    plt.text(ps1+ds,0,r'$\sigma_1$',va='bottom',ha='left',fontsize=fsz)
    plt.text(sigt+ds,0,r'$\sigma_t$',va='bottom',ha='left',fontsize=fsz)
    plt.plot([0],[cc],'o',color='#000000')
    plt.text(0.1,cc,r'$c$',va='bottom',ha='left',fontsize=fsz)

    ds=0.05*(ymax)
    ddl1=(cc-psm*np.tan(phi))*np.cos(phi)
    dds1=ddl1-rr
    ddl2=sigt-psm
    dds2=sigt-ps1
    dds=ds; dx=dds*np.cos(phi); dy=dds*np.sin(phi)
    x1=psm; y1=0; x2=psm+ddl1*np.sin(phi); y2=ddl1*np.cos(phi); barrow(x1,y1,x2,y2,r'$D_1$',ds,1)
    x1=psm+rr*np.sin(phi)+dx; y1=rr*np.cos(phi)-dy; x2=x2+dx; y2=y2-dy; barrow(x1,y1,x2,y2,r'$d_1$',ds,2)
    x1=psm-rr*np.cos(np.radians(20)); y1=rr*np.sin(np.radians(20)); x2=psm; y2=0; barrow(x1,y1,x2,y2,r'$\frac{\sigma_1-\sigma_2}{2}$',ds,1)
    dy=0.05*(ymax)
    x1=ps1; y1=-1*dy; x2=sigt; y2=y1; barrow(x1,y1,x2,y2,r'$d_2$',ds,2)
    x1=psm; y1=-3*dy; x2=sigt; y2=y1; barrow(x1,y1,x2,y2,r'$D_2$',ds,2)

    
    
def main():
    xmin=sigc-1; xmax=aa+1
    ymax=cc-xmin*np.tan(phi); ymin=-0.2*ymax
    fsz=20   # fontsize
    iw=10     # image width
    ih=iw*(ymax-ymin)/(xmax-xmin)
    plt.figure(figsize=(iw,ih),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.axis('off')
    plt.gca().set_aspect('equal',adjustable='box')

    drawfig(xmin,xmax,ymin,ymax,fsz)

    plt.tight_layout()
    fnameF='fig_mohr2.jpg'
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
    plt.show()


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

Thank you.

記事の先頭に行く

トンネル標準断面作図プログラム

両矢印を描く関数

片矢印を描く関数

指定した直線に対し文字列を描画する関数

プログラム全文

トンネル通水量計算・作図プログラム

貯水池水位容量曲線作図プログラム

洪水波形作図プログラム

Thank you.

Flood Routine 解析プログラム

Flood Routine 作図プログラム

Thank you.

概要

Conditions for Flood Routine Analysis

Hydrograph

Reservoir Storage Capacity Curve

Discharge Capacity Curve of Diversion Tunnel

Calculation of Discharge Capacity for Uniform Flow

Calculation of Discharge Capacity for Pressued Tunnel Flow

Result of Flood Routine Analysis

Thank you.

はじめに

ポイント

アノーテーション用矢印とテキスト描画

荷重用矢印

ひび割れ描画

塗りつぶし

ハッチングの場合

単色で塗りつぶす場合

円を描く

太字の描画

プログラム全文

Thank you.

内水圧を受けるRC圧力トンネル

出力

プログラム

外水圧を受けるRC圧力トンネル

出力

プログラム

Thank you.

概要

Design of RC Circular Pressure Tunnel

Basic Equations for Elastic Cylinder Model in Plane Strain State

Basic Equations for Isotropic Elastic Material in Plane Strain State

Stress - Strain Relationship

Strain - Displacement Relationship

Equilibrium Equation

General Expressions of Displacements and Stresses

Basic Equations for No-tension Material in Circumferential Direction

Model of RC Circular Tunnel under Internal Pressure

Double Reinforcement Section

Components of Model

Basic Equations for Each Material

Bedrock ()

Outer Cover Concrete ()

Outer Reinforcement ()

Middle Concrete ()

Inner Reinforcement ()

Inner Cover Concrete ()

Simultaneous linear equations

Matrix Expression of Simultaneous linear equations for Computer Programing

Single Reinforcement Section

Components of Model

Basic Equations for Each Material

Bedrock ()

Outer Cover Concrete ()

Reinforcement ()

Inner Cover Concrete ()

Simultaneous linear equations

Matrix Expression of Simultaneous linear equations for Computer Programing

Model of RC Circular Tunnel under External Pressure

Double Reinforcement Section

Components of Model

Basic Equations for Each Material

Outer Cover Concrete ()

Outer Reinforcement ()

Middle Concrete ()

Inner Reinforcement ()

Inner Cover Concrete ()

Simultaneous linear equations

Bedrock ( $r_6 \leqq r \leqq r_7$ )

Outer Cover Concrete ( $r_5 \leqq r \leqq r_6$ )

Outer Reinforcement ( $r_4 \leqq r \leqq r_5$ )

Middle Concrete ( $r_3 \leqq r \leqq r_4$ )

Inner Reinforcement ( $r_2 \leqq r \leqq r_3$ )

Inner Cover Concrete ( $r_1 \leqq r \leqq r_2$ )

Bedrock ( $r_4 \leqq r \leqq r_5$ )

Outer Cover Concrete ( $r_3 \leqq r \leqq r_4$ )

Reinforcement ( $r_2 \leqq r \leqq r_3$ )

Inner Cover Concrete ( $r_1 \leqq r \leqq r_2$ )

Outer Cover Concrete ( $r_5 \leqq r \leqq r_6$ )

Outer Reinforcement ( $r_4 \leqq r \leqq r_5$ )

Middle Concrete ( $r_3 \leqq r \leqq r_4$ )

Inner Reinforcement ( $r_2 \leqq r \leqq r_3$ )

Inner Cover Concrete ( $r_1 \leqq r \leqq r_2$ )

Outer Cover Concrete ( $r_3 \leqq r \leqq r_4$ )

Reinforcement ( $r_2 \leqq r \leqq r_3$ )

Inner Cover Concrete ( $r_1 \leqq r \leqq r_2$ )