Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Issues Replicating Simulink Model within MATLAB (2DOF)

  1. Mar 11, 2017 #1
    Hey all,

    I've been trying to debug a model I've created of a 4DOF geared system within Simulink by reproducing the results in MATLAB going from a 2DOF upwards to verify the dynamics are acting correctly.

    At the moment I am purely trying to model the stiffness between two gears, my equations are simply:
    From my reduced Simulink model, applying an initial displacement of the pinion of 1e-3m I get this response (sorry its hard to make out):
    Obviously from the figure above, you can see the equal amplitude oscillations with opposite direction (as these gears are not on the same shaft). I zoomed into the plot and verified that the initial condition of the gear is what I set it to be.

    However, I've written the following MATLAB code to try and recreate the above, but to no avail. I'm not sure where I am going wrong here..

    Code (Text):

    clear all;
    %% Defining global variables
    global r_p r_g I_p I_g k_mb
    % Assignment of ODE variable values:
    r_p = 0.1297;   % Pinion radius (m)
    r_g = 0.1297;   % Gear radius (m)
    I_p = 0.0025;   % Pinion inertia (kg.m^2)
    I_g = 0.0025;   % Gear inertia (kg.m^2)
    k_mb = 100;
    %% initial coditions for ode solver
    z0 = [1e-3,0,0,0]; % Initial Conditions for ODE45.
    t_span = [0 1];         % Timespan 0 - 1 second
    opts = odeset('Reltol', 1e-3, 'AbsTol', 1e-3);
    [t,z] = ode45(@ode2dof, t_span, z0, opts);
    %% Post Processing
    % Renaming the ode45 outputs for easier usage.
    theta_2 = z(:,1);
    theta_2d = z(:,2);
    theta_3 = z(:,3);
    theta_3d = z(:,4);

    % Plotting
    plot(t,theta_2,t,theta_3), xlabel('Time (s)'), ...
        ylabel('Angular Displacement (rad)'),title('Angular Displacements'),...
        legend('theta 2','theta 3','Location', 'northwest'),grid;
    plot(t,theta_2d,t,theta_3d), xlabel('Time (s)'), ...
        ylabel('Angular Velocity (rad/s)'), title('Angular Velocities'),...
        legend('theta 2d', 'theta 3d','Location', 'northwest'),grid;
    ODE Script
    Code (Text):

    function zdot = ode2dof(t,z)
    % defining global variables
    global r_p r_g I_p I_g k_mb
    %% S.S.E:
    zdot(2) = ((k_mb*r_p)/I_p)*(z(3)*r_g-z(1)*r_p);     % theta_2dd
    zdot(4) = ((k_mb*r_g)/I_g)*(z(1)*r_p-z(3)*r_g);     % theta_3dd

    z(2) = zdot(1);
    z(4) = zdot(3);
    % Outputs
    zdot = zdot';
    The output plot looks like this instead:
    Any help would be appreciated, as I'm lost at what's going on here even though its quite a simple system.
  2. jcsd
  3. Mar 16, 2017 #2
    Thanks for the thread! This is an automated courtesy bump. Sorry you aren't generating responses at the moment. Do you have any further information, come to any new conclusions or is it possible to reword the post? The more details the better.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted