八層框架的地震響應(yīng)計(jì)算和人工波生成的matlab實(shí)現(xiàn)及所需曲線圖的自動(dòng)存儲(chǔ)

1、一、 作業(yè)概況結(jié)構(gòu)根本參數(shù)：層間剪切型結(jié)構(gòu)，采用Rayleigh阻尼，第一、第二階阻尼比分別取3%、5%。圖1 結(jié)構(gòu)根本形狀表1 各層集中質(zhì)量 (105kg)層號(hào)12345678質(zhì)量3.403.403.203.202.802.802.702.60表2 各層層間剛度 (×108N/m)層號(hào)12345678層間剛度2.002.001.801.801.801.801.601.60二、 頻率及振型計(jì)算根據(jù)層間模型的假定，可以建立結(jié)構(gòu)的質(zhì)量矩陣以及剛度矩陣如下。 根據(jù)上面求得的質(zhì)量、剛度矩陣，即可求解特征方程： 1求解自振頻率以及陣型向量已經(jīng)演變成為典型的求解矩陣特征值以及特征向量的問題，特征

2、值即為圓頻率，特征向量即為振型向量。根據(jù)式1利用matlab編程計(jì)算，求解矩陣的特征值以及特征向量，進(jìn)而可以得到結(jié)構(gòu)前八階的自振圓頻率、自振頻率、自振周期如表2.1所示。表2.1 結(jié)構(gòu)前八階振型自振圓頻率及周期振型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而上述特征值問題求得的幅值向量實(shí)質(zhì)即為振型向量，在下文中記為，振型矩

3、陣記為。振型向量經(jīng)過(guò)歸一化處理之后如表2.2所示。表2.2 結(jié)構(gòu)前八階振型層號(hào)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 根據(jù)表2.2，可以繪制結(jié)構(gòu)的前八階振型圖，如圖2.1所示。 圖2.1 結(jié)構(gòu)前八階振型圖三、 采用振型分解法進(jìn)行地震時(shí)程計(jì)算多自由度體系結(jié)構(gòu)的位移反響可以表示為：2其中為廣義坐標(biāo)，為廣義坐標(biāo)向量。故線性結(jié)構(gòu)的動(dòng)力方程

5、可以表示為式3所示：3式中，為振型參與系數(shù)，滿足而阻尼矩陣采用瑞利阻尼假定，滿足：4其中利用振型的正交性，即：5在3式兩端左乘，并與4、5聯(lián)立，可得：6又因?yàn)橹g滿足：710、11聯(lián)立，兩邊同除，并記，可得：8至此，完成了對(duì)多自由度耦合的動(dòng)力方程的解耦，形成了假設(shè)干單自由度體系動(dòng)力方程，可以據(jù)此利用時(shí)程分析方法進(jìn)行各單自由度體系的時(shí)程計(jì)算。下面選用了紐馬克-法進(jìn)行單自由度體系的時(shí)程計(jì)算。在該方法中假定： 9其中、以及分別表示對(duì)應(yīng)第i個(gè)廣義坐標(biāo)在第k個(gè)時(shí)間點(diǎn)的相對(duì)位移、速度和加速度；而、以及那么分別表示對(duì)應(yīng)增量。和為參數(shù)，分別取1/2和1/6。又根據(jù)式8，其增量形式為：10所以將9帶入10中有：

