Simulate shock with spring ODE

Hello everyone !

I was given the project to compute and simulate a shock between two particles, using a spring differential equation. It's a 2 dimension problem, so there are eight equations to solve:

Where and are the initial velocities of the first particle (the second one is not moving at the beginning) and are constants. But when I run my code:

SystemInit=[0,0,30,30,2,2,0,0];
time=[0,100];
[t,y]=ode45(@(t,y) odefunc(t,y,a1,a2,K1x,K1y,K2x,K2y),time,SystemInit);

With my odefunc being:

function u=odefunc(t,y,a1,a2,K1x,K1y,K2x,K2y)
u=zeros(8,1);
u(1)=y(1);
u(2)=y(2);
u(3)=y(3);
u(4)=y(4);
u(5)=a1*y(5)+K1x;
u(6)=a1*y(6)+K1y;
u(7)=a2*y(7)+K2x;
u(8)=a2*y(8)+K2y;
end

I get a sort of exponential law for all the components of my y vector (excepted for the two firsts). So am I doing it right ? We've just started learning how to compute differential equations in Matlab, so there might be something I didn't manage to figure out.

I hope I've been clear enough.

[darova]

Everything looks correct. One thing:

% your case: u = y(2), du = y(1)
u(1)=y(1);
u(5)=a1*y(5)+K1x;
% usually: u = y(1), du = y(2)
u(1)=y(5);
u(5)=a1*y(1)+K1x;

You know what i mean?

[Ameer Hamza]

Actually, it is necessary to use the order as suggested by @darova, otherwise, you will get the wrong answer (as you mentioned exponential curve in your case).

Also, is the 2nd last equation correct? It again mentions .

[darova]

Indeed. Thanks Ameer Hamza

[Eugène PLANTEUR]

Yep, the 2nd last equation was wrong, just a typo, sorry !

I've tried switching the order as you recommanded but still got a wrong result :/

[Ameer Hamza]

If you can provide the values of a1,a2,K1x,K1y,K2x,K2y, then we can try to run the code and find the issue.

[darova]

Try to reduce time

time=[0,10];

[Eugène PLANTEUR]

Well, one of the final goals of the project is to study those constants (that depend on other constants) would change the phenomenon, and I've been told to first manage to solve the equations with matlab (so using any value for the constants) and then worry about them. But there one major problem, that I thought I'd solve after managing to set up my main code:

There are basically 3 steps in the problem. At the beginning, the first particle moves towards (or not) the second partcile. If, at any given time, the distance between them is lesser or equal to the sum of the radius of the two particles (i.e. they collide), then we have to solve the differential equations at that point in time. But, each K constant depends on the value of the position. For example:

Where k is the stifness constant of the spring, the mass of the first particle, the equilibrium length of the spring and d the distance bewtween the two particles. So I don't either know how I could compute that :/

[Ameer Hamza]

It appear from your equation tha is a function of , so you cannot ignore that.

[Ameer Hamza]

The values of parameters are also important. For example, if I run, (negative values of and ).

a1=-1;
a2=-1;
K1x=1;
K1y=1;
K2x=1;
K2y=1;

SystemInit=[0,0,30,30,2,2,0,0];
time=[0,100];
[t,y]=ode45(@(t,y) odefunc(t,y,a1,a2,K1x,K1y,K2x,K2y),time,SystemInit);

I don't get exponential curve.

The reason will become clear if you study analytical solution of 2nd order ODEs.

Remember:

function u=odefunc(t,y,a1,a2,K1x,K1y,K2x,K2y)
u=zeros(8,1);
u(1)=y(5);
u(2)=y(6);
u(3)=y(7);
u(4)=y(8);
u(5)=a1*y(1)+K1x;
u(6)=a1*y(2)+K1y;
u(7)=a2*y(3)+K2x;
u(8)=a2*y(4)+K2y;
end

[Eugène PLANTEUR]

That's what I thought at the geginning but our teacher just told us that it wasn't important at all and that we shouldn't care until we got the equations running, that's why I'm struggling, but thanks !

[Ameer Hamza]

Even in that case, you need to include the effect of , for example I changed the odefunc as follow

function u=odefunc(t,y,a1,a2,K1x,K1y,K2x,K2y)
u=zeros(8,1);
u(1)=y(5);
u(2)=y(6);
u(3)=y(7);
u(4)=y(8);
u(5)=a1*y(1)-K1x*y(3);
u(6)=a1*y(2)-K1y*y(4);
u(7)=a2*y(3)-K2x*y(1);
u(8)=a2*y(4)-K2y*y(2);
end

I lumped the effect of into K1x. I made the guess about the values of , and . In that case, the system can converge, even with the positive values of a1 and a2

a1=0.5;
a2=0.5;
K1x=1;
K1y=1;
K2x=1;
K2y=1;
SystemInit=[10,0,10,0,1,2,1,2];
time=[0,100];
[t,y]=ode45(@(t,y) odefunc(t,y,a1,a2,K1x,K1y,K2x,K2y),time,SystemInit);
plot(t,y)

Although it still diverges for some value of initial conditions, i guess that is actually the property of this system of differential equations.

[Eugène PLANTEUR]

Thanks ! But unfortunately I've tried re-implementing what could be some realistic values and I just got the same curve again (a1 and a2 would be around 15, and every K constant really depends on the position, so I don't know). I'll just ask my teacher tomorrow, hoping for more help than he gave. Thanks a lot anyway for trying to help me ! I really appreciate it !

[Ameer Hamza]

Glad to help.

Answer 1

[Eugène PLANTEUR]

So, I managed yo get something right ! Had to rearrange the code a bit, rewrite the equations in a much simpler form (I over-complicated myself), got some help from our teacher and here it is:

m1=1;m2=1;                 %Particles' masses
radius1=0.1;radius2=0.1;    %Particles' radius
k=10000;                    %Spring's stifness parameter
L0=rayon1+rayon2;           %Spring's equilibrium's length


SystemInit=[0,0,        %Initial position of the first particle
            1,0,        %Initial position of the second particle
            1,0.1,      %Initial speed of the first particle
            0,0];       %Initial speed of the second particle
time=[0,10];
[t,y]=ode45(@(t,y) odefonc(t,y,m1,m2,k,L0),time,SystemInit);

plot(y(:,1),y(:,2),'or',y(:,3),y(:,4),'db')
axis equal

And here's the odefunction:

function u=odefonc(t,y,m1,m2,k,L0)
u=zeros(size(y));
u(1)=y(5);              %Velocity of the first particle along the X axis
u(2)=y(6);              %Velocity of the first particle along the Y axis
u(3)=y(7);              %Velocity of the second particle along the X axis
u(4)=y(8);              %Velocity of the second particle along the Y axis

d=sqrt((y(3)-y(1))^2+((y(4)-y(2))^2)); %Distance between the two particles

a1=(k/d*m1)*(L0-d);
a2=(k/d*m2)*(L0-d);

if d<L0     %i.e if there is contact 
u(5)=a1*(y(1)-y(3));    %Acceleration of the first particle along the X axis
u(6)=a1*(y(2)-y(4));    %Acceleration of the first particle along the Y axis
u(7)=a2*(y(3)-y(1));    %Acceleration of the second particle along the X axis
u(8)=a2*(y(4)-y(2));    %Acceleration of the second particle along the Y axis
end
end

What I'm plotting is the trajectory of both particles.

Anyway, thanks to you all, you've helped a lot !