引言#
多刚体动力学模型是研究交通引发环境振动中车体系统的简化模型，是进行车体与道路（或轨道、桥梁）系统动力耦合分析的基础，相对于移动荷载更为真实。针对交通引发环境振动，常用的多刚体动力学模型包括四分之一车体模型、半模型和整车模型等，可以模拟汽车、卡车和铁路列车等的动力特性。其中，四分之一车体模型是最为简化的车体模型，对于一般的环境振动问题可以给出满意的结果。
本文给出了一个常用的铁路列车四分之一车体模型在简谐不平顺激励下不考虑轮轨耦合的动力方程，并给出了求解车体运动及轮轨接触力的程序。
内容#
图 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.