2
$\begingroup$

I am attempting to model the following SDOF system with a variable spring and having a sinusoidal input, $$m\ddot{x}+(k_0+k_1)x = m(2\pi f)^2\sin(2\pi ft)$$ where $m$ is the mass, $k_0$ the original spring and $k_1 = k_0$ for $x\dot{x}>0$ and $k_1 = 0$ otherwise.

I was able to code this in Matlab, but for low frequencies, say $f=0.05$, a chattering effect would happen in the plot of $x\dot{x}$ vs. $t$. The following is my code and plots of the result for low frequencies. 

f = 0.05;
t_f = 10;
tspan = [0 t_f];
z0 = [0 ; 0];

options = odeset( 'RelTol', 1e-13, 'AbsTol', 1e-13, 'Stats','on' );


[t, z] = ode23s(@(t,z) sys_13(t, z, f), tspan, z0, options);

figure(1)
set(gcf,'units','normalized','outerposition',[0 0 1 1])

subplot(2,1,1);
plot(t,z(:,1),'b');
title(['Time history of $y(t)\,\,(\ddot{u}_{gx}=-(2\pi f)^2\sin(2\pi f t)) \,\, f=\, $' num2str(f) ' $ $'],'interpreter','latex','fontsize',20)
ylabel('$y$','interpreter','latex','fontsize',15)
xlabel('$t$','interpreter','latex','fontsize',15)
xlim(tspan)

grid on
hold on

subplot(2,1,2);
plot(t,z(:,2),'b')
ylabel('$\dot{y}$','interpreter','latex','fontsize',15)
xlabel('$t$','interpreter','latex','fontsize',15)
xlim(tspan)

grid on
hold on

figure(2)
plot(t,z(:,1).*z(:,2),'b')

grid on
grid minor
hold on

k0 = 4;
k1 = 4;
ktot = k0+k1
k = calcKSmooth(z(:,1).*z(:,2), k0, ktot )*1e-5;

plot(t,k)


function dz = sys_13(t, z , f)
dz = zeros(2,1);

m = 1;
k0 = 4;
k1 = 4
% tol = 1e-6;
% if abs(z(1)*z(2)) > tol && z(1)*z(2) > 0
%     k1 = 4;
% else
%     k1 = 0;
% end

% if z(1)*z(2) > 0
%    k1 = 4;
% else
%    k1 = 0;
% end
ktot = k0 + k1;
k = calcKSmooth(z(1)*z(2) , k0 , ktot);

finput = ( -(2*pi*f)^2 )*sin( 2*pi*f*t );

dz(1) = z(2);
% dz(2) = - finput - ((k0 + k1)/m)*z(1);
dz(2) = - finput - (( k )/m)*z(1);
end

function k = calcKSmooth(xxdot, k0, ktot)
  k1 = k0; k2 = ktot; c = 0;
  r = 1e6;
  ecr = exp(c*r);
  erx = exp(r*xxdot);
  k = (k1*ecr + k2*erx)./(ecr+erx);
end

The plot of the time history of the position and velocity are enter image description here

The plot of $x\dot{x}$ vs. $t$ is enter image description here

If I zoom in close to where the plot is zero I see the following behavior enter image description here

I also tried applying a tolerance so to avoid the constantly switching on and off and instead the chattering occurs about the tolerance I set.

Now, I'm not sure what is causing these issues, but my best guess would be that the discontinuity in the state-space matrix $$ A=\begin{bmatrix}0&&1\\-\frac{k_0+k_1}{m}&&0\end{bmatrix}$$ when $k_1$ switches from $0$ to $k_0$ or vice versa is causing some sort of numerical errors in the solver. I was thinking about using an event location to find when $x\dot{x}$ passes through zero and start/stop the solver to avoid the discontinuity.

Thank you for any help or advice!

EDIT: Using the answer provided by Bill Greene, the chattering was removed, but the way stiffness is added to the system is gradual, yet physically the switching off or on should be more sudden, much like a step function.

The plot of $x\dot{x}$ and with $k=k_0+k_1$ superimposed is:

enter image description here

Improve this question
$\endgroup$
4
  • $\begingroup$ I don't know how to eliminate the chatter caused by the discontinuity in the k-coefficient. However, it is clear that $\partial f/\partial x$ is not defined at the discontinuity. So using an implicit method like ode23s is probably not a good choice because it will try to compute this Jacobian when solving the equations. An explicit solver (e.g. ode23) with a tight tolerance will likely produce a solution with less chatter. $\endgroup$
    – Bill Greene
    Commented Jan 19, 2017 at 18:44
  • $\begingroup$ Thank you for your response. I'll try what you suggested and edit the problem accordingly. What if I were to stop the solver (using an event function) when $x\dot{x}=0$, provide a new stiffness $k_1$ and continue. Would this work around the discontinuity in the k coefficient? $\endgroup$
    – S.D.
    Commented Jan 19, 2017 at 20:47
  • $\begingroup$ No, I don't believe that would be a good strategy. However, what I do believe might work , but forgot to mention in my comment (sorry) is to "smooth" the discontinuity around $x\dot x=0$ so that $f$ is continuous about that point. I think the transition could be quite steep as long as it is differentiable. If you search for smoothing of a step function, you can find suitable functions. $\endgroup$
    – Bill Greene
    Commented Jan 19, 2017 at 21:00
  • $\begingroup$ That makes sense. Something like a hard sigmoid function should work? Or a very steep, but smooth, sigmoid function? $\endgroup$
    – S.D.
    Commented Jan 19, 2017 at 21:06

1 Answer 1

Reset to default
2
$\begingroup$

This a followup to my comment about smoothing the expression for $k$ that is a function of $x$ and $\dot x$.

I converted your expression for $k$ into the following function to make experimentation easier. The input variable xxdot is just $x\dot x$

function k=calcK(xxdot)
  k = 4;
  if(xxdot>0)
   k = k + 4;
  end
end

Then I wrote a second version of this function that closely approximates the function above but is continuous and differentiable.

function k=calcKSmooth(xxdot)
  k1 = 4; k2 = 8; c = 0;
  r = 1e6;
  ecr = exp(c*r);
  erx = exp(r*xxdot);
  k = (k1*ecr + k2*erx)./(ecr+erx);
end

If you plot these two functions in the vicinity of $x\dot x=0$ they are almost indistinguishable. Using calcKSmooth in the ode calculation results in a solution that is very close to the one with calcK but doesn't exhibit the chattering. In addition, the execution time is much less because the algorithm doesn't need to cut the step size so drastically to maintain accuracy.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you for your answer. I used your function 'k=calcKSmooth(xxdot)', and yes the chattering is removed. Yet, when I plot $k$ vs. $t$, in the vicinity of where chattering occured, the graph of $k$ makes sudden changes or gradual decreases in stiffness. I went ahead and updated my question by including your code and displaying the plot of $x\dot{x}$ and the stiffness $k=k_0+k_1$. $\endgroup$
    – S.D.
    Commented Jan 22, 2017 at 11:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.