Asymmetric spring with two different stiffnesses

I'm trying to model a mass spring damper system and I want to describe the spring with two different values of k in bump and in rebound. I wrote the system made by the free body equations and what now what I want to do is to solve it with ode45, but I cannot find the way to compare y(t) and y(t-1) and to understand what stifness I have to use. I tried to define a vector in the function I wrote to solve the differential system, and the idea was to write in the i position the current value of the displacement but I don't know the lenght of the vector and at the beginning of every new iteration it goes zeros(n,1). Is there any solution to save and compare different values from different iterations in this case?

[Alan Stevens]

Upload what you have done so far.

[Paolo Alloisio]

function f = solve(t, y)
    load fileprova a
    m = 5;
    c = 1;
    k1 = 10;
    k2 = 1;
    f(1) = y(2);
    if abs(a)>=abs(y(2));
       f(2) = (-c*y(2) - k1*y(1))/m;
    else
       f(2) = (-c*y(2) - k2*y(1))/m;
    end
    f = [f(1); f(2)];
    a = f(2);
    save fileprova a
end

This is the code I wrote to solve the system made by differential equations. I'm sorry if it's not perfect, but I'm learning to use MATLAB by myself and sometimes I'm not able to see the problem from the right perspective. This morning I thought that would be smart to save y(2), which is the value I want to use to decide which k is the right one, in a variable which is external to the function, but I run the code and I don't know why a is always equal to 2 and not to y(2).

[darova]

Maybe you are looking for persistent

[Paolo Alloisio]

I didn't know the existence of this feature, so I really thank you. I'm having some troubles in comparing the right variables in the function used in the solution of the differential equations, but I hope I could find a way to model what I want to simulate.

[darova]

remember about clearing persistent variable before each run

clear functions     % clear persistent variables from functions

Answer 1

[Bjorn Gustavsson]

Yeah, I think that should be the correct test. In the modeling you do you have rather big difference in spring constant. When I test with this simplified version:

f_assyspring = @(t,y) [y(2);-y(1)*(1+1*double(sign(y(1))==sign(y(2))))];
t = linspace(0,50,3001);
y0 = [1,1];
[T2,Y2] = ode45(f_assyspring,t,y0);
subplot(2,1,1)
ph = plot(T2,Y2);                  
axis([0 20 -2 2])
grid on
legend(ph,'y-displacement','y-velocity')
subplot(2,1,2)
plot(T2,gradient(Y2(:,2),T2))
axis([0 20 -2 2])
legend('y-acceleration')
grid on

The solution behaves OK - in a mathematical sense.

However, from a physical point-of-view, this physical model decays rapidly - because it dissipates energy at the turning-points when there is a switch to the smaller spring-constant - at that point the simple Hooksian sping has a lot of stored energy (proportional to |ky^2| IIRC). By switching spring-constant there the difference in stored energy is magiced-away - which feels dodgy to me, for dampers I can see that it would be possible to have different leveles of dissipation in different directions.

[Paolo Alloisio]

I want to thank you once again, because you help me with the physical side of the problem while I was totally concerned with the code. To be completely honest, I didn't think about the energetical side of the matter, and so I couldn't understand where the problem was. But now that you help me to see things from the right point of view, I realize that my idea of starting to work with a spring instead of a damper was not so good. So, thank you so much Mr Gustavvson, your help was very precious to me. I'll keep in mind everything you told me when I'll be about to model the damper! But now I'm gonna start to work with a non linear spring, to understand the consequences of the hardening on the system.

Answer 2

[Bjorn Gustavsson]

ODE45 takes a function, f(t,y), describing the differential equation.

You have a differential equation something like this

where k has one value if the spring moves in one direction (dy/dy < 0) and another value when the spring moves in the other direction - how you connect those forces at the points where the velocity is zero I do not know.

ODE45 takes functions describing first-order odes, so the way to implement equations of motions is to convert the second-order ODE to 2 coupled first-order odes. That way we get something like this:

