地震工程学反应谱和地震时程波的相互转化matlab编程.docx
《地震工程学反应谱和地震时程波的相互转化matlab编程.docx》由会员分享，可在线阅读，更多相关《地震工程学反应谱和地震时程波的相互转化matlab编程.docx（17页珍藏版）》请在冰豆网上搜索。
地震工程学反应谱和地震时程波的相互转化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
ifT1arf_jg=lmt2*arfmax;
end
%%%%第三段T2~T3
ifT2arf_jg=（（Tg/T_jg）^sjzs）*lmt2*arfmax;
end
%%%%第四段T3～6.0
ifT3arf_jg=（lmt2*0.2^sjzs-lmt1*（T_jg-5*Tg））*arfmax;
end
%%%%第五段6.0～
if6.0arf_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人造地震波和初始反应谱
升级会员 