Car vibration

%       A simple model of car vibrations.
%       To produce hardcopies, delete the comment sign, %, before
%       the PRINT statements.

m = 1200;       % car mass, 4 passengers plus luggage, kg
k = 70000;      % spring rate, N/m
c = 5000;       % damping coefficient, Ns/m
omega = 5.7;    % excitation angular frequency, rad/s
Y = (k + c*omega*i)/(-m*omega^2 + k + c*omega*i);
amplification = abs(Y)          % real amplification factor
phase = angle(Y);               % response phase, rad
phase_d = 180*phase/pi          % response phase, deg
B = 0.1;        % excitation amplitude, m
A = amplification*B;   % response amplitude, m
Fs = [ 0 k ]*A;                 % spring force, N
Fd = Fs(2) + i*c*omega*A;       % damping force, N
Fi = Fd – m*omega^2*A;          % inertia force, N
output = [ Fs Fd Fi ];
Bs = [ 0 k ]*B*exp(-i*phase);    % spring component of excitation, N
Bd = Bs(2) + i*c*omega*B*exp(-i*phase);   % damping component of excitation, N
input = [ Bs Bd ];
% Plot now Figure 4.14
clg
plot(real(output), imag(output), real(input), imag(input))
axis(‘square’), axis([ 0 14000 0 14000 ])
title(‘Phasor diagram for omega = 5.7 rad/s’)
xlabel(‘Real, N’), ylabel(‘Imaginary, N’)
text(Fs(2)/2, 250, ‘k*A’)
text(11000, imag(Fd)/2, ‘i*c*omega*A’)
text((real(Fd) + real(Fi))/2, imag(Fi)+500, ‘m*omega^2*A’)
text(real(Bs(2))/2, imag(Bs(2))/2+400, ‘k*B’)
text((real(Bd) + real(Bs(2)))/2+150, imag(Bd)/2+950, ‘i*c*B’)
text(real(Bs(2))/2-500, imag(Bs(2))/4, ‘phase angle’)
% print
pause
zeta = c/(2*sqrt(k*m))
omega_n = sqrt(k*m);
r = 0:0.05:3;
s = 2*zeta*r*i;
Y = (1 + s)./(1 – r.^2 + s);
clg
% Plot now Figure 4.15
axis
axis(‘normal’)
subplot(2, 1, 1)
plot(r, abs(Y)), grid
title(‘Frequency response’)
ylabel(‘Amplification |A/B|’)
subplot(2, 1, 2)
plot(r, 180*angle(Y)/pi), grid
xlabel(‘Frequency ratio omega/omega_n’)
ylabel(‘Phase angle, deg’)
% print
pause
% begin BODE plot, Figure 4.16
clg
subplot(211)
semilogx(r, 20*log(abs(Y))), grid
title(‘Frequency response’)
ylabel(‘Amplification |A/B|, dB’)
subplot(212)
semilogx(r, 180*angle(Y)/pi), grid
xlabel(‘Frequency ratio omega/omega_n’)
ylabel(‘Phase angle, deg’)
% print