hi
I did this code to resolve the vibration beam using FEM by newmark integration but it didn't work. Can anyone tell me whats wrong with my code or recommend a matlab code?
Thank you.
NEWMARK TIME INTEGRATION
% Here for beta = 0.5 and alpha = 1/6, the following values are obtained. The MATLAB program is shown below:
alpha=-1/3;
beta=(1-alpha^2)/4;
Cdamp=alpha*Mtotal+beta*Ktotal; % Damping matrix; Mtotal:gmobal matrix; Ktotal :stifness matrix
F = zeros(sdof-4,1); % initialisation of force sdof (system degre of freedom)
F(sdof-6,1)=1000 ;
dt=0.005; T=2;
X0=zeros(sdof-4,1);; X0d=zeros(sdof-4,1);
X0dd=inv(Mtotal)*(F-Cdamp*X0d-Ktotal*X0);
beta=0.5;gamma=1/6;%0.25*(0.5+beta);
a0=1/(beta*dt^2);a1=gamma/(beta*dt);a2=1/(beta*dt);
a3=(1/2*beta)-1;a4=(gamma/beta-1);a5=0.5*(gamma/beta-2)*dt;
a6=dt*(1-beta);a7=beta*dt;
Kb=Ktotal+a0*Mtotal+a1*Cdamp;
i=1;
X(:,1)=X0;Xd(:,1)=X0d;Xdd(:,1)=X0dd;t=0;
fprintf('time(s)\t\tX1\t\tX2\n');
fprintf('%f\t%f\t%f\n',t,X(1,1),X(2,1));
for t=dt:dt:T
i=i+1;
F(sdof-6,1)=1000 ;
Fb=F+Mtotal*(a0*X(:,i-1)+a2*Xd(:,i-1)+a3*Xdd(:,i-1))+Cdamp*(a1*X(:,i-1)+a4*Xd(:,i-1)+a5*Xdd(:,i-1));
X(:,i)=inv(Kb)*Fb;
Xdd(:,i)=a0*(X(:,i)-X(:,i-1))-a2*Xd(:,i-1)-a3*Xdd(:,i-1);
Xd(:,i) =a1*(X(:,i)-X(:,i-1))-a4*Xd(:,i-1)-a5*Xdd(:,i-1);
end
figure(2)
t=[0:dt:T];
plot(t,X(sdof-4,:),'-pr')
xlabel('time(s)');
ylabel('displacement(m)');
legend('DOF–sdof-6');
grid on;