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

地震工程学作业课程名称: 地震工程学指导老师: ______ 翟永梅__________ 姓名: ________ 史先飞________ 学号: ______ 1232627 ________一、地震波生成反应谱1所取的地震波为Elce ntro地震波加速度曲线，如图1所示。

400系列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); % 相对位移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: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 A2-1)/(NatualFrequency A2*Dt);d_3t=Damp/(NatualFrequency A3*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*((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)=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)=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')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(' Z =0.02')Gridfigure %绘制速度反应谱plot(TA,MaxV(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('Absolute Acceleration')xlabel('Tn(s)')ylabel('absolute acceleration(m/s A2)')legend(' Z =0.02')Grid3运行的结果得到的反应谱OO MQS^tiL4i YHn£m g15 -l^gu■盒—H 图2位移反应谱RnlaJE一"UJ一r u'q :I.图3速度反应谱--丄4图4加速度反应谱、反应谱生成地震波1所取的反应谱为上海市设计反应谱设计反血［普图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%%%%地类别，设计地震分组选择if cd==1if fz==1Tg=0.25;endif fz==2Tg=0.30;endif fz==3Tg=0.35;endendif cd==2if fz==1Tg=0.35;end iffzE T g u o.4c?end=h f z s3TgM0.45- endendIFOO-HHWiffzMM 」TgM0.45-endiffzE TgHO.55-end=h f z s3Tg"0.65_ endend ifcdMM4iffzMM 」Tg"0.65_ endiffzE Tg"0.75_ end =h f z s 3T g u o .9c?end end %%%%%%%%% c e i s-M z n -%%%%<wm 民_2忒002+(005占毁5))00_if 3X 0-mMMO 八end _m R5巴 +(005占包5))二006+」.7*c e i s-)_if_m R5A 0.55-m R5M 0.55 八end sjzsMO.9+(oo5—ceis-)/(o.5+5*ceis-)八%%%<^^wT 」T2T3TglT3=5*Tg;T_jg=2*pi./pl;%%%第一段0〜T1if T_jg<=T1 arf」g=0.45*arfmax+(lmt2*arfmax-0.45*arfmax)/0.1*T」g;end%%%第二段T1〜T2if T1<T_jg &T_jgv=T2arf_jg=lmt2*arfmax;end%%%第三段T2~T3if T2<T_jg &T_jgv=T3arf_jg=((Tg/T_jg)A sjzs)*lmt2*arfmax;end%%%第四段T3〜6.0if T3<T_jg &T_jgv=6.0arf_jg=(lmt2*0.2Asjzs-lmt1*(T_jg-5*Tg))*arfmax;end%%%第五段6.0〜if 6.0<T_jgarf」g=(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.035for 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:40;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(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; %%%%%% 修改%%%%%%% solutionfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w A2;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=kk\pp;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);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=1;while max(rcrsp(:))>rc%%%%%!环体函数blxs=Syz./rspa;for X S X S A 5P -if 2*pi/p-(xsxsATyzu)b_xs 」(xsxs)Mb_xsu)八 end forSXSXM+nTyz 匕if(2*pi/p-(xsxs)VMTyz(sxsx))3(2*pi/p-(xsxsAMTyz(sxsxi))b_xs 」(xsxs)Mb_xs(sxsx)+(b_xs(sxsx+」)—b_xs(sxsx))*(2*pi/p-(XSXSTTyz(sxsx))/(Tyz(sxsx+」)—Tyz(sxsx))八end end if 2*pi/p-(xsxs)VTyz(nTyz)b_xs 」(xsxs)Mb_xs(nTyz=endend AWMAW.*b_xsn %%%%%%%%%%镇>鑒窗alrzerosmL 」)八aili-zeros(nL:)八for j "5ia i s A s *A W (i )*r e a -(e x p (s q r t s)*(p -(i )*s +f a i (i ))))八end a ll a i +a Fend %%%%%%**河m 请專电徊审裁百 r c m 5蚩4斓決%%%%%%%%%%%% response speara ofca三dar%%%%%%%parameter g u p o?XOMO 八VOMO 八WWM2*PL/TyN %%%%%%%%o a daguac%%%%%% 氢%%%%%%%socHonfory2nTyz Z U O O 37-WMWW(y)八CM2*z*wikuw>2p (i )p a g T 」)+a g eaoMm一i g(v c *v 9k *x o )八 k k M k +(d s2)一(6*m)+cn一(3*2ppup(i)+m*(cm(6*vo)+3*ao)+c*(3*vo+2v2:*ao))_ dxMkk 一 pp- d VMCm(3*dxT3*v92 一(di*ao)-x 」M X O +d ><XOMX 」八<M V O +d <V O M <八a s (T a pas(iHas(i)+ag(ix V S (T V PX S (T X Pendmaxas(y)Mmax(as=maxvs(y)Mmax(vs)八maxxs(y)Mmax(xs)八end for 耳5Tyzrspa(iHmaxas(一)八end%%%%% 黨<mfor耳5Tyz rcrsp(Tabs(rspa?syz(一))/max(syzc)=end jsnumMjsnumi max(rcrspc)) end%%%%%»凑罢河m请>涪酱请%%%%%%%%%%%%responsespearaofca 三dar%%%%%%%parameter %%Tjs"0.05o0rp? %%nTjsllengihcrjs)八g u p o?XOMO 八VOMO 八WWM2*PL/TyN %%%%%%%%o a dagMac%%%%%%%w 涔%%%%%%%socHonfor y=1:nTyzz=0.037;w=ww(y);c=2*z*w;k=w A2;for 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=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);endfor i=1:nTyzrspa(i)=maxas(i)/g;rspa_S(i)=r_s_1(2*pi/Tyz(i),ceita,8,2,1)/g;endsubplot(2,1,1);plot(t,at);subplot(2,1,2);plot(Tyz,rspa);hold on;plot(Tyz,rspa_S);4生成的人造地震波如图所示图6人造地震波和初始反应谱。