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

1 所取的地震波为Elcentro地震波加速度曲线，如图1所示。

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_3 t;

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/(NatualF requency^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 运行的结果得到的反应谱

1所取的反应谱为上海市设计反应谱

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;

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;

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

arf_jg=lmt2*arfmax;

end

%%%% 第三段 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)

f(i)=1;

end

if 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;

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)

blxs1(xsxs)=blxs(1);

end

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;

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;

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生成的人造地震波如图所示。