This MATLAB program simulates a full suspension mountain bike linkage. The simulation provides kinematics of all of the links including position, velocity, and acceleration. The simulation results show how the suspension will behave and how potential changes in the design will affect the suspension dynamics. The program makes use of algorithms for Jacobian matrix operations, vector loops to describe the kinematics of linkages, and programming a Newton-Rhapson algorithm to determine solutions to non-linear systems of equations.

Results are shown for a full suspension mountain bike through its full range of travel.

The first step in analyzing the VPP suspension design was to determine the geometry of the system. First, two vector loop equations were created from the links of the suspension. These vector loops were broken down into x and y components, leading to four position equations. The time derivates of these four position equations yields the four velocity and acceleration equations. Using trigonometry, the fixed lengths and angles of the suspension were calculated. After plugging the fixed geometries into the four position equations there are only four unknown variables left. Since there are four equations and four unknowns, the solution can be determined. However, this system of equations is non-linear so a Newton-Rhapson algorithm must be used. This algorithm makes an initial guess of the values of the unknowns and then calculates the next guess using the inverse Jacobian matrix of the four position equations multiplied by those four position equations. After a number of iterations, the initial guess converges very close to the true values of the unknowns. Now that all of the link angles and lengths are known, they are plugged into the four derived velocity equations. This leaves only four velocity unknowns. The same step is repeated for the acceleration equations to get the four acceleration unknowns. The math is done in MATLAB through matrix algebra, specifically, the inverse Jacobian matrix of the position equations multiplied by the four derived equations in a matrix excluding the unknown terms. This multiplication of matrices yields a third matrix of the solutions to the four unknowns.

• Plots of mountain bike suspension simulation results.

##### MATLAB Code
```        % 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
% vr5 = Velocity of linkage r5 (in/s)
% 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);  % V 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)
title('Spring Compression Relative to Input Angle');

% Figure 2
subplot(3,2,2)
plot(record_th2, record_vr5)
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)
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');
```