Introduction#
The multi-rigid-body dynamics model is a simplified model for studying the vehicle system in environmental vibrations caused by traffic, serving as the basis for dynamic coupling analysis between the vehicle and road (or track, bridge) systems, which is more realistic compared to moving loads. Commonly used multi-rigid-body dynamics models for traffic-induced environmental vibrations include the quarter-car model, half-car model, and full-car model, which can simulate the dynamic characteristics of cars, trucks, and railway trains. Among these, the quarter-car model is the most simplified vehicle model and can provide satisfactory results for general environmental vibration problems. This paper presents the dynamic equations of a commonly used quarter-car model for railway trains under harmonic unevenness excitation without considering wheel-rail coupling, and provides a program for solving vehicle motion and wheel-rail contact forces.
Content#
Figure 1 shows a quarter-car model representing a railway train with secondary suspension, where:
- mc, mb, mw, and uc, ub, uw represent the mass and displacement of the car body, bogies, and wheels (to be solved);
- kp and cp represent the stiffness and damping coefficient of the primary suspension;
- ks and cs represent the stiffness and damping coefficient of the secondary suspension;
- q(t) represents the contact force between the wheel and rail (to be solved);
- Ar represents the track unevenness (displacement excitation), assuming no detachment between the wheel and rail, then Ar=uw.
Figure 1
According to D'Alembert's principle, isolating each rigid body block, the motion equations of the vehicle model can be expressed as a system of differential equations:
References provide an explicit solution for the contact force through Fourier transform (including the effects of gravity and vehicle speed), and here we directly use ode45 to obtain the time-domain solution (see Figure 2). It is worth noting that whether the vehicle model parameters are representative needs further verification.
% Vehicle parameters definition (Bian et al., 2011)
m_c = 0.25 * 4240; % One quarter of the vehicle mass [kg]
m_b = 0.5 * 3400; % Half of the bogie mass [kg]
m_w = 2200; % Wheelset mass [kg]
c_p = 5e3; % Primary suspension damping [Ns/m]
k_p = 1.04e6; % Primary suspension stiffness [N/m]
c_s = 6e3; % Secondary suspension damping [Ns/m]
k_s = 4e5; % Secondary suspension stiffness [N/m]
% Road unevenness input, harmonic function form
A = 1e-3; % Road unevenness amplitude [m]
omega = 2*pi*5; % Road unevenness angular frequency [rad/s], determined by vehicle speed and individual unevenness wavelength
phi = pi/2; % Phase [rad]
y_r = @(t) A*sin(omega*t + phi); % Displacement input
dy_r = @(t) A*omega*cos(omega*t + phi); % Velocity input
d2y_r = @(t) -A*omega^2*sin(omega*t + phi); % Acceleration input
% Initial conditions
z0 = 0; % Vehicle displacement
dz0 = 0; % Vehicle velocity
z_s0 = 0; % Bogie displacement
dz_s0 = 0; % Bogie velocity
x0 = [z0; dz0; z_s0; dz_s0];
% Time range
tspan = [0, 10];
% Solve the differential equations
[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);
% Calculate excitation force
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;
% Plot results
figure;
subplot(3,1,1);
plot(t, x(:,1));
title('Vehicle Displacement z(t)');
xlabel('Time t [s]');
ylabel('Displacement z [m]');
subplot(3,1,2);
plot(t, x(:,3));
title('Bogie Displacement z_s(t)');
xlabel('Time t [s]');
ylabel('Displacement z_s [m]');
subplot(3,1,3);
plot(t, F);
title('Total Excitation Force F(t)');
xlabel('Time t [s]');
ylabel('Force F [N]');
% State derivative function
function dxdt = stateDerivative(t, x, m1, m2, m3, c1, k1, c2, k2, y_r, dy_r, d2y_r)
z = x(1); % Vehicle displacement
dz = x(2); % Vehicle velocity
z_s = x(3); % Bogie displacement
dz_s = x(4); % Bogie velocity
y_r_val = y_r(t); % Calculate the value of y_r(t)
dy_r_val = dy_r(t); % Calculate the value of dy_r(t)
d2y_r_val = d2y_r(t); % Calculate the value of d2y_r(t)
% Calculate vehicle acceleration
ddz = (-c1*(dz - dz_s) - k1*(z - z_s)) / m1;
% Calculate bogie acceleration
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;
% State variables
dxdt = [dz; ddz; dz_s; ddz_s];
end
Figure 2
References#
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.