はじめに

マレーシア時代に使っていた、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

$\begin{equation} \sigma_c= \begin{cases} \quad 0 & (\epsilon <0) \\ \quad \sigma_{cy}-\cfrac{\sigma_{cy}}{\epsilon_{cy}^2}(\epsilon - \epsilon_{cy})^2 & (0\leqq \epsilon \leqq \epsilon_{cy}) \\ \quad \sigma_{cy} & (\epsilon_{cy} < \epsilon \leqq \epsilon_{cu}) \\ \quad 0 & (\epsilon_{cu} < \epsilon) \end{cases} \end{equation}$

$\begin{gather} \sigma_{cy}=0.67\cdot \cfrac{f_{cu}}{\gamma_m} \qquad \epsilon_{cy}=2.4\times 10^{-4} \sqrt{\cfrac{f_{cu}}{\gamma_m}} \qquad \gamma_m=1.5 \qquad \epsilon_{cu}=0.0035 \end{gather}$

where, $f_{cu}$ means a compressive strength of concrete based on cubic specimen following BS. For example, C35 in BS means $f_{cu}=35$ (N/mm $^2$ ).

Reinforcement

$\begin{equation} \sigma_s= \begin{cases} \quad -\sigma_{sy} & (\epsilon < -\epsilon_{sy}) \\ \quad E_s \cdot \epsilon & (-\epsilon_{sy} \leqq \epsilon \leqq \epsilon_{sy}) \\ \quad \sigma_{sy} & (\epsilon_{sy} < \epsilon) \\ \end{cases} \end{equation}$

$\begin{gather} \sigma_{sy}=\cfrac{f_y}{\gamma_m} \qquad \epsilon_{sy}=\cfrac{f_y}{\gamma_m E_s} \qquad \gamma_m=1.05 \qquad E_s=2.0\times 10^5 \text{(N/mm$^2$)} \\ \end{gather}$

where, $f_y$ means a yield strength of rebar and usually $f_y=460$ (N/mm $^2$ ) 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.

$\begin{align} &\text{Axial force} &\qquad &N=\int_{0}^{h} \sigma \cdot dA=\sum_{i=1}^{m} \sigma_i \cdot b \cdot dy_i \\ &\text{Bending moment} &\qquad &M=\int_{0}^{h} \sigma \cdot (y-h/2) \cdot dA=\sum_{i=1}^{m} \sigma_i \cdot b \cdot dy_i \cdot (y_i-h/2) \end{align}$

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 $\epsilon_{cu}$
Set ultimate limit tensile strain of rebar as $-3 \epsilon_{sy}$ (conveniently assumed value)
Step-1: Fix the upper edge strain $\epsilon_1=\epsilon_{cu}$ . Vary the lower edge strain $\epsilon_2$ from $\epsilon_{cu}$ to $-3 \epsilon_{sy}$ and calculate the axial force and bending moment.
Step-2: Fix the lower edge strain $\epsilon_2=-3 \epsilon_{sy}$ . Vary the upper edge strain $\epsilon_1$ from $\epsilon_{cu}$ to $-3 \epsilon_{sy}$ and calculate the axial force and bending moment.

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).

f:id:damyarou:20190602175350j:plain

プログラム

自作関数のインポート

ここでは、曲げ耐荷力を計算する自作関数「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()