地震动反应谱特征周期计算地震荷载

选取同一类场地、震中距相近的20条地震动记录，地震动峰值均为0.7m/s2，单自由度结构的阻尼比为2%、5%、10%和15%，周期范围为0.1s~10s，计算位移反应谱、速度反应谱和伪速度反应谱、加速度反应谱和伪加速度反应谱，并分析比较速度反应谱和伪速度反应谱的区别，以及加速度反应谱和伪加速度反应谱的区别。

一.反应谱计算与绘图
反应谱的计算采用Newmark-β法计算，对于单自由度体系使用杜哈美积分来求解实际更为方便。

MATLAB的计算程序如下所示：
clc
clear
kesai=0.15; %阻尼比
m=1;
[acc,dt,N]=peer2acc('F:matlab-
learn','RSN3753_LANDERS_FVR135.AT2')
%peer2acc为处理原始地震动数据的程序
save('acc2','acc')
load('acc2.mat');
gama = 0.5;
beta = 0.25;
alpha0 = 1/beta/dt^2;
alpha1 = gama/beta/dt;
alpha2 = 1/beta/dt;
alpha3 = 1/2/beta - 1;
alpha4 = gama/beta - 1;
alpha5 = dt/2*(gama/beta-2);
alpha6 = dt*(1-gama);
alpha7 = gama*dt;
peak=9.8*max(abs(acc));
acc=acc*0.7/peak;
n=length(acc);
p=-m*9.8*acc;
j=0;
for T=0.1:0.01:10
j=j+1;
wn=2*pi/T;
k=m*wn^2;
c=kesai*2*m*wn;
Keq=k+ alpha0*m + alpha1*c;
wD=wn*(1-kesai^2)^0.5;
d=zeros(n,1);
v=zeros(n,1);
a=zeros(n,1);
for i=2:n
t=0.002*(i-1);
f=p(i) + m*(alpha0*d(i-1)+alpha2*v(i-1)+alpha3*a(i-1))+
c*(alpha1*d(i-1)+alpha4*v(i-1)+alpha5*a(i-1)); d(i) =f/Keq; %Newmark-β的计算程序
a(i) = alpha0*(d(i)-d(i-1))-alpha2*v(i-1)-alpha3*a(i-1);
v(i) = v(i-1) + alpha6*a(i-1) + alpha7*a(i);
end
sd(j)=max(abs(d)); %位移反应谱
sv(j)=max(abs(v)); %速度反应谱
sa(j)=max(abs(a)); %加速度反应谱
SA(j)=wn^2*sd(j); %伪加速度反应谱
SV(j)=wn*sd(j); %伪速度反应谱
end
选取的地震动记录如图
地震动记录一般在PEER网站下载。

下载只需要注册邮箱号即可免费
下载。

图1 地震动记录实例
位移反应谱
速度，伪速度反应谱
加速度，伪加速度反应谱
由上图可知，伪速度反应谱在周期较小时和速度反应谱接近，伪速度反应谱稍大于速度反应谱。

随着周期的增大，速度反应谱的值和伪速度反应谱之间的差别增大，速度反应谱远大于伪速度反应谱。

伪加速度反应谱和加速度反应谱的关系类似。

比较反应谱在阻尼不同时变化
位移反应谱
速度反应谱
伪速度反应谱
加速度反应谱
伪加速度反应谱
由上图可知，位移反应谱，速度反应谱，伪速度反应谱，加速度反应谱，伪加速度反应谱随着阻尼比的增大而减小，这是由于阻尼比越大，消耗的地震动能量越多，反应谱相应的减小。

二.设计谱
真实的设计谱是通过大量的地震动记录来计算得到的。

在这里我们通过20条地震动记录来大致计算设计谱。

设计谱计算程序如下所示。

下载的地震动记录为后缀.AT2的文件，程序里需要每次调用一个AT2文件，程序的开头部分即为实现此功能的程序。

通过取二十条地震动计算得到的伪加速度谱均值来得到设计谱。

设计谱在实质上是伪加速度谱。

%reponse spectrum calculate
clc
clear
kesai=0.15;
m=1;
cd('F:\matlab-learn');
MatFileList = dir('*.AT2'); %依次读取二十条地震动记录的程序
FileNum = length(MatFileList);
FileName = struct([]);VariableName = struct([]);datetimenow = struct([]);
stringcaselist =[];
for num = 1:FileNum
FileName{num} = MatFileList(num,1).name;
loadstring = [ '[acc,dt,N]=peer2acc(' '''F:matlab-learn''' ',''' FileName{num} ''');'];
eval(loadstring);
gama = 0.5;
beta = 0.25;
alpha0 = 1/beta/dt^2;
alpha1 = gama/beta/dt;
alpha2 = 1/beta/dt;
alpha3 = 1/2/beta - 1;
alpha4 = gama/beta - 1;
alpha5 = dt/2*(gama/beta-2);
alpha6 = dt*(1-gama);
alpha7 = gama*dt;
peak=9.8*max(abs(acc));
acc=acc*0.7/peak;
n=length(acc);
p=-m*9.8*acc;
j=0;
for T=0.1:0.01:10
j=j+1;
wn=2*pi/T;
k=m*wn^2;
c=kesai*2*m*wn;
Keq=k+ alpha0*m + alpha1*c;
wD=wn*(1-kesai^2)^0.5;
d=zeros(n,1);
v=zeros(n,1);
a=zeros(n,1);
for i=2:n
t=0.002*(i-1);
f=p(i) + m*(alpha0*d(i-1)+alpha2*v(i-1)+alpha3*a(i-1))+
c*(alpha1*d(i-1)+alpha4*v(i-1)+alpha5*a(i-1)); d(i) =f/Keq;
a(i) = alpha0*(d(i)-d(i-1))-alpha2*v(i-1)-alpha3*a(i-1);
v(i) = v(i-1) + alpha6*a(i-1) + alpha7*a(i);
end
sd(j)=max(abs(d));
% sv(j)=max(abs(v));
% sa(j)=max(abs(a));
SA(j)=wn^2*sd(j);
end
ReusltData(:,num)=SA';
end
T=0.1:0.01:10;
figure(1)
plot(T,ReusltData)
for m=1:991
alfai(m)=max(ReusltData(m,:))/9.8;
end
figure(2)
plot(T,alfai)
计算二十条地震动的伪加速度反应谱：
二十条伪加速度反应谱
求每个周期T处的平均的伪加速度来求设计谱。

设计谱。