MATLAB程序精确法求解反应谱

MATLAB程序精确法求解反应谱-2008-04-06（本文程序仅供参考，请勿直接抄袭）2010/01/21 09:52 .
1. 反应谱的概念
反应谱是在1932年由引入的，它是用来描述地面运动及其对结构的效应的一种实用工具。

现在，反应谱作为地震工程的核心概念，提供了一种方便的手段概括所有可能的线性单自由度体系对地面运动的某个特定分量的峰值反应。

它还提供了一种实用的方法，将结构动力学的知识应用于结构的设计以及建筑规范中侧向力条文的制定。

某个反应量的峰值作为体系的固有振动周期Tn,（或者循环频率fn）那样的相关参数的函数图形，称为该反应量的反应谱。

每一个这样的图形针对的是有一个具有固定阻尼比的单自由度体系，多个具有不同阻尼比的这类图形联合起来就能覆盖实际结构中遇到的阻尼值范围。

2. 反应谱的计算
反应谱数值计算方法
计算反应谱的方法有很多，又卷积计算法，傅立叶变换法，线性加速度法，中点加速度法，精确法等。

精确法
本文中采用精确法做计算，该方法是和于1969年提出的，此法的出发点是把地面运动的加速度记录相邻点间的值用分段线性差值表示，从而获得地面运动的连续表达式。

基于方程本身基础上进行，得到的结果全部采用精确的分析方法，没有任何的舍入误差，也不会产生任何的截断误差，所谓精确法就是指在这个意义上式精确的而然。

正因为这种方法不会引起数值计算的误差，所以它有较高的精度，只要进行较少的运算就可以达到采用其他方法需要较多次运算才能达到的精度。

由于在博客上的文章发表后，陆续有问参考文献的邮件，因此将参考文献版放上来供pdfsohu“大家学习、参考，请勿用于商业目的。

下载链接见强震观测与分析原理（精确法求解）地震动的谱分析入门%%%%%%%%%%%%%%%%%%%%% 反应谱精确法程序Begin
With %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
% ***********读入地震记录***********
fid = fopen('');
[Accelerate,count] = fscanf(fid,'%g'); %count 读入的记录的量
Accelerate=*Accelerate'; %单位统一为m和s
time=0::(count-1)*; %单位s
% ***********精确法计算各反应***********
初始化各储存向量%．
Displace=zeros(1,count); %相对位移
Velocity=zeros(1,count); %相对速度
AbsAcce=zeros(1,count); %绝对加速度
% ***********A,B矩阵***********
DampA=[0,,]; %三个阻尼比
TA=::6; %TA=::6; %结构周期
Dt=; %地震记录的步长
%记录计算得到的反应，MDis为某阻尼时最大相对位移，MVel为某阻尼
%时最大相对速度，MAcc某阻尼时最大绝对加速度，用于画图
MDis=zeros(3,length(TA));
MVel=zeros(3,length(TA));
MAcc=zeros(3,length(TA));
j=1; %在下一个循环中控制不同的阻尼比
for Damp=[0,,]
t=1; %在下一个循环中控制不同的结构自振周期
for T=::6
Frcy=2*pi/T ; %结构自振频率
DamFrcy=Frcy*sqrt(1-Damp*Damp); %计算公式化简
e_t=exp(-Damp*Frcy*Dt);
s=sin(DamFrcy*Dt);
c=cos(DamFrcy*Dt);
A=zeros(2,2);
A(1,1)=e_t*(s*Damp/sqrt(1-Damp*Damp)+c);
A(1,2)=e_t*s/DamFrcy;
A(2,1)=-Frcy*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)/(Frcy^2*Dt); %计算公式化简
d_3t=Damp/(Frcy^3*Dt);
B=zeros(2,2);
B(1,1)=e_t*((d_f+Damp/Frcy)*s/DamFrcy+(2*d_3t+1/Frcy^2)*c)-2*d_3t;
B(1,2)=-e_t*(d_f*s/DamFrcy+2*d_3t*c)-1/Frcy^2+2*d_3t;
B(2,1)=e_t*((d_f+Damp/Frcy)*(c-Damp/sqrt(1-Damp^2)*s)-(2*d_3t+1/Frcy^2)*(Dam Frcy*s+Damp*Frcy*c))+1/(Frcy^2*Dt);
B(2,2)=e_t*(1/(Frcy^2*Dt)*c+s*Damp/(Frcy*DamFrcy*Dt))-1/(Frcy^2*Dt);
for i=1:(count-1) %根据地震记录,计算不同的反应
Displace(i+1)=A(1,1)*Displace(i)+A(1,2)*Velocity(i)+B(1,1)*Accelerate(i)+B(1,2)*Accel erate(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*Frcy*Velocity(i+1)-Frcy^2*Displace(i+1);
end
MDis(j,t)=max(abs(Displace));
MVel(j,t)=max(abs(Velocity));
if T==
MAcc(j,t)=max(abs(Accelerate));
else
MAcc(j,t)=max(abs(AbsAcce));
end
Displace=zeros(1,count);%初始化各储存向量，避免下次不同周期计算时引用到前一个周期的结果
Velocity=zeros(1,count);
AbsAcce=zeros(1,count);
t=t+1;
end
j=j+1;
end
% ***********PLOT***********
close all
figure %绘制地震记录图
plot(time(:),Accelerate(:))
title('PEER STRONG MOTION DATABASE RECORD--CHI010')
xlabel('time(s)')
ylabel('acceleration(g)')
grid
figure %绘制位移反应谱
plot(TA,MDis(1,:),'',TA,MDis(2,:),'-r',TA,MDis(3,:),':k')
title('Displacement')
xlabel('Tn(s)')
ylabel('Displacement(m)')
legend('ζ=0','ζ=','ζ=')
grid
figure %绘制速度反应谱
plot(TA,MVel(1,:),'',TA,MVel(2,:),'-r',TA,MVel(3,:),':k')
title('Velocity')
xlabel('Tn(s)')
ylabel('velocity(m/s)')
legend('ζ=0','ζ=','ζ=')
grid
绘制绝对加速度反应谱figure %．
plot(TA,MAcc(1,:),'',TA,MAcc(2,:),'-r',TA,MAcc(3,:),':k')
title('Absolute Acceleration')
xlabel('Tn(s)')
ylabel('absolute acceleration(m/s^2)')
legend('ζ=0','ζ=','ζ=')
grid
figure %绘制标准加速度反应谱
M=max(abs(Accelerate)); %地震记录最大值
plot(TA,MAcc(1,:)/M,'',TA,MAcc(2,:)/M,'-r',TA,MAcc(3,:)/M,':k')
title('Normalized Absolute Acceleration')
xlabel('Tn(s)')
ylabel('Normalized absolute acceleration')
legend('ζ=0','ζ=','ζ=')
grid
%%%%%%%%%%%%%%%%%%%%%
End
With %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%。