MATLAB and mechanics (physics mostly) [part II]

Question

Resume form part 1

Basic scheme representing the system:

We assume that coupling of the two shafts is done with a clutch

Equations:

• J1*dw1/dt + Td(w12)+Ts(phi12) = T1;
• J2*dw2/dt - Td(w12) -Ts(phi12) = T2;
where w1 = dphi1/dt, 
    w2 = dphi2/dt, 
    phi12 = phi1 - phi2
    w12 = w1 - w2
    Td(w12) = c12 * w12
    Ts(phi12) = ks * phi12
c12 and ks are some coefficients
• dphi12/dt = w12
• dw12/dt = T1/J1 - T2/J2 - Td(w12)/Jeq - Ts(phi12)/Jeq

ccr = 2*Jeq*wn
wn = sqrt(ks/Jeq)
Jeq = (J1*J2)/(J1+J2)
T1(t) = T0*1(t), T0 = 1 T2(t) = 0
J1+J2 = 10 wn = 100 rad/s c12 = 0 Ts(phi12) = ks*phi12

My second task is to plot how nonlinear clutch affects the rotational inertia.

Given :

T1(t) = T0*1(t), T0 = 1, T2(t) = 0,
J1 = J2 = 5, c12 = 0, wn = 100,
Ts(phi12) = k1*phi12 + k2*(phi12)^3
k20 = 10e14, k2 = [10e-6 10e-5 10e-4 10e-3 10e-2 10e-1 10e0 10e1 10e2 10e3 10e4]

I need to plot the Tsmax/T0 = f(log(k2/k20));

So far my progress is:

MATLAB code

T2 = 0;
J1 = 5;
J2 = J1;
wn = 100;
Jz = J1.*J2/(J1+J2);
ckr = 2.*Jz*wn;
ks = (wn^2)*Jz;
c12 = 0;
k20 = 10e14;
k2 = [10e-6 10e-5 10e-4 10e-3 10e-2 10e-1 10e0 10e1 10e2 10e3 10e4];
Tsmax = zeros(size(k2));
k2a = 0;
plpl = zeros(size(k2));
k1 = 1;
sim ('model_2');

for ii=1:length(k2)
    k2a = k2(ii);
    sim ('model_2');
    Tsmax(ii) = max(Ts);
    plpl(ii) = log(k2a/k20);
end
figure()
    plot(plpl,Tsmax)  
   grid on

which yields

and Simulink model:

Well, I'm kind of sure my plot is not what it should be like.

Answer 1

Some points/suggestions:

• It seems that you parameterise your Simulink model with K2a (unless I can't see very well from the screenshot), but it's not defined anywhere. You define k2a in your code, but MATLAB is case-sensitive so it's not the same as as K2a.
• Your Fcn1 block doesn't do anything, you can get rid of it.

Apart from that, I don't really see anything wrong with what you have done, and we don't know what your plot should look like, so we can't really help you more than that I'm afraid.

EDIT:

I think that's what your plot is supposed to look like. I tried it in Octave, just using MATLAB code & the ode solver (no Simulink) with logspace so that I have 50 data points logarithmically space between 10^(-6) and 10^4, and got something very similar.

I think you are hitting numerical stiffness issues with larger values of k2 so I would suggest using ode15s and a max step size of 1e-4. This may take a bit longer to run, but you might avoid numerical issues when the system is numerically stiff.