地震工程学-反应谱和地震时程波的相互转化matlab编程

1、R 臨 +、fTONGJI UNIVERSITY地震工程学作业课程名称：地震工程学指导老师:翟永梅姓 名：史先飞学 号：1232627、地震波生成反应谱1所取的地震波为Elce ntro地震波加速度曲线，如图1所示。图1 Elce ntro地震波加速度曲线2所调用的Matlab程序为：% *读入地震记录 *ElCentro;Accelerate= EICentro(:,1)*9.8067;%单位统一为 m和 sN=length(Accelerate);%N读入的记录的量time=0:0.005:(N-1)*0.005; % 单位 s颊始化各储存向量Displace=zeros(1,N); %

2、相对位移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;for T=0.0:0.05:6Na

3、tualFrequency=2*pi/T ; % 结构自振频率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(

4、2,2)=e_t*(-s*Damp/sqrt(1-Damp*Damp)+c);d_f=(2*DampA2-1)/(NatualFrequencyA2*Dt);d_3t=Damp/(NatualFrequencyA3*Dt);B=zeros(2,2);B(1,1)=e_t*(d_f+Damp/NatualFrequency)*s/DampFrequency+(2*d_3t+1/NatualFrequencyA2)*c)-2*d_3 t;B(1,2)=-e_t*(d_f*s/DampFrequency+2*d_3t*c)-1/NatualFrequencyA2+2*d_3t;B(2,1)=e_t*

