設計 RC長方形断面の曲げ耐荷力計算
はじめに
マレーシア時代に使っていた、RC長方形断面の曲げ耐荷力計算プログラムを掲載しています。 部分安全係数、コンクリートおよび鉄筋の物性は BS 8110-1 に準拠したものとしています。
日本のコンクリート標準示方書での計算も、部分安全係数と物性値を変えればできるはずなので、後ほど試してみようと思います。
Ultimate Limit Strength Calculation of Rectangular RC member
Stress-strain curve
Stress-strain curves of concrete and rebar follow BS 8110-1.
The compressive stress and strain have positive sign for both of concrete and rebar.
Concrete
where, means a compressive strength of concrete based on cubic specimen following BS.
For example, C35 in BS means
(N/mm
).
Reinforcement
where, means a yield strength of rebar and usually
(N/mm
) is used.
| Stress-Strain Curves of Materials |
|---|
 |
Calculation of Axial Force and Bending Moment
Calculation formulas for axial force and bending moment are shown below. The moment is calculated with center point of the section.
| Model of Cross Section |
|---|
 |
Simplified Ultimate Limit State assuming a strain distribution in a section
At the ultimate limit state of a section, either compressive strain of concrete or tensile strain of rebar reachs its limit ultimate strain.
Considering above, simplified strain distribution at ultimate limit state of a section is assumed shown as above figure.
The calculation process is shown below:
- Set ultimate limit compressive strain of concrete as
- Set ultimate limit tensile strain of rebar as
(conveniently assumed value)
- Step-1: Fix the upper edge strain
. Vary the lower edge strain
from
to
- Step-2: Fix the lower edge strain
. Vary the upper edge strain
from
Obtained axial forces and bending moments from above process deam the N-M relationship at the ultimate limit state of a section.
| Concept of Calculation of Ultimate Strength |
|---|
 |
