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

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

《地震工程学反应谱和地震时程波的相互转化matlab编程.docx》由会员分享，可在线阅读，更多相关《地震工程学反应谱和地震时程波的相互转化matlab编程.docx（20页珍藏版）》请在冰豆网上搜索。

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

地震工程学作业

课程名称:

地震工程学______

指导老师:

_______翟永梅_________

姓名:

史先飞________

学　号:

1232627________

一、地震波生成反应谱

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

图1Elcentro地震波加速度曲线

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;

forT=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）;

fori=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））;

ifT==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***********

closeall

figure%绘制地震记录图

plot（time（:

）,Accelerate（:

））

title（'PEERSTRONGMOTIONDATABASERECORD'）

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（'AbsoluteAcceleration'）

xlabel（'Tn（s）'）

ylabel（'absoluteacceleration（m/s^2）'）

legend（'ζ=0.02'）

Grid

3运行的结果得到的反应谱

图2位移反应谱

图3速度反应谱

图4加速度反应谱

一、反应谱生成地震波

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

图5上海市设计反应谱

2反应谱取值程序为：

%%规范反应谱取值程序参照01年抗震规范

functionrs_z=r_s_1（pl,zn,ld,cd,fz）%%%pl圆频率,zn阻尼比,ld烈度,cd场地类型,场地分组fz

%%%%烈度选择

ifld==6

arfmax=0.11;

end

ifld==7

arfmax=0.23;

end

ifld==8

arfmax=0.45;

end

ifld==9

arfmax=0.90;

end

%%%%场地类别，设计地震分组选择

ifcd==1

iffz==1

Tg=0.25;

end

iffz==2

Tg=0.30;

end

iffz==3

Tg=0.35;

end

end

ifcd==2

iffz==1

Tg=0.35;

end

iffz==2

Tg=0.40;

end

iffz==3

Tg=0.45;

end

end

ifcd==3

iffz==1

Tg=0.45;

end

iffz==2

Tg=0.55;

end

iffz==3

Tg=0.65;

end

end

ifcd==4

iffz==1

Tg=0.65;

end

iffz==2

Tg=0.75;

end

iffz==3

Tg=0.90;

end

end

%%%%%%%%%

ceita=zn;%%%%%阻尼比

lmt1=0.02+（0.05-ceita）/8;

iflmt1<0

lmt1=0;

end

lmt2=1+（0.05-ceita）/（0.06+1.7*ceita）;

iflmt2<0.55

lmt2=0.55;

end

sjzs=0.9+（0.05-ceita）/（0.5+5*ceita）;

%%%%%分段位置T1T2T3

T1=0.1;

T2=Tg;

T3=5*Tg;

T_jg=2*pi./pl;

%%%%第一段0～T1

ifT_jg<=T1

arf_jg=0.45*arfmax+（lmt2*arfmax-0.45*arfmax）/0.1*T_jg;

end

%%%%第二段T1～T2

ifT1

arf_jg=lmt2*arfmax;

end

%%%%第三段T2~T3

ifT2

arf_jg=（（Tg/T_jg）^sjzs）*lmt2*arfmax;

end

%%%%第四段T3～6.0

ifT3

arf_jg=（lmt2*0.2^sjzs-lmt1*（T_jg-5*Tg））*arfmax;

end

%%%%第五段6.0～

if6.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

fori=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;%%%%衰减段参数

fori=1:

nt

ift（i）<=t1

f（i）=（t（i）/t1）^2;

end

ift（i）>t1&t（i）

f（i）=1;

end

ift（i）>=t2

f（i）=exp（-c*（t（i）-t2））;

end

end

%%%%%%%反应谱转换功率谱

fori=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）;

fori=1:

npl

fai（i）=rand

（1）*2*pi;

forj=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%的要求,需要时程动力分析

%%%%%%%%%%%%responsespectraofcallidar

%%%%%%%parameter

g=9.8;

m=1;

x0=0;

v0=0;

ww=2*pi./Tyz;

%%%%%%%%load

ag=at;%%%%%%%修改

%%%%%%%solution

fory=1:

nTyz

z=0.037;

w=ww（y）;

c=2*z*w;

k=w^2;

fori=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

fori=1:

nTyz

rspa（i）=maxas（i）;

end

%%%%%%%比较容差

fori=1:

nTyz

rcrsp（i）=abs（rspa（i）-Syz（i））/max（Syz（:

））;

end

jsnum=1;

whilemax（rcrsp（:

））>rc

%%%%%循环体函数

blxs=Syz./rspa;

forxsxs=1:

npl

if2*pi/pl（xsxs）

（1）

blxs1（xsxs）=blxs

（1）;

end

forsxsx=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

if2*pi/pl（xsxs）>Tyz（nTyz）

blxs1（xsxs）=blxs（nTyz）;

end

end

Aw=Aw.*blxs1;

%%%%%%%%%%%%%%合成地震动

at=zeros（nt,1）;

atj=zeros（nt,1）;

fori=1:

npl

forj=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%的要求

%%%%%%%%%%%%responsespectraofcallidar

%%%%%%%parameter

g=9.8;

m=1;

x0=0;

v0=0;

ww=2*pi./Tyz;

%%%%%%%%load

ag=at;%%%%%%%修改

%%%%%%%solution

fory=1:

nTyz

z=0.037;

w=ww（y）;

c=2*z*w;

k=w^2;

fori=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

fori=1:

nTyz

rspa（i）=maxas（i）;

end

%%%%%%%比较容差

fori=1:

nTyz

rcrsp（i）=abs（rspa（i）-Syz（i））/max（Syz（:

））;

end

jsnum=jsnum+1

max（rcrsp（:

））

end

%%%%%%%最终的反应谱与规范谱

%%%%%%%%%%%%responsespectraofcallidar

%%%%%%%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

fory=1:

nTyz

z=0.037;

w=ww（y）;

c=2*z*w;

k=w^2;

fori=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

fori=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）;

holdon;

plot（Tyz,rspa_S）;

4生成的人造地震波如图所示。

图6人造地震波和初始反应谱