内水圧を受ける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()