How to solve for the torsional degrees of freedom

I am trying to model the torsional degrees of freedom for a single stage helical gear system. My orignal code is here: https://uk.mathworks.com/matlabcentral/answers/1916150-keep-getting-a-blank-graph?s_tid=es_ans_avr_que_view

Please do not flag this question, I am trying to ask something different

For example-

I_a = moment of inertia of smaller driver gear

I_p = moment of inertia of bigger driven gear

R_a = radius of smaller driver gear

Cs = meshing damping

theta_a = torsional angle of smaller driver gear

K_s = meshing stiffness

R_p = radius of bigger driven gear

theta_a = torsional angle of bigger driven gear

F_y = meshing excitation force

How can I solve for theta_a and theta_p using ode45 in a similar way as we solve for x in a simple mass spring damper system

function Ydot = num_for(t, Y)
    Ydot = zeros(2,1);
    m  = 100; 
    k  = 1000; 
    c  = 160;
    w  = 5; 
    f  = 160; 
    F  = [0; (f/m)*cos(w*t)];	% input force vector
    A  = [0 1; -k/m -c/m];      % state matrix
    Ydot = A*Y + F;             % state-space model
end

function call-

Tspan  = [0 10]; 
y0     = [0.01; 0.1]; 
[t, y] = ode45(@num_for, Tspan, y0); 
figure 
plot(t, y(:,1), 'linewidth', 1.5) 
grid on 
xlabel('Time, t [sec]') 
ylabel({'displacement'}) 

Answer 1

[Sam Chak]

Hi @Siddharth Jain

You can try modifying the following code to fit your application.

tspan  = [0 0.1]; 
y0     = [0.1 0.2 0 0]; 
[t, y] = ode45(@heligear, tspan, y0); 

for j = 1:4
    figure(j)
    plot(t, y(:,j), 'linewidth', 1.5), grid on
    xlabel('Time, t [sec]'), 
    ystate = sprintf('y_{%d}', j);
    ylabel(ystate)
end

function ydot = heligear(t, y)
    ydot    = zeros(4, 1);
    
    % parameters
    m_a     = 0.5;
    m_p     = 0.5; 
    R_a     = 51.19e-3;
    R_p     = 139.763e-3;
    I_a     = 0.5*m_a*(R_a^2);
    I_p     = 0.5*m_p*(R_p^2);
    Freq    = 1000*20/60;
    Cs      = 0.1 + 0.01*sin(2*pi*Freq*t); 
    Ks      = 0.9e3 + 20*sin(2*pi*Freq*t);
    
    % ODEs
    % y1 is theta_a, y2 is theta_p
    % y3 is thetaDot_a, y4 is thetaDot_p
    ydot(1) = y(3);
    ydot(2) = y(4);
    ydot(3) = (- R_a*(Cs*(R_a*y(3) - R_p*y(4)) + Ks*(R_a*y(1) - R_p*y(2))))/I_a;
    ydot(4) = (  R_p*(Cs*(R_a*y(3) - R_p*y(4)) + Ks*(R_a*y(1) - R_p*y(2))))/I_p;
end