function dydtdvdt = your_ode(t,y)
k_p = 12; % Or whatever your spring-constants over mass are...
k_n = 10;

if y(2) > 0
  a = -k_p * y(1);
else
  a = -k_n * y(1);
end
dydtdvdt = [y(2);
            a];
end

Then you call ODE45 with this function the time-period of interest and the initial conditions (y(0) and v(0)), and then you get the numerical solution to your ODE. If you have a damper you'll have to add that force-term too.

HTH

[Paolo Alloisio]

The problem in zero displacement could be solved using an >= opeator, in that point the force of the spring is always zero. I have to say thanks to you, because this seems to be the best solution for my problem by now, the only thing that I have to solve is the condition

if y(2) > 0

because writing this implies that the software controls the velocity of the system, while I want a displacement control. I'm not sure I'm right, but running the script this is the plot I have

And I don't think this is correct, because the slope of the curve should change while it's shortening or stretching.

[Bjorn Gustavsson]

OK, then I misinterpreted your description. I read your description as if you had one spring-constant when the velocity, y(2), was in different directions. If you want the spring-constant different depending on whether the spring is extended or compressed you have to make the test for that - i.e. y(1).

[Paolo Alloisio]

And that's what I did. I follow your suggestion and I made the test for y(1), so the code was

function f = solve(t,x)
m = 1;
k1 = 10;
k2 = 1;
f(1) = x(2);
if x(1)>=0
    f(2) = -k1*x(1)/m;
else
    f(2) = -k2*x(1)/m;
end
f = [f(1); f(2)];
end

Then another problem arised, because this code describes a spring which is characterized by two different stifnesses, but f(2) it's influenced by the sign of x(1), and not by the fact that the spring is stretching or shortening. What I want to do is to compare x(1) from a iteration with the one from the previous one. So, if the difference is positive in absolute value I'm sure the spring is stretching even if displacement is negative and viceversa. Following the advice given by darova, I defined a variable called past using the feature

persistent

, and that's what I got

function f = solve(t,x)
persistent past
m = 1;
k1 = 10;
k2 = 1;
f(1) = x(2);
if abs(past) - abs(x(1))>=0
    f(2) = -k1*x(1)/m;
else
    f(2) = -k2*x(1)/m;
end
f = [f(1); f(2)];
past = x(1);
end

But now the system seems to diverge, and I couldn't really understand why

[Bjorn Gustavsson]

But then you want my initial solution. If the velocity is positive then the spring is stretching (possibly from a compressed state, but if the velocity is positive you know that the spring-length will increase), and the other way around.

You cannot use persistent here since there is no guarantee that ODE45 will only ever evaluate your ode-function in increasing order with constant step-size. Typically it might have to refine its step-size around points where the trajectory is "interesting".

[Paolo Alloisio]

Since I'm considering the system to have an harmonic response, I didn't consider the velocity as the right parameter to decide if the spring lenght is increasing or decreasing. And that's because the velocity has two different signs when the mass linked to the spring is oscillating from equilibrium position to the maximum positive distance from it( in this case the velocity is positive) and when it's getting away in the negative direction (while in this situation is negative). And in both cases the spring lenght is increasing, and so if I choose the velocity as the control parameter I'm going to describe the spring with two different stifnesses even if the behaviour would be the same. Forgive me if I'm not explaining too well, but I'm not so good with this software and I'm learning to use it by making some attempts. So are you telling me there's no way to compare two values from two different iteration of the algorithm ODE45?

[Bjorn Gustavsson]

There's two points you have to make clear. The first is what your force is supposed to be as a function of and . Write it clearly as equations, then you can think about how to implement them. The second is that you have to familiarise yourself with ODE-integrating methods, for example Runge-Kutta (both the basics and the adaptive section).

Small point to make is that you will only get "a harmonic response" if your acceleration term is of the form:

with k constant if you start to modify the constant with respect to sign of y or you will have a more complex response.

[Paolo Alloisio]

