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

1、同除为fTONGJI UNIVERSITY地震工程学作业课程名称:地震工程学指导老师:翟永梅姓 名:史先飞学 号:12326271 / 14word.一、地震波生成反应谱1所取的地震波为日centro地震波加速度曲线，如图 1所示。图1日centro地震波加速度曲线2所调用的Matlab程序为：*读入地震记录*ElCentro;Accelerate= ElCentro(:,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; %地震记录的步长MaxV为某阻尼最大相对速度，MaxA某阻尼时最大%己录计算得到的反应，MaxD为某阻尼时最大相对位移,绝对加速度，用于画图MaxD=zeros(3,length(TA);MaxV=zeros(3,length(TA);MaxA=zeros(3,length(TA);t=1;for T=0.0:0.0

3、5:6NatualFrequency=2*pi/T ; % 结构自振频率DampFrequency=NatualFrequency*sqrt(1-Damp*Damp); %计算公式化简2 / 14word. 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

4、(1-Damp*Damp);A(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*

5、d_3t;B(2,1)=e_t*(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) % 根据地震记录，计算不同的

6、反应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,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(Di

7、splace);MaxV(1,t)=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;End3 / 14word.*DI CT*close allfigure %绘制地震记录图plot(time(:),Accelerate(:)title(PEER STRONG MOTION

8、 DATABASE RECORD)xlabel(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(1=0.02)Gridfigure %绘制速度反应谱plot(TA,MaxV(1,:),-.b,TA,MaxV(2,:),-r,TA,MaxV(3,:),:k) title(Velocity)xlabel(Tn(s)ylab

9、el(velocity(m/s)legend(1=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/sA2)legend(1=0.02)Grid3运行的结果得到的反应谱4 / 14word.图2位移反应谱图3速度反应谱Ahi sbjlF- Nrr干FEtti图4加速度反应谱 45 / 14word.、反应谱生成地震波1所取的反应谱为上海市设计反

10、应谱图5上海市设计反应谱2反应谱取值程序为：%觑范反应谱取值程序参照01年抗震规范function 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%k类别，设计地震分组选择if cd=1if fz=1Tg=0.25;endif fz=2Tg=0.30;endif fz=3Tg=0.35;endendif cd=2if fz=1Tg=0.

11、35;6 / 14word.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;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 lmt10lmt1=0;endlmt2=1+(0.05-ceita)/(0.06+1.7*ceita);if lmt20.55lmt2=0.55;

12、endsjzs=0.9+(0.05-ceita)/(0.5+5*ceita);%位置 T1 T2 T3T1=0.1;T2=Tg;7 / 14word.T3=5*Tg;T_jg=2*pi./pl;%嗡一段0T1if T_jg=T1arf_jg=0.45*arfmax+(lmt2*arfmax-0.45*arfmax)/0.1*T_jg; end%嗡二段T1T2if T1T_jg&T_jg=T2arf_jg=lmt2*arfmax;end%嗡三段T2-T3if T2T_jg&T_jg=T3arf_jg=(Tg/T_jg)Asjzs)*lmt2*arfmax;end%嗡四段T36.0if T3T_j

13、g&T_jg=6.0arf_jg=(lmt2*0.2Asjzs-lmt1*(T_jg-5*Tg)*arfmax;end%嗡五段6.0if 6.0T_jgarf_jg=(lmt2*0.2Asjzs-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;% 阻尼比：

14、0.035 for 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*pi8 / 14word.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

15、.02;t=0:0.02:40;nt=length(t);%包络函数t1=8; %t 升段t2=8+24; %稳段；下降段则为 40 32 = 8sc=0.6; %短减段参数for i=1:ntif t(i)t1 & 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,1)A2/(-2*log(-1*pi*log(P)/(pl(i)*Td);Aw(i)=sqrt(4*Sw(i)*plc);end%动at=zeros(nt,1);atj=zeros(

16、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;9 / 14word.x0=0;v0=0;ww=2*pi./Tyz;% loadag=at; %脩改% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=wA2;fo

17、r i=1:nt-1p(i)=-ag(i+1)+ag(i);a0=m(-ag(i)-c*v0-k*x0);kk=k+(dtA2)(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:nTyzr

18、spa(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;10 / 14word.for xsxs=1:nplif 2*pi/pl(xsxs)=Tyz(sxsx) & (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

19、:ntatj(j)=f(j)*Aw(i)*real(exp(sqrt(-1)*(pl(i)*t(j)+fai(i);endat=at+atj;end%h算反应谱验证是否满足rc在5%勺要求% response spectra of callidar% parameterg=9.8;m=1;x0=0;v0=0;ww=2*pi./Tyz;% load ag=at; %改% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=wA2;for i=1:nt-111 / 14word.p(i)=-ag(i+1)+ag(i);a0=m(-ag(i)-c*v0-k*x

20、0);kk=k+(dtA2)(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%的反应谱与规范谱% response spectra of callidar% parameter% Tj