地震工程學-反應譜和地震時程波的相互轉化matlab編程

地震工程學作業(yè)課程名稱:地震工程學______指導老師:_______翟永梅_________姓名:史先飛________學號:1232627________一、地震波生成反應譜1所取的地震波為Elcentro地震波加速度曲線，如圖1所示。圖1Elcentro地震波加速度曲線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);%相對位移Velocity=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;%結構周期Dt=0.005;%地震記錄的步長%記錄計算得到的反應，MaxD為某阻尼時最大相對位移，MaxV為某阻尼最大相對速度，MaxA某阻尼時最大絕對加速度，用于畫圖MaxD=zeros(3,length(TA));MaxV=zeros(3,length(TA));MaxA=zeros(3,length(TA));t=1;forT=0.0:0.05:6NatualFrequency=2*pi/T;%結構自振頻率DampFrequency=NatualFrequency*sqrt(1-Damp*Damp);%計算公式化簡反應譜生成地震波所取的反應譜為上海市設計反應譜圖5上海市設計反應譜反應譜取值程序為：%%規(guī)范反應譜取值程序參照01年抗震規(guī)范functionrs_z=r_s_1(pl,zn,ld,cd,fz)%%%pl圓頻率,zn阻尼比,ld烈度,cd場地類型,場地分組fz%%%%烈度選擇ifld==6arfmax=0.11;endifld==7arfmax=0.23;endifld==8arfmax=0.45;endifld==9arfmax=0.90;end%%%%場地類別，設計地震分組選擇ifcd==1iffz==1Tg=0.25;endiffz==2Tg=0.30;endiffz==3Tg=0.35;endendifcd==2iffz==1Tg=0.35;endiffz==2Tg=0.40;endiffz==3Tg=0.45;endendifcd==3iffz==1Tg=0.45;endiffz==2Tg=0.55;endiffz==3Tg=0.65;endendifcd==4iffz==1Tg=0.65;endiffz==2Tg=0.75;endiffz==3Tg=0.90;endend%%%%%%%%%ceita=zn;%%%%%阻尼比lmt1=0.02+(0.05-ceita)/8;iflmt1<0lmt1=0;endlmt2=1+(0.05-ceita)/(0.06+1.7*ceita);iflmt2<0.55lmt2=0.55;endsjzs=0.9+(0.05-ceita)/(0.5+5*ceita);%%%%%分段位置T1T2T3T1=0.1;T2=Tg;T3=5*Tg;T_jg=2*pi./pl;%%%%第一段0～T1ifT_jg<=T1arf_jg=0.45*arfmax+(lmt2*arfmax-0.45*arfmax)/0.1*T_jg;end%%%%第二段T1～T2ifT1<T_jg&T_jg<=T2arf_jg=lmt2*arfmax;end%%%%第三段T2~T3ifT2<T_jg&T_jg<=T3arf_jg=((Tg/T_jg)^sjzs)*lmt2*arfmax;end%%%%第四段T3～6.0ifT3<T_jg&T_jg<=6.0arf_jg=(lmt2*0.2^sjzs-lmt1*(T_jg-5*Tg))*arfmax;end%%%%第五段6.0～if6.0<T_jgarf_jg=(lmt2*0.2^sjzs-lmt1*(6.0-5*Tg))*arfmax;end%%%%%%反應譜值擬加速度值rs_z=arf_jg*9.8;end生成人造地震波主程序：%%%主程序%%%%%%%%確定需要控制的反應譜Sa（T）(T=T1,...,TM)的坐標點數(shù)M，反應譜控制容差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.035fori=1:nTyzSyz(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1);%%%%8度，2類場地，第1地震分組end%%%%%%變換的頻率差：2*pi*0.005(可以保證長周期項5s附近有5項三角級數(shù))；%%%%頻率變化范圍N1=30,30*0.005*2*pi;N2=3000,5000*0.005*2*piplc=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);%%%%%%%衰減包絡函數(shù)t1=8;%%%%上升段t2=8+24;%%%%%平穩(wěn)段;下降段則為40－32＝8sc=0.6;%%%%衰減段參數(shù)fori=1:ntift(i)<=t1f(i)=(t(i)/t1)^2;endift(i)>t1&t(i)<t2f(i)=1;endift(i)>=t2f(i)=exp(-c*(t(i)-t2));endend%%%%%%%反應譜轉換功率譜fori=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(nt,1);atj=zeros(nt,1);fori=1:nplfai(i)=rand(1)*2*pi;forj=1:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i))));endat=at+atj;end%%%%%%%計算反應譜驗證是否滿足rc在5%的要求,需要時程動力分析%%%%%%%%%%%%responsespectraofcallidar%%%%%%%parameterg=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%%loadag=at;%%%%%%%修改%%%%%%%solutionfory=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;fori=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0));dx=kk\pp;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);endfori=1:nTyzrspa(i)=maxas(i);end%%%%%%%比較容差fori=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i))/max(Syz(:));endjsnum=1;whilemax(rcrsp(:))>rc%%%%%循環(huán)體函數(shù)blxs=Syz./rspa;forxsxs=1:nplif2*pi/pl(xsxs)<Tyz(1)blxs1(xsxs)=blxs(1);endforsxsx=1:nTyz-1if(2*pi/pl(xsxs)>=Tyz(sxsx))&(2*pi/pl(xsxs)<=Tyz(sxsx+1))blxs1(xsxs)=blxs(sxsx)+(blxs(sxsx+1)-blxs(sxsx))*(2*pi/pl(xsxs)-Tyz(sxsx))/(Tyz(sxsx+1)-Tyz(sxsx));endendif2*pi/pl(xsxs)>Tyz(nTyz)blxs1(xsxs)=blxs(nTyz);endendAw=Aw.*blxs1;%%%%%%%%%%%%%%合成地震動at=zeros(nt,1);atj=zeros(nt,1);fori=1:nplforj=1:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i))));endat=at+atj;end%%%%%%%計算反應譜驗證是否滿足rc在5%的要求%%%%%%%%%%%%responsespectraofcallidar%%%%%%%parameterg=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%%loadag=at;%%%%%%%修改%%%%%%%solutionfory=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;fori=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m\(-ag(i)-c*v0-k*x0);kk=k+(dt^2)\(6*m)+dt\(3*c);pp=p(i)+m*(dt\(6*v0)+3*a0)+c*(3*v0+2\(dt*a0));dx=kk\pp;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);endfori=1:nTyzrspa(i)=maxas(i);end%%%%%%%比較容差fori=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i))/max(Syz(:));endjsnum=jsnum+1max(rcrsp(:))end%%%%%%%最終的反應譜與規(guī)范譜%%%%%%%%%%%%responsespectraofcallidar%%%%%%%parameter%%Tjs=0.05:0.01:6;%%nTjs=length(Tjs);g=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%%loadag=at;%%%%%%%修改%%%%%%%solutionfory=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w^2;fori=1:nt-1p(i)=-ag(i+1)+ag(i);a0