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

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

2、city=zeros(1,N); %相對速度AbsAcce=zeros(1,N); %絕對加速度% *A,B矩陣 *Damp=0.02; %阻尼比 0.02TA=0.0:0.05:6; %TA=0.000001:0.02:6; %結(jié)構(gòu)周期Dt=0.005; % 地震記錄的步長%記錄計算得到的反應(yīng)，MaxD為某阻尼時最大相對位移，MaxV為某阻尼最大相對速度，MaxA某阻尼時最大絕對加速度，用于畫圖MaxD=zeros(3,length(TA);MaxV=zeros(3,length(TA);MaxA=zeros(3,length(TA);t=1;for T=0.0:0.05:6NatualFr

3、equency=2*pi/T ; %結(jié)構(gòu)自振頻率DampFrequency=NatualFrequency*sqrt(1-Damp*Damp); %計算公式化簡.精品文檔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

4、)=e_t*(-s*Damp/sqrt(1-Damp*Damp)+c);d_f=(2*Damp2-1)/(NatualFrequency2*Dt);d_3t=Damp/(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_3 t;B(1,2)=-e_t*(d_f*s/DampFrequency+2*d_3t*c)-1/NatualFrequency2+2*d_3t;B(2,1)=e_t*(d_f+Dam

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

6、ace(i)+A(1,2)*Velocity(i)+B(1,1)*Accelerate(i)+B(1,2)*Accelerate(i+1);Velocity(i+1)=A(2,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(Vel

7、ocity);if T=0.0MaxA(1,t)=max(abs(Accelerate);elseMaxA(1,t)=max(abs(AbsAcce);endDisplace=zeros(1,N);% 初始化各儲存向量，避免下次不同周期計算時引用到前一個周期的結(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')xlabe

8、l('time(s)')ylabel('acceleration(g)')gridfigure %繪制位移反應(yīng)譜plot(TA,MaxD(1,:),'-.b',TA,MaxD(2,:),'-r',TA,MaxD(3,:),':k')title('Displacement')xlabel('Tn(s)')ylabel('Displacement(m)')legend('=0.02')Gridfigure %繪制速度反應(yīng)譜plot(TA,MaxV(1,:)

9、,'-.b',TA,MaxV(2,:),'-r',TA,MaxV(3,:),':k')title('Velocity')xlabel('Tn(s)')ylabel('velocity(m/s)')legend('=0.02')Gridfigure %繪制絕對加速度反應(yīng)譜plot(TA,MaxA(1,:),'-.b',TA,MaxA(2,:),'-r',TA,MaxA(3,:),':k')title('Absolute Accel

10、eration')xlabel('Tn(s)')ylabel('absolute acceleration(m/s2)')legend('=0.02')Grid3 運(yùn)行的結(jié)果得到的反應(yīng)譜.精品文檔圖 2位移反應(yīng)譜圖 3 速度反應(yīng)譜圖 4加速度反應(yīng)譜.精品文檔一、反應(yīng)譜生成地震波1 所取的反應(yīng)譜為上海市設(shè)計反應(yīng)譜圖 5 上海市設(shè)計反應(yīng)譜2 反應(yīng)譜取值程序為：%規(guī)范反應(yīng)譜取值程序參照 01 年抗震規(guī)范function rs_z=r_s_1(pl,zn,ld,cd,fz) %pl圓頻率 ,zn 阻尼比 ,ld烈度 ,cd 場地類型 , 場地分組

11、fz%烈度選擇if ld=6arfmax=0.11;endif ld=7arfmax=0.23;endif ld=8arfmax=0.45;endif ld=9arfmax=0.90;end%場地類別，設(shè)計地震分組選擇if cd=1if fz=1Tg=0.25;endif fz=2Tg=0.30;endif fz=3Tg=0.35;endendif cd=2if fz=1Tg=0.35;.精品文檔endif fz=2Tg=0.40;endif fz=3Tg=0.45;endendif cd=3if fz=1Tg=0.45;endif fz=2Tg=0.55;endif fz=3Tg=0.65;

12、endendif cd=4if fz=1Tg=0.65;endif fz=2Tg=0.75;endif fz=3Tg=0.90;endend%ceita=zn; %阻尼比lmt1=0.02+(0.05-ceita)/8;if lmt1<0lmt1=0;endlmt2=1+(0.05-ceita)/(0.06+1.7*ceita);if lmt2<0.55lmt2=0.55;endsjzs=0.9+(0.05-ceita)/(0.5+5*ceita);%分段位置 T1 T2 T3T1=0.1;T2=Tg;.精品文檔T3=5*Tg;T_jg=2*pi./pl;%第一段 0 T1if T

13、_jg<=T1arf_jg=0.45*arfmax+(lmt2*arfmax-0.45*arfmax)/0.1*T_jg;end%第二段 T1 T2if T1<T_jg&T_jg<=T2arf_jg=lmt2*arfmax;end%第三段 T2T3if T2<T_jg&T_jg<=T3arf_jg=(Tg/T_jg)sjzs)*lmt2*arfmax;end%第四段 T3 6.0if T3<T_jg&T_jg<=6.0arf_jg=(lmt2*0.2sjzs-lmt1*(T_jg-5*Tg)*arfmax;end%第五段 6.0

14、 if 6.0<T_jgarf_jg=(lmt2*0.2sjzs-lmt1*(6.0-5*Tg)*arfmax;end%反應(yīng)譜值擬加速度值rs_z=arf_jg*9.8;end3 生成人造地震波主程序：%主程序 %確定需要控制的反應(yīng)譜Sa（ T）(T=T1,.,TM)的坐標(biāo)點(diǎn)數(shù)M， 反應(yīng)譜控制容差rcTyz=0.04:0.016:0.1,0.15:0.05:3.0,3.2:0.05:5.0;rc=0.06;nTyz=length(Tyz);ceita=0.035;% 阻尼比： 0.035for i=1:nTyzSyz(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1);

15、%8度， 2 類場地，第1 地震分組end%變換的頻率差：2*pi*0.005(可以保證長周期項5s 附近有 5 項三角級數(shù) ) ；%頻率變化范圍N1=30, 30*0.005*2*pi ;N2=3000, 5000*0.005*2*pi.精品文檔plc=2*pi*0.005;pl=30*0.005*2*pi:0.005*2*pi:10000*0.005*2*pi;npl=length(pl);P=0.9; %保證率%人造地震動持續(xù)時間40s，時間間隔： 0.02sTd=40;dt=0.02;t=0:0.02:40;nt=length(t);%衰減包絡(luò)函數(shù)t1=8; %上升段t2=8+24;

16、%平穩(wěn)段 ;下降段則為40328sc=0.6; %衰減段參數(shù)for i=1:ntif t(i)<=t1f(i)=(t(i)/t1)2;endif t(i)>t1 & t(i)<t2f(i)=1;endif t(i)>=t2f(i)=exp(-c*(t(i)-t2);endend%反應(yīng)譜轉(zhuǎn)換功率譜for i=1:nplSw(i)=(2*ceita/(pi*pl(i)*r_s_1(pl(i),ceita,8,2,1)2/(-2*log(-1*pi*log(P)/(pl(i)*Td);Aw(i)=sqrt(4*Sw(i)*plc);end%合成地震動at=zeros(

17、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%計算反應(yīng)譜驗證是否滿足rc 在 5%的要求 , 需要時程動力分析% response spectra of callidar% parameterg=9.8;m=1;.精品文檔x0=0;v0=0;ww=2*pi./Tyz;% loadag=at;%修改% solutionfor y=1:nTyzz=0.037;w=ww(y);c

18、=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)=x0;endmaxas(y)=max(as);maxvs(y)=max(vs);maxxs(y)=max(xs);end

19、for 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(1)blxs1(xsxs)=blxs(1);endfor sxsx=1:nTyz-1if (2*pi/pl(xsxs)>=Tyz(sxsx) & (2*pi/pl(xsxs)<=Tyz(sxsx+1)blx

20、s1(xsxs)=blxs(sxsx)+(blxs(sxsx+1)-blxs(sxsx)*(2*pi/pl(xsxs)-Tyz(sxsx)/(Tyz(sxsx+1)-Tyz(sxsx);endendif 2*pi/pl(xsxs)>Tyz(nTyz)blxs1(xsxs)=blxs(nTyz);endendAw=Aw.*blxs1;%合成地震動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%

21、計算反應(yīng)譜驗證是否滿足rc 在 5%的要求% response spectra of callidar% parameterg=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;% loadag=at;%修改% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w2;for i=1:nt-1.精品文檔p(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)-Sy