1D bar, fixed at one ends, Forced at the other end
1D bar, fixed at one ends, Forced at the other end
(OP)
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][u] = [F]
% =>[m][dv/dt]+[k][u]=[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)
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][u] = [F]
% =>[m][dv/dt]+[k][u]=[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)





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