2019/6/16

モンテカルロ法の例題

RC梁の耐力分布算定
信頼性工学特論
レポート課題例

2019/6/16

概要• 純曲げを受けるRC梁の終局曲げ耐力分布を求める• どうやって求めるか？…モンテカルロ法を利用• 材料パラメータ（強度等）と部材パラメータ（寸法）を乱数とする• RC梁の上縁部がコンクリート圧壊歪𝜀 に達する時点を終局曲げ破壊とする• ファイバー法による断面耐力の計算コードを作成する

2019/6/16

問題設定• RC梁断面を右のように定める

400• コンクリート• 主鉄筋D19✕5本（5-D19）• 1本の断面積：1.986(cm2)＝198.6(mm2)• 5-D19の断面積：1435(mm2)500• 材料パラメータ全てと主鉄筋の断面積が変動
60

300
22100

2019/6/16

材料パラメータと応力−歪曲線
𝜎=𝑓 +コンクリート𝑓 𝐸𝜀𝜀
1 − 𝑒𝑥𝑝 −𝜉−𝜉 𝐸𝜀𝜀

引張抵抗ゼロ𝜎=𝑓𝑓𝜎 = 𝐸𝜀𝜎 = −𝑓𝑓

𝜎=𝑓−𝜀塑性化歪：𝜀𝜀𝜀

𝜀+2𝜀−𝜀
圧壊歪：𝜀鋼材𝜀

圧縮強度：𝑓𝜀

ヤング率：𝐸

降伏強度：𝑓

降伏歪：𝜀 =硬化歪：𝜀

硬化係数：𝐸

硬化曲率：𝜉

2019/6/16

コンクリートの応力−歪曲線
strainがベクトルで与えられる場合を想定しIf文を避ける。
指定された範囲で1、それ以外で0を返す変数a、b、ｃを用意して利用する
Matlabでは 例えば
a=[1 2 3 4]に対してa>2は[0 0 1 1]を意味する

paramにコンクリートの材料パラメータを入れておく

引張域
𝜎=0a = strain <= 0; %-- tensile stress is zero

二次曲線

b = strain >= -ec;
s2 = param.fc * strain/ec.*(strain/ec+2);

塑性域

引張抵抗ゼロfunction [stress] = ss_concrete(strain, param)%SS_CONCRETE returns stress corresponding to the given strain
ec = param.ec;
eu = param.eu;

c = strain >= -eu & strain < -ec;
s3 = -param.fc;
stress = a.*(b.*s2 + c.*s3);
end

𝜎=𝑓

𝜎 = −𝑓𝑓
Matlabの掛け算は（*）は行列積（.*）は要素ごとの掛け算

塑性域−𝜀二次曲線−𝜀𝜀𝜀

𝜀+2𝜀2019/6/16

鋼材の応力−歪曲線
paramに鋼材の材料パラメータを入れておく

𝜎=𝑓 +function stress = ss_steal(strain, param)%SS_STEAL returns stress corresponding to the given strain
ey = param.fy/param.E;
弾性域
𝜎 = 𝐸𝜀a=abs(strain)<ey;
s1 = param.E * strain;

降伏棚
𝜎=𝑓b=abs(strain)<param.est & abs(strain)>=ey;s2 = param.fy*ones(size(strain));s2 = sign(strain).*s2;

硬化後strainがベクトルで与えられる場合を想定しIf文を避ける。
指定された範囲で1、それ以外で0を返す変数a、b、ｃを用意して利用する
sign関数は正負の符号を返すMatlab標準関数
strain./abs(strain)に相当するがゼロが分⺟に来た時エラーが⽣じるためsign関数を使うほうが便利

c=abs(strain)>=param.est;
alpha = param.fy / param.xi * param.Est / param.E;
s3 = param.fy + alpha * (1-exp(-param.xi*(abs(strain)/ey-param.est/ey)));s3 = sign(strain).*s3;
stress = a.*s1 + b.*s2 + c.*s3;
end𝑓 𝐸𝜀𝜀1 − 𝑒𝑥𝑝 −𝜉−𝜉 𝐸𝜀𝜀𝜎=𝑓𝑓

𝜎 = 𝐸𝜀𝜀𝜀2019/6/16

