地震工程学大作业代码(求ELCENTRO波的位移

M=10^5*diag([3.4 3.4 3.2 3.2 3.0 2.8 2.7 2.6]); %输入质量矩阵

K=10^8*(diag([4.4 4.2 4 3.8 3.6 3.4 3.2 1.6])+...

diag([-2.2 -2.0 -2.0 -1.8 -1.8 -1.6 -1.6],1)+...

diag([-2.2 -2.0 -2.0 -1.8 -1.8 -1.6 -1.6],-1)); %刚度矩阵

[fai,w]=eig(inv(K)*M); %特征值问题求解

omiga=1./diag(w.^0.5); %得到周期

for i=1:8 %归一化振型

[val,poi]=max(abs(fai(:,i)));

fai(:,i)=fai(:,i)/fai(poi,i);

end

gama=fai'*M*ones(8,1)./diag(fai'*M*fai); %计算振型参与系数

kesai=(0.1347+0.006306*omiga.^2)/2./omiga; %计算各阶阻尼

file_id='D:\地震动响应图\';

mkdir(file_id); %创建文件夹保存生成的所有图

A= textread('NS.txt','%n'); %读取地震波数据

N=length(A);

mint=0.01; %求时间间隔

T=(N-1)*mint;

B=xlsread('T');

allmotion=zeros(N,8); %用来存储所有时刻的位移

allacceleration=zeros(N,8); %用来存储所有时刻的加速度值

allacceleration(1,:)=A(1); %对0时刻加速度矩阵初始化

allvelocity=zeros(N,8); %用于存储所有时刻速度

vg(1)=0;

xg(1)=0;

for m=1:N-1 %计算场地速度、位移用于对加速度进行基线修正

vg(m+1)=vg(m)+A(m)*mint+(A(m+1)-A(m))*mint/2;

xg(m+1)=xg(m)+vg(m)*mint+A(m)*mint^2/2+mint^2*(A(m+1)-A(m))/6;

end

c1=28/13/T^2*(2*vg(N)-15/T^5*xg*(3*T*B(:,1).^2-2*B(:,1).^3)*mint); %计算系数c1

c0=(vg(N)-c1*T^2/2)/T; %计算系数c0

A=max(A)/max(abs(A-(c0+c1*B(:,1))))*(A-(c0+c1*B(:,1))); %修正后的地震波加速度

xg=xg'-(0.5*c0*B(:,1).^2+1/6*c1*B(:,1).^3); %修正后的场地位移

for i=1:8 %第三层循环用于对不同周期值计算反应

S=1+kesai(i)*omiga(i)*mint+((omiga(i)*mint)^2)/6;

for j=2:N %对每个周期利用NewMark方法进行数值积分得到对应值（具体请查看word文档）Q=2*omiga(i)*kesai(i)*allacceleration(j-1,i)*mint+...

(omiga(i)^2)*(allvelocity(j-1,i)*mint+0.5*allacceleration(j-1,i)*mint^2);

deltaa=-(gama(i)*(A(j)-A(j-1))+Q)/S;

allacceleration(j,i)=allacceleration(j-1,i)+deltaa;

allvelocity(j,i)= allvelocity(j-1,i)+allacceleration(j-1,i)*mint+deltaa*mint/2;

allmotion(j,i)= allmotion(j-1,i)+allvelocity(j-1,i)*mint+0.5*allacceleration(j-1,i)*mint^2+...

deltaa*(mint^2)/6;

end %积分计算反应并存放在前面的空矩阵里end

motion1=allmotion*fai'; %计算位移

acceleration=allacceleration*fai'; %计算加速度

velocity=allvelocity*fai'; %计算速度

for n=1:9

figure('color','white') %利用figure函数准确的控制画图窗口

if n~=9

axis([0,50,-0.03,0.03]) %控制坐标轴范围，便于比较

end

grid on

hold on

box off

%annotation('arrow',[0.132 0.132],[0.8 1]);

%annotation('arrow',[0.8 1],[0.108 0.108]); %产生坐标轴箭头（可能不够美观）if n==1 %计算相关数据和得到图名

name='底层层间位移时程';

motion(:,n)=motion1(:,n);

elseif n==9

name='顶层位移时程';

motion(:,9)=motion1(:,8)+xg;

else

motion(:,n)=motion1(:,n)-motion1(:,n-1);

name=strcat('第',num2str(n-1),'-',num2str(n),'层间位移时程');

end

plot(B(:,1),motion(:,n),'linewidth',2) %画相对位移时程

xlabel('时间（s）','FontName','宋体','FontSize',16); %x轴标注

ylabel('位移(m)','FontName','宋体','FontSize',16); %y轴标注

title(name,'FontName','宋体','FontSize',20) %标注图名

set(gcf,'position',get(0,'screensize')); %图形全屏，便于查看

shg

F=getframe(gcf);

imwrite(F.cdata,[file_id,strcat(name,'.png')]); %存储得到的时程图

end

GK=repmat(-[-2 -2 -1.8 -1.8 -1.8 -1.8 -1.6 -1.6]*10^8,N,1);%输入刚度矩阵

motion(:,8)=motion1(:,8)-motion1(:,7); %计算顶层的层间相对位移

Fq=motion(:,1:8).*GK; %计算剪力

figure('color','white') %生成图形窗口，背景白色plot(B(:,1),Fq(:,1)) %画底层剪力时程图

xlabel('时间（s）','FontName','宋体','FontSize',16); %x轴标注

ylabel('剪力（N）','FontName','宋体','FontSize',16); %y轴标注

title('底层剪力时程','FontName','宋体','FontSize',20); %标注图名

set(gcf,'position',get(0,'screensize')); %图形全屏，便于查看

F=getframe(gcf); %存储生成的图形imwrite(F.cdata,[file_id,strcat('底层剪力时程图','.png')]);

maxFq=zeros(1,16);

for f=1:8 %计算所需剪力和位移包络图数据，并得到便于画图的矩阵

maxmotion(2*f-1)=max(abs(motion(:,f)));

maxmotion(2*f)=maxmotion(:,2*f-1);

maxFq(2*f-1)=max(abs(Fq(:,f)));

maxFq(2*f)=max(abs(Fq(:,f)));

end

figure('color','white')

plot(maxmotion,[0 0.5 0.5 1 1 1.5 1.5 2 2 2.5 2.5 3 3 3.5 3.5 4],'k-');%生成层间位移包络图title('层间位移绝对值最大值包络图','FontName','宋体','FontSize',20); %标注图名

text(maxmotion(1:2:16),0.25:0.5:4,num2str(maxmotion(1:2:16)')) %标注位移数据set(gcf,'position',get(0,'screensize')); %图形全屏，便于查看

F=getframe(gcf);

imwrite(F.cdata,[file_id,strcat('层间位移绝对值最大值包络图','.png')]); %存储图形

figure('color','white')

plot(maxFq/1000,[0 0.5 0.5 1 1 1.5 1.5 2 2 2.5 2.5 3 3 3.5 3.5 4],'k-');

title('层间剪力绝对值最大值包络图','FontName','宋体','FontSize',20); %标注图名

text(maxFq(1:2:16)/1000,0.25:0.5:4,num2str(0.001*maxFq(1,1:2:16)'),'FontSize',10);

set(gcf,'position',get(0,'screensize')); %图形全屏，便于查看

F=getframe(gcf);

imwrite(F.cdata,[file_id,strcat('层间剪力绝对值最大值包络图','.png')]);