6、11在利用matlab編制程序時(shí)，記，所以有：12即當(dāng)?shù)卣鸺铀俣葧r(shí)程以及上一時(shí)刻單自由度體系加速度、速度以及位移時(shí)，利用式1以及5即可求得體系在下一時(shí)刻的加速度、速度以及位移。不斷遞推計(jì)算最終可以求得體系在時(shí)域的地震響應(yīng)。計(jì)算前，假定在初始0時(shí)刻，有：13利用matlab編制程序程序見附，求得了八個(gè)廣義坐標(biāo)的地震反響時(shí)程。采用此方法計(jì)算得到的位移時(shí)程為相對(duì)位移時(shí)程。因此為了得到結(jié)構(gòu)頂部的絕對(duì)位移時(shí)程，還需要計(jì)算出場(chǎng)地位移時(shí)程后，與頂部相對(duì)位移時(shí)程相加。由于地震儀在記錄地震動(dòng)時(shí)記錄紙的蛇行運(yùn)動(dòng)和放大器的不穩(wěn)定等，記錄的零線會(huì)產(chǎn)生很小的搖擺錯(cuò)位，因此在時(shí)程計(jì)算之前需要對(duì)地震波進(jìn)行基線校正。采用線性

7、修正的方法，地震波位移、速度、加速度的修正值如式14所示。14其中，系數(shù)由下面的式子計(jì)算：15經(jīng)過(guò)上述過(guò)程編程計(jì)算后，結(jié)構(gòu)頂部位移時(shí)程如圖3.1所示。圖3.1 EL-Centro波作用下結(jié)構(gòu)頂層位移時(shí)程各層層間位移時(shí)程如圖3.2所示。圖3.2 EL-Centro波作用下結(jié)構(gòu)各層層間位移時(shí)程底層層間剪力時(shí)程如圖3.3所示。圖3.3 EL-Centro波作用下結(jié)構(gòu)底層層間剪力時(shí)程各層的層間位移與層間剪力絕對(duì)值包絡(luò)圖如圖3.4所示。(a)層間位移絕對(duì)值包絡(luò)圖m b層間剪力絕對(duì)值包絡(luò)圖kN圖3.4 EL-Centro波作用下結(jié)構(gòu)各層層間位移與層間剪力絕對(duì)值包絡(luò)圖下面對(duì)上述計(jì)算程序及結(jié)果的正確性進(jìn)行簡(jiǎn)要

8、的驗(yàn)證。驗(yàn)證采用的方法為振型分解反響譜法，即利用該方法得到結(jié)構(gòu)在EL-Centro波作用下各層層間剪力的最大值，與上文中時(shí)程計(jì)算得到的剪力絕對(duì)值包絡(luò)圖進(jìn)行比照。根據(jù)相關(guān)標(biāo)準(zhǔn)，在罕遇地震作用下結(jié)構(gòu)阻尼取5%計(jì)算得到EL-Centro波NS向的加速度反響譜，反響譜圖如圖3.5所示。圖3.5 EL-Centro波在5%阻尼比下絕對(duì)加速度反響譜將結(jié)構(gòu)前八階振型的周期域反響譜進(jìn)行比對(duì)利用線性內(nèi)插法，可以得到相應(yīng)的結(jié)構(gòu)最大絕對(duì)加速度。結(jié)構(gòu)的前八階振型以及振型參與系數(shù)、最大絕對(duì)加速度通過(guò)計(jì)算列于表3.1中。表3.1 振型分解反響譜法所需結(jié)構(gòu)參數(shù)振型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中的數(shù)據(jù)，利用下式可以求得各振型各質(zhì)點(diǎn)的最大地震作用其中j表示振型編號(hào)，i表示由下自上的質(zhì)點(diǎn)號(hào)。16列表如3.2所示。表3.2 各振型各質(zhì)點(diǎn)的最大地震作用單位：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 因此，通過(guò)表3.2可以得到各個(gè)振型下各層間剪力的最大值，如表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 由于各振型下層間剪力的最大值并不一定同時(shí)出現(xiàn)，因此采用平方和開方的方式，求得地震波作用下層間剪力絕對(duì)值的最大值，如表3.4所示。表3.4 各層間剪力的最大值單位：kN編號(hào)123456785166.0 4824.7 4423.4 4124.5 3894.5 3

