MATLAB Answers

Asymmetric spring with two different stiffnesses

조회 수: 16(최근 30일)
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?

  댓글 수: 5

표시 이전 댓글 수: 2
darova
darova 6 Jul 2020
Maybe you are looking for persistent
Paolo Alloisio
Paolo Alloisio 6 Jul 2020
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
darova 6 Jul 2020
remember about clearing persistent variable before each run
clear functions % clear persistent variables from functions

Sign in to comment.

채택된 답변

Bjorn Gustavsson
Bjorn Gustavsson 8 Jul 2020
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.

  댓글 수: 1

Paolo Alloisio
Paolo Alloisio 8 Jul 2020
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.

Sign in to comment.

추가 답변(1개)

Bjorn Gustavsson
Bjorn Gustavsson 6 Jul 2020
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

  댓글 수: 9

표시 이전 댓글 수: 6
Paolo Alloisio
Paolo Alloisio 7 Jul 2020
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
Bjorn Gustavsson 8 Jul 2020
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
Paolo Alloisio 8 Jul 2020
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.

Sign in to comment.


Translated by