検証①｜適切に動作するか確認
コンクリートの材料パラメータをparamにまとめる
実行スクリプト%-- Setting Parameters
fc = 30; %-- Compression Strength
ec= 0.002; %-- Platcization Strain
eu=0.0035; %-- Ultmate Strain
param_concrete = init_concrete(fc, ec, eu);
strain = [-4000:100:1000]/1000000;
stress = ss_concrete(strain ,param_concrete);figure(3)plot(strain, stress, 'r-')ylim([-40 5])

正しく応力を返すかあらゆる歪を与え、グラフ化して確認する

function [param] = init_concrete(fc, ec, eu)%INIT_CONCRETE sets the parameters of steel material
param.fc = fc;
param.ec = ec;param.eu = eu;
end

歪を与えると応力を返す関数｜ss_concrete
function [stress] = ss_concrete(strain, param)%SS_CONCRETE returns stress corresponding to the given strain
ec = param.ec;
eu = param.eu;
a = strain <= 0; %-- tensile stress is zerob = strain >= -ec;
s2 = param.fc * strain/ec.*(strain/ec+2);
c = strain >= -eu & strain < -ec;
s3 = -param.fc;
stress = a.*(b.*s2 + c.*s3);end2019/6/16

検証②｜適切に動作するか確認
鋼材の材料パラメータをparamにまとめる
実行スクリプト%-- Setting Parameters
E = 200000; %-- Young Modulus 200GPa
fy= 300; %-- Yield Strength
xi=0.050; %-- Hardening Curvature
Est=4500; %-- Hardening Coefficient
est=0.020000; %-- Hardening Strain (20000u)
param_steel = init_steal(E,fy,xi,Est,est);
strain = [0:100:2000 4000:2000:40000 50000:10000:300000]'/1000000;
strain = [flipud(-strain);strain];stress = ss_steal(strain,param_steel);figure(1)plot(strain, stress,'r-');figure(2)plot(strain,stress,'r-');xlim([0 60000]/1000000)

正しく応力を返すかあらゆる歪を与え、グラフ化して確認する

function [param] = init_steal(E,fy,xi,Est,est)%INIT_STEAL sets the parameters of steel material
param.E = E;param.fy = fy;param.xi = xi;
param.Est = Est;
param.est = est;
end

歪を与えると応力を返す関数｜ss_steal
function stress = ss_steal(strain, param)%SS_STEAL returns stress corresponding to the given strain
ey = param.fy/param.E;
a=abs(strain)<ey;
s1 = param.E * strain;b=abs(strain)<param.est & abs(strain)>=ey;s2 = param.fy*ones(size(strain));s2 = sign(strain).*s2;
c=abs(strain)>=param.est;
alpha = param.fy / param.xi * param.Est / param.E;
s3 = param.fy + alpha * (1-exp(-param.xi*(abs(strain)/ey-param.est/ey)));s3 = sign(strain).*s3;
stress = a.*s1 + b.*s2 + c.*s3;
end2019/6/16

検証③｜適切に動作するか確認
sample_01を実行%-- Setting Parameters
fc = 30; %-- Compression Strength
ec= 0.002; %-- Platcization Strain
eu=0.0035; %-- Ultmate Strain
param_concrete = init_concrete(fc, ec, eu);>> sample_01

strain = [-4000:100:1000]/1000000;
stress = ss_concrete(strain ,param_concrete);figure(3)plot(strain, stress, 'r-')ylim([-40 5])%-- Setting Parameters
E = 200000; %-- Young Modulus 200GPa
fy= 300; %-- Yield Strength
xi=0.050; %-- Hardening Curvature
Est=4500; %-- Hardening Coefficient
est=0.020000; %-- Hardening Strain (20000u)
param_steel = init_steal(E,fy,xi,Est,est);

ss_concreteとss_stealが正しく動作することを確認

strain = [0:100:2000 4000:2000:40000 50000:10000:300000]'/1000000;
strain = [flipud(-strain);strain];stress = ss_steal(strain,param_steel);figure(1)plot(strain, stress,'r-');figure(2)plot(strain,stress,'r-');xlim([0 60000]/1000000)

実行結果2019/6/16

