地震动反应谱分析的几种积分方法

2012-05-27 14:27阅读：

http://blog.sina.cn/dpool/blog/u/2279091464

第一次写博客，最近老师布置的地震工程学的大作业，要求编写几种不同的数值积分方法求解地震反应谱，到目前只做了前两个小题，以下是用matlab写的几个不同积分法的函数，可以比较几种积分法的结果，与大家分享，不足之处请读者慷慨指出。

（1）精细积分法
function [Umax,Vmax,Amax]= JXJF(t,g,ksi,T0)
dt=0.02;
du0=0;
u0=0;
Y=[du0;u0];
%YY(1,1)=Y(1,1)
%YY(2,1)=Y(2,1)
du(1)=du0;
u(1)=u0;
LT=length(T0);
lt=length(t)-1;
for j=1:LT
w(j)=6.28/T0(j);
%计算迭代参数
%首先定义矩阵A,注意力向量Fe（t）是变化的，不在循环外定义。
A=[-2*ksi*w(j),-w(j)^2;1,0];
invA=inv(A);
T=expm(A*dt);
fo

r i=1:lt
Ftk=[-g(i);0];
rk=[(-g(i+1)+g(i))/dt;0];
Y=T*[Y+invA*(Ftk+invA*rk)]-invA*[Ftk+rk*dt+invA*rk];
du(i+1)=Y(1,1);
u(i+1)=Y(2,1);
end
Umax(j)=abs(max(u));
%Vmax(j)=abs(max(du)); 相对速度
%Amax(j)=abs(max(ddu)); 相对加速度
Vmax(j)=Umax(j)*w(j); %伪速度，反应谱用
Amax(j)=Umax(j)*w(j)*w(j); %伪加速度，反应谱用
end
plot(T0,Umax)
legend('阻尼比0.02');
xlabel('周期Tn/s'),ylabel('U/m');
title('精细积分法相对位移反应谱');
grid on
figure
plot(T0,Vmax)
legend('阻尼比0.02');
xlabel('周期Tn/s'),ylabel('V/（m/s）');
title('精细积分法伪速度反应谱');
grid on
figure
plot(T0,Amax)
legend('阻尼比0.02');
xlabel('周期Tn/s'),ylabel('A/(m/s^2)');
title('精细积分法法伪加速度反应谱');
grid on
end
输入el-centrol-NS分量：

得到以下反应谱：

（2）Newmark-b法
function [Umax,Vmax,Amax]= FYP(t,g,ksi,T0,dt,u0,du0)
%用线性加速度法
LT=length(T0)
lt=length(t)-1
for j=1:LT
w(j)=6.28/T0(j);
%首先定义初始参数,初始加速度
ddu0=-g(1)-2*ksi*w(j)*du0-w(j)^2*u0;
u(1)=u0;
du(1)=du0;
ddu(1)=ddu0;
%计算迭代参数
kba=w(j)^2+3/dt*2*ksi*w(j)+6/dt^2;
a=6/dt+3*2*ksi*w(j);
b=3+dt/2*2*ksi*w(j);
for i=1:lt
deltaP(i)=-g(i+1)+g(i)+a*du(i)+b*ddu(i);
deltaU(i)=deltaP(i)/kba;
deltaDU(i)=3/dt*deltaU(i)-3*du(i)-dt/2*ddu(i);
deltaDDU(i)=6/dt^2*(deltaU(i)-dt*du(i))-3*ddu(i);
u(i+1)=u(i)+deltaU(i);
du(i+1)=du(i)+deltaDU(i);
ddu(i+1)=ddu(i)+deltaDDU(i);
end
Umax(j)=abs(max(u));
%Vmax(j)=abs(max(du)); 相对速度
%Amax(j)=abs(max(ddu)); 相对加速度
Vmax(j)=Umax(j)*w(j); %伪速度，反应谱用
Amax(j)=Umax(j)*w(j)*w(j); %伪加速度，反应谱用
end
plot(T0,Umax)
legend('阻尼比0');
xlabel('周期Tn/s'),ylabel('U/m');
title('Newmark-b法相对位移反应谱');
grid on
figure
plot(T0,Vmax)
legend('阻尼比0');
xlabel('周期Tn/s'),ylabel('V/（m/s）');
title('Newmark-b法伪速度反应谱');
grid on
figure
plot(T0,Amax)
legend('阻尼比0');
xlabel('周期Tn/s'),ylabel('A/(m/s^2)');
title('Newmark-b法伪加速度反应谱');
grid on
end
输入el-centrol-NS分量，可得以下反应谱：

（3）广义a方法
function [Umax,Vmax,Amax]= GuangyiA(t,g,ksi,T0)
dt=0.02;
du0=0;
u0=0;
LT=length(T0);
lt=length(t)-1;
ro=0.5;
am=(2*ro-1)/(1+ro);
af=ro/(1+ro);
beta=0.25*(1-am+af)^2;
gama=0.5-am+af;
for j=1:LT
w(j)=6.28/T0(j);
%首先定义初始参数,初始加速度
ddu0=-g(1)-2*ksi*w(j)*du0-w(j)^2*u0;
u(1)=u0;
du(1)=du0;
ddu(1)=ddu0;
%计算迭代参数
%kba=w(j)^2*(1-af)+3/dt*2*ksi*w(j)+6/dt^2;
Ma=1-am;
Ca=2*ksi*w(j)*(1-af);
Ka=w(j)^2*(1-af);
kba=1/(beta*dt*dt)*Ma+gama/(beta*dt)*Ca+Ka;
%a=6/dt+3*2*ksi*w(j);
%b=3+dt/2*2*ksi*w(j);
a=1/(beta*dt)*Ma+gama/beta*Ca;
b=1/(2*beta)*Ma-(2*beta-gama)/(2*beta)*dt*Ca;
for i=1:lt
Rk=-g(i)-(ddu(i)+2*ksi*w(j)*du(i)+w(j)^2*u(i));
deltaFa=(1-af)*(-g(i+1)+g(i))+Rk;
deltaP(i)=deltaFa+a*du(i)+b*ddu(i);
deltaU(i)=deltaP(i)/kba;
deltaDU(i)=gama/(beta*dt)*deltaU(i)-gama/beta*du(i)+dt*(1-gama/(2*beta))*ddu(i);
deltaDDU(i)=1/(beta*dt^2)*(deltaU(i)-dt*du(i))-1/(2*beta)*ddu(i);
u(i+1)=u(i)+deltaU(i);
du(i+1)=du(i)+deltaDU(i);
ddu(i+1)=ddu(i)+deltaDDU(i);
end
Umax(j)=abs(max(u));
%Vmax(j)=abs(max(du)); 相对速度
%Amax(j)=abs(max(ddu)); 相对加速度
Vmax(j)=Umax(j)*w(j); %伪速度，反应谱用
Amax(j)=Umax(j)*w(j)*w(j); %伪加速度，反应谱用
end
plot(T0,Umax)
legend('阻尼比0.02');
xlabel('周期Tn/s'),ylabel('U/m');
title('广义a法相对位移反应谱');
grid on
figure
plot(T0,Vmax)
legend('阻尼比0.02');
xlabel('周期Tn/s'),ylabel('V/（m/s）');
title('广义a法伪速度反应谱');
grid on
figure
plot(T0,Amax)
legend('阻尼比0.02');
xlabel('周期Tn/s'),ylabel('A/(m/s^2)');
title('广义a法伪加速度反应谱');
grid on
end