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