地震工程學(xué)反應(yīng)譜和地震時(shí)程波的相互轉(zhuǎn)化編程

1、地震工程學(xué)作業(yè)課程名稱：地震工程學(xué)指導(dǎo)老師:翟永梅姓名:史先飛學(xué)號(hào)：1232627一、地震波生成反應(yīng)譜1所取的地震波為Elcentro地震波加速度曲線，如圖1所示。圖1 Elcentro地震波加速度曲線2所調(diào)用的Mat lab程序?yàn)?% * 讀入地震記錄 *ElCentro;Accelerate ElCentro(:,1)*;%單位統(tǒng)一為 m 和 sN=length(Accelerate);%N 讀入的記錄的量time=0:(N-1)*; %單位 s%初始化各儲(chǔ)存向量Displace=zeros(1,N); %相對(duì)位移Velocity=zeros(1,N); %相對(duì)速度AbsAcce=zero

2、s(1,N); %絕對(duì)加速度% *A,B 矩陣*Damp=; %阻尼比TA=:6; %TA=:6; %結(jié)構(gòu)周期Dt=; %地震記錄的步長(zhǎng)%記錄計(jì)算得到的反應(yīng)，MaxD為某阻尼時(shí)最大相對(duì)位移，MaxV為某阻尼最大相對(duì)速度，MaxA某阻尼時(shí)最大 絕對(duì)加速度，用于畫圖MaxD=zeros(3,length(TA);MaxV=zeros(3,length(TA);MaxA=zeros(3,length(TA);t=1;for T=:6NatualFrequency=2*pi/T ; %結(jié)構(gòu)自振頻率DampFrequency=NatualFrequency*sqrt(1-Damp*Damp); %計(jì)算公

3、式化簡(jiǎn)e_t=exp(-Damp*NatualFrequency*Dt);s=sin(DampFrequency*Dt);c=cos(DampFrequency*Dt);A=zeros(2,2);A(1,1)=e_t*(s*Damp/sqrt(1-Damp*Damp)+c);A(1,2)=e_t*s/DampFrequency;A(2,1)=-NatualFrequency*e_t*s/sqrt(1-Damp*Damp);A(2,2)=e_t*(-s*Damp/sqrt(1-Damp*Damp)+c);d_f=(2*Damp2T)/(NatualFrequency2*Dt);d_3t=Damp

4、/(NatualFrequency3*Dt);B=zeros(2,2);B(1,1)=e_t*(d_f+Damp/NatualFrequency)*s/DampFrequency+(2*d_3t+1/NatualFrequency2)*c)-2*d_3t;B(1,2)=-e_t*(d_f*s/DampFrequency+2*d_3t*c)T/NatualFrequency2+2*d_3t;B(2,1)=e_t*(d_f+Damp/NatualFrequency)*(c-Damp/sqrt(1-Damp2)*s)-(2*d_3t+1/NatualFrequency”2)*(DampFrequen

5、cy*s+Damp*NatualFrequency*c)+1/(NatualFrequency2*Dt);B(2,2)=e_t*(1/(NatualFrequency2*Dt)*c+s*Damp/(NatualFrequency*DampFrequency*Dt)T/(NatualF requency2*Dt);for i=1:(N-1) %根據(jù)地震記錄，計(jì)算不同的反應(yīng)Displace(i+1)=A(1,1)*Displace(i)+A(1,2)*Velocity(i)+B(1,1)*Accelerate(i)+B(1,2)*Accelerate(i+1);Velocity(i+1)=A(2,

6、1)*Displace(i)+A(2,2)*Velocity(i)+B(2,1)*Accelerate(i)+B(2,2)*Accelerate(i+1);AbsAcce(i+1)=-2*Damp*NatualFrequency*Velocity(i+1)-NatualFrequency2*Displace(i+1);endMaxD(1,t)=max(abs(Displace);MaxV(1,t)=max(abs(Velocity);if T=MaxA(1,t)=max(abs(Accelerate);elseMaxA(1,t)=max(abs(AbsAcce);endDisplace=zer

7、os(1,N);%初始化各儲(chǔ)存向量，避免下次不同周期計(jì)算時(shí)引用到前一個(gè)周期的結(jié)果Velocity=zeros(1,N);AbsAcce=zeros(1,N);t=t+1;End% *PLOT*close allfigure %繪制地震記錄圖plot(time(:),Accelerate(:)title(PEER STRONG MOTION DATABASE RECORD)xlabel(time(s)ylabel(acceleration(g)gridfigure %繪制位移反應(yīng)譜plot(TA,MaxD(1,:),TA,MaxD(2,:),-r,TA,MaxD(3,:),:k) title(D

8、isplacement)xlabel(Tn(s)ylabel(Displacement(m)legend(Z=)Grid figure %繪制速度反應(yīng)譜plot(TA,MaxV(1,:),TA,MaxV(2,:),-r,TA,MaxV(3,:),:k)title(Velocity)xlabel(Tn(s)ylabel(velocity(m/s)legend(Z=)Gridfigure %繪制絕對(duì)加速度反應(yīng)譜plot(TA,MaxA(1,:),TA,MaxA(2,:),-r,TA,MaxA(3,:),:k)title(Absolute Acceleration)xlabel(Tn(s)ylabe

9、l(absolute acceleration(m/s2) legend(Z=)Grid3運(yùn)行的結(jié)果得到的反應(yīng)譜圖2位移反應(yīng)譜圖3速度反應(yīng)譜圖4加速度反應(yīng)譜一、反應(yīng)譜生成地震波1所取的反應(yīng)譜為上海市設(shè)計(jì)反應(yīng)譜圖5上海市設(shè)計(jì)反應(yīng)譜2反應(yīng)譜取值程序?yàn)椋?規(guī)范反應(yīng)譜取值程序參照01年抗震規(guī)范function rs_z=r_s_1(pl,zn,ld,cd,fz) %pl 圓頻率,zn 阻尼比，ld 烈度,cd 場(chǎng)地類型，場(chǎng)地分組 fz%烈度選擇if ld=6arfmax=;endif ld=7arfmax=;endif ld=8arfmax=;endif ld=9arfmax=;end%場(chǎng)地類別，設(shè)計(jì)

10、地震分組選擇if cd=1if fz=1Tg=;endif fz=2Tg=;endif fz=3Tg=;endendif cd=2if fz=1Tg=;endif fz=2Tg=;endif fz=3Tg=;endendif cd=3if fz=1Tg=;endif fz=2Tg=;endif fz=3Tg=;endendif cd=4if fz=1Tg=;end if fz=2Tg=;end if fz=3Tg=；endend%ceita=zn; % 阻尼比lmt1=+/8;if lmt10lmt1=0;endlmt2=1+/+*ceita);if lmt2lmt2=;endsjzs=+/+

11、5*ceita);%分段位置T1 T2 T3T1 = ;T2=Tg;T3=5*Tg;T_jg=2*pi./pl;%第一段0T1if T_jg=T1arf_jg=*arfmax+(lmt2*arfmax)/*T_jg;end%第二段T1T2if T1T_jg&T_jg=T2arf_jg=lmt2*arfmax;end%第三段T2T3if T2T_jg&T_jg=T3arf_jg=(Tg/T_jg)sjzs)*lmt2*arfmax;end%第四段T3if T3T_jg&T_jg=arf_jg=(lmt2*sjzs-lmt1*(T_jg-5*Tg)*arfmax;end%第五段if T_jgarf

12、_jg=(lmt2*sjzs-lmt1*Tg)*arfmax;end%反應(yīng)譜值擬加速度值rs_z=arf_jg*;end3生成人造地震波主程序：%主程序%確定需要控制的反應(yīng)譜Sa(T)(T=T1TM)的坐標(biāo)點(diǎn)數(shù)M，反應(yīng)譜控制容差rcTyz=:,:,:;rc=;nTyz=length(Tyz);ceita=;% 阻尼比：for i=1:nTyzSyz =r_s_1(2*pi/Tyz(i),ceita,8,2,1); %8 度，2 類場(chǎng)地，第 1 地震分組end%變換的頻率差：2*pi*(可以保證長(zhǎng)周期項(xiàng)5s附近有5項(xiàng)三角級(jí)數(shù))；% 頻率變化范圍 N1=30, 30*2*pi ;N2=3000,

13、5000*2*piplc=2*pi*;pl=30*2*pi:*2*pi:10000*2*pi;npl = length(pl);P=; %保證率%人造地震動(dòng)持續(xù)時(shí)間40s，時(shí)間間隔：Td=40;dt=;t=0:40;nt=length(t);%衰減包絡(luò)函數(shù)t1=8; %上升段t2=8+24; %平穩(wěn)段；下降段則為40 32 = 8sc=; %衰減段參數(shù)for i=1:ntif t(i)t1 & t(i)=t2f(i)=exp(-c*(t -t2);endend%反應(yīng)譜轉(zhuǎn)換功率譜for i=1:nplSw(i) = (2*ceita/(pi*pl(i)*r_s_1(pl(i),ceita,8,2

14、,1)2/(-2*log(T*pi*log(P)/(pl(i)*Td);Aw=sqrt(4*Sw (i) *plc);end%合成地震動(dòng)at=zeros(nt,1);atj=zeros(nt,1);for i=1:nplfai(i)=rand (1)*2*pi;for j=1:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i);endat=at+atj;end%計(jì)算反應(yīng)譜驗(yàn)證是否滿足rc在5%的要求，需要時(shí)程動(dòng)力分析% response spectra of callidar% parameterg=;m=1;x0=0;v0=0;w

15、w=2*pi./Tyz;% loadag=at;% 修改% solutionfor y=1:nTyzz=;w=ww(y);c=2*z*w;k=w2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m(-ag(i)-c*v0-k*x0);kk=k+(dt2)(6*m)+dt(3*c);pp=p(i)+m*(dt(6*v0)+3*a0)+c*(3*v0+2(dt*a0);dx=kkpp;dv=dt(3*dx)-3*v0-2(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i

16、)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i);end%比較容差for i=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i)/max(Syz(:);endjsnum=1;while max(rcrsp(:)rc%循環(huán)體函數(shù)blxs=Syz./rspa;for xsxs=1:nplif 2*pi/pl(xsxs)=Tyz(sxsx) & (2*pi/pl(xsxs)Tyz(nTyz)blxs1(xsxs)=blxs(nTyz);endendAw=

17、Aw.*blxs1;%合成地震動(dòng)at=zeros(nt,1);atj=zeros(nt,1);for i=1:nplfor j=1:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i);endat=at+atj;end%計(jì)算反應(yīng)譜驗(yàn)證是否滿足rc在5%的要求% response spectra of callidar% parameterg=;m=1;x0=0;v0=0;ww=2*pi./Tyz;% loadag=at;% 修改% solutionfor y=1:nTyzz=;w=ww(y);c=2*z*w;k=w2;for i=1:

18、nt-1p(i)=-ag(i+1)+ag(i);a0=m(-ag(i)-c*v0-k*x0);kk=k+(dt2)(6*m)+dt(3*c);pp=p(i)+m*(dt(6*v0)+3*a0)+c*(3*v0+2(dt*a0);dx=kkpp;dv=dt(3*dx)-3*v0-2(dt*a0);x1=x0+dx;x0=x1;v1=v0+dv;v0=v1;as(i)=a0;as(i)=as(i)+ag(i);vs(i)=v0;xs(i)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);endfor i=1:nTyzrspa(i)=maxas(i);end%比較容差for i=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i)/max(Syz(:);endjsnum=jsnum+1max(rcrsp(:)end%最終的反應(yīng)譜與規(guī)范譜% response spectra of