# Issues Replicating Simulink Model within MATLAB (2DOF)

Tags:
1. Mar 11, 2017

### Sirsh

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..

Mainscript
Code (Text):

clear all;
clc;
%% 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
subplot(2,1,1)
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;
subplot(2,1,2)
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. Mar 16, 2017

### PF_Help_Bot

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.