２軸曲げを受けるRC部材の終局強度 (2022.02.02)

軸力およびモーメントの計算

軸力 $N$ 、x軸回りの曲げモーメント $M_x$ 、y軸回りの曲げモーメント $M_y$ を受けるRC部材を考える。座標軸は、下図に示すように、断面の図心を原点とし、強軸をx軸、弱軸をy軸と定義する。

f:id:damyarou:20220202173319j:plain

上図を参照して、積荷点の座標 $(e_x, e_y)$ は以下の通り計算できる。

$\begin{equation} e_x=\cfrac{M_y}{N} \qquad e_y=\cfrac{M_x}{N} \end{equation}$

ここで、軸力・モーメントと部材寸法の単位はあわせる必要があることに注意する。例えば、軸力の単位が、kN、モーメントの単位がkN-m、部材寸法の単位が mm であれば、上式で計算した $e_x$ および $e_y$ は1000倍する必要がある。

耐荷力計算においては、モーメントを計算する軸が載荷点上を通るよう回転させ、この回転させた軸に対し、モーメントの計算を行う。軸の回転角 $\theta$ は、以下により計算できる。

$\begin{equation} \theta=\tan^{-1}\left(\cfrac{e_x}{e_y}\right) \end{equation}$

ここで、コンクリート断面を、もとの座標上で予め小さな領域に分割し、すべての鉄筋についても、その位置と断面積を元の x-y 座標上で指定し、それら一つ一つを $(x_i, y_i)$ とする。もとの x-y 座標系においての座標 $(x_i, y_i)$ は、回転させたあとの x' 軸からの距離 $r$ によりモーメントを計算することになり、 $r$ の値は、以下により計算できる。

$\begin{equation} r_i=y_i\cdot \cos\theta+x_i\cdot\sin\theta \end{equation}$

以上より、新しいx' 軸回りの軸力およびモーメントは、以下のように求めることができる。

$\begin{gather} N = \int \sigma_i dA = \sum \left(\sigma_i \cdot\Delta_i \right) \\ M = \int \sigma_i \cdot r_i dA = \sum \left(\sigma_i \cdot r_i \cdot\Delta_i \right) \\ \end{gather}$

ここに、 $\Delta_i$ は、座標 $(x_i, y_i)$ に代表される微小領域の面積を示す。

材料の応力-ひずみ曲線

f:id:damyarou:20220202173404j:plain

コンクリートの応力-ひずみ曲線は、コンクリート圧縮ひずみが0.002までの範囲で、以下の式で表す。

$\begin{equation} \sigma_c=k_1\cdot fc\times\cfrac{\epsilon_c}{0.002}\times\left(2-\cfrac{\epsilon_c}{0.002}\right) \end{equation}$

コンクリート圧縮ひずみが0.002を超えてからは一定強度 $k_1\cdot f_c$ を保ち、終局ひずみは $\epsilon_{cu}=0.0035$ とする。

鉄筋の応力-ひずみ曲線は、圧縮降伏点から引張降伏点までは弾性係数 $E_s$ をもつ直線とし、降伏後は降伏点強度 $f_y$ を保つものとする。鉄筋の伸びはコンクリートの終局ひずみ 0.0035 に比べ十分に大きいため、伸びの最大・最小は規定しない。

なお、プログラムでの計算上、ひずみ、コンクリート応力・鉄筋応力共に、圧縮を正とする。

計算例

以下の諸元の断面・荷重に対するM-N線図を作成する。

Item	Symbol	Value
Concrete strength	fc	21 MPa
Yield strength of rebar	fy	390 MPa
Axial force	N	5000 kN
Moment around x-axis	Mx	4000 kN-m
Moment around y-axis	My	3000 kN-m
Width of section	b	1000 mm
Height of section	h	2000 mm
Number of division of width	mb	100
Number of division of height	mh	200
Number of reinforcements	ms	20 nos

プログラムの中で値を指定している変数は以下の通り。

k1=0.85    # reducing coefficient of concrete strength
ecu=0.0035 # ultimate strain of concrete
Es=2.0e5   # (MPa) Elastic modulus of steel bar

入力データ・フォーマットは以下の通り。

fc fy N Mx My b h mb mh ms
xs[i] ys[i] area[i]
    xs   : x-coordinate of rebar
    ys   : y-coordinate of rebar
    area : section area of rebar
..... i = 1 ~ ms .....

以下に入力データ事例を示す。

21 390 5000 4000 3000 1000 2000 100 200 20
-400  900 794
-200  900 794
   0  900 794
 200  900 794
 400  900 794
-400 -900 794
-200 -900 794
   0 -900 794
 200 -900 794
 400 -900 794
-400  700 794
-200  700 794
   0  700 794
 200  700 794
 400  700 794
-400 -700 794
-200 -700 794
   0 -700 794
 200 -700 794
 400 -700 794

以下に計算結果図を示す。

f:id:damyarou:20220202173449j:plain

計算プログラム

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as tick


def inpdata(fnameR):
    f=open(fnameR,'r')
    text=f.readline()
    text=text.strip()
    text=text.split()
    fc =float(text[0]) # concrete strength
    fy =float(text[1]) # concrete strength
    FN =float(text[2]) # axial force
    MX =float(text[3]) # bending moment around x-axis
    MY =float(text[4]) # bending moment around y-axis
    bb =float(text[5]) # width of rectangle section
    hh =float(text[6]) # height of rectangle section
    mb =int(text[7]) # number of division of width
    mh =int(text[8]) # number of division of height
    ms =int(text[9]) # number of reinforcement bar
    xs  =np.zeros(ms,dtype=np.float64) # x-coordinate of rebar
    ys  =np.zeros(ms,dtype=np.float64) # y-coordinate of rebar
    area=np.zeros(ms,dtype=np.float64) # section area of rebar
    for i in range(0,ms):
        text=f.readline()
        text=text.strip()
        text=text.split()
        xs[i]=float(text[0])
        ys[i]=float(text[1])
        area[i]=float(text[2])
    f.close()
    return fc,fy,FN,MX,MY,bb,hh,mb,mh,ms,xs,ys,area


def drawfig(theta0,an_0,mo_0,theta1,an_1,mo_1,theta2,an_2,mo_2,fax,mox):
    fsz=10
    xmin=0 ; xmax=5; dx=1
    ymin=-10; ymax=30; dy=5
    iw,ih=4,4
    plt.figure(figsize=(iw,ih),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='monospace'
    tstr='M-N diagram'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Moment M/b/h**2 (MPa)')
    plt.ylabel('Axial force N/b/h (MPa)')
    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)
    s0=r'$\theta$={0:4.0f} deg'.format(np.degrees(theta0))
    s1=r'$\theta$={0:4.0f} deg'.format(np.degrees(theta1))
    s2=r'$\theta$={0:4.0f} deg'.format(np.degrees(theta2))
    sp=r'Loaded point'
    plt.plot(mo_0,an_0,'-',color='#000080',lw=1,label=s0)
    plt.plot(mo_1,an_1,'-.',color='#333333',lw=1,label=s1)
    plt.plot(mo_2,an_2,'--',color='#333333',lw=1,label=s2)
    plt.plot([mox],[fax],'o',ms=4,color='#000080',label=sp)
    plt.legend(loc='upper right')
    plt.tight_layout()
    fnameF='fig_2mn.jpg'
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
    #plt.show()