軸力の計算方法
平面保持仮定が成り立つと仮定−𝜀（歪分布は直線）
歪：𝜀

応力：𝜎

終局限界状態では𝑦 = 0で𝜀 = −𝜀
未知数は曲率：𝜙のみ
𝜀 = 𝜙𝑦 − 𝜀

平面保持仮定𝜀 = 𝜙𝑦 − 𝜀𝜙𝜎は𝜀から求められる
軸力𝑁はゼロなので
𝑁=𝜎 d𝐴 =𝜎𝐵d𝑦 + 𝜎 𝐴 = 0

を満たす曲率：𝜙を求めれば良い
位置：𝑦2019/6/16

ニュートン法による曲率の計算
𝜙 の値を最大値と最小値で適当に決め、それぞれの場合で軸力𝑁を求める

−𝜀𝑁

𝑁

歪：𝜀𝑁
𝜙𝜙𝜙𝜙
位置：𝑦𝜙
𝑁
くり返すとゼロに収束するはず

𝜙𝜙

2つの点を結ぶ直線を 𝑁 図の推定値として、ゼロを与える𝜙を求める
実際の𝑁図は直線ではないので、ここでの𝜙に対する𝑁値はゼロとはならない
𝑁がゼロより大きいときは最大値を、ゼロより小さいときは最小値を更新する

2019/6/16

軸力の計算

𝑁𝑁

itr_maxは繰り返し上限値
収束しなかった場合、永遠に終わらないことのないように打ち切る回数を決めておくparam=mk_model_02();
itr_max = 500;
err_max = 0.001; %-- it is assumed to be zeroニュートン法では、𝑁はゼロに近づくだけでゼロにはならない
ので、err_maxでゼロと見なす最大値を定めておく%-- fibery = param.member.y;
eu = param.concrete.eu;
hs = param.member.hs;%-- calculating stress distribution
ec = -eu + phi * y;er = -eu + phi * hs;
sc = ss_concrete(ec, param.concrete);sr = ss_steal(er, param.steel);
N = axial_force(y, sc, sr, param);%-- min/max curvaturephi1 = 0;
phi2 = 0.1;

𝜀分布を𝜙と𝜙それぞれで求めるの予め作っておいた
ss_concrete、ss_stealで応力分布を求める軸力を求める%-- strain distribution
ec1 = -eu + phi1 * y;er1 = -eu + phi1 * hs;
ec2 = -eu + phi2 * y;er2 = -eu + phi2 * hs;

ecはコンクリートの歪分布（ベクトル）erは鉄筋の歪（スカラー）

scはコンクリートの応力分布（ベクトル）%-- stress distribution
sc1 = ss_concrete(ec1, param.concrete);
sr1 = ss_steal(er1, param.steel);
sc2 = ss_concrete(ec2, param.concrete);sr2 = ss_steal(er2, param.steel);%-- axial force
N1=axial_force(y, sc1, sr1, param);
N2=axial_force(y, sc2, sr2, param);

axial_forceという関数を新たに作る%-- Newton method
for kk = 1:itr_max%-- calculating next phi
a=(N2-N1)/(phi2-phi1);
phi=phi2-N2/a;srは鉄筋の応力（スカラー）

𝜙𝜙𝜙𝑁𝜙と𝑁 、𝜙と𝑁から𝜙と𝑁を求める%-- checking axial force
if N<=err_max
break;
axial_force関数はここでも使う
elseif N<0phi1 = phi;
N1 = N;
軸力𝑁がゼロ（err_max以下）ならbreak（ループを出る）
elseif N>0
phi2 = phi;
𝜙か𝜙を𝜙で更新
N2 = N;
else%-- not usedbreak;
endend
軸力ゼロのときの𝜙（から求められる応力分布）で曲げモーメントを求める
y0=eu/phi;
M = bending_moment(y,y0,sc, sr, param);2019/6/16

軸力の計算
function [N] = axial_force(y, sc, sr, param)%AXIAL_FORCE returns axial force of the RC beam
Nc = 0; Ncはコンクリートが分担する軸力
for ii=1:(length(y)-1)
dy = y(ii+1)-y(ii);bw = param.member.bw;
Δ𝑦Nc = Nc + (sc(ii)+sc(ii+1)) * dy / 2 * bw;
end
応力分布を台形近似しながら積分する
As = param.member.As;
Ns = sr * As; Nsは鉄筋が分担する軸力
N = Nc + Ns;
end