13、537.0 2854.8 1658.6 從表3.4可以得到利用振型分解反響譜法求得的在El Centro波NS向作用下各層間剪力絕對(duì)值的包絡(luò)圖，如圖3.6a所示，3.6b為上文中求得的時(shí)程計(jì)算得到的各層間剪力絕對(duì)值的包絡(luò)圖。圖3.6 兩種方法計(jì)算EL-Centro波作用下各層層間剪力絕對(duì)值包絡(luò)圖通過(guò)比照可以發(fā)現(xiàn)時(shí)程計(jì)算的結(jié)果與振型分解反響譜法的計(jì)算結(jié)果在底層吻合的很好，在結(jié)構(gòu)高程的差距相對(duì)較大，但整體上兩者的結(jié)果是在同一數(shù)量級(jí)，且差值合理的，最大差值為26.9%?？紤]到振型分解反響譜法本身上是一種為設(shè)計(jì)效勞求得最大地震作用的近似參考方法，所以通過(guò)與本作業(yè)中編程進(jìn)行時(shí)程計(jì)算的結(jié)果相比，在一定程度

14、上從側(cè)面印證了前文時(shí)程計(jì)算結(jié)果的合理性與正確性。四、 人工地震波合成本文采用經(jīng)典譜表達(dá)方法來(lái)進(jìn)行強(qiáng)震加速度模擬。經(jīng)典譜表達(dá)方法如式17所示。17式中：為時(shí)間包絡(luò)函數(shù)，表達(dá)式如下所示：18相關(guān)參數(shù)選擇按照表4.1選?。罕?.1 時(shí)間包絡(luò)函數(shù)的參數(shù)值及持時(shí)為主振平穩(wěn)段段首時(shí)間，對(duì)設(shè)計(jì)的兩條人工波假定分別對(duì)應(yīng)二類及三類場(chǎng)地的第一組地震波，可分別取0.8s以及1.2s；為主振平穩(wěn)段段尾時(shí)間，對(duì)兩條人工波分別取7s以及9s；為衰減系數(shù)，對(duì)兩條人工波分別取0.35s-1以及0.25 s-1。式17中的可由下式計(jì)算。19即為頻域樣本點(diǎn)，為對(duì)應(yīng)樣本點(diǎn)的相位，為所選用的功率譜。滿足：20其中為截?cái)囝l率，為頻率下

15、界，為頻域的樣本點(diǎn)數(shù)目。所以，在上式的根底上，需要功率譜密度函數(shù)模型以及相位譜，才能夠生成人工地震波。而根據(jù)作業(yè)要求， 選用功率譜密度函數(shù)模型Clough-Penzien譜模型，其表達(dá)式為：21式中相關(guān)參數(shù)選擇按照表4.2選?。罕?.2 Clough-Penzien模型場(chǎng)地土參數(shù)設(shè)計(jì)值為場(chǎng)地卓越周期，對(duì)設(shè)計(jì)的兩條人工波假定分別對(duì)應(yīng)二類及三類場(chǎng)地的第一組地震波，故分別取20.94與13.96；為二次過(guò)濾參數(shù)，取為0.15一般可取0.1-0.2；為場(chǎng)地地基阻尼，對(duì)兩條地震波分別取為0.72以及0.8；為二次過(guò)濾參數(shù)，取=；為譜強(qiáng)度因子，在本文中，為使兩條人工地震波的加速度峰值在0.4g左右，經(jīng)過(guò)調(diào)

16、整試算，分別取0.025、0.032。地震動(dòng)持時(shí)采用相對(duì)分?jǐn)?shù)持時(shí)，即地震動(dòng)加速度記錄最后一次到達(dá)峰值加速度0.05倍時(shí)所對(duì)應(yīng)的時(shí)間。上述參數(shù)的選擇參考文獻(xiàn)：張猛,張哲,李天. 與標(biāo)準(zhǔn)反響譜相對(duì)應(yīng)Clough-Penzien模型參數(shù)研究J. 世界地震工程,2007,01:56-60.兩條地震波對(duì)應(yīng)的Clough-Penzien模型如圖4.1所示。(a)人工波1(b)人工波2圖4.1人工波Clough-Penzien模型功率譜圖相位譜那么采用隨機(jī)方式生成隨機(jī)相位譜。故綜上利用matlab編制人工地震波合成程序程序見附，可以得到兩條人工地震波地震波數(shù)據(jù)見附，兩條波的加速度時(shí)程如圖4.2所示。 (a)

