地震工程学反应谱和地震时程波的相互转化matlab编程Word文档格式.docx
《地震工程学反应谱和地震时程波的相互转化matlab编程Word文档格式.docx》由会员分享，可在线阅读，更多相关《地震工程学反应谱和地震时程波的相互转化matlab编程Word文档格式.docx（17页珍藏版）》请在冰豆网上搜索。
%绝对加速度
%***********A,B矩阵***********
Damp=;
%阻尼比
TA=:
6;
%TA=:
%结构周期
Dt=;
%地震记录的步长
%记录计算得到的反应，MaxD为某阻尼时最大相对位移，MaxV为某阻尼最大相对速度，MaxA某阻尼时最大绝对加速度，用于画图
MaxD=zeros（3,length（TA））;
MaxV=zeros（3,length（TA））;
MaxA=zeros（3,length（TA））;
t=1;
forT=:
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==
MaxA（1,t）=max（abs（Accelerate））;
else
MaxA（1,t）=max（abs（AbsAcce））;
%初始化各储存向量，避免下次不同周期计算时引用到前一个周期的结果
t=t+1;
End
%***********PLOT***********
closeall
figure%绘制地震记录图
plot（time（:
）,Accelerate（:
））
title（'
PEERSTRONGMOTIONDATABASERECORD'
）
xlabel（'
time（s）'
ylabel（'
acceleration（g）'
）
grid
figure%绘制位移反应谱
plot（TA,MaxD（1,:
）,'
'
TA,MaxD（2,:
-r'
TA,MaxD（3,:
k'
Displacement'
Tn（s）'
Displacement（m）'
legend（'
ζ='
Grid
figure%绘制速度反应谱
plot（TA,MaxV（1,:
TA,MaxV（2,:
TA,MaxV（3,:
Velocity'
velocity（m/s）'
figure%绘制绝对加速度反应谱
plot（TA,MaxA（1,:
TA,MaxA（2,:
TA,MaxA（3,:
AbsoluteAcceleration'
absoluteacceleration（m/s^2）'
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=;
ifld==7
ifld==8
ifld==9
%%%%场地类别，设计地震分组选择
ifcd==1
iffz==1
Tg=;
end
iffz==2
iffz==3
ifcd==2
ifcd==3
ifcd==4
%%%%%%%%%
ceita=zn;
%%%%%阻尼比
lmt1=+/8;
iflmt1<
lmt1=0;
lmt2=1+/+*ceita）;
iflmt2<
lmt2=;
sjzs=+/+5*ceita）;
%%%%%分段位置T1T2T3
T1=;
T2=Tg;
T3=5*Tg;
T_jg=2*pi./pl;
%%%%第一段0～T1
ifT_jg<
=T1
arf_jg=*arfmax+（lmt2**arfmax）/*T_jg;
%%%%第二段T1～T2
ifT1<
T_jg&
T_jg<
=T2
arf_jg=lmt2*arfmax;
%%%%第三段T2~T3
ifT2<
=T3
arf_jg=（（Tg/T_jg）^sjzs）*lmt2*arfmax;
%%%%第四段T3～
ifT3<
=
arf_jg=（lmt2*^sjzs-lmt1*（T_jg-5*Tg））*arfmax;
%%%%第五段～
if<
T_jg
arf_jg=（lmt2*^sjzs-lmt1**Tg））*arfmax;
%%%%%%反应谱值拟加速度值
rs_z=arf_jg*;
3生成人造地震波主程序：
%%%主程序%%%%
%%%%确定需要控制的反应谱Sa（T）（T=T1,...,TM）的坐标点数M，反应谱控制容差rc
Tyz=[:
:
];
rc=;
nTyz=length（Tyz）;
ceita=;
%%%阻尼比：
nTyz
Syz（i）=r_s_1（2*pi/Tyz（i）,ceita,8,2,1）;
%%%%8度，2类场地，第1地震分组
%%%%%%变换的频率差：
2*pi*（可以保证长周期项5s附近有5项三角级数）；
%%%%频率变化范围N1=30,30**2*pi;
N2=3000,5000**2*pi
plc=2*pi*;
pl=30**2*pi:
*2*pi:
10000**2*pi;
npl=length（pl）;
P=;
%%%保证率
%%%%%%人造地震动持续时间40s，时间间隔：
Td=40;
dt=;
t=0:
40;
nt=length（t）;
%%%%%%%衰减包络函数
t1=8;
%%%%上升段
t2=8+24;
%%%%%平稳段;
下降段则为40－32＝8s
c=;
%%%%衰减段参数
nt
ift（i）<
=t1
f（i）=（t（i）/t1）^2;
ift（i）>
t1&
t（i）<
t2
f（i）=1;
end
=t2
f（i）=exp（-c*（t（i）-t2））;
%%%%%%%反应谱转换功率谱
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）;
%%%%%%%%%%%%%%合成地震动
at=zeros（nt,1）;
atj=zeros（nt,1）;
npl
fai（i）=rand
（1）*2*pi;
forj=1:
atj（j）=f（j）*Aw（i）*real（exp（sqrt（-1）*（pl（i）*t（j）+fai（i））））;
at=at+atj;
%%%%%%%计算反应谱验证是否满足rc在5%的要求,需要时程动力分析
%%%%%%%%%%%%responsespectraofcallidar
%%%%%%%parameter
g=;
m=1;
x0=0;
v0=0;
ww=2*pi./Tyz;
%%%%%%%%load
ag=at;
%%%%%%%修改
%%%%%%%solution
fory=1:
z=;
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;
maxas（y）=max（as）;
maxvs（y）=max（vs）;
maxxs（y）=max（xs）;
rspa（i）=maxas（i）;
%%%%%%%比较容差
rcrsp（i）=abs（rspa（i）-Syz（i））/max（Syz（:
））;
jsnum=1;
whilemax（rcrsp（:
））>
rc
%%%%%循环体函数
blxs=Syz./rspa;
forxsxs=1:
if2*pi/pl（xsxs）<
Tyz
（1）
blxs1（xsxs）=blxs
（1）;
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））;
if2*pi/pl（xsxs）>
Tyz（nTyz）
blxs1（xsxs）=blxs（nTyz）;
Aw=Aw.*blxs1;
%%%%%%%%%%%%%%合成地震动
at=zeros（nt,1）;
atj=zeros（nt,1）;
%%%%%%%计算反应谱验证是否满足rc在5%的要求
%%%%%%%%%%%%responsespectraofcallidar
%%%%%%%parameter
g=;
jsnum=jsnum+1
max（rcrsp（:
））
%%%%%%%最终的反应谱与规范谱
%%Tjs=:
%%nTjs=length（Tjs）;
m=1;
x0=0;
v0=0;
ww=2*pi./Tyz;
%%%%%%%%load
ag=at;
%%%%%%%solution
fory=1:
rspa（i）=maxas（i）/g;
rspa_S（i）=r_s_1（2*pi/Tyz（i）,ceita,8,2,1）/g;
subplot（2,1,1）;
plot（t,at）;
subplot（2,1,2）;
plot（Tyz,rspa）;
holdon;
plot（Tyz,rspa_S）;
4生成的人造地震波如图所示。
图6人造地震波和初始反应谱
升级会员 