mountain-bike-suspension.m

% 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');