17、人工波1(b)人工波2圖4.2 人工波加速度時(shí)程五、 人工波時(shí)程計(jì)算結(jié)果將四中人工生成的兩條地震波用前文的方法進(jìn)行時(shí)程計(jì)算，得到了響應(yīng)的計(jì)算結(jié)果。人工波1作用下，結(jié)構(gòu)頂層位移時(shí)程如圖5.1所示。圖5.1 人工波1作用下結(jié)構(gòu)頂層位移時(shí)程人工波1作用下，結(jié)構(gòu)各層層間位移時(shí)程如圖5.2所示。圖5.2 人工波1作用下結(jié)構(gòu)各層層間位移時(shí)程在人工波1作用下，結(jié)構(gòu)底層剪力時(shí)程如圖5.3所示。圖5.3 人工波1作用下結(jié)構(gòu)底層剪力時(shí)程在人工波1作用下，各層的層間位移與層間剪力絕對(duì)值包絡(luò)圖如圖5.4所示。(a)層間位移絕對(duì)值包絡(luò)圖b層間剪力絕對(duì)值包絡(luò)圖圖5.4 人工波1作用下結(jié)構(gòu)各層層間位移與層間剪力絕對(duì)值包絡(luò)圖

18、人工波2作用下，結(jié)構(gòu)頂層位移時(shí)程如圖5.5所示。圖5.5 人工波2作用下結(jié)構(gòu)頂層位移時(shí)程人工波2作用下，結(jié)構(gòu)各層層間位移時(shí)程如圖5.6所示。圖5.6 人工波2作用下結(jié)構(gòu)各層層間位移時(shí)程在人工波2作用下，結(jié)構(gòu)底層剪力時(shí)程如圖5.7所示。圖5.7 人工波2作用下結(jié)構(gòu)底層剪力時(shí)程在人工波2作用下，各層的層間位移與層間剪力絕對(duì)值包絡(luò)圖如圖5.8所示。(a)層間位移絕對(duì)值包絡(luò)圖 b層間剪力絕對(duì)值包絡(luò)圖圖5.8 人工波2作用下結(jié)構(gòu)各層層間位移與層間剪力絕對(duì)值包絡(luò)圖1、程序在2021Ra上可以正常運(yùn)行2、運(yùn)行時(shí)需要把ThreeKindsWaves復(fù)制到當(dāng)前工作路徑下3、文件的存儲(chǔ)地址可以根據(jù)自身情況進(jìn)行修改

19、，或者不需要可以直接注釋%利用解特征方程的方法得到結(jié)構(gòu)的8階自振周期和振型，并歸一化，然后利用循環(huán)產(chǎn)生陣型圖，保存在文件陣型圖內(nèi)tic %計(jì)時(shí)開始clear all;%去除變量digits(4) %設(shè)定精度M=105*diag(3.4 3.4 3.2 3.2 2.8 2.8 2.7 2.6);%輸入質(zhì)量矩陣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內(nèi)置函數(shù)解特征方程，周期平方存為omiga，振型存為faiomiga=omiga.0.5; %得到周期transfer=fai(:,8);fai(:,8)=fai(:,7);fai(:,7)=transfer;file_id='C:UserslenovoDesktop第二次作業(yè)用圖振型圖'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');%利用內(nèi)置插值函數(shù)使連線光滑 plot(pp,yi,'k-') plot(fai(:,i),0.5:0.5:4,'o','MarkerFaceColor','k','

