Python RC部材の設計
RC部材の設計をしたので流れを記録しておく。
リンク
- 平面骨組解析プログラムはこちら（Qiita）。
- 配筋計算プログラムの自作関数はこちら（はてな）。
出力
モデル図
断面力図
配鉄設計
プログラム
モデル図
import matplotlib.pyplot as plt import numpy as np from scipy import interpolate def barrow(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 drawfig(): tstr='Model of Tailrace Outlet Channel' fsz=12 xmin=-11 xmax=11 ymin=-8 ymax=2 iw=12 ih=(ymax-ymin)/(xmax-xmin)*iw 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.xlabel('x_direction (m)') plt.ylabel('y_direction (m)') plt.gca().spines['right'].set_visible(False) plt.gca().spines['top'].set_visible(False) plt.gca().yaxis.set_ticks_position('left') plt.gca().xaxis.set_ticks_position('bottom') plt.gca().set_aspect('equal',adjustable='box') plt.title(tstr,loc='left',fontsize=fsz) x=np.array([-3.75,3.75,3.75,2.75,2.75,-2.75,-2.75,-3.75,-3.75]) y=np.array([-5.0,-5.0,0.0,0.0,-4.0,-4.0,0.0,0.0,-5.0]) plt.fill(x,y,color='#eeeeee') plt.plot(x,y,'-',color='#000000',lw=1) x1=-3.75; y1=-1.0; x2=-2.75; y2=y1; barrow(x1,y1,x2,y2) xs=0.5*(x1+x2); ys=y1; plt.text(xs,ys,'{0:.3f}'.format(x2-x1),va='bottom',ha='center',fontsize=fsz) x1=-2.75; y1=-1.0; x2= 2.75; y2=y1; barrow(x1,y1,x2,y2) xs=0.5*(x1+x2); ys=y1; plt.text(xs,ys,'{0:.3f}'.format(x2-x1),va='bottom',ha='center',fontsize=fsz) x1= 2.75; y1=-1.0; x2= 3.75; y2=y1; barrow(x1,y1,x2,y2) xs=0.5*(x1+x2); ys=y1; plt.text(xs,ys,'{0:.3f}'.format(x2-x1),va='bottom',ha='center',fontsize=fsz) x1=-1.5; y1=-5.0; x2=x1; y2=-4.0; barrow(x1,y1,x2,y2) xs=x1; ys=0.5*(y1+y2); plt.text(xs,ys,'{0:.3f}'.format(y2-y1),va='center',ha='right',rotation=90,fontsize=fsz) x1=-1.5; y1=-4.0; x2=x1; y2=0.0; barrow(x1,y1,x2,y2) xs=x1; ys=0.5*(y1+y2); plt.text(xs,ys,'{0:.3f}'.format(y2-y1),va='center',ha='right',rotation=90,fontsize=fsz) x1=-4.5; y1=-1.3; x2=x1; y2=0.0; barrow(x1,y1,x2,y2) xs=x1; ys=0.5*(y1+y2); plt.text(xs,ys,'{0:.3f}'.format(y2-y1),va='center',ha='right',rotation=90,fontsize=fsz) # lateral water pressure x1=-3.75; p1=3.7*0.5 x=np.array([x1,x1,x1-p1,x1]) y=np.array([-5.0,-1.3,-5.0,-5.0]) plt.fill(x,y,color='#00ffff',alpha=0.4) plt.fill(-x,y,color='#00ffff',alpha=0.4) yy=np.array([-5.0,-4.5,-4.0,-3.5,-3.0,-2.5,-2.0]) f1 = interpolate.interp1d([-5.0,-1.3], [3.7,0]) uu = f1(yy)*0.5 xx=x1-uu vv=np.zeros(len(yy),dtype=np.float64) plt.quiver(xx,yy,uu,vv,angles='xy',scale_units='xy',scale=1,width=0.003) plt.quiver(-xx,yy,-uu,vv,angles='xy',scale_units='xy',scale=1,width=0.003) xs=-3.75; ys=-5.0; ss='{0:.1f}kN/m$^2$'.format(37) plt.text(xs,ys,ss,va='top',ha='right',fontsize=fsz) xs=3.75; ys=-5.0; ss='{0:.1f}kN/m$^2$'.format(37) plt.text(xs,ys,ss,va='top',ha='left',fontsize=fsz) xs=-3.75; ys=0; ss='Lateral\nwater\npressure' plt.text(xs,ys,ss,va='bottom',ha='right',fontsize=fsz) xs=-x1; ys=0 plt.text(xs,ys,ss,va='bottom',ha='left',fontsize=fsz) # uplift x=np.array([-3.75,3.75,3.75,-3.75,-3.75]) y=np.array([-5.0-p1,-5.0-p1,-5.0,-5.0,-5.0-p1]) plt.fill(x,y,color='#00ffff',alpha=0.4) xx=np.linspace(-3.75,3.75,16) yy=np.zeros(len(xx),dtype=np.float64)-5.0-3.7/2 uu=np.zeros(len(xx),dtype=np.float64) vv=np.zeros(len(xx),dtype=np.float64)+3.7/2 plt.quiver(xx,yy,uu,vv,angles='xy',scale_units='xy',scale=1,width=0.003) xs=0; ys=-5.0-p1; ss='Uplift {0:.1f}kN/m$^2$'.format(37) plt.text(xs,ys,ss,va='top',ha='center',fontsize=fsz) # lateral earth pressure k0=0.5; qq=1.0 db=5.0; pb=k0*(1.3*db+qq)*0.5 di=1.3; pi=k0*(2.3*di+qq)*0.5 dt=0.0; pt=k0*(2.3*dt+qq)*0.5 x1=-6.0 x=np.array([x1,x1,x1-pt,x1-pi,x1-pb,x1]) y=np.array([-5.0,0.0,0.0,-di,-5.0,-5.0]) plt.fill(x,y,color='#00ff00',alpha=0.4) plt.fill(-x,y,color='#00ff00',alpha=0.4) yy=np.array([-5.0,-4.5,-4.0,-3.5,-3.0,-2.5,-2.0,-1.5,-1.0,-0.5]) f1 = interpolate.interp1d([-5.0,-1.3,0], [pb,pi,pt]) uu = f1(yy) xx=x1-uu vv=np.zeros(len(yy),dtype=np.float64) plt.quiver(xx,yy,uu,vv,angles='xy',scale_units='xy',scale=1,width=0.003) plt.quiver(-xx,yy,-uu,vv,angles='xy',scale_units='xy',scale=1,width=0.003) xs=x1-pt; ys=dt; ss='{0:.1f}kN/m$^2$'.format(pt*2*10) plt.text(xs,ys,ss,va='center',ha='right',fontsize=fsz) xs=x1-pi; ys=-di; ss='{0:.1f}kN/m$^2$'.format(pi*2*10) plt.text(xs,ys,ss,va='center',ha='right',fontsize=fsz) xs=x1-pb; ys=-db; ss='{0:.1f}kN/m$^2$'.format(pb*2*10) plt.text(xs,ys,ss,va='center',ha='right',fontsize=fsz) xs=-(x1-pt); ys=dt; ss='{0:.1f}kN/m$^2$'.format(pt*2*10) plt.text(xs,ys,ss,va='center',ha='left',fontsize=fsz) xs=-(x1-pi); ys=-di; ss='{0:.1f}kN/m$^2$'.format(pi*2*10) plt.text(xs,ys,ss,va='center',ha='left',fontsize=fsz) xs=-(x1-pb); ys=-db; ss='{0:.1f}kN/m$^2$'.format(pb*2*10) plt.text(xs,ys,ss,va='center',ha='left',fontsize=fsz) xs=x1; ys=-5.5; ss='Lateral\nearth\npressure' plt.text(xs,ys,ss,va='top',ha='right',fontsize=fsz) xs=-x1; ys=-5.5 plt.text(xs,ys,ss,va='top',ha='left',fontsize=fsz) plt.tight_layout() fnameF='fig_model.jpg' plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1) plt.show() def main(): drawfig() if __name__ == '__main__': main()
断面力図
import matplotlib.pyplot as plt import numpy as np import sys def drawfig(disx,ax1,sh1,mo1,disy,ax2,sh2,mo2,x,y,node,nele): nmax=np.max([np.max(np.abs(ax1)),np.max(np.abs(ax2))]) smax=np.max([np.max(np.abs(sh1)),np.max(np.abs(sh2))]) mmax=np.max([np.max(np.abs(mo1)),np.max(np.abs(mo2))]) dmax=np.max([np.max(np.abs(disx)),np.max(np.abs(disy))]) fsz=10 xmin=-7 xmax=7 ymin=-7 ymax=4 scl_dis=1.0 scl_axi=1.0 scl_she=1.0 scl_mom=2.0 iw=10 ih=(ymax-ymin)/(xmax-xmin)*iw plt.figure(figsize=(iw,ih),facecolor='w') plt.rcParams['font.size']=fsz plt.rcParams['font.family']='sans-serif' for nnn in [1,2,3,4]: nplot=220+nnn plt.subplot(nplot) plt.xlim([xmin,xmax]) plt.ylim([ymin,ymax]) plt.xlabel('x_direction (m)') plt.ylabel('y_direction (m)') plt.gca().spines['right'].set_visible(False) plt.gca().spines['top'].set_visible(False) plt.gca().yaxis.set_ticks_position('left') plt.gca().xaxis.set_ticks_position('bottom') plt.gca().set_aspect('equal',adjustable='box') if nnn==1: # displacement tstr='Displacement mode' ls1='disx_max={0:8.3f}mm'.format(np.max(disx)) ls2='disx_min={0:8.3f}mm'.format(np.min(disx)) ls3='disy_max={0:8.3f}mm'.format(np.max(disy)) ls4='disy_min={0:8.3f}mm'.format(np.min(disy)) dx=x+disx/dmax*scl_dis dy=y+disy/dmax*scl_dis for ne in range(0,nele): n1=node[0,ne]-1 n2=node[1,ne]-1 plt.plot([x[n1],x[n2]],[y[n1],y[n2]],color='gray',linewidth=0.5) plt.plot([dx[n1],dx[n2]],[dy[n1],dy[n2]],color='black',linewidth=1) if nnn==2: # axial force diagram tstr='Axial force diagram' ls1='N_max={0:10.3f}kN'.format(np.max([np.max(ax1),np.max(ax2)])) ls2='N_min={0:10.3f}kN'.format(np.min([np.min(ax1),np.min(ax2)])) ls3='' ls4='' d1=ax1/nmax*scl_axi d2=ax2/nmax*scl_axi for ne in range(0,nele): x1,x2,x3,x4,y1,y2,y3,y4=calc(ne,node,x,y,d1,d2) if d1[ne]<=0.0: # compression plt.fill([x1,x2,x3,x4],[y1,y2,y3,y4],color='black',alpha=0.1) else: # tension plt.fill([x1,x2,x3,x4],[y1,y2,y3,y4],color='black',alpha=0.2) for ne in range(0,nele): n1=node[0,ne]-1 n2=node[1,ne]-1 plt.plot([x[n1],x[n2]],[y[n1],y[n2]],color='black',linewidth=0.5) if nnn==3: # shearing force tstr='Shear force diagram' ls1='S_max={0:10.3f}kN'.format(np.max([np.max(-sh1),np.max(-sh2)])) ls2='S_min={0:10.3f}kN'.format(np.min([np.min(-sh1),np.min(-sh2)])) ls3='' ls4='' d1=sh1/smax*scl_she d2=sh2/smax*scl_she for ne in range(0,nele): x1,x2,x3,x4,y1,y2,y3,y4=calc(ne,node,x,y,d1,d2) plt.fill([x1,x2,x3,x4],[y1,y2,y3,y4],color='black',alpha=0.1) for ne in range(0,nele): n1=node[0,ne]-1 n2=node[1,ne]-1 plt.plot([x[n1],x[n2]],[y[n1],y[n2]],color='black',linewidth=0.5) if nnn==4: # moment tstr='Bending moment diagram' ls1='M_max={0:10.3f}kN-m'.format(np.max([np.max(-mo1),np.max(-mo2)])) ls2='M_min={0:10.3f}kN-m'.format(np.min([np.min(-mo1),np.min(-mo2)])) ls3='' ls4='' d1=mo1/mmax*scl_mom d2=mo2/mmax*scl_mom for ne in range(0,nele): x1,x2,x3,x4,y1,y2,y3,y4=calc(ne,node,x,y,d1,d2) plt.fill([x1,x2,x3,x4],[y1,y2,y3,y4],color='black',alpha=0.1) for ne in range(0,nele): n1=node[0,ne]-1 n2=node[1,ne]-1 plt.plot([x[n1],x[n2]],[y[n1],y[n2]],color='black',linewidth=0.5) plt.plot(xmin,ymin,'.',label=ls1) plt.plot(xmin,ymin,'.',label=ls2) plt.plot(xmin,ymin,'.',label=ls3) plt.plot(xmin,ymin,'.',label=ls4) plt.legend(loc='upper right',numpoints=1,markerscale=0, frameon=False,prop={'family':'monospace','size':12}) plt.title(tstr,loc='left',fontsize=fsz) plt.tight_layout() fnameF='fig_force.jpg' plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1) plt.show() def calc(ne,node,x,y,d1,d2): i=node[0,ne]-1 j=node[1,ne]-1 x1=x[i] y1=y[i] x2=x[j] y2=y[j] al=np.sqrt((x2-x1)**2+(y2-y1)**2) theta=np.arccos((x2-x1)/al) x4=x1-d1[ne]*np.sin(theta) y4=y1+d1[ne]*np.cos(theta) x3=x2-d2[ne]*np.sin(theta) y3=y2+d2[ne]*np.cos(theta) return x1,x2,x3,x4,y1,y2,y3,y4 def main(): # data input #args = sys.argv #fnameR=args[1] # input data file fnameR='out_fem_outlet.txt' f=open(fnameR,'r') text=f.readline() text=f.readline() text=text.strip() text=text.split() npoin=int(text[0]) # Number of nodes nele =int(text[1]) # Number of elements nsec =int(text[2]) # Number of sections npfix=int(text[3]) # Number of restricted nodes nlod =int(text[4]) # Number of loaded nodes x =np.zeros(npoin,dtype=np.float64) # Coordinates of nodes y =np.zeros(npoin,dtype=np.float64) # Coordinates of nodes node=np.zeros([3,nele],dtype=np.int) # Node-element relationship disx=np.zeros(npoin,dtype=np.float64) # Coordinates of nodes disy=np.zeros(npoin,dtype=np.float64) # Coordinates of nodes ax1 =np.zeros(nele,dtype=np.float64) # Section force vector sh1 =np.zeros(nele,dtype=np.float64) # Section force vector mo1 =np.zeros(nele,dtype=np.float64) # Section force vector ax2 =np.zeros(nele,dtype=np.float64) # Section force vector sh2 =np.zeros(nele,dtype=np.float64) # Section force vector mo2 =np.zeros(nele,dtype=np.float64) # Section force vector text=f.readline() for i in range(0,nsec): text=f.readline() text=f.readline() for i in range(0,npoin): text=f.readline() text=text.strip() text=text.split() x[i]=float(text[1]) # x-coordinate y[i]=float(text[2]) # y-coordinate text=f.readline() for i in range(0,npfix): text=f.readline() text=f.readline() for i in range(nele): text=f.readline() text=text.strip() text=text.split() node[0,i]=int(text[1]) #node_1 node[1,i]=int(text[2]) #node_2 node[2,i]=int(text[2]) #mat-no text=f.readline() for i in range(0,npoin): text=f.readline() text=text.strip() text=text.split() disx[i]=float(text[1]) # displacement in x-direction disy[i]=float(text[2]) # displacement in y-direction text=f.readline() for i in range(0,nele): text=f.readline() text=text.strip() text=text.split() ax1[i]=-float(text[1]) # axial force at node-1 sh1[i]= float(text[2]) # shear force at node-1 mo1[i]= float(text[3]) # moment at node-1 ax2[i]= float(text[4]) # axial force at node-2 sh2[i]=-float(text[5]) # shear force at node-2 mo2[i]=-float(text[6]) # moment at node-2 f.close() disx=disx*1000 disy=disy*1000 ax1=ax1*10 sh1=sh1*10 mo1=mo1*10 ax2=ax2*10 sh2=sh2*10 mo2=mo2*10 drawfig(disx,ax1,sh1,mo1,disy,ax2,sh2,mo2,x,y,node,nele) if __name__ == '__main__': main()
配筋設計
プログラム中で用いている、自作関数 im_rebar.py については、 ここ：https://damyarou.hatenablog.com/entry/2019/06/02/172405 に紹介しています。
#============================================= # 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/data_kk/pypro') import im_rebar def drawfig(fnp,fmp,fm0,nreq,mreq,strt,fnameF): fsz=16 fig=plt.figure(figsize=(6,6),facecolor='w') plt.rcParams["font.size"] = fsz plt.rcParams['font.family'] ='sans-serif' plt.xlabel('Bending moment $M\: / \:bh^2$ (MPa)') plt.ylabel('Axial force $N\: / \:bh$ (MPa)') plt.grid(which='major',color='#777777',linestyle='solid') plt.plot(fmp,fnp,'-',color='#000080',label='Ultimate strength') plt.plot(mreq,nreq,'o',color='#ff0000',label='required') xs=mreq[0]; ys=nreq[0]+0.5; ss='{0:.3f}'.format(xs) plt.text(xs,ys,ss,va='bottom',ha='center',fontsize=fsz,rotation=90) plt.legend(loc='upper right') plt.title(strt,loc='left',fontsize=fsz-2) plt.tight_layout() plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1) plt.show() def main(): start=time.time() a32=32*32*3.14/4 a25=25*25*3.14/4 a20=20*20*3.14/4 fcu=35.0 # design strength of concrete for iii in [1]: if iii==1: bb=1000 hh=1000 As=np.array([a20*5,a20*5]) ys=np.array([100,hh-100]) strt='Outlet channel B1000xH1000,C35-T20-200(D)' fnameF='fig_rc_channel.png' mh=int(hh/10) nreq=np.array([0]) # kN mreq=np.array([257.561*1.4]) # kN-m nreq=nreq*1000/bb/hh mreq=mreq*1000*1000/bb/hh/hh fnp,fmp,fm0=im_rebar.calmain(mh,fcu,bb,hh,As,ys) drawfig(fnp,fmp,fm0,nreq,mreq,strt,fnameF) print(time.time()-start) #============== # Execution #============== if __name__ == '__main__': main()
以 上
