damyarou

python, GMT などのプログラム

Ultimate limit strength of rectangular RC section (2022.03.18)

  • A method of calculation of ultimate limit strength of rectangular RC member using rectangular stress block is described.
  • In this document, the positive sign of stress and strain for concretee and compressive reinforcement mean compression, and the positive sign of stress and asrein for tensile reinforcement means tension.
  • The calculation program is basically for symmetric rebar arrangement.

General formulas

A rectangular concrete cross section with double reinforcement is considered.


\begin{align}
&b & & \text{width of section} \\
&h & & \text{height of section} \\
\end{align}

From Bernoulli-Euler hypothsis, the strains of tensile and compressive re-bars can be expressed as followings.


\begin{align}
& \epsilon_t=\cfrac{d-x}{x}\cdot \epsilon_{cu} & & +: \text{tensile}\\
& \epsilon_c=\cfrac{x-d'}{x}\cdot \epsilon_{cu} & & +: \text{compressive}
\end{align}

\begin{align}
&\epsilon_t & & \text{strain of tensile re-bar} \\
&\epsilon_c & & \text{strain of compressive re-bar} \\
&\epsilon_{cu} & & \text{ultimate strain of concrete} \\
&x & & \text{location of neutral axis from compressive edge}\\
&d & & \text{effective depth of tensile reinforcement} \\
&d' & & \text{cover of compressive reinforcement}
\end{align}

Using above, the stress of re-bar can be calculated as follows.


\begin{align}
&-f_y/\gamma_m \leq \sigma_s =E_s\cdot \epsilon_t \leq f_y/\gamma_m & & +:\text{tensille}\\
&-f_y/\gamma_m \leq \sigma_s'=E_s\cdot \epsilon_c \leq f_y/\gamma_m & & +:\text{compressive}
\end{align}

\begin{align}
&\sigma_s & & \text{sstress of tensile re-bar} \\
&\sigma_s' & & \text{stress of compressive re-bar} \\
&E_s & & \text{elastic modulus of reinforcement} \\
&f_y & & \text{yield strength of reinforcement} \\
&\gamma_m & & \text{partial safety factor for material}
\end{align}

The height of compressive stress block of cncrete of a is defined as follow using the location of neutral axis of x.


\begin{equation}
a=\beta_1\cdot x \leq h
\end{equation}

\begin{align}
&a & & \text{height of concrete stress block}\\
&\beta_1 & & \text{coefficient for height of concrete stress block}
\end{align}

The axsial force of N and the bending moment of M can be expressed as followings.


\begin{equation}
N=\beta_2\cdot f_{cu}\cdot b\cdot a - A_s\cdot \sigma_s + A_s'\cdot \sigma_s'
\end{equation}

\begin{equation}
M=N\cdot e=\beta_2\cdot f_{cu}\cdot b\cdot d^2 \left(\cfrac{a}{d}\right)
\left\{1-\cfrac{1}{2}\left(\cfrac{a}{d}\right)\right\} + A_s'\cdot \sigma_s'\cdot (d-d')
-N\cdot c
\end{equation}

\begin{align}
&N & & \text{axial force}\\
&M & & \text{bending moment}\\
&e & & \text{distance from centroid to loading point of $N$} \\
&c & & \text{distance from centroid to tensile re-bar}\\
&f_{cu} & & \text{concrete strength} \\
&\beta_2 & & \text{coefficient for width of concrete stress block} \\
&A_s & & \text{section area of tensile reinforcement} \\
&A_s' & & \text{section area of compressive reinforcement} 
\end{align}

f:id:damyarou:20220318143524j:plain

Location of neutral axis

Balanced failure

In the case of balanced failure, following conditions are applied.

  • the concrete strain reaches the ultimate strain of \epsilon_{cu}.
  • the same time, the tensile re-bar stress reaches the yield stress.

From above, the location of neutral axis of x can be obtained as follow.


\begin{gather} 
\epsilon_t=\cfrac{f_y}{\gamma_m\cdot E_s}=\cfrac{d-x}{x}\cdot \epsilon_{cu} \\
x=\cfrac{d}{1+\cfrac{f_y}{\gamma_m\cdot E_s\cdot \epsilon_{cu}}}
\end{gather}

When above x is re-defined as x_0, the axial force of N_0 for balanced failure can be calculated using x_0.

Pure bending

In case of pure bending, following conditions are applied.

  • the axial force is zero (N=0).
  • the concrete strain has reached the ultimate strain of \epsilon_{cu}.
  • if 0 \geq N_0, the failure mode is tensile failure. This means the tensile re-bar stress has reached the yield stress.
  • if N_0 \lt 0, the failure mode is compressive failure. This means the compressive re-bar stress has reached the yield stress.

From above, the location of neutral axis of x can be obtained as a root of following quadratic eqiation which is derived from calculation formula of N.

For tensile failure (N=0)


\begin{equation}
\beta_2\cdot f_{cu}\cdot b\cdot \beta_1\cdot x^2
+\left(A_s'\cdot E_s\cdot \epsilon_{cu}-A_s\cdot \cfrac{f_y}{\gamma_m}-N\right)\cdot x
-A_s'\cdot E_s\cdot d'\cdot \epsilon_{cu}=0  
\end{equation}

\begin{align}
&c_1=\beta_2\cdot f_{cu}\cdot b\cdot \beta_1 \\
&c_2=A_s'\cdot E_s\cdot \epsilon_{cu}-A_s\cdot \cfrac{f_y}{\gamma_m}-N \\
&c_3=-A_s'\cdot E_s\cdot d'\cdot \epsilon_{cu}
\end{align}

For compressive failure (N=0)


\begin{equation}
\beta_2\cdot f_{cu}\cdot b\cdot \beta_1\cdot x^2
+\left(A_s\cdot E_s\cdot \epsilon_{cu}+A_s'\cdot \cfrac{f_y}{\gamma_m}-N\right)\cdot x
-A_s\cdot E_s\cdot d\cdot \epsilon_{cu}=0  
\end{equation}

\begin{align}
&c_1=\beta_2\cdot f_{cu}\cdot b\cdot \beta_1 \\
&c_2=A_s\cdot E_s\cdot \epsilon_{cu}+A_s'\cdot \cfrac{f_y}{\gamma_m}-N \\
&c_3=-A_s\cdot E_s\cdot d\cdot \epsilon_{cu}
\end{align}

Root of quadratic equation


\begin{equation}
x=\cfrac{-c_2+\sqrt{c_2{}^2-4\cdot c_1\cdot c_3}}{2\cdot c_1}
\end{equation}

Calculation of ultimate strength

Constants based on BS 8110-1


\begin{align}
&\gamma_m=1.5 & & \text{for concrete} \\
&\gamma_m=1.05 & & \text{for reinforcement} \\
&\beta_1=0.9 & & \text{coefficient for height of concrete stress block} \\
&\beta_2=0.67 / \gamma_m & & =0.447 \quad (\gamma_m=1.5) \\
&E_s=200,000 & & \text{in MPa} \\
&\epsilon_{cu}=0.0035 & & \text{ultimate strain of concrete} \\
\end{align}

f:id:damyarou:20220318144300j:plain

Calculation steps

Step Description
1 Inpput parameters
3 Calculate neutral axis for balanced failure and pure bending
2 Specify neutral axis
3 Calculate axial force and moment
4 Draw M-N diagram

Program

import numpy as np
import matplotlib.pyplot as plt


def drawfig(a_ax,a_mo,tstr):
    fsz=10
    xmin=0 ; xmax=5; dx=1
    ymin=-10; ymax=30; dy=5
    iw,ih=5,5
    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('Bending moment M/b/h**2 (MPa)')
    plt.ylabel('Axialforce 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)
    plt.plot(a_mo,a_ax,'o-',ms=4,clip_on=False)
    plt.tight_layout()
    fnameF='fig_mn_rebar1.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    #plt.show()


def pure_bend(fcu,beta1,beta2,ecu,fy,rms,Es,ast,asc,dd,dc,bb,ax0):
    ax=0 # axial force
    if 0<ax0: # tensile failure
        c1=beta2*fcu*bb*beta1
        c2=asc*Es*ecu-ast*fy/rms-ax
        c3=-asc*Es*dc*ecu
    else: # compressive failure
        c1=beta2*fcu*bb*beta1
        c2=ast*Es*ecu+asc*fy/rms-ax
        c3=-ast*Es*dd*ecu
    x=(-c2+np.sqrt(c2**2-4*c1*c3))/2/c1 # neutral axis
    return x


def cal_mn(bb,hh,fcu,x,dd,dc,cc,beta1,beta2,ecu,Es,fy,rms,ast,asc):
    et=(dd-x)/x*ecu
    ec=(x-dc)/x*ecu
    sigt=Es*et
    sigc=Es*ec
    if fy/rms < sigt: sigt=fy/rms
    if sigt < -fy/rms: sigt=-fy/rms
    if fy/rms < sigc: sigc=fy/rms      
    if sigc < -fy/rms: sigc=-fy/rms
    aa=beta1*x
    if hh<aa: aa=hh
    ax=beta2*fcu*bb*aa-ast*sigt+asc*sigc
    mo=beta2*fcu*bb*dd**2*(aa/dd)*(1-1/2*(aa/dd))+asc*sigc*(dd-dc)-ax*cc
    return ax,mo


def calc(fcu,fy,bb,hh,ast,asc,dt,dc):
    # constants
    rmc=1.5 # partial safety factor for concrete
    rms=1.05 # partial safety factor for reinforcement
    beta1=0.9 # coefficient for height of concrete stress block
    beta2=0.67/rmc # coefficient for width of concrete stress block 
    Es=200000 # elastic modulus of reinforcement
    ecu=0.0035 # ultimate strain of concrete
    dd=hh-dt
    cc=hh/2-dt
    # Balanced failure
    x0=dd/(1+fy/rms/Es/ecu)
    ax0,mo0=cal_mn(bb,hh,fcu,x0,dd,dc,cc,beta1,beta2,ecu,Es,fy,rms,ast,asc)
    # Pure bending
    xb=pure_bend(fcu,beta1,beta2,ecu,fy,rms,Es,ast,asc,dd,dc,bb,ax0)
    # Max and Min axial force
    ax_max=beta2*fcu*bb*hh+ast*fy/rms+asc*fy/rms
    ax_min=-(ast*fy/rms+asc*fy/rms)
    # Definition of neutral axis for calculation
    xxx=np.linspace(1,1.2*hh,20)
    xxx=np.concatenate([xxx,[x0,xb]])
    xxx=np.sort(xxx)
    # calculation
    a_ax=np.zeros(len(xxx),dtype=np.float64)
    a_mo=np.zeros(len(xxx),dtype=np.float64)
    for i,x in enumerate(xxx):
        ax,mo=cal_mn(bb,hh,fcu,x,dd,dc,cc,beta1,beta2,ecu,Es,fy,rms,ast,asc)
        a_ax[i]=ax
        a_mo[i]=mo
    a_ax=a_ax/bb/hh
    a_mo=a_mo/bb/hh**2
    a_ax=np.concatenate([a_ax,[ax_max/bb/hh,ax_min/bb/hh]])
    a_mo=np.concatenate([a_mo,[0,0]])
    ii=np.argsort(a_ax)
    a_ax=a_ax[ii]
    a_mo=a_mo[ii]
    return a_ax,a_mo


def main():
    # parameteers
    fcu=30 # concrete strength
    fy=460 # yield strength of reinforcement
    bb=2000 # width of member
    hh=2000 # height of member
    ast=np.pi*32**2/4*20 # section area of tensile reinforcement
    asc=np.pi*32**2/4*20 # section area of compressive reinforcement
    dt=200 # cover of tensile reinforcement
    dc=200 # cover of compressive reinforcement
    a_ax,a_mo=calc(fcu,fy,bb,hh,ast,asc,dt,dc)    
    tstr='2T32ctc200x2layer,double,bxh=2000x2000'
    drawfig(a_ax,a_mo,tstr)
    
    
#---------------
# Execute
#---------------
if __name__ == '__main__': main()
Example of analysis result

f:id:damyarou:20220318143559j:plain

Thank you.

MATECH 100W 充電器購入(2022.03.10)

2022年3月10日、昨日注文した MATECH 100W 充電器(Sonicharge 100W Pro)が届いた。¥6,980也。100W 充電ケーブルも付属している。ちなみに現在保有している CIO 100W 充電器(2021年7月12購入)は、¥6,578也。付属品はなし。

現在 CIO の充電器は2台同時接続時は 45W+45W で給電されるため、iPad Pro 12.9 と iPad mini 6 の充電に使っているのだが、今後出張を考えると MacBook Pro 14 用にもう1台 100W 充電器がほしかったため、今回 MATECH のものを選定した。これは2台以上接続時でも、1 port 65Wの出力が保証されるため、MacBook Pro 14 の充電に最適と考えたため選定したものである。各ポートからの出力は以下の通りである。

MATECH 100W 充電器出力

Port C1 C2 A C1+C2 C1+A C2+A C1+C2+A
USB-C1 100W --- --- 65W 65W --- 65W
USB-C2 --- 65W --- 30W --- (18W) (18W)
USB-A --- --- 30W --- 30W (18W) (18W)

CIO 100W 充電器出力

Port C1 C2 C1+C2
USB-C1 100W --- 45W
USB-C2 --- 100W 45W

MATECH 充電器写真

f:id:damyarou:20220310122735j:plain f:id:damyarou:20220310122754j:plain f:id:damyarou:20220310122813j:plain

CIO 充電器との比較写真

f:id:damyarou:20220310122831j:plain f:id:damyarou:20220310122847j:plain f:id:damyarou:20220310122902j:plain

以 上

iPad mini 6の充電テスト(2022.03.03)

2022年3月2日、先日購入した Anker のモバイルバッテリーで iPad mini 6 の充電を行ってみた。iPhone は Lightning 端子なのに対し、iPad mini 6 は USB-C 端子なのでケーブルはモバイルバッテリー付属の USB-C<=>USB-C ケーブルを使用。バッテリー残量 13% からはじめての結果は以下の通り。1時間半の充電で92%まで回復。満タンでほぼ連続使用しても7〜8時間は持つので1時間充電すれば約60%回復するので3〜4時間は使えそうな感じだ。なお1時間半充電時点で、モバイルバッテリー側の青いランプは4つのうち2つが点灯。80%充電するのに50〜75%使った勘定になる(詳細は不明)。モバイルバッテリー側の残量が正確ではないが、概算、4000mAh(5100 x 80% )を満たすのにバッテリー側は 5000〜7500mAh 使っていることになるので、効率は 80% から 53% の間ということになる。YouTube の情報によれば効率は一般的には 60〜70% らしい。 Amazon の製品紹介によれば、モバイルバッテリーは 18W 出力の PD 対応充電器および USB-C & USB-C ケーブルを使用した場合、所要充電時間は 2.8 時間となっている。充電には iPad 付属の 20W 充電器を使っているので、約3時間でフル充電できるということになる?(未検証)

時間 (分) 0 30 60 90
充電 (%) 13 46 74 92

(参考)ネットで調べたバッテリー容量

Device capacity
iPhone SE 1,812 mAh
iPad mini 6 5,124 mAh
iPad Pro 12.9 7,600 mAh

f:id:damyarou:20220303070132j:plain

以 上

10cm充電ケーブル購入(2022.03.01)

今日は会社の健康診断で千駄ヶ谷に行きそのあと出社。午後半休をとって12:00に退社。12:40の新幹線で帰路につく。帰りがけにヤマダ電機で充電用ケーブル(USB-C<=>Lightning)を購入。¥1,233也。長さは10cm。ELECOM製で60Wまで送れるとのこと。iPhone SEと充電器を並べるとケーブル長は10cmで丁度いい感じである。

以 上

モバイルバッテリー他購入(2022.02.28)

2022年2月27日、前日にAmazonで注文した品物が届いた。 ロシアのウクライナ侵攻もあり、色々品不足の時代が来るかもしれない。主要なもの(MacBook pro 14"、iPad Pro 12.9", iPad mini 6)は、昨年更新してある。ということで、Amazonで3品を購入。Amazonポイントが384円分あったので活用。購入品は以下の通り。

Anker USB-C & USB-C Thunderbolt 3 ケーブル (0.7m ブラック)【100W出力 / 40Gbps / 高速データ転送 / 4K対応 / 5K対応】MacBook iPad Pro 他対応

MOSISO ラップトップ スリーブケース 耐衝撃 対応MacBook Pro 14 インチ 2021 2022 M1 Pro/M1 Max A2442、対応MacBook Air/Pro Retina、13-13.3 インチ Notebook 撥水 ポリエステル 保護ケース ポケット付き(ラベンダー グレー)

Anker PowerCore 10000 PD Redux 25W(モバイルバッテリー 10000mAh 大容量 )【USB Power Delivery対応/PPS規格対応/PowerIQ搭載/PSE技術基準適合】 iPad iPhone Galaxy Android スマートフォン タブレット 各種 その他機器対応

購入品 金額
Anker Thunderbolt 3 ケーブル ¥2,797
MOSISO ラップトップスリーブケース ¥1,899
Anker PowerCore 10000 PD Redux 25W ¥3192
小計 ¥7,888
Amazon ポイント -¥384
合計 ¥7,504

Anker モバイルバッテリーは CIO の 30W のものと迷ったが、結構価格差があったし、充電対象は基本的に iPhone SE なので、お安い Anker の 25W にした。モバイルバッテリーへの充電は、製品付属の USB-C<=>USB-C ケーブルで、iPad に付属してきた 20W 充電器で行う。iPhone SE への充電は、モバイルバッテリーと、iPhone SE 付属の Lightning<=>USB-C ケーブルを用いて行う。ここの操作が面倒ではあるが、しょうがない。充電テストを行った結果は以下の通り。スタートが 27% とちょっと高めであったが問題なく動いている。充電率が100%に近くなると速度が遅くなるのはしょうがない。なお、iPhone SE 付属の Lightning<=>USB-C ケーブルはモバイルバッテリーからの充電用としては長いので、短いものを購入したほうが良いかもしれない。

時間 (分) 0 30 60 90 120 150
充電 (%) 27 59 83 91 95 97

iPad Pro 12.9" を持ち運ぶためにケースを探していたのだが、YouYube で見つけて、MOSISO のラップトップケースを購入することに。色はラベンダーグレー。MacBook Pro 14" も問題なく入るが、MacBook 用には7年来使い続けているケースがあり、特に不具合もないのでそのまま使い続けることにしている。今回購入したケースは、メインの収納には iPad Pro 12.9" をいれるが、ポケットがついており、iPad mini 6 がちょうど収まりいい感じである。出張時はこれに、iPad Pro 12.9" と iPad mini 6 を入れて歩くことにするつもりである。

Ankerのサンダーボルト3ケーブルは、予備のために購入。100W充電・データ転送・映像出力に対応している。

f:id:damyarou:20220228084305j:plain

f:id:damyarou:20220228084332j:plain

f:id:damyarou:20220228084356j:plain

f:id:damyarou:20220228084414j:plain

f:id:damyarou:20220228084434j:plain

f:id:damyarou:20220228084451j:plain

以 上

2軸曲げを受ける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()

以上

ワイヤレスポータブルスピーカー購入(2022.01.23)

昨日、近所のビックカメラでワイヤレスポータブルスピーカー(SONY SRS-XB23)を購入。これには充電器は同梱されておらず、スピーカー側USB-C、充電器側USB-Aのケーブルのみが同梱されているので、USB-C充電器を使うため、変換アダプタ(USB-Aメス<=>USB-C:USB3-AFCMADBK)も同時に購入。 スピーカーはiPhoneBluetooth接続しどこにでも持ち歩ける。充電はiPad付属の20W充電器(USB-C)に接続し実施。

購入品 購入日 購入場所 購入金額
SRS-XB23 (SONY) 2022.01.22 ビックカメラ ¥11,200
USB3-AFCMADBK (ELECOM) 2022.01.22 ビックカメラ ¥1,290

f:id:damyarou:20220123100145j:plain

以 上