Hello everyone i am solving multi degree freedom of system by using Newmarks beta method in which loading function F(x,v). When i am ruing the code, the result i
17 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
% dynamic analysis using direct integration method
% % the porbleme is: Mx''+Cx'+Kx=(C1x'+K1x) *eps*cos(Omeg*t)+K2x* eps^2cos(Omeg*t)^2
clc;
clear all
format short e
%close all
m=[18438.6 0 0;0 13761 0;0 0 9174];
disp(' mass matrix')
m;
[ns,ms]=size(m);
% if forces are acting at degrees of freedom
m=[40000 0 0;0 20000 0;0 0 20000];
c0=[0 0 0;0 -8000 8000;0 -8000 8000];
c1=[0 0 0;0 -4000 8000;0 -8000 4000];
k0=[30000 -10000 0;-10000 20000 -10000;0 -10000 10000];
k1=[30000 -10000 0;-10000 50000 -10000;0 -10000 50000];
k2=[30000 -10000 0;-10000 20000 -40000;0 -40000 10000];
% disp(' force at various degrees of freedom');
% f;
%if base ground acceleration is given
% dis='disp.dat'
% di=load(dis)
% % convert to equivalent nodal loads
% for i=1:ns
% f(:,i)=-di*m(i,i)
% end
disp(' damping matrix')
c0;
disp(' stiffness matrix')
k0;
format long;
kim=inv(k0)*m;
[evec,ev]=eig(kim);
for i=1:ns
omega(i)=1/sqrt(ev(i,i));
end
disp(' natural frequencies')
omega;
% give gamma=0.5 and beta=0.25 for Newmark average accln method
%gama=0.5;
%beta=0.25;
%give gamma=0.5 and beta=0.1667 for Newmark linear accln method
gama=0.5;
beta=0.167;
%give initial conditions for displacements
u0=[0 0.01 0.05];
disp(' initial displacements')
u0;
%give initial condition for velocities
v0=[0. 0. 0.];
%y0=[0;0.01;0.05;0;0;0];
disp(' initial velocities')
v0;
om=5; eps=0.01;
IM=inv(m);
X=u0'; Z=v0';
tt=0;
% S1=-IM*k0+IM*k1*eps.*cos(om*tt)+IM*k2*eps^2.*cos(om*tt)^2;
% S2=-IM*c0+IM*c1*eps.*cos(om*tt);
dt=0.02;
S=k2*dt^2*beta*eps^2*cos(om*tt)^2+k1*dt^2*beta*eps*cos(om*tt)+...
c1*gama*dt*eps*cos(om*tt)+k0*dt^2*beta+c0*gama*dt+m;
S1=k1*eps.*cos(om*tt)+k2*eps^2.*cos(om*tt)^2;
S2=c1*eps.*cos(om*tt);
f(1,:)=(S1*X+S2*Z);
%for i=1:ns
a0=-inv(m)*(f(1,:)'+c0*v0'+k0*u0');
%end
kba=k0+(gama/(beta*dt))*c0+(1/(beta*dt*dt))*m;
kin=inv(kba);
aa=(1/(beta*dt))*m+(gama/beta)*c0;
bb=(1/(2.0*beta))*m+dt*(gama/(2.0*beta)-1)*c0;
u(1,:)=u0;
v(1,:)=v0;
a(1,:)=a0;
t=linspace(0,5,251);
for i=2:10%251
X=u(i-1,:)'; Z=v(i-1,:)';
S1=k1*eps.*cos(om*tt)+k2*eps^2.*cos(om*tt)^2;
S2=c1*eps.*cos(om*tt);
f(i,:)=(S1*X+S2*Z); %%%%% ??????
df=f(i,:)-f(i-1,:);
dfb(i,:)=df+v(i-1,:)*aa'+a(i-1,:)*bb';
du(i,:)=dfb(i,:)*kin;
dv(i,:)=(gama/(beta*dt))*du(i,:)-(gama/beta)*v(i-1,:)+dt*...
(1-gama/(2.0*beta))*a(i-1,:);
da(i,:)=(1/(beta*dt^2))*du(i,:)-(1/(beta*dt))*v(i-1,:)-(1/(2.0*beta))*a(i-1,:);
u(i,:)=u(i-1,:)+du(i,:);
v(i,:)=v(i-1,:)+dv(i,:);
a(i,:)=a(i-1,:)+da(i,:);
end
%figure(1);
hold on
plot(tt,u(:,1),'k');
xlabel(' time in secs');
ylabel(' roof displacement');
title(' displacement response of the roof');
%figure(2);
hold on
plot(tt,u(:,2),'k');
xlabel(' time in secs');
ylabel(' roof velocity');
title('velocity response of the roof');
%figure(3);
hold on
plot(tt,u(:,3),'k');
xlabel(' time in secs');
ylabel(' roof acceleration');
title(' acceleration response of the roof')
0 commentaires
Réponses (0)
Voir également
Catégories
En savoir plus sur Lighting, Transparency, and Shading dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!