CLASSICAL STUDY OF THE NONLINEAR MASS-DAMPER-SPRING STABILITY

In this part, a numerical study of the nonlinear mass-damper-spring has been investigated. In the first step, the spring is taken as a linear one and the damper is omitted. Then, the spring is implemented with its nonlinear characteristics and the damper is introduced.

In order to solve diffrential equations system, a MATLAB program based on the forth order Runge Kutta scheme is written. Results are compared to a preimplemented MATLAB function called ODE45 which directly solves diffrential equations. So, good agreement is found between program solutions and ODE45 ones.

MATLAB  MAIN PROGRAM

close all;
clear all;
echo on

%%%%%%%%%%%%%%%%%%%%%Etude du ressort a raideur polynomiale%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%                                      %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%              K0*x + K1*x^3           %%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%              Première Partie         %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Cette premiere partie traite un ressort de raideur k(x) = k0.x sans et un systeme
% sans amortissement (frottement).
% l'equation est : -k0.x + F(t) = mx". cette equation se traduit par un systeme :
% x' = y et y' = F(t)/m -k0.x/m.

% Un shema inconditionnelement stable = Runge Kutta qui s'ecrit pour l'ordre 4 :
% x(i+1)=x(i) + h/6*(f1+2f2+2f3+f4);
% y(i+1)=y(i) + h/6*(g1+2g2+2g3+g4)

% f1(i) = f(t(i),x(i), y(i));
% f2(i) = f(t(i)+h/2,x(i)+h/2*f1, y(i)+h/2*g1);
% f3(i) = f(t(i)+h/2,x(i)+h/2*f2, y(i)+h/2*g2);
% f4(i) = f(t(i)+h,x(i)+h*f3, y(i)+h*g3));

echo off

%%%%%%%%%%%%%%%%%%%%% Variables %%%%%%%%%%%%%%%%%%%%%%%

global h Tmin Tmax m k0 k1 k2 a b w;

m = input('Donner la masse du corps, m = ');
h = input('Donner le pas de temps du calcul, h = ');
Tmin = input('Donner la borne temporelle inferieure Tmin = ');
Tmax = input('Donner la borne temporelle superieure Tmax = ');
a = input('Donner l"amplitude de la force de teste, a = ');
w = input('Donner la pulsation de la force de teste, w = ');
k0 = input('Donner la raideur linéaire du ressort, k0 = ');

%%%%%%%%%%%%%%%%%%%%% Programme %%%%%%%%%%%%%%%%%%%%%%%

t = Tmin:h:Tmax;
n = length(t) ;
F = a*sin(w*t) ;

x(1) = 0;
y(1) = 0;
k1 = 0;
k2 = 0;
b = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Solution avec une fonction integree de matlab : ode45
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[t,Y] = ode45('odefile',t,[0 0]);
plot(t,Y(:,1),'r');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Methode de Runge Kutta%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

k1 = 0;
k2 = 0;
b = 0;

for i = 1:n-1,
f1 = f(t(i),x(i), y(i));
g1 = g(t(i),x(i), y(i));
f2 = f(t(i)+h/2,x(i)+h/2*f1, y(i)+h/2*g1);
g2 = g(t(i)+h/2,x(i)+h/2*f1, y(i)+h/2*g1);
f3 = f(t(i)+h/2,x(i)+h/2*f2, y(i)+h/2*g2);
g3 = g(t(i)+h/2,x(i)+h/2*f2, y(i)+h/2*g2);
f4 = f(t(i)+h,x(i)+h*f3, y(i)+h*g3);
g4 = g(t(i)+h,x(i)+h*f3, y(i)+h*g3);

x(i+1) = x(i) + h/6*(f1 + 2*f2 + 2*f3 + f4);
y(i+1) = y(i) + h/6*(g1 + 2*g2 + 2*g3 + g4);
end,

subplot(2,1,1)
title('position x en fonction du temps')
plot(t,x,'r');

subplot(2,1,2)
title('vitesse x' en fonction du temps')
plot(t,y,'b');

pause;

echo on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%             Deuxième Partie          %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Cette deuxieme partie traite un ressort de raideur k(x) = k0.x + k1.x^3 avec
% amortissement.
% l'equation est : -k0.x + k1.x^3 + F(t) = mx". cette equation se traduit par le
% système : x' = y et y' = F(t)/m -k0.x/m - k1/m.x^3.

echo off

k1 = input('Donner la raideur K1 du ressort d"ordre 3, k1 = ' );
b  = input('Donner la valeur de l"amortissement, b = ');

for i = 1:n-1,
f1 = f(t(i),x(i), y(i));
g1 = g(t(i),x(i), y(i));
f2 = f(t(i)+h/2,x(i)+h/2*f1, y(i)+h/2*g1);
g2 = g(t(i)+h/2,x(i)+h/2*f1, y(i)+h/2*g1);
f3 = f(t(i)+h/2,x(i)+h/2*f2, y(i)+h/2*g2);
g3 = g(t(i)+h/2,x(i)+h/2*f2, y(i)+h/2*g2);
f4 = f(t(i)+h,x(i)+h*f3, y(i)+h*g3);
g4 = g(t(i)+h,x(i)+h*f3, y(i)+h*g3);

x(i+1) = x(i) + h/6*(f1 + 2*f2 + 2*f3 + f4);
y(i+1) = y(i) + h/6*(g1 + 2*g2 + 2*g3 + g4);

end,

subplot(2,1,1)
plot(t,x)
title('position en fonction du temps')

subplot(2,1,2)
title('vitesse en fonction du temps')
plot(t,y,'b');

end

MATLAB SUBPROGRAMS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%               f function             %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function s = f(t,x,y)

s = y;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%               g function             %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

functiond = g(t,x,y)

globala k0 k1 b m w;

d = a*sin(w*t)/m - k0*x/m - k1/m*x^3 - b/m*y*abs(y);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%           odefile function           %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ode = odefile(t,Y),

global k0 k1 m a w b k2;

ode = [Y(2) ; a/m*sin(w*t)- k0/m*Y(1)- k1/m*(Y(1))^3 - b/m*Y(1)*abs(Y(1))];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%       Phase portrait computing       %%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all;
clear all;

global h Tmin Tmax m k0 k1 a b w;

m    = input('Donner la masse du corps, m = ');
h    = input('Donner le pas de calcul, h = ');
Tmin = input('Donner la borne temporelle inferieure Tmin = ');
Tmax = input('Donner la borne temporelle superieure Tmax = ');
a    = input('Donner l"amplitude de la force de teste, a = ');
w    = input('Donner la pulsation de la force de teste, w = ');
k0   = input('Donner la raideur K0 du ressort, k0 = ');
k1   = input('Donner la raideur K1 du ressort, k1 = ');
b    = input('Donner la valeur de l"amortissement, b = ');

%%%%%%%%%%%%%%%%%%%%%%%%             programme            %%%%%%%%%%%%%%%%%%%%%%%%

t= Tmin:h:Tmax;
n = length(t);
F = a*sin(w*t);

x(1) = 0;
y(1) = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%  Solution avec la fonction integrée de matlab : ode45 %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[t,Y] = ode45('odefile',t,[0 0]);
plot(Y(:,1),Y(:,2),'r');

RESULTS

Corrolation between Program resultes and preimplemented MATLAB function (ODE45) is controlled by plotting masse center gravity postion (x) versus time (t) whith both methods on the same graph. Curves are found to be the same.

When nonlinear characteristic of the spring (k1) is introduced, Frequency Instability appeares by changing k1 value.

k1= 0 to k1 = 1 by step of  0.1