Equation of Motion is Blowing Up

Hello all, I am trying to add a torsional spring to my equation of motion for a pendulum, which is the last column of 'dy' seen below. I only want the spring to kick in if the angle is greater than 70 degrees on either side. This ends up blowing up though and I'm not sure why. You can see in the picture below the damping that occurs past -70 degrees as expected but then when the pendulum reaches +70 degrees it immediately blows up and continues increasing indefinitely. Is there something wrong with the logic of the code or possibly the way I implemented the equation for a torsional spring? Any help would be greatly appreciated.

Note: y(3) is the angle in radians

if abs(y(3)) > 70*pi/180 %converted to rads
    k = 70; %spring constant 
else
    k = 0;
end

dy = [y(2)-c0*y(1);
      -k0*y(1);
      y(4);
      -alpha1*sin(y(3))-ee*gamma1*l*y(4)+sqrt(ee)*lambda*sin(y(3))*A*sin(omega_ex*t) + k*(abs(y(3))-70*pi/180)];

[darova]

Can you show something more? Equations? Pictures?

Answer 1

[David Goodmanson]

HI Samuel,

I won't have access to Matlab for a few days so this needs to be verified, but I believe that in the last line of the code for dy, the last term should not be

k*(abs(y(3))-70*pi/180)

but rather

-k*y(3).

That is, the 'if check' does not belong in the kinematics itself.

This assumes a constant-k spring, as opposed to something more complicated such as a spring constant that increases the more you exceed +-70 degrees.

Is it true that the spring constant k = 70 is coincidentally the same as the turn-on-the spring angle of +-70 degrees?

[RPS19]

Hi David, thanks for the help. Yes it is just a coincidence that they both equal 70. The reason that I included that bit in the last term of dy is to account for the delta theta. (The 70 in the dy is for the angle just to be clear). So it would be a constant k (coincidentally also 70) multiplied by the change in angle.

[RPS19]

I was able to fix it by changing that portion of the code to

- k*((y(3)/abs(y(3))) * (abs(y(3))-70*pi/180))

that way it accounts for the restorative factor of a spring with the first negative, like you mentioned, and normalizing the vector with y(3)/abs(y(3)) retains the +/- of the input angle. Thanks a lot for your help.

Answer 2

[Jim Riggs]

It's seems pretty clear from the plot of angle vs. speed that the system is gaining energy. You really need to draw a picture of the system in order to clearly define the angles and sign conventions of all of the terms. We can't tell how any of your terms are defined. It is likely that that there is a wrong sign on a term and/or your integration method is causing a numerical error.

[RPS19]

Hi Jim, it did turn out to be a matter of sign convention, thanks for the help.