地震工程学-反应谱和地震时程波的相互转化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

%%%% 第三段 T2~T3 if T2

arf_jg=((Tg/T_jg)^sjzs)*lmt2*arfmax; end

%%%% 第四段 T3～6.0 if T3

arf_jg=(lmt2*0.2^sjzs-lmt1*(T_jg-5*Tg))*arfmax; end

%%%% 第五段 6.0～ if 6.0

arf_jg=(lmt2*0.2^sjzs-lmt1*(6.0-5*Tg))*arfmax; end

%%%%%%反应谱值 拟加速度值 rs_z=arf_jg*9.8; end

3 生成人造地震波主程序：

%%%主程序%%%%

%%%%确定需要控制的反应谱Sa（T）(T=T1,...,TM)的坐标点数M， 反应谱控制容差rc Tyz=[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.035 for i=1: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.02s Td=40; dt=0.02; t=0:0.02:40; nt=length(t); %%%%%%% 衰减包络函数 t1=8; %%%%上升段

t2=8+24; %%%%%平稳段; 下降段则为40－32＝8s c=0.6; %%%%衰减段参数 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)

for sxsx=1:nTyz-1

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 人造地震波和初始反应谱