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

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

2、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;

3、for T=0.0:0.05:6NatualFrequency=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/sq

4、rt(1-Damp*Damp); A(2,2)=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_3t;B(1,2)=-e_t*(d_f*s/DampFrequency+2*d_3t*c)-1/NatualFrequency2+2*d

5、_3t;B(2,1)=e_t*(d_f+Damp/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/(NatualFrequency2*Dt); for i=1:(N-1) %根据地震记录,计算不同的反应Displa

6、ce(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,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);end MaxD(1,t)=max(abs(Displace);

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

8、RD)xlabel(time(s)ylabel(acceleration(g) grid figure %绘制位移反应谱plot(TA,MaxD(1,:),-.b,TA,MaxD(2,:),-r,TA,MaxD(3,:),:k) title(Displacement) xlabel(Tn(s)ylabel(Displacement(m)legend(=0.02) Grid figure %绘制速度反应谱plot(TA,MaxV(1,:),-.b,TA,MaxV(2,:),-r,TA,MaxV(3,:),:k)title(Velocity)xlabel(Tn(s)ylabel(velocity(

9、m/s)legend(=0.02)Gridfigure %绘制绝对加速度反应谱plot(TA,MaxA(1,:),-.b,TA,MaxA(2,:),-r,TA,MaxA(3,:),:k)title(Absolute Acceleration)xlabel(Tn(s)ylabel(absolute acceleration(m/s2)legend(=0.02)Grid3 运行的结果得到的反应谱 图2 位移反应谱图3 速度反应谱图4 加速度反应谱一、 反应谱生成地震波1 所取的反应谱为上海市设计反应谱图5 上海市设计反应谱2 反应谱取值程序为：%规范反应谱取值程序 参照01年抗震规范functio

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

11、=3 Tg=0.45; endendif cd=3 if fz=1 Tg=0.45; end if fz=2 Tg=0.55; end if fz=3 Tg=0.65; endendif cd=4 if fz=1 Tg=0.65; end if fz=2 Tg=0.75; end if fz=3 Tg=0.90; endend%ceita=zn; %阻尼比lmt1=0.02+(0.05-ceita)/8;if lmt10 lmt1=0;endlmt2=1+(0.05-ceita)/(0.06+1.7*ceita);if lmt20.55 lmt2=0.55;endsjzs=0.9+(0.05-

12、ceita)/(0.5+5*ceita);%分段位置 T1 T2 T3T1=0.1;T2=Tg;T3=5*Tg;T_jg=2*pi./pl;% 第一段 0T1if T_jg=T1 arf_jg=0.45*arfmax+(lmt2*arfmax-0.45*arfmax)/0.1*T_jg;end% 第二段 T1T2if T1T_jg&T_jg=T2 arf_jg=lmt2*arfmax;end% 第三段 T2T3if T2T_jg&T_jg=T3 arf_jg=(Tg/T_jg)sjzs)*lmt2*arfmax;end% 第四段 T36.0if T3T_jg&T_jg=6.0 arf_jg=(

13、lmt2*0.2sjzs-lmt1*(T_jg-5*Tg)*arfmax;end% 第五段 6.0if 6.0T_jg arf_jg=(lmt2*0.2sjzs-lmt1*(6.0-5*Tg)*arfmax;end%反应谱值 拟加速度值rs_z=arf_jg*9.8;end3 生成人造地震波主程序：%主程序%确定需要控制的反应谱Sa（T）(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.035for i=1

14、:nTyz Syz(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*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; %保证率%人造地震动持续时间40s，时间间隔：0.02sTd=40; dt=0.02;t=0:0.02:40;nt=l

15、ength(t);% 衰减包络函数t1=8; %上升段t2=8+24; %平稳段; 下降段则为40328sc=0.6; %衰减段参数for i=1:nt if t(i)t1 & t(i)=t2 f(i)=exp(-c*(t(i)-t2); end end% 反应谱转换功率谱for i=1:npl Sw(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);fo

16、r i=1:npl fai(i)=rand(1)*2*pi; for j=1:nt atj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i); end at=at+atj;end% 计算反应谱验证是否满足rc在5%的要求,需要时程动力分析% response spectra of callidar% parameterg=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;% loadag=at; %修改% solutionfor y=1:nTyz z=0.037; w=ww(y); c=2*z*w; k=w2; for i=1:n

17、t-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; end maxas(y)=max(as); maxvs(y)=max(vs); maxxs(y)=max(xs);endfor

18、 i=1:nTyz rspa(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; for xsxs=1:npl if 2*pi/pl(xsxs)=Tyz(sxsx) & (2*pi/pl(xsxs)Tyz(nTyz) blxs1(xsxs)=blxs(nTyz); end end Aw=Aw.*blxs1; % 合成地震动 at=zeros(nt,1); atj=zeros(nt,1); fo

19、r i=1:npl for j=1:nt atj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i); end at=at+atj; end % 计算反应谱验证是否满足rc在5%的要求 % response spectra of callidar % parameter g=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;% loadag=at; %修改% solutionfor y=1:nTyz z=0.037; w=ww(y); c=2*z*w; k=w2; for i=1:nt-1 p(i)=-ag(i+1)+ag(i);

20、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; end maxas(y)=max(as); maxvs(y)=max(vs); maxxs(y)=max(xs);endfor i=1:nTyz rspa(i)=maxas(i

21、);end% 比较容差for i=1:nTyzrcrsp(i)=abs(rspa(i)-Syz(i)/max(Syz(:); end jsnum=jsnum+1 max(rcrsp(:)end% 最终的反应谱 与规范谱 % response spectra of callidar % parameter% Tjs=0.05:0.01:6; % nTjs=length(Tjs); g=9.8; m=1; x0=0; v0=0; ww=2*pi./Tyz; % load ag=at; %修改 % solution for y=1:nTyz z=0.037; w=ww(y); c=2*z*w; k=

22、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; end maxas(y)=max(as); maxvs(y)=max(vs); maxxs(y)=max(xs); end for i=1:nTyz rspa(i)=maxas(i)/g; rspa_S(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1)/g; end subplot(2,1,1);plot(t,at);subplot(2,1,2);plot(Tyz,rspa);hold on;plot(Tyz,rspa_S);4 生成的人造地震波如图所示。图6 人造地震波和初始反应谱