5、(d_f+Damp/NatualFrequency)*(c-Damp/sqrt(1-DampA2)*s)-(2*d_3t+1/NatualFrequencyA2 )*(DampFrequency*s+Damp*NatualFrequency*c)+1/(NatualFrequencyA2*Dt);B(2,2)=e_t*(1/(NatualFrequencyA2*Dt)*c+s*Damp/(NatualFrequency*DampFrequency*Dt)-1/(NatualF requencyA2*Dt);for i=1:(N-1) % 根据地震记录，计算不同的反应Displace(i+1)=

6、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,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)-NatualFrequencyA2*Displace(i+1);endMaxD(1,t)=max(abs(Displace);MaxV(1,t

7、)=max(abs(Velocity);if T=0.0MaxA(1,t)=max(abs(Accelerate);elseMaxA(1,t)=max(abs(AbsAcce);endDisplace=zeros(1,N);%初始化各储存向量，避免下次不同周期计算时引用到前一个周期的结果Velocity=zeros(1,N);AbsAcce=zeros(1,N);t=t+1;End*PLOT*close allfigure %绘制地震记录图plot(time(:),Accelerate(:)title('PEER STRONG MOTION DATABASE RECORD')x

8、label('time(s)')ylabel('acceleration(g)')gridfigure %绘制位移反应谱plot(TA,MaxD(1,:),'-.b',TA,MaxD(2,:),'-r',TA,MaxD(3,:),':k')title('Displacement')xlabel('Tn(s)')ylabel('Displacement(m)')legend(' Z =0.02')Gridfigure %绘制速度反应谱plot(TA,Ma

9、xV(1,:),'-.b',TA,MaxV(2,:),'-r',TA,MaxV(3,:),':k') title('Velocity')xlabel('Tn(s)')ylabel('velocity(m/s)')legend(' Z =0.02')Gridfigure %绘制绝对加速度反应谱plot(TA,MaxA(1,:),'-.b',TA,MaxA(2,:),'-r',TA,MaxA(3,:),':k')title('Abs

10、olute Acceleration')xlabel('Tn(s)')ylabel('absolute acceleration(m/sA2)')legend(' Z =0.02')Grid3运行的结果得到的反应谱LAQStiL4i,YHn£图2位移反应谱OOMmg15 -lgu 盒H-aa图3速度反应谱111-1I I I I I 一 9 s 4vnHJEE-图4加速度反应谱一、反应谱生成地震波1所取的反应谱为上海市设计反应谱设计反血普图5上海市设计反应谱2反应谱取值程序为：%规范反应谱取值程序参照01年抗震规范functio

11、n rs_z=r_s_1(pl,zn,ld,cd,fz) %pl圆频率，zn 阻尼比，ld 烈度，cd 场地类型，场地分组 fz%度选择if ld=6arfmax=0.11;endif ld=7arfmax=0.23;endif ld=8arfmax=0.45;endif ld=9arfmax=0.90;end%地类别，设计地震分组选择if cd=1if fz=1Tg=0.25;endif fz=2Tg=0.30;endif fz=3Tg=0.35;endendif cd=2if fz=1Tg=0.35;*3nH«endiffzETguo.4c? end =hfzs3TgM0.45-

12、end end if oo-hhwiffzMMTgM0.45- end iffzETgHO.55- end =hfzs3Tg"0.65_end end if cdMM4iffzMMTg"0.65_endiffzETg"0.75_end=hfzs3Tguo.9c?endend%ceis-Mzn- %<wm 民_2忒002+(005占毁5)00_if 3X0-mMMO 八end_mR5巴 +(005占包5)二006+.7*ceis-)-if _mR5A0.55-mR5M0.55 八endsjzsMO.9+(oo5ceis-)/(o.5+5*ceis-)八%*再&

13、#187;TT2 T3TglT2MTF7 泮S°TBT3=5*Tg;T_jg=2*pi./pl;%第一段0T1if T_jg<=T1arfg=0.45*arfmax+(lmt2*arfmax-0.45*arfmax)/0.1*Tg;end%第二段 T1 T2if T1<T_jg &T_jgv=T2arf_jg=lmt2*arfmax; end%第三段 T2T3if T2<T_jg &T_jgv=T3arf_jg=(Tg/T_jg)Asjzs)*lmt2*arfmax;end%第四段 T36.0if T3<T_jg &T_jgv=6.0a

14、rf_jg=(lmt2*0.2Asjzs-lmt1*(T_jg-5*Tg)*arfmax; end%第五段6.0if 6.0<T_jgarfg=(lmt2*0.2Asjzs-lmt1*(6.0-5*Tg)*aifmax; end%反应谱值拟加速度值rs_z=arf_jg*9.8;end3生成人造地震波主程序：%主 程序 %确定需要控制的反应谱Sa(T)仃=T1,TM)的坐标点数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.035fo

15、r i=1:nTyzSyz(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1); %8度，2 类场地，第 1 地震分组end%变换的频率差：2*pi*0.005(可以保证长周期项 5s附近有5项三角级数)；%频率变化范围 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; %保 证率%人造地震动持续时间40s，时间间隔：0.02sTd=40;dt=0.02;t=0:0.02:4

16、0;nt=length(t);%衰减包络函数t1=8; %t升段t2=8+24; %平稳段;下降段则为 40- 32= 8sc=0.6; %瘵减段参数for i=1:ntif t(i)v=t1f(i)=(t(i)/t1F2;endif t(i)>t1 & t(i)<t2f(i)=1;endif t(i)>=t2f(i)=exp(-c*(t(i)-t2);endend%反应谱转换功率谱for i=1:nplSw(i)=(2*ceita/(pi*pl(i)*r_s_1(pl(i),ceita,8,2,1F2心2*log(-1*pi*log(P)/(pl(i)*Td);Aw

17、(i)=sqrt(4*Sw(i)*plc);end%合成地震动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%计算反应谱验证是否满足rc在5%的要求,需要时程动力分析% response spectra of callidar% parameterg=9.8;m=1;xO=O;vO=O;ww=2*pi./Tyz;% loadag=at; % 修改% solu

18、tionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=wA2;for i=1:nt-1p(i)=-ag(i+1)+ag(i);aO=m(-ag(i)-c*vO-k*xO);kk=k+(dtA2)(6*m)+dt(3*c);pp=p(i)+m*(dt(6*vO)+3*aO)+c*(3*vO+2(dt*aO);dx=kkpp;dv=dt(3*dx)-3*vO-2(dt*aO);x仁 xO+dx;xO=x1;v1= vO+dv;vO=v1;as(i)=aO;as(i)=as(i)+ag(i);vs(i)=vO;xs(i)=xO;endmaxas(y)=max(as);m

19、axvs(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%!环体函数blxs=Syz./rspa;*3nH«for XSXSA5P-if 2*pi/p-(xsxsATyzu) b_xs(xsxs)Mb_xsu)八 endfor SXSXM+nTyz 匕 if(2*pi/p-(xsxs)VMTyz(sxsx) 3 (2*pi/p

20、-(xsxsAMTyz(sxsxi)b_xs(xsxs)Mb_xs(sxsx)+(b_xs(sxsx+)b_xs(sxsx)*(2*pi/p-(XSXSTTyz(sxsx)/(Tyz(sxsx+)Tyz(sxsx)八endendif 2*pi/p-(xsxs)VTyz(nTyz) b_xs(xsxs)Mb_xs(nTyz= endendAWMAW.*b_xsn%镇>鑒窗alrzerosmL)八aili-zeros(nL:)八for j"5iaisAs*AW(i)*rea-(exp(sqrts)*(p-(i)*s+fai(i)八endallai+aFend%*河m请專电徊审裁百

21、rcm5蚩4斓決 % response speara of ca 三 dar % parametergupo?XOMO 八VOMO 八WWM2*PL/TyN%oadaguac% 氢% socHonfor y2nTyzZUOO37-WMWW(y)八CM2*z*wikuw>2for im 匕二沐syB*3nH« p(i)pagT)+age aoMm 一ig(vc*v9k*xo)八 kkMk+(ds2)一(6*m)+cn一(3*2 ppup(i)+m*(cm(6*vo)+3*ao)+c*(3*vo+2v2:*ao)- dxMkk 一 pp- d VMCm(3*dxT3*v92 一(di*ao)- xMXO+d><XOMX八<MVO+d<VOM<八 as(Tap as(iHas(i)+ag(ixVS(TVPXS(TXPend maxas(y)Mmax(as= maxvs(y)Mmax(vs)八 maxxs(y)Mmax(xs)八 endfor 耳5Tyz rspa(iHmaxas(一)八 end% 黨<mfor 耳5Tyzrcrsp(Tabs(rspa?syz(一)/max(syzc)=endjsnumMjsnumimax(rcrspc)end%»