22、MarkerSize',20)%控制點(diǎn)的形狀和大小代表集中質(zhì)量 text(fai(:,i)+0.2,0.6:0.5:4.1,num2str(fai(:,i)%標(biāo)注振型數(shù)值 for j=1:8 plot(0 fai(j,i),0.5*j 0.5*j,'k-')%畫相關(guān)的輔助線 end name=strcat('第',num2str(9-i),'階振型圖'); title(name,'FontSize',16,'position',0 -0.3)%標(biāo)注圖名 F=getframe(figure(i); imwri

23、te(F.cdata,file_id,strcat(name,'.png'); %保存振型圖片 endf=f'T=1./f;%輸出周期和頻率Toc%通過(guò)解特征方程得到振型和周期，并進(jìn)行歸一化，為下面進(jìn)行時(shí)程分析做鋪墊。讀取Excel表格存儲(chǔ)的地震動(dòng)數(shù)據(jù)，%利用紐馬克方法進(jìn)行時(shí)程分析得到地面速度和地面位移，利用相關(guān)公式對(duì)數(shù)據(jù)進(jìn)行基線校正得到更為合理的數(shù)據(jù)。%利用紐馬克方法進(jìn)行時(shí)程分析得到結(jié)構(gòu)各階振型的響應(yīng)，然后通過(guò)公式X=A*u得到結(jié)構(gòu)總的位移響應(yīng)，然后利用%結(jié)構(gòu)力學(xué)相關(guān)知識(shí)得到剪力和剪力包絡(luò)圖。tic;clear all ;clf ; %計(jì)時(shí)開始，并清楚其他無(wú)關(guān)變量M=

24、105*diag(3.4 3.4 3.2 3.2 2.8 2.8 2.7 2.6); %輸入質(zhì)量矩陣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); %計(jì)算振型參與系數(shù)kesai=(0.1329+0.0067*omiga.2)/2./omiga; %計(jì)算各階阻尼file_id='C:UserslenovoDesktop第二次作業(yè)用圖地震動(dòng)響應(yīng)圖'mkdir(file_id); %創(chuàng)立文件夾保存生成的所有圖A=xlsread('ThreeKindsWaves','ELcentro-NS'); %讀取地震波數(shù)據(jù)N=length(A(:,1); mint=A(2,1)-

26、A(1,1); %求時(shí)間間隔T=(N-1)*mint;allmotion=zeros(N,8); %用來(lái)存儲(chǔ)所有時(shí)刻的位移allacceleration=zeros(N,8); %用來(lái)存儲(chǔ)所有時(shí)刻的加速度值allacceleration(1,:)=A(1,2); %對(duì)0時(shí)刻加速度矩陣初始化allvelocity=zeros(N,8); %用于存儲(chǔ)所有時(shí)刻速度vg(1)=0;xg(1)=0;for m=1:N-1 %計(jì)算場(chǎng)地速度、位移用于對(duì)加速度進(jìn)行基線修正 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); %計(jì)算系數(shù)c1 c0=(vg(N)-c1*T2/2)/T; %計(jì)算系數(shù)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); %修正后的場(chǎng)地位移 for i=1:8

28、%第三層循環(huán)用于對(duì)不同周期值計(jì)算反響S=1+kesai(i)*omiga(i)*mint+(omiga(i)*mint)2)/6; for j=2:N %對(duì)每個(gè)周期利用NewMark方法進(jìn)行數(shù)值積分得到對(duì)應(yīng)值具體請(qǐng)查看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 %積分計(jì)算反響并存放在前面的空矩陣?yán)?endmotion1=allmotion*fai' %計(jì)算位移acceleration=allacce

30、leration*fai' %計(jì)算加速度velocity=allvelocity*fai' %計(jì)算速度f(wàn)or n=1:9figure('color','white') %利用figure函數(shù)準(zhǔn)確的控制畫圖窗口if n=9axis(0,50,-0.03,0.03) %控制坐標(biāo)軸范圍，便于比較endgrid on hold onbox off%annotation('arrow',0.132 0.132,0.8 1);%annotation('arrow',0.8 1,0.108 0.108); %產(chǎn)生坐標(biāo)軸箭頭可能不

