八层框架的地震响应计算和人工波生成的matlab实现及所需曲线图的自动存储

1、一、 作业概况结构根本参数：层间剪切型结构，采用Rayleigh阻尼，第一、第二阶阻尼比分别取3%、5%。图1 结构根本形状表1 各层集中质量 (105kg)层号12345678质量3.403.403.203.202.802.802.702.60表2 各层层间刚度 (×108N/m)层号12345678层间刚度2.002.001.801.801.801.801.601.60二、 频率及振型计算根据层间模型的假定，可以建立结构的质量矩阵以及刚度矩阵如下。 根据上面求得的质量、刚度矩阵，即可求解特征方程： 1求解自振频率以及阵型向量已经演变成为典型的求解矩阵特征值以及特征向量的问题，特征

2、值即为圆频率，特征向量即为振型向量。根据式1利用matlab编程计算，求解矩阵的特征值以及特征向量，进而可以得到结构前八阶的自振圆频率、自振频率、自振周期如表2.1所示。表2.1 结构前八阶振型自振圆频率及周期振型14.76050.75771.3198213.36202.12660.4702321.56383.43200.2914429.12494.63540.2157536.16175.75530.1738641.24946.56500.1523744.87367.14180.1400848.03557.64510.1308而上述特征值问题求得的幅值向量实质即为振型向量，在下文中记为，振型矩

3、阵记为。振型向量经过归一化处理之后如表2.2所示。表2.2 结构前八阶振型层号1阶振型2阶振型3阶振型4阶振型5阶振型6阶振型7阶振型8阶振型10.175 -0.482 -0.784 0.830 -0.970 0.714 -0.703 0.144 20.344 -0.817 -0.948 0.463 0.216 -0.637 1.000 -0.277 30.516 -0.914 -0.298 -0.687 1.000 -0.091 -0.912 0.463 40.668 -0.721 0.599 -0.801 -0.541 0.730 0.440 -0.696 50.793 -0.299 1.

4、000 0.293 -0.824 -0.658 0.216 1.000 60.890 0.206 0.678 1.000 0.569 -0.305 -0.685 -0.893 70.963 0.710 -0.236 0.311 0.834 1.000 0.715 0.584 81.000 1.000 -0.965 -0.823 -0.742 -0.567 -0.315 -0.212 根据表2.2，可以绘制结构的前八阶振型图，如图2.1所示。 图2.1 结构前八阶振型图三、 采用振型分解法进行地震时程计算多自由度体系结构的位移反响可以表示为：2其中为广义坐标，为广义坐标向量。故线性结构的动力方程

5、可以表示为式3所示：3式中，为振型参与系数，满足而阻尼矩阵采用瑞利阻尼假定，满足：4其中利用振型的正交性，即：5在3式两端左乘，并与4、5联立，可得：6又因为之间满足：710、11联立，两边同除，并记，可得：8至此，完成了对多自由度耦合的动力方程的解耦，形成了假设干单自由度体系动力方程，可以据此利用时程分析方法进行各单自由度体系的时程计算。下面选用了纽马克-法进行单自由度体系的时程计算。在该方法中假定： 9其中、以及分别表示对应第i个广义坐标在第k个时间点的相对位移、速度和加速度；而、以及那么分别表示对应增量。和为参数，分别取1/2和1/6。又根据式8，其增量形式为：10所以将9带入10中有：

6、11在利用matlab编制程序时，记，所以有：12即当地震加速度时程以及上一时刻单自由度体系加速度、速度以及位移时，利用式1以及5即可求得体系在下一时刻的加速度、速度以及位移。不断递推计算最终可以求得体系在时域的地震响应。计算前，假定在初始0时刻，有：13利用matlab编制程序程序见附，求得了八个广义坐标的地震反响时程。采用此方法计算得到的位移时程为相对位移时程。因此为了得到结构顶部的绝对位移时程，还需要计算出场地位移时程后，与顶部相对位移时程相加。由于地震仪在记录地震动时记录纸的蛇行运动和放大器的不稳定等，记录的零线会产生很小的摇摆错位，因此在时程计算之前需要对地震波进行基线校正。采用线性

7、修正的方法，地震波位移、速度、加速度的修正值如式14所示。14其中，系数由下面的式子计算：15经过上述过程编程计算后，结构顶部位移时程如图3.1所示。图3.1 EL-Centro波作用下结构顶层位移时程各层层间位移时程如图3.2所示。图3.2 EL-Centro波作用下结构各层层间位移时程底层层间剪力时程如图3.3所示。图3.3 EL-Centro波作用下结构底层层间剪力时程各层的层间位移与层间剪力绝对值包络图如图3.4所示。(a)层间位移绝对值包络图m b层间剪力绝对值包络图kN图3.4 EL-Centro波作用下结构各层层间位移与层间剪力绝对值包络图下面对上述计算程序及结果的正确性进行简要

8、的验证。验证采用的方法为振型分解反响谱法，即利用该方法得到结构在EL-Centro波作用下各层层间剪力的最大值，与上文中时程计算得到的剪力绝对值包络图进行比照。根据相关标准，在罕遇地震作用下结构阻尼取5%计算得到EL-Centro波NS向的加速度反响谱，反响谱图如图3.5所示。图3.5 EL-Centro波在5%阻尼比下绝对加速度反响谱将结构前八阶振型的周期域反响谱进行比对利用线性内插法，可以得到相应的结构最大绝对加速度。结构的前八阶振型以及振型参与系数、最大绝对加速度通过计算列于表3.1中。表3.1 振型分解反响谱法所需结构参数振型11.31981.3033 2.05420.4702-0.4

9、600 7.92630.2914-0.2548 6.69640.21570.1673 6.34350.17380.1090 6.86460.15230.0852 5.22970.1400-0.06106.56180.13080.0145 7.334所以利用表3.1中的数据，利用下式可以求得各振型各质点的最大地震作用其中j表示振型编号，i表示由下自上的质点号。16列表如3.2所示。表3.2 各振型各质点的最大地震作用单位：kNji123456781159.5597.2 427.8 281.9 203.1 89.1 75.9 4.0 2312.81013.1 517.4 157.3 -45.3 -

10、79.5 -108.1 -7.7 3469.81133.6 162.5 -233.2 -209.4 -11.3 98.5 12.8 4607.9894.2 -326.7 -272.0 113.3 91.1 -47.6 -19.2 5721.4371.1 -545.9 99.4 172.6 -82.1 -23.4 27.6 6809.6-255.2 -370.2 339.7 -119.2 -38.0 74.0 -24.7 7876.6-880.0 128.7 105.8 -174.7 124.7 -77.3 16.1 8910.1-1239.6 526.6 -279.5 155.3 -70.7

11、34.0 -5.9 因此，通过表3.2可以得到各个振型下各层间剪力的最大值，如表3.3所示。表3.3 各振型各层间剪力的最大值单位：kNji1234567814867.81634.4 520.3 199.3 95.7 23.4 26.2 3.1 24708.3 1037.2 92.5 -82.7 -107.4 -65.7 -49.7 -0.9 34395.5 24.1 -425.0 -240.0 -62.1 13.8 58.3 6.8 43925.7 -1109.5 -587.5 -6.8 147.3 25.1 -40.2 -6.0 53317.8 -2003.7 -260.8 265.3 3

12、4.0 -66.0 7.4 13.2 62596.4 -2374.8 285.1 165.9 -138.6 16.1 30.8 -14.4 71786.8 -2119.6 655.3 -173.7 -19.4 54.1 -43.3 10.3 8910.1 -1239.6 526.6 -279.5 155.3 -70.7 34.0 -5.9 由于各振型下层间剪力的最大值并不一定同时出现，因此采用平方和开方的方式，求得地震波作用下层间剪力绝对值的最大值，如表3.4所示。表3.4 各层间剪力的最大值单位：kN编号123456785166.0 4824.7 4423.4 4124.5 3894.5 3

13、537.0 2854.8 1658.6 从表3.4可以得到利用振型分解反响谱法求得的在El Centro波NS向作用下各层间剪力绝对值的包络图，如图3.6a所示，3.6b为上文中求得的时程计算得到的各层间剪力绝对值的包络图。图3.6 两种方法计算EL-Centro波作用下各层层间剪力绝对值包络图通过比照可以发现时程计算的结果与振型分解反响谱法的计算结果在底层吻合的很好，在结构高程的差距相对较大，但整体上两者的结果是在同一数量级，且差值合理的，最大差值为26.9%。考虑到振型分解反响谱法本身上是一种为设计效劳求得最大地震作用的近似参考方法，所以通过与本作业中编程进行时程计算的结果相比，在一定程度

14、上从侧面印证了前文时程计算结果的合理性与正确性。四、 人工地震波合成本文采用经典谱表达方法来进行强震加速度模拟。经典谱表达方法如式17所示。17式中：为时间包络函数，表达式如下所示：18相关参数选择按照表4.1选取：表4.1 时间包络函数的参数值及持时为主振平稳段段首时间，对设计的两条人工波假定分别对应二类及三类场地的第一组地震波，可分别取0.8s以及1.2s；为主振平稳段段尾时间，对两条人工波分别取7s以及9s；为衰减系数，对两条人工波分别取0.35s-1以及0.25 s-1。式17中的可由下式计算。19即为频域样本点，为对应样本点的相位，为所选用的功率谱。满足：20其中为截断频率，为频率下

15、界，为频域的样本点数目。所以，在上式的根底上，需要功率谱密度函数模型以及相位谱，才能够生成人工地震波。而根据作业要求， 选用功率谱密度函数模型Clough-Penzien谱模型，其表达式为：21式中相关参数选择按照表4.2选取：表4.2 Clough-Penzien模型场地土参数设计值为场地卓越周期，对设计的两条人工波假定分别对应二类及三类场地的第一组地震波，故分别取20.94与13.96；为二次过滤参数，取为0.15一般可取0.1-0.2；为场地地基阻尼，对两条地震波分别取为0.72以及0.8；为二次过滤参数，取=；为谱强度因子，在本文中，为使两条人工地震波的加速度峰值在0.4g左右，经过调

16、整试算，分别取0.025、0.032。地震动持时采用相对分数持时，即地震动加速度记录最后一次到达峰值加速度0.05倍时所对应的时间。上述参数的选择参考文献：张猛,张哲,李天. 与标准反响谱相对应Clough-Penzien模型参数研究J. 世界地震工程,2007,01:56-60.两条地震波对应的Clough-Penzien模型如图4.1所示。(a)人工波1(b)人工波2图4.1人工波Clough-Penzien模型功率谱图相位谱那么采用随机方式生成随机相位谱。故综上利用matlab编制人工地震波合成程序程序见附，可以得到两条人工地震波地震波数据见附，两条波的加速度时程如图4.2所示。 (a)

17、人工波1(b)人工波2图4.2 人工波加速度时程五、 人工波时程计算结果将四中人工生成的两条地震波用前文的方法进行时程计算，得到了响应的计算结果。人工波1作用下，结构顶层位移时程如图5.1所示。图5.1 人工波1作用下结构顶层位移时程人工波1作用下，结构各层层间位移时程如图5.2所示。图5.2 人工波1作用下结构各层层间位移时程在人工波1作用下，结构底层剪力时程如图5.3所示。图5.3 人工波1作用下结构底层剪力时程在人工波1作用下，各层的层间位移与层间剪力绝对值包络图如图5.4所示。(a)层间位移绝对值包络图b层间剪力绝对值包络图图5.4 人工波1作用下结构各层层间位移与层间剪力绝对值包络图

18、人工波2作用下，结构顶层位移时程如图5.5所示。图5.5 人工波2作用下结构顶层位移时程人工波2作用下，结构各层层间位移时程如图5.6所示。图5.6 人工波2作用下结构各层层间位移时程在人工波2作用下，结构底层剪力时程如图5.7所示。图5.7 人工波2作用下结构底层剪力时程在人工波2作用下，各层的层间位移与层间剪力绝对值包络图如图5.8所示。(a)层间位移绝对值包络图 b层间剪力绝对值包络图图5.8 人工波2作用下结构各层层间位移与层间剪力绝对值包络图1、程序在2021Ra上可以正常运行2、运行时需要把ThreeKindsWaves复制到当前工作路径下3、文件的存储地址可以根据自身情况进行修改

19、，或者不需要可以直接注释%利用解特征方程的方法得到结构的8阶自振周期和振型，并归一化，然后利用循环产生阵型图，保存在文件阵型图内tic %计时开始clear all;%去除变量digits(4) %设定精度M=105*diag(3.4 3.4 3.2 3.2 2.8 2.8 2.7 2.6);%输入质量矩阵K=108*(diag(4 3.8 3.6 3.6 3.6 3.4 3.2 1.6)+diag(-2 -1.8 -1.8 -1.8 -1.8 -1.6 -1.6,1)+. diag(-2 -1.8 -1.8 -1.8 -1.8 -1.6 -1.6,-1);%输入刚度矩阵fai,omiga=e

20、ig(M-1*K);%利用matlab内置函数解特征方程，周期平方存为omiga，振型存为faiomiga=omiga.0.5; %得到周期transfer=fai(:,8);fai(:,8)=fai(:,7);fai(:,7)=transfer;file_id='C:UserslenovoDesktop第二次作业用图振型图'mkdir(file_id);for i=1:8 %归一化振型并画图 val,poi=max(abs(fai(:,i); fai(:,i)=fai(:,i)/fai(poi,i); f(i)=omiga(i,i)/(2*pi); figure('c

21、olor','white');%建立图形窗口，背景白色 axis equal axis off hold on plot(0 0,0,4,'k-') x=0 fai(:,i)' y=linspace(0,4,9); yi=0:0.001:4; pp=interp1(y,x,yi,'cubic');%利用内置插值函数使连线光滑 plot(pp,yi,'k-') plot(fai(:,i),0.5:0.5:4,'o','MarkerFaceColor','k','

22、MarkerSize',20)%控制点的形状和大小代表集中质量 text(fai(:,i)+0.2,0.6:0.5:4.1,num2str(fai(:,i)%标注振型数值 for j=1:8 plot(0 fai(j,i),0.5*j 0.5*j,'k-')%画相关的辅助线 end name=strcat('第',num2str(9-i),'阶振型图'); title(name,'FontSize',16,'position',0 -0.3)%标注图名 F=getframe(figure(i); imwri

23、te(F.cdata,file_id,strcat(name,'.png'); %保存振型图片 endf=f'T=1./f;%输出周期和频率Toc%通过解特征方程得到振型和周期，并进行归一化，为下面进行时程分析做铺垫。读取Excel表格存储的地震动数据，%利用纽马克方法进行时程分析得到地面速度和地面位移，利用相关公式对数据进行基线校正得到更为合理的数据。%利用纽马克方法进行时程分析得到结构各阶振型的响应，然后通过公式X=A*u得到结构总的位移响应，然后利用%结构力学相关知识得到剪力和剪力包络图。tic;clear all ;clf ; %计时开始，并清楚其他无关变量M=

24、105*diag(3.4 3.4 3.2 3.2 2.8 2.8 2.7 2.6); %输入质量矩阵K=108*(diag(4 3.8 3.6 3.6 3.6 3.4 3.2 1.6)+diag(-2 -1.8 -1.8 -1.8 -1.8 -1.6 -1.6,1)+. diag(-2 -1.8 -1.8 -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(:,

25、i)/fai(poi,i);endgama=fai'*M*ones(8,1)./diag(fai'*M*fai); %计算振型参与系数kesai=(0.1329+0.0067*omiga.2)/2./omiga; %计算各阶阻尼file_id='C:UserslenovoDesktop第二次作业用图地震动响应图'mkdir(file_id); %创立文件夹保存生成的所有图A=xlsread('ThreeKindsWaves','ELcentro-NS'); %读取地震波数据N=length(A(:,1); mint=A(2,1)-

26、A(1,1); %求时间间隔T=(N-1)*mint;allmotion=zeros(N,8); %用来存储所有时刻的位移allacceleration=zeros(N,8); %用来存储所有时刻的加速度值allacceleration(1,:)=A(1,2); %对0时刻加速度矩阵初始化allvelocity=zeros(N,8); %用于存储所有时刻速度vg(1)=0;xg(1)=0;for m=1:N-1 %计算场地速度、位移用于对加速度进行基线修正 vg(m+1)=vg(m)+A(m,2)*mint+(A(m+1,2)-A(m,2)*mint/2; xg(m+1)=xg(m)+vg(m

27、)*mint+A(m,2)*mint2/2+mint2*(A(m+1,2)-A(m,2)/6;end c1=28/13/T2*(2*vg(N)-15/T5*xg*(3*T*A(:,1).2-2*A(:,1).3)*mint); %计算系数c1 c0=(vg(N)-c1*T2/2)/T; %计算系数c0 A(:,2)=max(A(:,2)/max(abs(A(:,2)-(c0+c1*A(:,1)*(A(:,2)-(c0+c1*A(:,1); %修正后的地震波加速度 xg=xg'-(0.5*c0*A(:,1).2+1/6*c1*A(:,1).3); %修正后的场地位移 for i=1:8

28、%第三层循环用于对不同周期值计算反响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)*mint2); deltaa=-(gama(i)*(A(j,2)-A(j-1,2)+Q)/S; allacceleration(j,i)

29、=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)*mint2+. deltaa*(mint2)/6; end %积分计算反响并存放在前面的空矩阵里 endmotion1=allmotion*fai' %计算位移acceleration=allacce

30、leration*fai' %计算加速度velocity=allvelocity*fai' %计算速度for n=1:9figure('color','white') %利用figure函数准确的控制画图窗口if n=9axis(0,50,-0.03,0.03) %控制坐标轴范围，便于比较endgrid on hold onbox off%annotation('arrow',0.132 0.132,0.8 1);%annotation('arrow',0.8 1,0.108 0.108); %产生坐标轴箭头可能不

31、够美观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),'层间位移时程');endplot(A(:,1),motion(:,n),'linewidth&#

32、39;,2) %画相对位移时程xlabel('时间s','FontName','宋体','FontSize',16); %x轴标注ylabel('位移(m)','FontName','宋体','FontSize',16); %y轴标注title(name,'FontName','宋体','FontSize',20) %标注图名set(gcf,'position',get(0,'screensize

33、'); %图形全屏，便于查看shg F=getframe(gcf); imwrite(F.cdata,file_id,strcat(name,'.png'); %存储得到的时程图 endGK=repmat(-2 -2 -1.8 -1.8 -1.8 -1.8 -1.6 -1.6*108,N,1);%输入刚度矩阵motion(:,8)=motion1(:,8)-motion1(:,7); %计算顶层的层间相对位移Fq=motion(:,1:8).*GK; %计算剪力figure('color','white') %生成图形窗口，背景白色plo

34、t(A(:,1),Fq(:,1) %画底层剪力时程图xlabel('时间s','FontName','宋体','FontSize',16); %x轴标注ylabel('剪力N','FontName','宋体','FontSize',16); %y轴标注title('底层剪力时程','FontName','宋体','FontSize',20); %标注图名set(gcf,'position'

35、,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

36、)=max(abs(Fq(:,f);endfigure('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)

37、9;) %标注位移数据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('

38、层间剪力绝对值最大值包络图','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('层间剪力绝对

39、值最大值包络图','.png'); toc %计时结束%利用地震波的经典表达方式和Clough-Penzien功率谱表达式计算对应功率以及振幅和功率谱的关系得到振幅，%利用随机函数得到相位，进而得到各类场地地震动数据。并把生成的地震动数据和相关图形保存function Artificial_wavetic %计时开始clear all; %去除无关变量 w=linspace(0,100,10000);%生成周期矩阵kesai=0.72 0.8; %输入阻尼矩阵w0=20.94 13.96; %输入周期矩阵w1=0.15*w0;kesai1=kesai;t=linspac

40、e(0,40,2001);%生成时间矩阵，便于后期画图t1=0.8 1.2; %输入强度函数相关参数t2=7 9;c=0.35 0.25;file_id='C:UserslenovoDesktop第二次作业用图人造波'mkdir(file_id); %创立文件夹保存生成的所有图和地震动数据for i=1:2 %利用循环以此计算两类地震动数据 for j=1:10000 %利用循环计算功率谱 a=(1+4*kesai(1)2*(w(j)/w0(i)2)*(w(j)/w1(i)4; b=(1-(w(j)/w0(i)2)2+4*kesai(i)2*(w(j)/w0(i)2; d=(1

41、-(w(j)/w1(i)2)2+4*kesai1(i)2*(w(j)/w1(i)2; S(j)=a*0.0252/(b*d); A(j)=sqrt(2*S(j)*(w(2)-w(1)/pi); %计算对应振幅 end %功率谱计算结束 fai=repmat(2*pi*rand(10000,1)-pi,1,2001);%生成相位矩阵 for n=1:2001 %计算强度函数值 gt(n)=g(t(n),t1(i),t2(i),c(i); end acceleration=A*cos(w'*t+fai).*gt; %计算地震动数据 N= find(abs(acceleration)>

42、=0.02*max(abs(acceleration);%寻找截断点 figure(2*i-1); %利用figure函数控制图形窗口 plot(w,S); %得到功率谱图 name=strcat(num2str(i+1),'类场地人工波功率谱'); title(name); set(gcf,'position',get(0,'screensize'); %图形全屏，便于查看 F=getframe(gcf); imwrite(F.cdata,file_id,strcat(name,'.png'); %存储功率谱图形 figure(

43、2*i); plot(t(1:N(end),acceleration(1:N(end); %画地震动时程图 name=strcat(num2str(i+1),'类场地人工波加速度时程'); title(name); set(gcf,'position',get(0,'screensize'); %图形全屏，便于查看 F=getframe(gcf); imwrite(F.cdata,file_id,strcat(name,'.png'); %存储地震动时程图 xlswrite(strcat(file_id,'第',n

44、um2str(i+1),'类场地人工波.xls'),t(1:N(end);acceleration(1:N(end)');%存储地震动数据endendfunction gt=g(t,t1,t2,c)%强度包络函数计算 if t<t1 gt=(t/t1)2; elseif t>t2 gt=exp(-c*(t-t2); else gt=1; endend%该代码与Earthquake_Response.m根本一致故不进行详细标注tic;clear all ;clf ; %计时开始M=105*diag(3.4 3.4 3.2 3.2 2.8 2.8 2.7 2.6

45、);K=108*(diag(4 3.8 3.6 3.6 3.6 3.4 3.2 1.6)+diag(-2 -1.8 -1.8 -1.8 -1.8 -1.6 -1.6,1)+. diag(-2 -1.8 -1.8 -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);endgama=fai'*M*ones(8,1)./diag(fai'*M*fai);kes

46、ai=(0.1329+0.0067*omiga.2)/2./omiga;file_id='C:UserslenovoDesktop第二次作业用图人工地震动响应图'mkdir(file_id); %创立文件夹保存生成的所有图B='第2类场地人工波','第3类场地人工波'for a=1:2clear A N mint TA=xlsread('ThreeKindsWaves',Ba); %读取地震波数据A(:,2)=-A(:,2);N=length(A(:,1); mint=A(2,1)-A(1,1); %求时间间隔T=(N-1)*mi

47、nt;allmotion=zeros(N,8); %用来存储所有时刻的位移allacceleration=zeros(N,8); %用来存储所有时刻的加速度值allacceleration(1,:)=A(1,2); %对0时刻加速度矩阵初始化allvelocity=zeros(N,8); %用于存储所有时刻速度motion1=zeros(N,8);motion=zeros(N,8);vg=zeros(1,N);xg=zeros(1,N);for m=1:N-1 %计算场地速度、位移用于对加速度进行基线修正 vg(m+1)=vg(m)+A(m,2)*mint+(A(m+1,2)-A(m,2)*m

48、int/2; xg(m+1)=xg(m)+vg(m)*mint+A(m,2)*mint2/2+mint2*(A(m+1,2)-A(m,2)/6;end for i=1:8 %第三层循环用于对不同周期值计算反响S=1+kesai(i)*omiga(i)*mint+(omiga(i)*mint)2)/6; for j=2:N %对每个周期利用NewMark方法进行数值积分得到对应值 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)*mint2); deltaa=-(gam