1D bar, fixed at one ends, Forced at the other end

[mina1363]

Hi All,
I am having difficulty understanding a set of results I have obtained by solving the dynamic equations of a 1D bar, fixed on the left hand side, and excited by an external force of F0*sin(w*t)on the right hand side. I have integrated in time using an explicit Euler scheme. I suspect for a simple case, that should not cause instabilities.
I have also tried to work out the natural frequencies of the forced system analytically.
I was expecting the results to be similar to a free vibration results, when you have the amplitude waves that either damp over time (provided frequency not close to the natural frequency). But here I get waves that are not based around a straight line at (0,0) but the fluctuate around a line that has a gradient. Not sure if I am clear here, so I have attached my code which produces the results. The velocity results seem more reasonable.
I would appreciate it if you could run my code and guide me through this. I have spent quite a bit of time on it. I want to understand this fully. Are there ways to check my answers? How can I use the eigenfrequency and eigenvectors to check my results?
I would appreciate any guidance ASAP. Thank you very much.

clear all
clc

%%
% MATERIAL PROPERTIES
nn = 4;
E = 30e6;
A =1 ;
L0 = 100*(nn-1);
le0 = (L0/nn);
irhos0 = 0.00073*386.22;

k1 = [1 -1 ; -1 1]; k2 = k1; k3 = k1;
m1= [2 1;1 2]; m2 = m1; m3 = m1;
K = ((E*A)/le0)*[k1(1,1) k1(1,2) 0 0;
k1(2,1) k1(2,2)+k2(1,1) k2(1,2) 0;
0 k2(2,1) k2(2,2)+k3(1,1) k3(1,2);
0 0 k3(2,1) k3(2,2)];
M = ((irhos0*A*le0)/6)*[m1(1,1) m1(1,2) 0 0 ;
m1(2,1) m1(2,2)+m2(1,1) m2(1,2) 0;
0 m2(2,1) m2(2,2)+m3(1,1) m3(1,2);
0 0 m3(2,1) m3(2,2)];

F0 = 1000;
%%
% ANALYTICAL, CALCULATES THE NATURAL FREQUENCIES OF THE FORCED SYSTEM
% Let x(t) = Xsin(wt), xdot = wXcos(wt), xddot = -(w^2)X(sin(xt)
% Let f(t) = Fsin(wt)
% sub above in [M]{xddot}+[k]{x}={f}
% => ([k]-[m](w^2)){X}={F}
% => {X} = inv[[k]-[m](w^2)]*{F}
% => CHOOSE F function, here F = F0*w^2

% EIGENVALUES
X = zeros(nn,1); % AMPLITUDE OF VIBRATIONS
b = 1;
w1 = 0:3*180;
freq = w1/(2*pi);
Fex1 = 0; % EXTERNAL FORCE
F1 = zeros(nn,1);
for j=1:length(w1)
Fex1(j) = F0*(w1(j)^2);
F1,j) = [zeros(nn-1,1);Fex1(j)];
A = K-M*(w1(j)^2);
X,j) = inv(A)*F1,j);
end
figure
subplot(4,1,1)
plot(w1,abs(X(1,),'-r')
subplot(4,1,2)
plot(w1,abs(X(2,),'-r')
subplot(4,1,3)
plot(w1,abs(X(3,),'-r')
subplot(4,1,4)
plot(w1,abs(X(4,),'-r')
xlabel('\omega rads/sec')
ylabel('Displacement')
% close all


% Eigenvectors
% ??

%%
% FEA - EXPLICIT EULER TIME SCHEME
m = 5000; % NUMBER OF TIME STEPS
wn = sqrt(E/irhos0); %[E*A/L/rho]=
w = wn/100;
tp = (2*pi)/w;
tf = 2*tp;
dt = tf/m;

% [m][d2u/dt2]+[k] = [F]
% =>[m][dv/dt]+[k]=[F]
% => dv/dt = inv[m]*([F]-[k]{u}) EQ1
% du/dt = v EQ2
u = zeros(nn,1);
v = zeros(nn,1);
u,1) = 0;
v,1) = 0;
u(1, = 0;
v(1, = 0;
Fex2 = 0;
F2,1) = zeros(nn,1);
t = 0;
ww = 0;
for i = 1:m
t(i+1) = t(i)+dt;
Fex(i+1) = F0*sin(t(i)*w);
F2,i) = [zeros(nn-1,1);Fex(i)];
v,i+1)= v,i)+dt*inv(M)*(F2,i)-K*u,i));
u,i+1)= u,i)+dt*v,i);
t(i+1) = t(i)+dt;
end

% PLOT DISPLACEMENT WITH TIME
figure
subplot(4,1,1)
plot(t,(u(1:nn,))
subplot(4,1,2)
plot(t,(u(2,))
subplot(4,1,3)
plot(t,(u(3,))
subplot(4,1,4)
plot(t,(u(4,))

% PLOT OF VELOCITY WITH TIME
figure
subplot(4,1,1)
plot(t,(v(1,))
subplot(4,1,2)
plot(t,(v(2,))
subplot(4,1,3)
plot(t,(v(3,))
subplot(4,1,4)
plot(t,(v(4,))

% PLOT OF FORCE WITH TIME
figure
plot(t(1:end-1),Fex)

 

[mina1363]

I just found out that I have not really implemented the boundary condition in the code. Have done that now and results seem more reasonable. Having said that, I would appreciate some guidance on how I can evaluate the results, and how to work out the eigenvectors of a forced system.