Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
% Program:
% mountain-bike-suspension.m
% Mountain bike suspension linkage dynamics simulation.
%
% Description:
% Determines the positions, velocities, and accelerations of a full
% mountain bike suspension linkage through the use of vector loops Newton
% Raphson iteration and the Jacobian matrix. The graphs plot system
% dynamics of key components of the suspension as the shock is fully
% compressed.
%
% Variable List:
% r1 = Linkage length 1
% r2 = Linkage length 2
% r3 = Linkage length 3
% r4 = Linkage length 4
% r5 = Linkage length 5
% r7 = Linkage length 7
% th1 = Linkage angle 1
% th2 = Linkage angle 2
% th3 = Linkage angle 3
% th4 = Linkage angle 4
% th5 = Linkage angle 5
% th7 = Linkage angle 7
% vth2 = Angular velocity of linkage angle 1 (rad/s)
% vth3 = Angular velocity of linkage angle 3 (rad/s)
% vth4 = Angular velocity of linkage angle 4 (rad/s)
% vth5 = Angular velocity of linkage angle 5 (rad/s)
% vr5 = Velocity of linkage r5 (in/s)
% ath2 = Angular acceleration of linkage angle 2 (rad/s^2)
% ath3 = Angular acceleration of linkage angle 3 (rad/s^2)
% ath4 = Angular acceleration of linkage angle 4 (rad/s^2)
% ath5 = Angular acceleration of linkage angle 5 (rad/s^2)
% ar5 = Acceleration of linkage r5 (in/s)
% F1 = Vector loop equation 1
% F2 = Vector loop equation 2
% F3 = Vector loop equation 3
% F4 = Vector loop equation 4
% F = Vector loop equations array
% veq1 = Velocity equation 1
% veq2 = Velocity equation 2
% veq3 = Velocity equation 3
% veq4 = Velocity equation 4
% veqans = Array of velocity equations
% aeq1 = Acceleration equation 1
% aeq2 = Acceleration equation 2
% aeq3 = Acceleration equation 3
% aeq4 = Acceleration equation 4
% aeqans = Acceleration equations array
% x = Positions array
% v = Velocity array
% a = Acceleration array
% unknown = Unknown variables array
% m = Counter
% n = Counter
% Tolerance = Newton Raphson error tolerance
% G = Jacobian matrix
% change = New array for defining the sum of delta x
% check = Total change value
% delta_x = Delta x used in Newton Rhapson method
% step_size = Step size in radians
% stopping_point = Point to stop running simulation
% record_th3 = Record th3 for plotting
% record_th4 = Record th4 for plotting
% record_th5 = Record th5 for plotting
% record_r5 = Record r5 for plotting
% record_vth3 = Record vth3 for plotting
% record_vth4 = Record vth4 for plotting
% record_vth5 = Record vth5 for plotting
% record_vr5 = Record vr5 for plotting
% record_ath3 = Record ath3 for plotting
% record_ath4 = Record ath4 for plotting
% record_ath5 = Record ath5 for plotting
% record_ar5 = Record ar5 for plotting
% record_th2 = Record th2 for plotting
% xtravel= Back axle x travel relative to frame
% ytravel = Back axle y travel relative to frame
% velvert = Vertical velocity of back axle
clear, clc % Clear command window and workspace
% Define angle sympols for symbolic functions
th1 = sym('th1');
th2 = sym('th2');
th3 = sym('th3');
th4 = sym('th4');
th5 = sym('th5');
th7 = sym('th7');
% Define length sympols for symbolic functions
r1 = sym('r1');
r2 = sym('r2');
r3 = sym('r3');
r4 = sym('r4');
r5 = sym('r5');
r7 = sym('r7');
% Define angular velocity sympols for symbolic functions
vth2 = sym('vth2');
vth3 = sym('vth3');
vth4 = sym('vth4');
vth5 = sym('vth5');
vr5 = sym('vr5');
% Define angular acceleration sympols for symbolic functions
ath2 = sym('ath2');
ath3 = sym('ath3');
ath4 = sym('ath4');
ath5 = sym('ath5');
ar5 = sym('ar5');
% Input lengths and angles of suspension system
Tolerance = 10^-6; % Newton Raphson error tolerance
r3 = 9.0; % Length of link #3
r2 = 2.236; % Length of link #2
r4 = 2.0616; % Length of link #4
r7 = 7.3824; % Length of link #7
r1 = 9.124; % Length of link #1
th1 = 1.4056; % Angle of link #1
th7 = 5.789; % Angle of link #7
vth2 = 1; % Velocity of link #2
ath2 = 0; % Acceleration of link #2
% Position equations derived from vector loops 1 and 2 in x and y components
F1 = r4*cos(th4) - r7*cos(th7) + r5*cos(th5); % 1, x
F2 = r4*sin(th4) - r7*sin(th7) + r5*sin(th5); % 1, y
F3 = r2*cos(th2) + r3*cos(th3) + r4*cos(th4) - r1*cos(th1); % 2, x
F4 = r2*sin(th2) + r3*sin(th3) + r4*sin(th4) - r1*sin(th1); % 2, y
F = [F1; F2; F3; F4]; % Vector loop equations array
unknown = [th3;th4;th5;r5]; % Unknown variable array
G = jacobian(F,unknown); % Compute Jacobian matrix with symbols
delta_x = -inv(G)*F; % Delta x for Newton Rhapson method
step_size = 0.01; % Step size in radians
% Velocity equations derived from position equations
veq1 = 0;
veq2 = 0;
veq3 = r2*vth2*sin(th2);
veq4 = -r2*vth2*cos(th2);
veqans = [veq1; veq2; veq3; veq4];
% Acceleration equations derived from position equations
aeq1 = -r4*vth4^2*cos(th4) - vr5*sin(th5)*vth5 - r5*vth5^2*cos(th5);
aeq2 = -r4*vth4^2*sin(th4) + vr5*cos(th5)*vth5 - r5*vth5^2*sin(th5) + vr5*cos(th5)*vth5;
aeq3 = -r2*ath2*sin(th2) - r2*vth2^2*cos(th2) - r3*vth3^2*cos(th3) - r4*vth4^2*cos(th4);
aeq4 = r2*ath2*cos(th2) - r2*vth2^2*sin(th2) - r3*vth3^2*sin(th3) - r4*vth4^2*sin(th4);
aeqans = [aeq1; aeq2; aeq3; aeq4];
x = [1.6; 1.3; 5.6; 8.1]; % Initial guess of unknown positions
r5 = x(4,1); % Initial guess of r5
th2 = 5.176; % Initial th2
m = 1; % Initialize counter
stopping_point = 5.8885; % r5 = 8.1385 - 2.25 = 5.8885
while r5 > stopping_point % Run until stop point
n = 1; % Initialize counter
check = 1; % Initialize check
while check > Tolerance % Stops iterations when delta x is small
% Call initial guess values from x array above
th3 = x(1,1);
th4 = x(2,1);
th5 = x(3,1);
r5 = x(4,1);
x = x + eval(delta_x); % Newton Rhapson increment
change = eval(delta_x); % New array for sum of delta x
% Sum of delta x used for while loop
check = abs(change(1,1) + change(2,1) + change(3,1) + change(4,1));
n = n + 1; % Increment counter
end
% Set unknown positions to last iteration of Newton Rhapson
th3 = x(1,1);
th4 = x(2,1);
th5 = x(3,1);
r5 = x(4,1);
% Compute unknown velocities using vel equation array and positions
% solved from Newton Rhapson
v = eval(inv(G)*veqans);
% Set unknown velocities variables to the corresponding answer from
% solution array 'v' above
vth3 = v(1,1);
vth4 = v(2,1);
vth5 = v(3,1);
vr5 = v(4,1);
% Compute accel unkowns from accel equations and calculated velocities
% and positions
a = eval(inv(G)*aeqans);
% Set unknown accel variables to the corresponding answer from solution
% array 'a' above
ath3 = a(1,1);
ath4 = a(2,1);
ath5 = a(3,1);
ar5 = a(4,1);
% Record unknown positions at given current input angle (for plotting)
record_th3(m) = x(1,1);
record_th4(m) = x(2,1);
record_th5(m) = x(3,1);
record_r5(m) = x(4,1);
% Record unknown velocities at given current input angle (for plotting)
record_vth3(m) = v(1,1);
record_vth4(m) = v(2,1);
record_vth5(m) = v(3,1);
record_vr5(m) = v(4,1);
% Record unknown accelerations at given current input angle (for
% plotting)
record_ath3(m) = a(1,1);
record_ath4(m) = a(2,1);
record_ath5(m) = a(3,1);
record_ar5(m) = a(4,1);
record_th2(m) = th2; % Record current th2 for this configuration
xtravel(m) = r2*cos(record_th2(m)) + 15*cos(record_th3(m) - 1.5708); % Back axle x travel
ytravel(m) = r2*sin(record_th2(m)) + 15*sin(record_th3(m) - 1.5708); % Back axle y travel
velvert(m) = vth2*r2*cos(th2) + record_vth3(m)*15*cos(record_th3(m) - 1.5708); % Vertical vel of back axle
fprintf('Run = %3.3i \n', m);
th2 = th2 + step_size; % Increment th2
m = m + 1; % Increment counter
end
% Figure 1
subplot(3,2,1)
plot(record_th2, record_r5)
xlabel('Input Angle (radians)');
ylabel('Length of Spring Link (inches)');
title('Spring Compression Relative to Input Angle');
% Figure 2
subplot(3,2,2)
plot(record_th2, record_vr5)
xlabel('Input Angle (radians)');
ylabel('Velocity of Spring (in/s)');
title('Spring Compression Rate Relative to Input Angle');
% Figure 3
subplot(3,2,3)
plot(record_th2, record_ar5)
xlabel('Input Angle (radians)');
ylabel('Acceleration of Spring (in/s^2)');
title('Acceleration of Spring Compression Relative to Input Angle');
% Figure 4
subplot(3,2,4)
plot(xtravel, ytravel)
xlabel('Position in X direction (in)');
ylabel('Position in Y direction (in)');
title('X-Y Motion of Center of Rear Wheel Relative to Frame at Input Link');
% Figure 5
subplot(3,2,5)
plot(velvert, record_vr5)
xlabel('Vertical Component of Velocity of Rear Wheel (in/s)');
ylabel('Velocity of Spring (in/s)');
title('Spring Compression Rate Relative to Vertical Velocity of Rear Axle');
% Figure 6
subplot(3,2,6)
plot(ytravel, record_r5)
xlabel('Y motion of Rear Wheel (in/s)');
ylabel('Compression of Spring (in/s)');
title('Spring Compression Relative to Y Motion of Rear Axle');
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment