地震工程学反应谱和地震时程波的相互转化matlab编程.docx
《地震工程学反应谱和地震时程波的相互转化matlab编程.docx》由会员分享,可在线阅读,更多相关《地震工程学反应谱和地震时程波的相互转化matlab编程.docx(18页珍藏版)》请在冰豆网上搜索。
![地震工程学反应谱和地震时程波的相互转化matlab编程.docx](https://file1.bdocx.com/fileroot1/2023-1/12/045e39a9-5c63-4db2-8570-fc1b8b5fbf73/045e39a9-5c63-4db2-8570-fc1b8b5fbf731.gif)
地震工程学反应谱和地震时程波的相互转化matlab编程
地震工程学-反应谱和地震时程波的相互转化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人造地震波和初始反应谱