为了正常的体验网站,请在浏览器设置里面开启Javascript功能!

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

2019-02-22 28页 doc 55KB 127阅读

用户头像

is_266065

暂无简介

举报
地震工程学-反应谱和地震时程波的相互转化matlab编程 地震工程学作业 课程名称:      地震工程学______ 指导老师:_______翟永梅_________ 姓    名:        史先飞________       学   号:      1232627________ 一、地震波生成反应谱 1 所取的地震波为Elcentro地震波加速度曲线,如图1所示。 图1 Elcentro地震波加速度曲线 2 所调用的Matlab程序为: % ***********读入地震记录*********** ElCentro; Accelerate= ElCentro(:,1...
地震工程学-反应谱和地震时程波的相互转化matlab编程
地震工程学作业 课程名称:      地震工程学______ 指导老师:_______翟永梅_________ 姓    名:        史先飞________       学   号:      1232627________ 一、地震波生成反应谱 1 所取的地震波为Elcentro地震波加速度曲线,如图1所示。 图1 Elcentro地震波加速度曲线 2 所调用的Matlab程序为: % ***********读入地震记录*********** ElCentro; Accelerate= ElCentro(:,1)*9.8067;%单位统一为m和s N=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.02 TA=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:6 NatualFrequency=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); 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)-NatualFrequency^2*Displace(i+1); end MaxD(1,t)=max(abs(Displace)); MaxV(1,t)=max(abs(Velocity)); if T==0.0 MaxA(1,t)=max(abs(Accelerate)); else MaxA(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 RECORD') 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(m/s)') legend('ζ=0.02') Grid figure %绘制绝对加速度反应谱 plot(TA,MaxA(1,:),'-.b',TA,MaxA(2,:),'-r',TA,MaxA(3,:),':k') title('Absolute Acceleration') xlabel('Tn(s)') ylabel('absolute acceleration(m/s^2)') legend('ζ=0.02') Grid 3 运行的结果得到的反应谱 图2  位移反应谱 图3 速度反应谱 图4  加速度反应谱 一、 反应谱生成地震波 1 所取的反应谱为上海市反应谱 图5 上海市设计反应谱 2 反应谱取值程序为: %%反应谱取值程序 参照01年抗震规范 function rs_z=r_s_1(pl,zn,ld,cd,fz) %%%pl 圆频率,zn阻尼比,ld烈度,cd场地类型,场地分组fz %%%%烈度选择 if ld==6 arfmax=0.11; end if ld==7 arfmax=0.23; end if ld==8 arfmax=0.45; end if 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; end end if cd==2 if fz==1 Tg=0.35; end if fz==2        Tg=0.40; end if fz==3 Tg=0.45; end end if cd==3 if fz==1 Tg=0.45; end if fz==2        Tg=0.55; end if fz==3 Tg=0.65; end end if cd==4 if fz==1 Tg=0.65; end if fz==2        Tg=0.75; end if fz==3 Tg=0.90; end end %%%%%%%%% ceita=zn; %%%%%阻尼比 lmt1=0.02+(0.05-ceita)/8; if lmt1<0 lmt1=0; end lmt2=1+(0.05-ceita)/(0.06+1.7*ceita); if lmt2<0.55 lmt2=0.55; end sjzs=0.9+(0.05-ceita)/(0.5+5*ceita); %%%%%分段位置 T1 T2 T3 T1=0.1; T2=Tg; T3=5*Tg; T_jg=2*pi./pl; %%%% 第一段 0~T1 if T_jg<=T1 arf_jg=0.45*arfmax+(lmt2*arfmax-0.45*arfmax)/0.1*T_jg; end %%%% 第二段 T1~T2 if T1参数
for i=1:nt if t(i)<=t1 f(i)=(t(i)/t1)^2; end 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); for 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 %%%%%%% parameter 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=w^2; for i=1:nt-1 p(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; end maxas(y)=max(as); maxvs(y)=max(vs); maxxs(y)=max(xs); end for i=1:nTyz rspa(i)=maxas(i); end %%%%%%% 比较容差 for i=1:nTyz rcrsp(i)=abs(rspa(i)-Syz(i))/max(Syz(:));  end jsnum=1; while max(rcrsp(:))>rc %%%%%循环体函数 blxs=Syz./rspa; for xsxs=1:npl if 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)); end end if 2*pi/pl(xsxs)>Tyz(nTyz) blxs1(xsxs)=blxs(nTyz); end end Aw=Aw.*blxs1; %%%%%%%%%%%%%% 合成地震动 at=zeros(nt,1); atj=zeros(nt,1); for 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; %%%%%%%% load ag=at;  %%%%%%%修改 %%%%%%% solution for y=1:nTyz z=0.037; w=ww(y); c=2*z*w; k=w^2; for i=1:nt-1 p(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; end maxas(y)=max(as); maxvs(y)=max(vs); maxxs(y)=max(xs); end for i=1:nTyz rspa(i)=maxas(i); end %%%%%%% 比较容差 for i=1:nTyz rcrsp(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=w^2; for i=1:nt-1 p(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; 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 人造地震波和初始反应谱
/
本文档为【地震工程学-反应谱和地震时程波的相互转化matlab编程】,请使用软件OFFICE或WPS软件打开。作品中的文字与图均可以修改和编辑, 图片更改请在作品中右键图片并更换,文字修改请直接点击文字进行修改,也可以新增和删除文档中的内容。
[版权声明] 本站所有资料为用户分享产生,若发现您的权利被侵害,请联系客服邮件isharekefu@iask.cn,我们尽快处理。 本作品所展示的图片、画像、字体、音乐的版权可能需版权方额外授权,请谨慎使用。 网站提供的党政主题相关内容(国旗、国徽、党徽..)目的在于配合国家政策宣传,仅限个人学习分享使用,禁止用于任何广告和商用目的。

历史搜索

    清空历史搜索