𝜎 𝑦𝜎 𝑦2019/6/16

曲げモーメントの計算−𝜀歪：𝜀応力：𝜎𝑀==

平面保持仮定𝜀 = 𝜙𝑦 − 𝜀𝜙

位置：𝑦𝜎 𝑦−𝑦d𝐴

𝜎 𝑦 − 𝑦 𝐵d𝑦 + 𝜎 ℎ − 𝑦 𝐴

応力に腕の⻑さ（中立軸からの距離）を掛け算するとその応力のモーメントが求められる2019/6/16

曲げモーメントの計算
y0は中立軸位置

function [M] = bending_moment(y, y0, sc, sr, param)%BENDING_MOMENT returns bending moment of the RC beam
Mc = 0;for ii=1:(length(y)-1)
応力分布を台形近似しながら積分する
dy = y(ii+1)-y(ii);bw = param.member.bw;
Mc = Mc + (sc(ii)*(y(ii)-y0)+sc(ii+1)*(y(ii+1)-y0)) * dy / 2 * bw;
end
As = param.member.As;
hs = param.member.hs;
Ms = sr * As * (hs-y0);
M = Mc + Ms;
end

中立軸位置からの腕の⻑さを掛け算して応力モーメントを求めてから積分する

2019/6/16

検証④｜正しく動作するか確認sample_02を実行

param=mk_model_02();
itr_max = 500;
err_max = 0.001; %-- it is assumed to be zero%-- fibery = param.member.y;
eu = param.concrete.eu;
hs = param.member.hs;%-- Newton method
for kk = 1:itr_max%-- calculating next phi
a=(N2-N1)/(phi2-phi1);
phi=phi2-N2/a;%-- calculating stress distribution
ec = -eu + phi * y;er = -eu + phi * hs;
sc = ss_concrete(ec, param.concrete);sr = ss_steal(er, param.steel);
N = axial_force(y, sc, sr, param);%-- min/max curvaturephi1 = 0;
phi2 = 0.1;%-- strain distribution
ec1 = -eu + phi1 * y;er1 = -eu + phi1 * hs;
ec2 = -eu + phi2 * y;er2 = -eu + phi2 * hs;%-- checking axial force
if N<=err_max
break;
elseif N<0phi1 = phi;
N1 = N;elseif N>0
phi2 = phi;
N2 = N;
else%-- not usedbreak;
end%-- stress distribution
sc1 = ss_concrete(ec1, param.concrete);
sr1 = ss_steal(er1, param.steel);
sc2 = ss_concrete(ec2, param.concrete);sr2 = ss_steal(er2, param.steel);%-- axial force
N1=axial_force(y, sc1, sr1, param);
N2=axial_force(y, sc2, sr2, param);

さらに、>> fprintf([‘phi=‘ num2str(phi,’%f’) ‘ Mu=‘ num2str(M,’%f’) ‘¥n’])を実行すれば収束後のphiとMの値が確認できる

end
y0=eu/phi;
M = bending_moment(y,y0,sc, sr, param);2019/6/16

曲げモーメントの繰り返し計算の準備
function [M, phi] = capacity_calc(param)%CAPACITY_CALC returns capacity of RC cross section from given parameters
itr_max = 500;
err_max = 0.001; %-- it is assumed to be zero%-- fibery = param.member.y;
eu = param.concrete.eu;
hs = param.member.hs;%-- min/max curvaturephi1 = 0;
phi2 = 0.1;%-- Newton method
for kk = 1:itr_max%-- calculating next phi
a=(N2-N1)/(phi2-phi1);
phi=phi2-N2/a;%-- calculating stress distribution
ec = -eu + phi * y;er = -eu + phi * hs;
sc = ss_concrete(ec, param.concrete);sr = ss_steal(er, param.steel);
N = axial_force(y, sc, sr, param);

