ENGINEERING

Home Engineering
  • Engineering
  • Senior Project
  • Projects
  • 3 Doors Game Simulation
  • Woodworking: Bedstand
  • Woodworking: Desk
  • Numbers Game
  • Path Planning
  • Snake Video Game AI
  • Minimum Spanning Tree
  • Heat Map
  • MATLAB: Rocket Simulation
  • MATLAB: Bike Suspension
  • MATLAB: Thermal FEA SS
  • MATLAB: Thermal FEA Tr
  • Notes
  • Statistics
  • Steoscopic Estimates
  • Cubed Roots
  • Other
  • Free Fall Distance
  • Speed of Sound
  • Travel
  • Solar System Scale
  • Videos
  • Yellowstone Geysers
  • Geysers and Waterfalls
  • Yellowstone and Glacier
  • Grinnell Glacier
  • Yellowstone Geysers
  • Photography

    Subjects

    MATLAB Mountain Bike Full Suspension Simulation



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


    The first step in analyzing the VPP suspension design was to determine the geometry of the system. First of all, two vector loop equations were created from the links of the suspension (See appendix C for all hand calculations regarding kinematic equations). These vector loops were broken down into x and y components, thus generating 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 (Appendix D). After plugging in the fixed geometries into the four position equations there were only four unknown variables left; Since there are four equations and four unknowns, we can solve for the unknowns. However, this system of equations is not 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 change of the 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 4 velocity unknowns; The same step is repeated for the acceleration equations to get the four acceleration unknowns; The “plugging in and solving” 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.


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