clear;
%用Newmark方法计算结构动力反应(无阻尼)
m=[5.0e5 0;
0 2.0e4];
k=[6.8e6 -1.8e6;
-1.8e6 1.8e6];
beta=1/6; %定义系数beta
n=length(m);
tt=13.1;
pbc=xlsread('mydata.xlsx');
pp1=-1.8e6*pbc(:,1);
pp2=1.8e6*pbc(:,1);
p(:,:)=[pp1';pp2'];
y(:,1)=zeros(n,1); %初始位移向量
y1(:,1)=zeros(n,1); %初始速度向量
y2(:,1)=inv(m)*(p(:,1)-k*y(:,1)); %初始加速度向量
detat=0.02; %计算取用时间间隔
a0=1/(beta*detat^2);
a1=0.5/(beta*detat);
a2=1/(beta*detat);
a3=1/(2*beta)-1;
a4=0.5/beta-1;
a5=detat/2*(0.5/beta-2);
a6=detat*(1-0.5);
a7=0.5*detat;
k2=k+a0*m; %等效刚度
for t=detat:detat:tt
j=round(t/detat)+1;
p(:,j)=[pp1(j,1);pp2(j,1)];
%计算j时刻位移、加速度和速度向量
y(:,j)=inv(k2)*(p(:,j)+m*(a0*y(:,j-1)+a2*y1(:,j-1)+a3*y2(:,j-1)));
y2(:,j)=a0*(y(:,j)-y(:,j-1))-a2*y1(:,j-1)-a3*y2(:,j-1);
y1(:,j)=y1(:,j-1)+a6*y2(:,j-1)+a7*y2(:,j);
end
tt2=0:detat:tt;
subplot(2,2,1),plot(tt2,y(1,:))
title('The Newmark Method')
xlabel('t')
ylabel('X1')
subplot(2,2,2),plot(tt2,y(2,:))
title('The Newmark Method')
xlabel('t')
ylabel('X2')
fid=fopen('shuju2.m','w+');
fprintf(fid,'\nt x1 x2\n');
for i=1:j
fprintf(fid,'\n% +7.5f %+12.4e %12.4e',tt2(i),y(1,i),y(2,i));
end
fclose(fid);
评论0