sample_02は元々スクリプトだが、paramを入力引数とする関数へ変更%-- checking axial force
if N<=err_max
break;
elseif N<0phi1 = phi;
N1 = N;elseif N>0
phi2 = phi;
N2 = N;
else%-- not usedbreak;
end%-- strain distribution
ec1 = -eu + phi1 * y;er1 = -eu + phi1 * hs;
ec2 = -eu + phi2 * y;er2 = -eu + phi2 * hs;%-- stress distribution
sc1 = ss_concrete(ec1, param.concrete);
sr1 = ss_steal(er1, param.steel);
sc2 = ss_concrete(ec2, param.concrete);sr2 = ss_steal(er2, param.steel);end%-- axial force
N1=axial_force(y, sc1, sr1, param);
N2=axial_force(y, sc2, sr2, param);

Mを返す

y0=eu/phi;
M = bending_moment(y,y0,sc, sr, param);end2019/6/16

モンテカルロ法の準備
Matlabでは対数正規分布を返す関数LognormalDistributionを使用するにはStatistics and Machine Learning Toolboxを購入する必要がある

平均と分散を与えると、対数正規分布のパラメータ𝜇と𝜎を求めて対数正規乱数を発⽣させる関数を自作した
function [x] = lognorm(n,mean,sd)%LOGNORM returns n random values obeying lognorm distribution
mu = log(mean^2/sqrt(mean^2+sd^2));
sigma = sqrt(log(sd^2/mean^2+1));
x = exp(sigma*randn(n,1)+mu);
end

randnは標準正規分布乱数を発⽣させる標準関数

信頼性工学特論
話題提供#02の資料参照
𝑋が対数正規分布（平均：𝜇 、分散：𝜎 ）なら、𝑌（= ln 𝑋）は正規分布（平均：𝜇 、分散𝜎 ）𝑌 = ln 𝑋
𝜇 = ln𝜇𝜇 +𝜎

𝜎 = ln𝜎+1𝜇2019/6/16

曲げモーメント使用
n=1000;%-- material
fc = lognorm(n, 30, 5); %-- Compression Strength
ec = lognorm(n, 0.002000, 0.000200); %-- Platcization Strain
eu = lognorm(n, 0.003500, 0.000200); %-- Ultmate Strain
E = lognorm(n, 200000, 1000); %-- Young Modulus 200GPa
fy= lognorm(n, 300, 10); %-- Yield Strength
xi=lognorm(n, 0.050, 0.010); %-- Hardening Curvature
Est=lognorm(n, 4500, 1000); %-- Hardening Coefficient
est=lognorm(n, 0.020000, 0.003000); %-- Hardening Strain (20000u)

sample_03
材料・部材パラメータの値を乱数として発⽣させる
たとえばfc（･･･𝑓 ：コンクリート圧縮強度）は平均：30[N/mm2]、分散：5[N/mm2]とした%-- member
H = lognorm(n, 600,1); %-- height
bw=lognorm(n, 400, 3); %-- width [mm]As=lognorm(n, 1435, 1); %-- Area of Riber
hs=lognorm(n, 500, 1); %-- height of Riber
M=zeros(n,1);
n個発⽣させた乱数の組み合わせの内、phi=zeros(n,1);
ii番目の値の組み合わせ（粒子という）を用いて
for ii = 1:n
終局曲げ耐力を計算
y=0:H(ii)/120:H(ii); %-- fiber height
concrete = init_concrete(fc(ii), ec(ii), eu(ii));
steel = init_steal(E(ii),fy(ii),xi(ii),Est(ii),est(ii));
member = init_member(y, bw(ii), As(ii), hs(ii));
param.concrete = concrete;param.steel = steel;param.member = member;[M(ii), phi(ii)] = capacity_calc(param);
fprintf([num2str(ii,'%0.4d') ': phi=' num2str(phi(ii),'%f')])fprintf([' Mu=' num2str(M(ii),'%f') '¥n'])end
histogram(M)ニュートン法で曲げモーメントを算出
sample_02を関数化した
capacity_calc関数を利用

以下を実行する>> sample_03

2019/6/16

実行結果
材料・部材パラメータを対数正規乱数とした時のRC梁の終局曲げ耐力分布が得られた>>histogram(fc)コンクリートの圧縮強度>>histogram(M)RC梁の終局曲げ耐力分布