引言#
多剛體動力學模型是研究交通引發環境振動中車體系統的簡化模型，是進行車體與道路（或軌道、橋樑）系統動力耦合分析的基礎，相對於移動荷載更為真實。針對交通引發環境振動，常用的多剛體動力學模型包括四分之一車體模型、半模型和整車模型等，可以模擬汽車、卡車和鐵路列車等的動力特性。其中，四分之一車體模型是最為簡化的車體模型，對於一般的環境振動問題可以給出滿意的結果。
本文給出了一個常用的鐵路列車四分之一車體模型在簡諧不平順激勵下不考慮輪軌耦合的動力方程，並給出了求解車體運動及輪軌接觸力的程序。
內容#
圖 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.