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;