31、夠美觀if n=1 %計(jì)算相關(guān)數(shù)據(jù)和得到圖名 name='底層層間位移時(shí)程' motion(:,n)=motion1(:,n);elseif n=9 name='頂層位移時(shí)程' motion(:,9)=motion1(:,8)+xg;else motion(:,n)=motion1(:,n)-motion1(:,n-1); name=strcat('第',num2str(n-1),'-',num2str(n),'層間位移時(shí)程');endplot(A(:,1),motion(:,n),'linewidth&#

32、39;,2) %畫相對(duì)位移時(shí)程xlabel('時(shí)間s','FontName','宋體','FontSize',16); %x軸標(biāo)注ylabel('位移(m)','FontName','宋體','FontSize',16); %y軸標(biāo)注title(name,'FontName','宋體','FontSize',20) %標(biāo)注圖名set(gcf,'position',get(0,'screensize

33、'); %圖形全屏，便于查看shg F=getframe(gcf); imwrite(F.cdata,file_id,strcat(name,'.png'); %存儲(chǔ)得到的時(shí)程圖 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); %計(jì)算頂層的層間相對(duì)位移Fq=motion(:,1:8).*GK; %計(jì)算剪力figure('color','white') %生成圖形窗口，背景白色plo

34、t(A(:,1),Fq(:,1) %畫底層剪力時(shí)程圖xlabel('時(shí)間s','FontName','宋體','FontSize',16); %x軸標(biāo)注ylabel('剪力N','FontName','宋體','FontSize',16); %y軸標(biāo)注title('底層剪力時(shí)程','FontName','宋體','FontSize',20); %標(biāo)注圖名set(gcf,'position'

35、,get(0,'screensize'); %圖形全屏，便于查看 F=getframe(gcf); %存儲(chǔ)生成的圖形 imwrite(F.cdata,file_id,strcat('底層剪力時(shí)程圖','.png'); maxFq=zeros(1,16);for f=1:8 %計(jì)算所需剪力和位移包絡(luò)圖數(shù)據(jù)，并得到便于畫圖的矩陣 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-');%生成層間位移包絡(luò)圖title('層間位移絕對(duì)值最大值包絡(luò)圖','FontName','宋體','FontSize',20); %標(biāo)注圖名text(maxmotion(1:2:16),0.25:0.5:4,num2str(maxmotion(1:2:16)

37、9;) %標(biāo)注位移數(shù)據(jù)set(gcf,'position',get(0,'screensize'); %圖形全屏，便于查看F=getframe(gcf); imwrite(F.cdata,file_id,strcat('層間位移絕對(duì)值最大值包絡(luò)圖','.png'); %存儲(chǔ)圖形 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、層間剪力絕對(duì)值最大值包絡(luò)圖','FontName','宋體','FontSize',20); %標(biāo)注圖名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('層間剪力絕對(duì)

39、值最大值包絡(luò)圖','.png'); toc %計(jì)時(shí)結(jié)束%利用地震波的經(jīng)典表達(dá)方式和Clough-Penzien功率譜表達(dá)式計(jì)算對(duì)應(yīng)功率以及振幅和功率譜的關(guān)系得到振幅，%利用隨機(jī)函數(shù)得到相位，進(jìn)而得到各類場(chǎng)地地震動(dòng)數(shù)據(jù)。并把生成的地震動(dòng)數(shù)據(jù)和相關(guān)圖形保存function Artificial_wavetic %計(jì)時(shí)開始clear all; %去除無(wú)關(guān)變量 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);%生成時(shí)間矩陣，便于后期畫圖t1=0.8 1.2; %輸入強(qiáng)度函數(shù)相關(guān)參數(shù)t2=7 9;c=0.35 0.25;file_id='C:UserslenovoDesktop第二次作業(yè)用圖人造波'mkdir(file_id); %創(chuàng)立文件夾保存生成的所有圖和地震動(dòng)數(shù)據(jù)for i=1:2 %利用循環(huán)以此計(jì)算兩類地震動(dòng)數(shù)據(jù) for j=1:10000 %利用循環(huán)計(jì)算功率譜 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); %計(jì)算對(duì)應(yīng)振幅 end %功率譜計(jì)算結(jié)束 fai=repmat(2*pi*rand(10000,1)-pi,1,2001);%生成相位矩陣 for n=1:2001 %計(jì)算強(qiáng)度函數(shù)值 gt(n)=g(t(n),t1(i),t2(i),c(i); end acceleration=A*cos(w'*t+fai).*gt; %計(jì)算地震動(dòng)數(shù)據(jù) N= find(abs(acceleration)>

42、=0.02*max(abs(acceleration);%尋找截?cái)帱c(diǎn) figure(2*i-1); %利用figure函數(shù)控制圖形窗口 plot(w,S); %得到功率譜圖 name=strcat(num2str(i+1),'類場(chǎng)地人工波功率譜'); title(name); set(gcf,'position',get(0,'screensize'); %圖形全屏，便于查看 F=getframe(gcf); imwrite(F.cdata,file_id,strcat(name,'.png'); %存儲(chǔ)功率譜圖形 figure(

43、2*i); plot(t(1:N(end),acceleration(1:N(end); %畫地震動(dòng)時(shí)程圖 name=strcat(num2str(i+1),'類場(chǎng)地人工波加速度時(shí)程'); title(name); set(gcf,'position',get(0,'screensize'); %圖形全屏，便于查看 F=getframe(gcf); imwrite(F.cdata,file_id,strcat(name,'.png'); %存儲(chǔ)地震動(dòng)時(shí)程圖 xlswrite(strcat(file_id,'第',n

44、um2str(i+1),'類場(chǎng)地人工波.xls'),t(1:N(end);acceleration(1:N(end)');%存儲(chǔ)地震動(dòng)數(shù)據(jù)endendfunction gt=g(t,t1,t2,c)%強(qiáng)度包絡(luò)函數(shù)計(jì)算 if t<t1 gt=(t/t1)2; elseif t>t2 gt=exp(-c*(t-t2); else gt=1; endend%該代碼與Earthquake_Response.m根本一致故不進(jìn)行詳細(xì)標(biāo)注tic;clear all ;clf ; %計(jì)時(shí)開始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第二次作業(yè)用圖人工地震動(dòng)響應(yīng)圖'mkdir(file_id); %創(chuàng)立文件夾保存生成的所有圖B='第2類場(chǎng)地人工波','第3類場(chǎng)地人工波'for a=1:2clear A N mint TA=xlsread('ThreeKindsWaves',Ba); %讀取地震波數(shù)據(jù)A(:,2)=-A(:,2);N=length(A(:,1); mint=A(2,1)-A(1,1); %求時(shí)間間隔T=(N-1)*mi

47、nt;allmotion=zeros(N,8); %用來(lái)存儲(chǔ)所有時(shí)刻的位移allacceleration=zeros(N,8); %用來(lái)存儲(chǔ)所有時(shí)刻的加速度值allacceleration(1,:)=A(1,2); %對(duì)0時(shí)刻加速度矩陣初始化allvelocity=zeros(N,8); %用于存儲(chǔ)所有時(shí)刻速度motion1=zeros(N,8);motion=zeros(N,8);vg=zeros(1,N);xg=zeros(1,N);for m=1:N-1 %計(jì)算場(chǎng)地速度、位移用于對(duì)加速度進(jìn)行基線修正 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 %第三層循環(huán)用于對(duì)不同周期值計(jì)算反響S=1+kesai(i)*omiga(i)*mint+(omiga(i)*mint)2)/6; for j=2:N %對(duì)每個(gè)周期利用NewMark方法進(jìn)行數(shù)值積分得到對(duì)應(yīng)值 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