EE341.01: MATLAB M-FILE FOR SUSPENSION SIMULATION

This example shows a MATLAB M-file for simulating a suspension system, animating the system, and plotting the input and output positions for a sinusoidal input.

%
% Filename: example8.m
%
% Description: M-file to simulate front bike suspension for sinusoidal
%              inputs.

clear; 					% clear memory and figure
clf;

d = 0.9;				% equil distance: bike-to-ground in m
m = 40;					% mass of 1/2(bike+person) in kg
b = 100;				% damping coefficient
k = 79000;				% spring coefficient
g = 9.81;				% acceleration of gravity

t   = 0:0.0001:0.02;			% define range of times
wo  = 2*pi/0.006;			% T=6ms for 4cm bump spacing at 15mph
A   = 0.02;                      	% 2cm bumps
xg  = A*cos(wo*t);			% ground position (input)

xb1 = (d - m*g/k)*ones(size(t));	% bike position due to gravity input
H = (j*b*wo+k)/(-m*wo*wo+j*b*wo+k);	% FT of impulse response
xb2 = A*abs(H)*cos(wo*t+angle(H));      % bike position due to xg input

xb = xb1 + xb2;                         % sum outputs together by linearity

%*************** Plot ground position ************************************
figure(1);

subplot(2,2,1);
plot(t,xg);
xlabel('t (seconds)')
ylabel('xg(t)')
title('Ground Position');

%*************** Plot bike position **************************************
subplot(2,2,2);				
plot(t,xb);
xlabel('t (seconds)')
ylabel('xb(t)')
title('Bike Position');

%*************** Plot magnitude spectrum of H(w) *************************
w = -500:1:500;
H = (j*b*w+k)./(-m*w.*w+j*b*w+k);
subplot(2,2,3);
plot(w,abs(H));
xlabel('w (rad/sec)');
ylabel('|H(w)|');
title('Magnitude Spectrum of H(w)');

%*************** Plot phase spectrum of H(w) *****************************
subplot(2,2,4);
plot(w,angle(H)*180/pi);
xlabel('w (rad/sec)');
ylabel('/_ H(w)');
title('Phase Spectrum of H(w)');

%*************** Call drawsusp.m to animate suspension *******************
figure(2);
drawsusp(t, xg, xb);

function drawsusp(t, xg, xb)
% 
% Filename: drawsusp.m
%
% Description: This M-file contains a function to animate the motion
%              of a bicycle suspension system.  It needs to be sent
%              the time vector t, the ground position xg, and the bike
%              position xb.

fll = 0.4;				% fixed length of lower fork
wr = 0.33;                              % wheel radius of 13"
theta = 83*pi/180;                      % fork rake angle of 83 degrees

arcstep = 1;				% angle increment (in degrees) 
                                        % to sweep out arc lengths when
                                        % creating sides of wheel

lt       = length(t);			% find length of t
wheelctr = xg + wr;			% find wheel center

j = 0:arcstep:(360-arcstep);     	% create wheel sides using arclengths
arcx = wr * cos((j+arcstep) * pi/180); 	% x coordinates of wheel sides
arcy = wr * sin((j+arcstep) * pi/180); 	% y coordinates of wheel sides

				        % construct initial fork with wide line
					% LineWidth of 0.5 in normal
FL = plot([0 fll*cos(theta)], [wheelctr(1) (wheelctr(1)+fll*sin(theta))], ...
          'k', 'EraseMode', 'xor', 'LineWidth', [15.0]);
hold;

FU = plot([fll*cos(theta) (xb(1)-wheelctr(1))/tan(theta)], ...
          [(wheelctr(1)+fll*sin(theta)) xb(1)], ...
          'r', 'EraseMode', 'xor','LineWidth', [13.0]);

axis([-0.75 0.75 -0.25 1.25]);		% constuct square axes
axis('square')
					% construct initial wheel
W = plot(arcx, arcy+wheelctr(1), 'k', 'EraseMode', 'xor', 'LineWidth',[24.0]); 

hold;

'Press a key to continue'

pause

for i = 2:lt,				% loop over data values

  set(FL, 'XData', [0 fll*cos(theta)]);	% change x coordinates of low fork ends
  					% change y coordinates of low fork ends
  set(FL, 'YData', [wheelctr(i) (wheelctr(i)+fll*sin(theta))]);
					% change x coordinates of hi fork ends
  set(FU, 'XData', [fll*cos(theta) (xb(i)-wheelctr(i))/tan(theta)]);
					% change y coordinates of hi fork ends
  set(FU, 'YData', [(wheelctr(i)+fll*sin(theta)) xb(i)]);
  set(W, 'XData', arcx);		% change x coordinate of wheel
  set(W, 'YData', arcy+wheelctr(i));	% change y coordinate of wheel
  drawnow;
  
  for j = 1:50,                   	% delay loop
    cos(j.*cos(j));
  end;
end;