Test Running
Used properties of RC cross section and parameters are shown below.
| b=1000 (mm) | section width |
| h=1000 (mm) | section height |
| dt=100 (mm) | cover of tensile rebar (from tensile edge to rebar center) |
| dc=100 (mm) | cover of compressive rebar (from compressive edge to rebar center) |
| p1 (%) | ratio of tensile reinforcement sres to RC cross section area |
| p2 (%) | ratio of compressive reinforcement sres to RC cross section area |
| Concrete grade | C25 |
| Yield strength of reinforcement | 460 MPa |
In the following program, the strain difference between $\epsilon{cu}$ and $-3 \epsilon{sy}$ is divided by 20 (nn=20).
プログラム
自作関数のインポート
ここでは、曲げ耐荷力を計算する自作関数「im_rebar.py」を作成し、メインプログラムから呼ぶようにしています。 Jupyter 上で作業する場合、計算ケースが変わる場合は違うセルを用いますが、同じ部分を何度もコピーして使っていると、意に反して変更してしまい、トラブルの原因にもなるので、この方法を用いています。
自作関数が、Jupyter を実行するのと同じディレクトリにあれば、以下のように関数をインポートすれば使えます。
import im_rebar
自作関数が、Jupyter を実行するのと違うディレクトリにあれば、以下のように関数が格納されているディレクトリのフルパスを import sys で加えた後に関数をインポートします。
import sys; sys.path.append('/Users/kk/pypro') import im_rebar
負の曲げモーメント
引張鉄筋量と圧縮鉄筋量が等しい場合は負のモーメントは発生しませんが、これらが異なる場合、圧縮応力あるいは引張応力が大きい場合、鉄筋量の差により負のモーメントが発生します。この問題の処理として、このプログラムでは負のモーメントは無視しています。すなわち断面耐荷力の有効範囲は、計算上の正のモーメントの領域のみとしています。
自作関数（im_rebar.py）
# im_rebar.py import numpy as np def calmain(mh,fcu,bb,hh,asr,ys): # Properties of concrete rmc=1.5 # partial safety factor for concrete ecu=0.0035 # ultimate limit strain of concrete # properties of rebar fsy=460.0 # yield strength of rebar rms=1.05 # partial safety factor for rebar Es=2.0e5 # elastic modulus of steel # Division of section height yg=hh/2 # location of gravity center y=np.zeros(mh,dtype=np.float64) dy=hh/float(mh) y=np.linspace(0.5*dy,hh-0.5*dy,mh) # distance from tensile edge to each strip # Definition of strain for calculation ec_lim=1*ecu # assumed limit compressive strain et_lim=-3*fsy/rms/Es # assumed limit tensile strain nn=20 # number of division between ec_lim and et_lim _e=np.linspace(et_lim,ec_lim,nn+1) epsilon=_e[::-1] # strain values for calculation fn=np.zeros(2*(nn+1),dtype=np.float64) # memory of axial force fm=np.zeros(2*(nn+1),dtype=np.float64) # memory of bending moment for (k,ee) in enumerate(epsilon): _fn1,_fm1=cal_nm(ec_lim,ee,bb,hh,y,dy,asr,ys,yg,ecu,fcu,rmc,fsy,rms,Es) # Step1 (compressive strain fixed) _fn2,_fm2=cal_nm(ee,et_lim,bb,hh,y,dy,asr,ys,yg,ecu,fcu,rmc,fsy,rms,Es) # Step2 (tensile strain fixed) fn[k]=_fn1 fm[k]=_fm1 fn[k+nn+1]=_fn2 fm[k+nn+1]=_fm2 # Arrangement for output fnn=fn/bb/hh # axial force for plotting in stress unit fmm=fm/bb/hh**2 # bending moment for plotting in stress unit fnc,fnt,fm0=src(nn,fnn,fmm) # search the limit axial forces and pure bending moment _i=np.where(0<=fmm) fnnp=np.zeros(len(fnn[_i])+2,dtype=np.float64) fmmp=np.zeros(len(fmm[_i])+2,dtype=np.float64) fnnp[0]=fnc fmmp[0]=0 fnnp[1:len(fnn[_i])+1]=fnn[_i] fmmp[1:len(fmm[_i])+1]=fmm[_i] fnnp[len(fnn[_i])+1]=fnt fmmp[len(fmm[_i])+1]=0 return fnnp,fmmp,fm0 def sigc(ee,ecu,fcu,rmc): # Stress-strain curve of concrete scy=0.67*fcu/rmc ecy=2.4e-4*np.sqrt(fcu/rmc) sig=scy-scy/ecy**2*(ee-ecy)**2 if ee<=0.0: sig=0.0 if ecu<=ee: sig=0.0 if ecy<ee and ee<=ecu: sig=scy return sig def sigs(ee,fsy,rms,Es): # Stress-strain curve of Steel rebar ssy=fsy/rms esy=fsy/rms/Es sig=Es*ee if ee<-esy: sig=-ssy if esy<ee: sig=ssy return sig def cal_nm(e1,e2,bb,hh,y,dy,asr,ys,yg,ecu,fcu,rmc,fsy,rms,Es): # Calculation of axial force and bending moment fn=0.0 fm=0.0 sc=np.zeros(len(y),dtype=np.float64) # concrete stress eps=e1-(e1-e2)/hh*(hh-y) for i,ee in enumerate(eps): sc[i]=sigc(ee,ecu,fcu,rmc) ss=np.zeros(len(ys),dtype=np.float64) # rebar stress eps=e1-(e1-e2)/hh*(hh-ys) for i,ee in enumerate(eps): ss[i]=sigs(ee,fsy,rms,Es)-sigc(ee,ecu,fcu,rmc) fn=np.sum(sc*dy*bb)+np.sum(ss*asr) fm=np.sum(sc*dy*bb*(y-yg))+np.sum(ss*asr*(ys-yg)) return fn,fm def src(nn,fnn,fmm): # search limit axial forces and pure bending moment by interpolation fnc=fnn[0] fnt=fnn[-1] for k in range(0,2*(nn+1)-1): if fmm[k]<0.0 and 0.0<fmm[k+1]: fnc=fnn[k]+(fnn[k+1]-fnn[k])/(fmm[k+1]-fmm[k])*(-fmm[k]) if fmm[k+1]<0.0 and 0.0<fmm[k]: fnt=fnn[k]-(fnn[k]-fnn[k+1])/(fmm[k]-fmm[k+1])*(fmm[k]) if 0.0<fnn[k] and fnn[k+1]<0.0: # estimate ultimate pure bending moment fm0=fmm[k]-(fmm[k]-fmm[k+1])/(fnn[k]-fnn[k+1])*fnn[k] return fnc,fnt,fm0
メインプログラム
#============================================= # Ultimate strehgth of Rectangular RC section # (use module: im_rebar) #============================================= import numpy as np import matplotlib.pyplot as plt import time #import sys; sys.path.append('/Users/kk/pypro') import im_rebar def sarrow(x1,y1,x2,y2,ss,fsz): sv=0 plt.annotate( ss, xy=(x1,y1), xycoords='data', xytext=(x2,y2), textcoords='data', fontsize=fsz, ha='left', va='center', arrowprops=dict(shrink=sv,width=0.1,headwidth=5,headlength=10, connectionstyle='arc3',facecolor='#777777',edgecolor='#777777')) def main(): start=time.time() bb=1000 hh=1000 mh=int(hh/20) fcu=25.0 # design strength of concrete lfile=['fig_nm0.jpg','fig_nm1.jpg','fig_nm2.jpg','fig_nm3.jpg'] lco=['#cccccc','#bbbbbb','#aaaaaa','#999999','#888888','#777777','#666666','#444444','#222222','#000000"'] pp1=np.array([0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0]) # tensile bar pp2=np.array([0.0, 1.0, 2.0, 3.0]) # compressive bar fsz=16 xmin=0; xmax=20; dx=5 ymin=-40; ymax=40; dy=10 for fnameF,p2 in zip(lfile,pp2): tstr='B x H = 1000 x 1000 (C{0:.0f}, p2 = {1:.1f} %)'.format(fcu,p2) fig=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.xlabel('Bending moment $M\: / \:bh^2$ (MPa)') plt.ylabel('Axial force $N\: / \:bh$ (MPa)') plt.xticks(np.arange(xmin,xmax+dx,dx)) plt.yticks(np.arange(ymin,ymax+dy,dy)) plt.grid(which='major',color='#777777',linestyle='dashed') for p1,col in zip(pp1,lco): asr=np.array([p1*bb*hh*0.01,p2*bb*hh*0.01]) ys =np.array([100,hh-100]) fnp,fmp,fm0=im_rebar.calmain(mh,fcu,bb,hh,asr,ys) plt.plot(fmp,fnp,'-',color=col,label='{0:.1f}'.format(p1)) if p1==0.5: plt.fill_betweenx(fnp, 0, fmp, color='#ffffcc') ss='Available area\nfor p1=0.5%' sarrow(np.mean(fmp),np.mean(fnp),xmin+0.40*(xmax-xmin),ymin+0.85*(ymax-ymin),ss,fsz-2) plt.legend(loc='lower right',fontsize=fsz-2,title='p1 (%)').get_title().set_fontsize(fsz-2) plt.title(tstr,loc='left',fontsize=fsz-2) plt.tight_layout() plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1) plt.show() print(time.time()-start) #============== # Execution #============== if __name__ == '__main__': main()
Thank you.
設計 確率雨量の算定（４）追補
平方根指数型最大値分布による確率雨量計算
最近ある方から、昔作った文書http://civilyarou.web.fc2.com/WANtaroHP_html5_win/f90_ENGI/dir_HFA/suimon.pdfの平方根指数型最大値分布に関連して、メールを頂いた。 そのころは、C あるいは Fortran90 を使ってプログラムを作っていたのだが、最近使っている Python で書いてみたくなった。 このため、計算理論は焼き直しであるが、あたらしく Python でプログラムを作ってみた。 ２変数確率分布であるが、未知パラメータを算定するには工夫が必要である。
プログラムとしては、観測データからプロッティングポジションを通さずに直接未知パラメータを計算する従来法と、プロッティングポジションを通して確率分布関数に含まれる未知パラメータを、Scipy の非線形最小二乗法で求める方法の、２種類を紹介している。
計算理論
SQRT exponential-type distribution of maximum (SQRT-ET)
Probabulity Density Function
Cumulative Distribution Function
Quantile 
for non-exceedance probability 
In above, the conditions of and
shall be satisfied, because the terms in logarithmic function and square-root have to be positive values.
Calculation Method of 
Since the relationship between and
is non-linear, Newton-Raphson Method is used to obtain the value of
at non-exceedance probability of
.
Following functions are considered, where
is the differential of function
.
To obtain the value of at
, iterative calculation will be carried out using following equation.
As an initial value of ,
can be used.
Estimation of Parameters (Maximum Likelihood Method)
The parameters and
are determined with a condition that a log-likelihoog function
shown below has the maximum value.
From the conditions that a partial differential of by
becomes
and a partial differential of
by
becomes
,
following relations can be obtained.
Since a function has the maximum value at the condition of
, the parameter
can be obtained using Bisection Method considering following equation.
It shall be noted that a relation of is everytime true, but
has to satisfy the following relation to ensure the condition of
.
Above relation can be used as the lowest limit of in Bisection Method.
Regarding the highest limit of
for Bisection Method, it is no information.
Therefore, following method in which searching range of solution is moved with an interval of 0.5 is applied for the calculation.
# Sample Code for Bisection method by Python # (xx : vector of observed values) # (bb : unknown scalar parameter to be solved) def func(bb,xx): # function for Bisection Method nn=len(xx) a1=(np.sum(np.sqrt(bb*xx))-2*nn)/np.sum((bb*xx)*np.exp(-np.sqrt(bb*xx))) a2=nn/(np.sum(np.exp(-np.sqrt(bb*xx)))+np.sum(np.sqrt(bb*xx)*np.exp(-np.sqrt(bb*xx)))) return a1-a2 imax=100 # maximum number of iterations err=1.0e-10 # allowable error of the zero value x1=(2*len(xx)/np.sum(np.sqrt(xx)))**2 x2=x1+0.5 for i in range(0,imax): xi=0.5*(x1+x2) f1=func(x1,xx) f2=func(x2,xx) fi=func(xi,xx) if 0<f1*f2: x1=x2 x2=x1+0.5 elif abs(x2-x1) < err: break # converged else: if f1*fi < 0.0: x2=xi if fi*f2 < 0.0: x1=xi return xi
結果出力
用いたデータは、「設計 確率雨量の推定（３）」で用いた１つである、宮崎の日最大雨量データである。
| 従来法 | Scipy 非線形回帰 |
|---|---|
 |
 |
