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

1、M=10A5*diag(3.4 3.4 3.2 3.2 3.0 2.8 2.7 2.6);K=10A8*(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);%输入质量矩阵喇I度矩阵fai,w=eig(inv(K)*M);omiga=1./diag(w.A0.5);for i=1:8%寺征值问题求解%导到周期%0一化振型val,poi=max(abs(fai(:,i);fai(:,i)=fai(:,

2、i)/fai(poi,i);endgama=fai'*M*ones(8,1)./diag(fai'*M*fai);kesai=(0.1347+0.006306*omiga.A2)/2./omiga;%十算振型参与系数%十算各阶阻尼地震动响应图'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);

3、allacceleration(1,:)=A(1);allvelocity=zeros(N,8);vg(1)=0;xg(1)=0;for m=1:N-1烟取地震波数据%求时间间隔% 创立文件夹保存生成的所%用来存储所有时刻的位移%用来存储所有时刻的加速度值%寸0时刻加速度矩阵初始化%用于存储所有时刻速度计算场地速度、位移用于对加速度进行基线修正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)*mintA2/2+mintA2*(A(m+1)-A(m)/6;%十算endc1=28/13/TA2*(2*

4、vg(N)-15/TA5*xg*(3*T*B(:,1).A2-2*B(:,1).A3)*mint);系数c1c0=(vg(N)-c1*TA2/2)/T;%分正后的地震波加速度%啰正计算系数C0A=max(A)/max(abs(A-(c0+c1*B(:,1)*(A-(c0+c1*B(:,1);xg=xg'-(0.5*c0*B(:,1).A2+1/6*c1*B(:,1).A3);后的场地位移for i=1:8第三层循环用于对不同周期值计算反响S=1+kesai(i)*omiga(i)*mint+(omiga(i)*mint)A2)/6;for j=2:N对每个周期利用NewMark方法进行

5、数值积分得到对应值(具体请查看word文档)Q=2*omiga(i)*kesai(i)*allacceleration(j-1,i)*mint+.(omiga(i)A2)*(allvelocity(j-1,i)*mint+0.5*allacceleration(j-1,i)*mintA2);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+del

6、taa*mint/2;allmotion(j,i)= allmotion(j-1,i)+allvelocity(j-1,i)*mint+0.5*allacceleration(j-1,i)*mintA2+.deltaa*(mintA2)/6;end%积分计算反响并存放在前面的空矩阵里%计算位移% 十算加速度 %计算速度%RJ用figure函数准确的限制画图窗口魅制坐标轴范围,便于比拟end motion1=allmotion*fai'acceleration=allacceleration*fai'velocity=allvelocity*fai'for n=1:9 f

7、igure('color','white')if n=9axis(0,50,-0.03,0.03) endgrid onhold 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=9name='顶层位移时程,;motion(:

8、,9)=motion1(:,8)+xg;elsemotion(:,n)=motion1(:,n)-motion1(:,n-1);name=strcat('第',num2str(n-1),'-',num2str(n),' 层间位移时程,);%®相对位移时程%x轴标注%y轴标注%标注图名洞形全屏,便于查看end plot(B(:,1),motion(:,n),'linewidth',2)xlabel('时间(s) ','FontName','宋体,'FontSize',16);

9、ylabel('位移(m)','FontName','宋体','FontSize',16);title(name,'FontName','宋体','FontSize',20)set(gcf,'position',get(0,'screensize');shgF=getframe(gcf);imwrite(F.cdata,file_id,strcat(name,'.png');筛储得到的时程图end%计算顶层的层间相对位移%计算剪力GK=

10、repmat(-2 -2 -1.8 -1.8 -1.8 -1.8 -1.6 -1.6*10A8,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);ylabel('剪力(N) ','FontName',&

11、#39;宋体','FontSize',16);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%生成图形窗口,背景白色%画底层剪力时程图%x轴

12、标注%y轴标注晰注图名%图形全屏,便于查看%存储生成的图形%计算所需剪力和位移包络图数据,并得到便于画图的矩阵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-'

13、;);% 生成层间位移包络图 title('层间位移绝对值最大值包络图','FontName','宋体','FontSize',20);两注图名text(maxmotion(1:2:16),0.25:0.5:4,num2str(maxmotion(1:2:16)')%B 注位移数据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.2