In Projects

MATLAB Mountain Bike Full Suspension Simulation

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.

  • #
    Figure 1
    Mountain Bike Suspension Motion

    Plots of mountain bike suspension simulation results.

MATLAB Code Download

Download the MATLAB program.

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