1自由度系非線形地震応答をMATLABで計算するサンプルコード
バイリニア型履歴復元力特性を有する1自由度系の非線形地震応答解析をするMATLABのコードです． 改変等，適宜修正していただいて，ご自由にお使いください．
% バイリニア型1自由度系の非線形地震応答解析
% 2021-11-22 伊津野和行
% 必要なファイル eq.txt 地震加速度のみがgal=cm/s/s単位で入っているファイル
% 結果は results.xlsx に時刻歴応答が出力される
% 積分は酒井らの方法 https://doi.org/10.2208/jscej.1995.507_137
% 以下のコメントで式(?)とあるのは，同上論文の式番号．
clear;
% パラメータの入力 ここから
m = 100; % 質量 t（あるいはkg）
k1 = 10000; % 初期剛性 kN/m（質量がkgならN/m）
k2 = 1000; % 二次剛性 kN/m（質量がkgならN/m）
yy = 500; % 降伏耐力 kN（質量がkgならN）
h = 0.05; % 減衰定数 無次元
deltaT = 0.02; % 地震波の時間刻み
% パラメータの入力 ここまで
omega = ( k1 / m )^0.5; % 固有円振動数 rad/s
T = 2 * pi / omega; % 固有周期 s
msg = ['固有周期 ' num2str(T) ' (s)'];
disp(msg)
c = 2 * h * omega * m; % 粘性減衰係数 kN･s/m
bun = 10; % 計算時間間隔分割数
dt = deltaT / bun; % 計算用時間刻み
eq = readmatrix('eq.txt'); % 地震波の読み込み gal=cm/s/s
eq = 0.01 * eq; % 地震加速度の単位変換 m/s/s
n = size(eq,1); % 地震波データ数
dt2 = dt * dt;
tmp1 = k1 + 2 * c / dt + 4 * m / dt2; % 計算途中で使う係数
tmp2 = m + dt * c / 2;
output = zeros(n,6); % 出力用の変数
% 初期値
x = 0; % 変位
v = 0; % 速度
a = 0; % 加速度
f = 0; % 前ステップにおける復元力
Q = 0; % 調整外力（非線形と線形との差）
z0 = 0;
for i = 1:n-1
z1 = eq(i + 1);
dz = (z1 - z0) / bun;
for j = 1:bun
z0 = z0 + dz;
P = -m * z0; % 地震外力
x2 = (P + Q + m * (4 * x / dt2 + 4 * v / dt + a) ...
+ c * (2 * x / dt + v)) / tmp1; % 式(26)
f2 = bilinear( x2, x, yy, k1, k2, f ); % 変位x2に対するバイリニアの復元力
Q2 = k1 * x2 - f; % 式(3)
dq = Q2 - Q;
a = -a - 4 * v / dt + 4 * (x2 - x) / dt2 + dq / tmp2; % 式(28)
v = -v + 2 * (x2 - x) / dt + dq * dt / 2 / tmp2; % 式(27)
% 式(28)を式(27)より先に計算しないとvが更新されてしまうので注意
x = x2; % 値を更新
f = f2;
Q = Q2;
end % for j
% ファイルへの結果出力
output(i+1,1) = i * deltaT;
output(i+1,2) = eq(i+1);
output(i+1,3) = x;
output(i+1,4) = v;
output(i+1,5) = a + eq(i+1);
output(i+1,6) = f;
end % for i
cells = ["時間 (s)" "地震加速度 (m/s/s)" "変位 (m)" ...
"速度 (m/s)" "加速度 (m/s/s)" "復元力 (kN)"];
% 質量の単位をkgにした場合は，復元力 (kN)を復元力 (N)に変更．
% 次の2行でエクセルファイルに出力．
writematrix(cells,'results.xlsx','WriteMode','overwritesheet');
writematrix(output,'results.xlsx','Range','A2');
disp('NORMAL END')
function y = bilinear( x, xold, yy, k1, k2, yold )
% 変位xに対するバイリニア型（移動硬化型バイリニア）履歴復元力を計算
% x:変位 xold:前ステップの変位 yy:降伏耐力
% k1:初期剛性 k2:降伏後の二次剛性 yold:前ステップの復元力
% 考え方 https://www.ritsumei.ac.jp/se/rv/izuno/programming.html
xy = yy / k1; % 降伏変位
yLinear = k1 * (x - xold) + yold;
% 前ステップの変位から初期剛性で変化した場合の復元力
yPlus = k2 * (x - xy) + yy; % 第二剛性で決まる上側の骨格曲線上の復元力
yMinus = k2 * (x + xy) - yy; % 第二剛性で決まる下側の骨格曲線上の復元力
y = median( [yLinear yPlus yMinus] );
% 3つの復元力の中央値が，次ステップの復元力になる
end
必要な地震波ファイル eq.txt は，次のようなファイルです．
0 0.09 -60.1 -69.97 -13.78 24.19
![[ページのトップへ戻る]](images/totop.png)