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

地震工程学-反应谱和地震时程波的相互转化matlab编程地震工程学-反应谱和地震时程波的相互转化matlab编程地震工程学-反应谱和地震时程波的相互转化matlab编程xxx公司地震工程学-反应谱和地震时程波的相互转化matlab编程文件编号：文件日期：修订次数：第1.0次更改批准审核制定方案设计，管理制度地震工程学作业课程名称:地震工程学______指导老师:_______翟永梅_________姓名:史先飞________学号:1232627________一、地震波生成反应谱1所取的地震波为Elcentro地震波加速度曲线，如图1所示。图1Elcentro地震波加速度曲线2所调用的Matlab程序为：%***********读入地震记录***********ElCentro;Accelerate=ElCentro(:,1)*;%单位统一为m和sN=length(Accelerate);%N读入的记录的量time=0::(N-1)*;%单位s%初始化各储存向量Displace=zeros(1,N);%相对位移Velocity=zeros(1,N);%相对速度AbsAcce=zeros(1,N);%绝对加速度%***********A,B矩阵***********Damp=;%阻尼比TA=::6;%TA=::6;%结构周期Dt=;%地震记录的步长%记录计算得到的反应，MaxD为某阻尼时最大相对位移，MaxV为某阻尼最大相对速度，MaxA某阻尼时最大绝对加速度，用于画图MaxD=zeros(3,length(TA));MaxV=zeros(3,length(TA));MaxA=zeros(3,length(TA));t=1;forT=::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/sqrt(1-Damp*Damp);A(2,2)=e_t*(-s*Damp/sqrt(1-Damp*Damp)+c);d_f=(2*Damp^2-1)/(NatualFrequency^2*Dt);d_3t=Damp/(NatualFrequency^3*Dt);B=zeros(2,2);B(1,1)=e_t*((d_f+Damp/NatualFrequency)*s/DampFrequency+(2*d_3t+1/NatualFrequency^2)*c)-2*d_3t;B(1,2)=-e_t*(d_f*s/DampFrequency+2*d_3t*c)-1/NatualFrequency^2+2*d_3t;B(2,1)=e_t*((d_f+Damp/NatualFrequency)*(c-Damp/sqrt(1-Damp^2)*s)-(2*d_3t+1/NatualFrequency^2)*(DampFrequency*s+Damp*NatualFrequency*c))+1/(NatualFrequency^2*Dt);B(2,2)=e_t*(1/(NatualFrequency^2*Dt)*c+s*Damp/(NatualFrequency*DampFrequency*Dt))-1/(NatualFrequency^2*Dt);fori=1:(N-1)%根据地震记录,计算不同的反应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)-NatualFrequency^2*Displace(i+1);endMaxD(1,t)=max(abs(Displace));MaxV(1,t)=max(abs(Velocity));ifT==MaxA(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***********closeallfigure%绘制地震记录图plot(time(:),Accelerate(:))title('PEERSTRONGMOTIONDATABASERECORD')xlabel('time(s)')ylabel('acceleration(g)')gridfigure%绘制位移反应谱plot(TA,MaxD(1,:),'',TA,MaxD(2,:),'-r',TA,MaxD(3,:),':k')title('Displacement')xlabel('Tn(s)')ylabel('Displacement(m)')legend('ζ=')Gridfigure%绘制速度反应谱plot(TA,MaxV(1,:),'',TA,MaxV(2,:),'-r',TA,MaxV(3,:),':k')title('Velocity')xlabel('Tn(s)')ylabel('velocity(m/s)')legend('ζ=')Gridfigure%绘制绝对加速度反应谱plot(TA,MaxA(1,:),'',TA,MaxA(2,:),'-r',TA,MaxA(3,:),':k')title('AbsoluteAcceleration')xlabel('Tn(s)')ylabel('absoluteacceleration(m/s^2)')legend('ζ=')Grid3运行的结果得到的反应谱图2位移反应谱图3速度反应谱图4加速度反应谱反应谱生成地震波所取的反应谱为上海市设计反应谱图5上海市设计反应谱反应谱取值程序为：%%规范反应谱取值程序参照01年抗震规范functionrs_z=r_s_1(pl,zn,ld,cd,fz)%%%pl圆频率,zn阻尼比,ld烈度,cd场地类型,场地分组fz%%%%烈度选择ifld==6arfmax=;endifld==7arfmax=;endifld==8arfmax=;endifld==9arfmax=;end%%%%场地类别，设计地震分组选择ifcd==1iffz==1Tg=;endiffz==2Tg=;endiffz==3Tg=;endendifcd==2iffz==1Tg=;endiffz==2Tg=;endiffz==3Tg=;endendifcd==3iffz==1Tg=;endiffz==2Tg=;endiffz==3Tg=;endendifcd==4iffz==1Tg=;endiffz==2Tg=;endiffz==3Tg=;endend%%%%%%%%%ceita=zn;%%%%%阻尼比lmt1=+/8;iflmt1<0lmt1=0;endlmt2=1+/+*ceita);iflmt2<lmt2=;endsjzs=+/+5*ceita);%%%%%分段位置T1T2T3T1=;T2=Tg;T3=5*Tg;T_jg=2*pi./pl;%%%%第一段0～T1ifT_jg<=T1arf_jg=*arfmax+(lmt2**arfmax)/*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～ifT3<T_jg&T_jg<=arf_jg=(lmt2*^sjzs-lmt1*(T_jg-5*Tg))*arfmax;end%%%%第五段～if<T_jgarf_jg=(lmt2*^sjzs-lmt1**Tg))*arfmax;end%%%%%%反应谱值拟加速度值rs_z=arf_jg*;end生成人造地震波主程序：%%%主程序%%%%%%%%确定需要控制的反应谱Sa（T）(T=T1,...,TM)的坐标点数M，反应谱控制容差rcTyz=[::,::,::];rc=;nTyz=length(Tyz);ceita=;%%%阻尼比：fori=1:nTyzSyz(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1);%%%%8度，2类场地，第1地震分组end%%%%%%变换的频率差：2*pi*(可以保证长周期项5s附近有5项三角级数)；%%%%频率变化范围N1=30,30**2*pi;N2=3000,5000**2*piplc=2*pi*;pl=30**2*pi:*2*pi:10000**2*pi;npl=length(pl);P=;%%%保证率%%%%%%人造地震动持续时间40s，时间间隔：Td=40;dt=;t=0::40;nt=length(t);%%%%%%%衰减包络函数t1=8;%%%%上升段t2=8+24;%%%%%平稳段;下降段则为40－32＝8sc=;%%%%衰减段参数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=;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%%loadag=at;%%%%%%%修改%%%%%%%solutionfory=1:nTyzz=;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%%%%%循环体函数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=;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%%loadag=at;%%%%%%%修改%%%%%%%solutionfory=1:nTyzz=;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%%%%%%%最终的反应谱与规范谱%%%%%%%%%%%%responsespectraofcallidar%%%%%%%parameter%%Tjs=::6;%%nTjs=length(Tjs);g=;m=1;x0=0;v0=0;ww=2*pi./Tyz;%%%%%%%%