def calc(xc,yc,da,xs,ys,area,theta,mb,mh,ms,k1,fc,ecu,Es,fy):
    def sig_con(k1,fc,ec,ecu):
        # concrete stress
        sigc=k1*fc*ec/0.002*(2-ec/0.002)
        if 0.002<=ec: sigc=k1*fc
        if ec<0: sigc=0.0
        if ecu*1.0001<ec: sigc=0.0
        return sigc
    def sig_bar(Es,fy,eb):
        # steel bar stress
        sigs=Es*eb    
        if fy<sigs: sigs=fy
        if sigs<-fy: sigs=-fy
        return sigs
    def cal_mn(e1,e2,k1,fc,ecu,rc,da,Es,fy,rs,area,l_an,l_mo):
        # calculation of axial force and moment
        rmin=np.min(rc)
        rmax=np.max(rc)
        a=(e2-e1)/(rmax-rmin)
        b=(rmax*e1-rmin*e2)/(rmax-rmin)
        e_con=a*rc+b
        e_bar=a*rs+b
        sig_c=np.zeros(len(e_con),dtype=np.float64)
        sig_b=np.zeros(len(e_bar),dtype=np.float64)
        sig_r=np.zeros(len(e_bar),dtype=np.float64)
        for i,ec in enumerate(e_con):
            sig_c[i]=sig_con(k1,fc,ec,ecu)
        for i,eb in enumerate(e_bar): 
            sig_b[i]=sig_bar(Es,fy,eb)
            sig_r[i]=sig_con(k1,fc,eb,ecu)
        an=np.sum(sig_c*da)+np.sum(sig_b*area)-np.sum(sig_r*area)
        mo=np.sum(sig_c*da*rc)+np.sum(sig_b*area*rs)-np.sum(sig_r*area*rs)
        l_an=l_an+[an]
        l_mo=l_mo+[mo]
        return l_an, l_mo
    def mo_nega(a_an,a_mo):
        # treatment of negative moment
        if 0<=a_mo[0]:
            a_an=np.append(a_an[0]-1e-6,a_an)
            a_mo=np.append(0,a_mo)
        else:
            for i in range(0,len(a_mo)-1):
                if a_mo[i]<=0 and 0<a_mo[i+1]:
                    x1,y1=a_mo[i],a_an[i]
                    x2,y2=a_mo[i+1],a_an[i+1]
                    b=(x2*y1-x1*y2)/(x2-x1)
                    break
            a_an=np.append(b,a_an)
            a_mo=np.append(0,a_mo)
        if 0<=a_mo[-1]:
            a_an=np.append(a_an[-1]+1e-6,a_an)
            a_mo=np.append(0,a_mo)
        else:
            for i in reversed(range(1,len(a_mo))):
                if a_mo[i]<=0 and 0<a_mo[i-1]:
                    x1,y1=a_mo[i],a_an[i]
                    x2,y2=a_mo[i-1],a_an[i-1]
                    b=(x2*y1-x1*y2)/(x2-x1)
                    break
            a_an=np.append(a_an,b)
            a_mo=np.append(a_mo,0)
        ii=np.argsort(a_an)
        _a_an=a_an[ii]
        _a_mo=a_mo[ii]
        jj=np.where(-1e-6<=_a_mo)
        a_an=_a_an[jj]
        a_mo=_a_mo[jj]
        return a_an, a_mo

    # main calculation
    rc=yc*np.cos(theta)+xc*np.sin(theta)
    rs=ys*np.cos(theta)+xs*np.sin(theta)
    
    esy=fy/Es
    nn=20
    num=3
    l_an=[] # axial force
    l_mo=[] # bending moment
    eps=np.linspace(-num*esy,ecu,nn)
    e1=-num*esy
    for e2 in eps:
        l_an,l_mo=cal_mn(e1,e2,k1,fc,ecu,rc,da,Es,fy,rs,area,l_an,l_mo)
    e2=ecu
    for e1 in eps:
        l_an,l_mo=cal_mn(e1,e2,k1,fc,ecu,rc,da,Es,fy,rs,area,l_an,l_mo)
    a_an=np.array(l_an)
    a_mo=np.array(l_mo)
    a_an,a_mo=mo_nega(a_an,a_mo) # treatment of negative moment
    return a_an, a_mo


def main():   
    fnameR='inp_rebar.txt'
    fc,fy,FN,MX,MY,bb,hh,mb,mh,ms,xs,ys,area=inpdata(fnameR)
    xc=np.zeros(mb*mh,dtype=np.float64)
    yc=np.zeros(mb*mh,dtype=np.float64)

    k1=0.85 # reducing coefficient of concrete strength
    ecu=0.0035 # ultimate strain of concrete
    Es=2.0e5 # (MPa) Elastic modulus of steel bar

    #eccentricity and angle of rotation
    if 1e-6 < np.abs(FN):
        ecy=MX/FN
        ecx=MY/FN
        theta0=np.arctan2(ecx,ecy)
    else:
        if np.abs(MX)<1e-6: theta0=np.pi/2
        if np.abs(MY)<1e-6: theta0=0.0
    # coordinates of concrete element and rebars
    dx=bb/float(mb)
    dy=hh/float(mh)
    da=dx*dy
    k=-1
    for i in range(0,mh):
        y=0.5*hh-0.5*dy-dy*float(i)
        for j in range(0,mb):
            x=-0.5*bb+0.5*dx+dx*float(j)
            k=k+1
            xc[k]=x
            yc[k]=y

    nmin=np.sum(fy*area)/bb/hh
    nmax=np.sum(fy*area)/bb/hh+k1*fc
    print('Max.tension, Max.conmresion')    
    print('{0:10.3f}{1:10.3f}'.format(nmin,nmax))
    theta1, theta2 = 0, np.pi/2
    for kkk,theta in enumerate([theta0, theta1, theta2]):
        a_an,a_mo=calc(xc,yc,da,xs,ys,area,theta,mb,mh,ms,k1,fc,ecu,Es,fy)
        a_an=a_an/bb/hh
        a_mo=a_mo/bb/hh**2
        print('****',np.degrees(theta))
        for i in range(len(a_an)):
            print('{0:10.3f}{1:10.3f}'.format(a_an[i],a_mo[i]))    
        if kkk==0:
            an_0=a_an
            mo_0=a_mo
        if kkk==1:
            an_1=a_an
            mo_1=a_mo
        if kkk==2:
            an_2=a_an
            mo_2=a_mo
    fax=FN*1e3/bb/hh
    mox=np.sqrt(MX**2+MY**2)*1e6/bb/hh**2
    drawfig(theta0,an_0,mo_0,theta1,an_1,mo_1,theta2,an_2,mo_2,fax,mox)
    

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

以上