| Year | 2 | 5 | 10 | 20 | 50 | 100 | 200 | 500 | 1000 | 2000 | 5000 | 10000 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Conventional | 166.7 | 230.3 | 277.4 | 326.2 | 394.8 | 450.0 | 508.4 | 590.4 | 656.1 | 725.1 | 821.2 | 897.6 |
| Scipy leasrsq | 166.7 | 237.5 | 290.4 | 345.6 | 423.5 | 486.5 | 553.3 | 647.4 | 723.1 | 802.6 | 913.7 | 1002.2 |
従来法によるプログラム
上述した、二分法により未知パラメータ および
を計算する方法によるプログラムを以下に示す。
プロッティング・ポジション公式は、Hazen 公式としている。
import numpy as np import matplotlib.pyplot as plt def qqplot(xx,ye,ss): fsz=12 xmin=0; xmax=700; dx=100 ymin=0; ymax=700; dy=100 plt.figure(figsize=(4,4),facecolor='w') plt.rcParams['font.size']=fsz plt.rcParams['font.family']='sans-serif' plt.xlim([xmin,xmax]) plt.ylim([ymin,ymax]) plt.xlabel('Observed (mm/day)') plt.ylabel('Predicted (mm/day)') plt.xticks(np.arange(xmin,xmax+dx,dx)) plt.yticks(np.arange(ymin,ymax+dy,dy)) plt.grid(color='#999999',linestyle='solid') plt.gca().set_aspect('equal',adjustable='box') plt.plot(xx,ye,'o',clip_on=False) plt.plot([xmin,xmax],[ymin,ymax],'-',color='#ff0000',lw=1) xs=xmin+0.05*(xmax-xmin) ys=ymin+0.95*(ymax-ymin) plt.text(xs,ys,ss,va='top',ha='left',fontsize=fsz,fontname='Ricty Diminished') plt.tight_layout() fnameF='fig_qq_M_org.jpg' plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1) plt.show() def sqrtet(xx,pp,pyear): def est_sqrtet(aa,bb,pp,x0): # estimation of quantile by Newton-Raphson Method def gx(x,aa,p): ggx=np.log(1+x)-x-np.log(-1/aa*np.log(p)) gdx=1/(1+x)-1 return ggx,gdx imax=100 eps=1.0e-10 tp=np.zeros(len(pp),dtype=np.float64) x1=x0 for i,p in enumerate(pp): for j in range(imax): ggx,gdx=gx(x1,aa,p) x2=x1-ggx/gdx if np.abs(x2-x1)<eps: break x1=x2 tp[i]=x2 return tp**2/bb def bisection(xx): # estimation of coefficient bb by Bisection Method def func(bb,xx): # function for Bisection Method nn=len(xx) a1=(np.sum(np.sqrt(bb*xx))-2*nn)/np.sum((bb*xx)*np.exp(-np.sqrt(bb*xx))) a2=nn/(np.sum(np.exp(-np.sqrt(bb*xx)))+np.sum(np.sqrt(bb*xx)*np.exp(-np.sqrt(bb*xx)))) return a1-a2 imax=100 # maximum number of iterations err=1.0e-10 # allowable error of the zero value x1=(2*len(xx)/np.sum(np.sqrt(xx)))**2 x2=x1+0.5 for i in range(0,imax): xi=0.5*(x1+x2) f1=func(x1,xx) f2=func(x2,xx) fi=func(xi,xx) if 0<f1*f2: x1=x2 x2=x1+0.5 elif abs(x2-x1) < err: break # converged else: if f1*fi < 0.0: x2=xi if fi*f2 < 0.0: x1=xi return xi bb=bisection(xx) aa=(np.sum(np.sqrt(bb*xx))-2*len(xx))/np.sum((bb*xx)*np.exp(-np.sqrt(bb*xx))) # estimation x0=np.sqrt(bb*xx[-1]) ye=est_sqrtet(aa,bb,pp,x0) rr=np.corrcoef(xx, ye)[0][1] # probable value prob=1.0-1.0/pyear yp=est_sqrtet(aa,bb,prob,x0) res=np.array([aa,bb,rr]) ss='SQRT-ET\na={0:8.3f}\nb={1:8.3f}\nr={2:8.3f}\nn={3:4.0f}'.format(res[0],res[1],res[2],len(xx)) return res,ye,yp,ss def main(): pyear=np.array([2,5,10,20,50,100,200,500,1000,2000,5000,10000]) xx=np.array([490.2, 175.7, 285.6, 152.9, 199.1, 98.4, 135.3, 188.3, 143.1, 157.2, 165.1, 119.5, 175.9, 198.7, 293.0, 122.4, 104.3, 154.8, 112.0, 305.8, 108.6, 178.8, 135.3, 136.5, 142.9, 175.4, 197.2, 165.0, 232.5, 97.4, 328.9, 316.7, 125.1, 194.9, 232.2, 182.8, 298.1, 115.9, 175.4, 180.1, 99.8, 321.7, 146.7, 156.7, 191.5, 163.2, 83.6, 217.9, 430.4, 143.7, 137.6, 132.8, 202.9, 587.2, 157.5, 170.1, 213.4, 301.9, 112.9, 205.2, 168.0, 129.5, 175.8, 155.6, 221.4, 341.9, 171.5, 181.1, 173.0, 151.5, 116.3, 145.4, 124.6, 172.7, 135.0, 197.0, 124.8, 172.6, 207.7, 131.1, 387.3, 98.9, 242.0, 112.5, 100.0, 200.0, 186.5, 108.0, 132.5, 162.5, 121.0, 155.5, 93.5, 286.5, 167.0, 103.5, 141.5, 315.5, 102.5, 207.5, 86.0, 233.5, 146.5, 228.5, 437.5, 108.5, 120.0, 317.5, 139.5, 203.0, 132.5, 234.0, 144.5, 217.5, 133.0, 217.0, 91.0, 135.5, 188.0, 310.0, 155.0, 199.5, 280.5, 137.5, 125.5, 150.0, 247.5, 119.0, 200.0, 124.5, 222.0, 268.0]) # Miyazaki # coefficient for plotting position formula #alpha=0.0 # Weibull #alpha=0.4 # Cunnane alpha=0.5 # Hazen xx=np.sort(xx) # setting protting position n=len(xx) ii=np.arange(n)+1 pp=(ii-alpha)/(n+1-2*alpha) res,ye,yp,ss=sqrtet(xx,pp,pyear) qqplot(xx,ye,ss) s0='{0:4s}'.format('Year') s1='{0:4s}'.format('SQRT') for i in range(len(pyear)): s0=s0+',{0:6.0f}'.format(pyear[i]) s1=s1+',{0:6.1f}'.format(yp[i]) print(s0) print(s1) #============== # Execution #============== if __name__ == '__main__': main()
Scipyの非線形最小二乗法による回帰プログラム
未知パラメータ および
を、Scipy の非線形最小二乗法による回帰計算で求めるプログラムを示す。
ここで、非線形回帰を行うために
および
の初期値を入力する必要があるが、ここでは以下のように初期値を定めた。
の初期値
の初期値としては小さめの値を与えたいため、二分法の解の検索区間の最小値である以下の値を
の初期値とした。
の初期値
クオンタイル を求めるための
と
および
の関係式より、以下の関係を導ける。
上式において、大き目の を初期値として与えたいため、
、
としたときの
を初期値とした。すなわち、
プログラム全文
プログラム全文を以下に示す。プロッティング・ポジション公式は、Hazen 公式としている。
import numpy as np from scipy import optimize # for use of leastsq import matplotlib.pyplot as plt def qqplot(xx,ye,ss): fsz=12 xmin=0; xmax=700; dx=100 ymin=0; ymax=700; dy=100 plt.figure(figsize=(4,4),facecolor='w') plt.rcParams['font.size']=fsz plt.rcParams['font.family']='sans-serif' plt.xlim([xmin,xmax]) plt.ylim([ymin,ymax]) plt.xlabel('Observed (mm/day)') plt.ylabel('Predicted (mm/day)') plt.xticks(np.arange(xmin,xmax+dx,dx)) plt.yticks(np.arange(ymin,ymax+dy,dy)) plt.grid(color='#999999',linestyle='solid') plt.gca().set_aspect('equal',adjustable='box') plt.plot(xx,ye,'o',clip_on=False) plt.plot([xmin,xmax],[ymin,ymax],'-',color='#ff0000',lw=1) xs=xmin+0.05*(xmax-xmin) ys=ymin+0.95*(ymax-ymin) plt.text(xs,ys,ss,va='top',ha='left',fontsize=fsz,fontname='Ricty Diminished') plt.tight_layout() fnameF='fig_qq_M_reg.jpg' plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1) plt.show() def sqrtet(xx,pp,pyear): def est_sqrtet(aa,bb,pp,x0): # estimation of quantile by Newton-Raphson Method def fx(tp,aa,p): return np.log(1+tp)-tp-np.log(-1/aa*np.log(p)) tp=np.zeros(len(pp),dtype=np.float64) for i,p in enumerate(pp): root = optimize.newton(fx,x0,args=(aa,p)) tp[i]=root return tp**2/bb def func_sqrtet(parameter,xx,pp): # function for leastsq aa=parameter[0] bb=parameter[1] x0=np.sqrt(bb*xx[-1]) est=est_sqrtet(aa,bb,pp,x0) residual=xx-est return residual # parameter calculation by scipy.optimize.leastsq b0=(2*len(xx)/np.sum(np.sqrt(xx)))**2 t0=np.sqrt(b0*xx[-1]) a0=-np.log(1e-6)/np.exp(np.log(1+t0)-t0) npara=2; parameter0 = np.array([a0,b0]) result = optimize.leastsq(func_sqrtet,parameter0,args=(xx,pp)) aa=result[0][0] bb=result[0][1] # estimation x0=np.sqrt(bb*xx[-1]) ye=est_sqrtet(aa,bb,pp,x0) rr=np.corrcoef(xx, ye)[0][1] # probable value prob=1.0-1.0/pyear yp=est_sqrtet(aa,bb,prob,x0) res=np.array([aa,bb,rr]) ss='SQRT-ET\na={0:8.3f}\nb={1:8.3f}\nr={2:8.3f}\nn={3:4.0f}'.format(res[0],res[1],res[2],len(xx)) return res,ye,yp,ss def main(): pyear=np.array([2,5,10,20,50,100,200,500,1000,2000,5000,10000]) xx=np.array([490.2, 175.7, 285.6, 152.9, 199.1, 98.4, 135.3, 188.3, 143.1, 157.2, 165.1, 119.5, 175.9, 198.7, 293.0, 122.4, 104.3, 154.8, 112.0, 305.8, 108.6, 178.8, 135.3, 136.5, 142.9, 175.4, 197.2, 165.0, 232.5, 97.4, 328.9, 316.7, 125.1, 194.9, 232.2, 182.8, 298.1, 115.9, 175.4, 180.1, 99.8, 321.7, 146.7, 156.7, 191.5, 163.2, 83.6, 217.9, 430.4, 143.7, 137.6, 132.8, 202.9, 587.2, 157.5, 170.1, 213.4, 301.9, 112.9, 205.2, 168.0, 129.5, 175.8, 155.6, 221.4, 341.9, 171.5, 181.1, 173.0, 151.5, 116.3, 145.4, 124.6, 172.7, 135.0, 197.0, 124.8, 172.6, 207.7, 131.1, 387.3, 98.9, 242.0, 112.5, 100.0, 200.0, 186.5, 108.0, 132.5, 162.5, 121.0, 155.5, 93.5, 286.5, 167.0, 103.5, 141.5, 315.5, 102.5, 207.5, 86.0, 233.5, 146.5, 228.5, 437.5, 108.5, 120.0, 317.5, 139.5, 203.0, 132.5, 234.0, 144.5, 217.5, 133.0, 217.0, 91.0, 135.5, 188.0, 310.0, 155.0, 199.5, 280.5, 137.5, 125.5, 150.0, 247.5, 119.0, 200.0, 124.5, 222.0, 268.0]) # Miyazaki # coefficient for plotting position formula #alpha=0.0 # Weibull #alpha=0.4 # Cunnane alpha=0.5 # Hazen xx=np.sort(xx) # setting protting position n=len(xx) ii=np.arange(n)+1 pp=(ii-alpha)/(n+1-2*alpha) res,ye,yp,ss=sqrtet(xx,pp,pyear) qqplot(xx,ye,ss) s0='{0:4s}'.format('Year') s1='{0:4s}'.format('SQRT') for i in range(len(pyear)): s0=s0+',{0:6.0f}'.format(pyear[i]) s1=s1+',{0:6.1f}'.format(yp[i]) print(s0) print(s1) #============== # Execution #============== if __name__ == '__main__': main()
Thank you.
設計 仮排水路設計における Flood Routine の活用（3）
ここでは、「設計 仮排水路設計における Flood Routine の活用（２）」で述べた事項に関連する計算・作図プログラムのうち、以下の4つを掲載する。
- トンネル標準断面作図プログラム
- トンネル通水量計算・作図プログラム
- 貯水池水位容量曲線作図プログラム
- 洪水波形作図プログラム
トンネル標準断面作図プログラム
このような説明図を描くプログラムが意外に難しいというか面倒くさい。
断面形は、幅と高さが等しい複合円である。
この断面の場合、トンネル内幅を とすれば、上半円の半径は
、下半円の半径は
となり、インバート幅は
となる。
道路トンネルの場合、インバート幅をきっちりした数字とする場合が多く、下半円の半径が中途半端な数字となることが多いが、水路トンネルの場合のインバートの制約は、施工用重機が往来できればいい。このため、水路トンネルとしてのこの断面形状は、幾何学的にきれいで結構気に入っている。
トンネル形状は、中心を原点として、水平右の点から角度を指定して (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.
設計 仮排水路設計における 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.
設計 仮排水路設計における 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.
Reservoir Storage Capacity Curve
Following reservoir capacity curve is used as a given condition.
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.
| 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: Bending angle: |
Bending radius: Bending angle: |
| Curve-2 | Bending radius: Bending angle: |
Bending radius: Bending angle: |
| The distance between tunnels are set to 30m (center to center) | ||
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
| discharge | |
| flow section area | |
| average velocity | |
| manning's roughness coefficient (=0.015) | |
| hydraulic radius | |
| invert gradient |
Calculation of Discharge Capacity for Pressued Tunnel Flow
| discharge | |
| section area of diversion tunnel | |
| average velocity | |
| difference of water head between upstream and downstream of diversion tunnel | |
| head loss coefficient of entrance (=0.25) | |
| head loss coefficient of exit (=1.00) | |
| head loss coefficient of bending for curve-1 | |
| head loss coefficient of bending for curve-2 | |
| friction head loss coefficient | |
| Manning's roughness coefficient (=0.015) | |
| length of tunnel | |
| hydraulic radius of tunnel | |
| internal diameter of tunnel | |
| radius of curvature of bending | |
| inter angle of bending | |
| gravity acceleration (=9.8 m/s |
The calculated discharge capacity of diversion tunnels are shown below.
Result of Flood Routine Analysis
The result of flood routine analysis with parameters of tunnel diameter from 5m to 12m is shown below.
| 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 |
Thank you.
matplotlib RC円形圧力トンネルモデル図
はじめに
以下に示す、「設計 RC円形圧力トンネルの配筋設計（１）」で示した図を作成するプログラムである。
ポイント
アノーテーション用矢印とテキスト描画
- 説明用のボックス内テキストと矢印を描画する。
- 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)
ひび割れ描画
コンクリート内のひび割れは、サインカーブで定義。
まず、座標 -
間でサインカーブを定義し、これを座標変換しながら
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 で、 から
を分割して角度を指定して
座標を求め、これを折れ線で結んでいく。
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' で太字を描画できる。
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.
設計 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()
