引言#
多刚体动力学モデルは、交通によって引き起こされる環境振動における車体システムの簡略化モデルであり、車体と道路（または軌道、橋梁）システムの動的結合分析を行うための基礎であり、移動荷重に対してより現実的です。交通によって引き起こされる環境振動に対して、一般的に使用される多刚体动力学モデルには四分の一車体モデル、半モデル、全車モデルなどが含まれ、これらは自動車、トラック、鉄道列車などの動的特性をシミュレートできます。その中で、四分の一車体モデルは最も簡略化された車体モデルであり、一般的な環境振動問題に対して満足のいく結果を提供できます。
本稿では、一般的な鉄道列車の四分の一車体モデルが簡諧不平順励起の下で輪軌結合を考慮せずに動的方程式を示し、車体の運動および輪軌接触力を求めるプログラムを提供します。
内容#
図 1 は、二系悬挂を持つ鉄道列車の四分の一車体モデルを示しています。ここで：
- mc、mb、mw および uc、ub、uw はそれぞれ車両（car bodys）、台車（bogies）、および車輪（wheels）の質量と変位 （求解量） を表します；
- kp および cp はそれぞれ一系悬挂（the primary suspension）の剛性と減衰係数を表します；
- ks および cs はそれぞれ二系悬挂（the secondary suspension）の剛性と減衰係数を表します；
- q (t) は輪軌間の接触力 （求解量） を表します；
- Ar は軌道の不平順（変位励起）を表し、輪軌が脱空しないと仮定すると、Ar=uw となります。
図 1
ダランベールの原理に基づいて各剛体ブロックを隔離体分析し、車体モデルの運動方程式は微分方程式系として表されます：
参考文献では、フーリエ変換を用いて接触力の明示的な解答を直接示しています（重力と車速の影響を含む）。ここでは ode45 を用いて時域解を求めます（図 2 参照）。なお、車体モデルのパラメータが代表性を持つかどうかはさらに検証が必要です。
% 車体パラメータ定義 (Bian et al., 2011)
m_c = 0.25 * 4240; % 車体質量 [kg]の四分の一
m_b = 0.5 * 3400; % 台車質量 [kg]の二分の一
m_w = 2200; % 車輪対質量 [kg]
c_p = 5e3; % 一系悬挂減衰 [Ns/m]
k_p = 1.04e6; % 一系悬挂剛性 [N/m]
c_s = 6e3; % 二系悬挂減衰 [Ns/m]
k_s = 4e5; % 二系悬挂剛性 [N/m]
% 道路不平順入力、簡諧関数形式
A = 1e-3; % 道路不平順振幅 [m]
omega = 2*pi*5; % 道路不平順角周波数 [rad/s]、車速と単一の不平順波長によって決定される
phi = pi/2; % 位相 [rad]
y_r = @(t) A*sin(omega*t + phi); % 変位入力
dy_r = @(t) A*omega*cos(omega*t + phi); % 速度入力
d2y_r = @(t) -A*omega^2*sin(omega*t + phi); % 加速度入力
% 初期条件
z0 = 0; % 車体変位
dz0 = 0; % 車体速度
z_s0 = 0; % 台車変位
dz_s0 = 0; % 台車速度
x0 = [z0; dz0; z_s0; dz_s0];
% 時間範囲
tspan = [0, 10];
% 微分方程式を解く
[t, x] = ode45(@(t,x) stateDerivative(t, x, m_c, m_b, m_w, c_s, k_s, c_p, k_p, y_r, dy_r, d2y_r), tspan, x0);
% 激振力を計算
F_c = m_w.*d2y_r(t) + c_p.*(dy_r(t) - x(:,4)) + k_p.*(y_r(t) - x(:,3));
F = (m_c + m_b + m_w) * 9.81 - F_c;
% 結果をプロット
figure;
subplot(3,1,1);
plot(t, x(:,1));
title('車体変位 z(t)');
xlabel('時間 t [s]');
ylabel('変位 z [m]');
subplot(3,1,2);
plot(t, x(:,3));
title('台車変位 z_s(t)');
xlabel('時間 t [s]');
ylabel('変位 z_s [m]');
subplot(3,1,3);
plot(t, F);
title('総激振力 F(t)');
xlabel('時間 t [s]');
ylabel('力 F [N]');
% 状態導出関数
function dxdt = stateDerivative(t, x, m1, m2, m3, c1, k1, c2, k2, y_r, dy_r, d2y_r)
z = x(1); % 車体変位
dz = x(2); % 車体速度
z_s = x(3); % 台車変位
dz_s = x(4); % 台車速度
y_r_val = y_r(t); % y_r(t) の値を計算
dy_r_val = dy_r(t); % dy_r(t) の値を計算
d2y_r_val = d2y_r(t); % d2y_r(t) の値を計算
% 車体加速度を計算
ddz = (-c1*(dz - dz_s) - k1*(z - z_s)) / m1;
% 台車加速度を計算
ddz_s = (-c1*dz_s + c1*dz - k1*z_s + k1*z - 2*c2*dz_s + 2*c2*dy_r_val - 2*k2*y_r_val + 2*k2*y_r(t) + m3*d2y_r_val) / m2;
% 状態変数
dxdt = [dz; ddz; dz_s; ddz_s];
end
図 2
参考文献#
Bian, Xue-cheng, et al. "A 2.5 D finite element approach for predicting ground vibrations generated by vertical track irregularities." Journal of Zhejiang University-Science A 12.12 (2011): 885-894.