Yes, probably you are right. It's possible that I thought how to write the code to solve the problem without thinking enough about the problem itself. I'm working on a thesis about the vertical dynamics of the vehicle and I'm learning to use MATLAB by trying and making mistakes. I started my study writing some code to analize the classic quarter car model; in that case I supposed all the springs and the damper to have linear behaviour, and so I solved the problem with the superposition of modal vector while there was no force applied at the mass, and in the space state. But my objective is now to define a more complex model and so I started to think how to describe a spring which have two different stifness in stretching and shortening, while it's moving around the equilibrium position. The system I asked about is the more simple I could think about and is made by a mass linked to the ground by a spring. I choose this example because my aim is to understand how to use MATLAB to solve this simple problem before analyzing the 2 dof system. So the free body equation of the simple sistem I wrote is , which is the same you wrote before. Thinking about what you said, I would say that the stifness k is a function of the position, and I would define as k_stretch if , while it's k_short if . And if I have to describe what my problem is, I would say that I'm not able to verify the previous relationship to understand what stifness is the right one on every iteration. I also have to say that you are right when you say that I have to better understand the Runge-Kutta algorithms, and so I want to ask you if the Wikipedia page you linked me could help me to understand how the MATLAB functions works, or I need to find some other material because the situation is more complex.

[Bjorn Gustavsson]

So that characteristic you should be able to implement as a combination of checks on the velocity, and the sign of y, you get four cases to check for.

However, this model will have an assymetry around the turning-points where changes sign. That is where the |ky| will be at its maximum. How you see this connected to the physics of the springs is up to you, but to me it looks peculiar.

ODE45 uses a 4-5 version of Runge-Kutta, so the Wikipedia-page should be a good introduction. If you need to get into the details the function is implemented in plain .m-code so you can read it, there are also references to publications in the ode45.m file.

[Paolo Alloisio]

I admit I quitted the idea this morning, because you are right when you say this is a peculiar spring. I was thinking about how to produce this kind of elastic element and I really don't know how to do this materially speaking, because I guess the use of a metallic material processed with the common production techniques gives a symmetric element. I had this idea while I was thinking about some structural materials as concrete, which has a different behaviour while subject to tensile and compressive forces, and it exhibits two different elastic modules. What I'm gonna try now is to describe a spring with a non costant stiffness, but which is still symmetric in traction and compression, because if I'm not wrong it would be desirable to have a stiffness which grows with the displacement to avoid the suspension to reach the end of its travel. But I still want to ask you one thing, because you had really helped me so far. I also know that what I have tried to model it's not really common in the field of elastic elements, but it's something often done in the tailoring of the dampers. I read that hydraulic elements which are part of the suspension are dual rate with a normal three to one ratio between jounce and rebound damping, and so I will probably need to describe the damper in the same way I tried to do with the spring. This morning I've been thinking about what you told me about the check of velocity and position, and I wrote this

function f = solve(t,x)
m = 1;
k_stretch = 100;
k_short = 1;
f(1) = x(2);
if x(1)>=0 && x(2)>=0
    f(2) = -k_stretch*x(1)/m;
elseif x(1)>=0 && x(2)<0
    f(2) = -k_short*x(1)/m;
elseif x(1)<0 && x(2)<0
    f(2) = -k_stretch*x(1)/m;
else x(1)<0 && x(2)>=0;
    f(2) = -k_short*x(1)/m;
end
f = [f(1); f(2)];
end

The idea was to choose the right value of k analyzing the sign of the displacement and the velocity. I think that the result will be the same if I write the check condition as

if sign(x(1))==sign(x(2))
    f(2) = -k_stretch*x(1)/m;
else
    f(2) = -k_short*x(1)/m;
end

And running the scritp this is what I get

I plotted the displacement of two springs with constant and to analyse better the behaviour of the element I was trying to model and I'm still non sure of what I got. What I want to ask you is if this seems to have sense to you, because I want to undestand if I could use this way of thinking when I'll have to model an asymmetric damper. My main concern is what happen when I change the two stifness, because if I define and the system still diverges.