-
-
Save KevinKParsons/b4968adbb65659ccf3648de40c6f2aee to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
% 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