MATLAB程序精確法求解反應譜-2008-04-06（本文程序僅供參考，請勿直接抄襲）2010/01/21 09:52 P.M.1. 反應譜的概念 反應譜是在1932年由M.A.Biot引入的，它是用來描述地面運動及其對結構的效應的一種實用工具。現在，反應譜作為地震工程的核心概念，提供了一種方便的手段概括所有可能的線性單自由度體系對地面運動的某個特定分量的峰值反應。它還提供了一種實用的方法，將結構動力學的知識應用于結構的設計以及建筑規范中側向力條文的制定。 某個反應量的峰值作為體系的固有振動周期Tn,（或者循環頻率fn）那樣的相關參數的函數圖形，稱為該反應量的反應譜。每一個這樣的圖形針對的是有一個具有固定阻尼比的單自由度體系，多個具有不同阻尼比的這類圖形聯合起來就能覆蓋實際結構中遇到的阻尼值范圍。2. 反應譜的計算2.1反應譜數值計算方法 計算反應譜的方法有很多，又卷積計算法，傅立葉變換法，線性加速度法，中點加速度法，精確法等。2.2精確法 本文中采用精確法做計算，該方法是N.C.Nigam和P.C.Jennings于1969年提出的，此法的出發點是把地面運動的加速度記錄相鄰點間的值用分段線性差值表示，從而獲得地面運動的連續表達式。基于方程本身基礎上進行，得到的結果全部采用精確的分析方法，沒有任何的舍入誤差，也不會產生任何的截斷誤差，所謂精確法就是指在這個意義上式精確的而然。正因為這種方法不會引起數值計算的誤差，所以它有較高的精度，只要進行較少的運算就可以達到采用其他方法需要較多次運算才能達到的精度。“由于在sohu博客上的文章發表后，陸續有問參考文獻的郵件，因此將參考文獻pdf版放上來供大家學習、參考，請勿用于商業目的。下載鏈接見地震動的譜分析入門 /linshibin/blog/item/3b03ce3482f6aed6a3cc2b15.html強震觀測與分析原理 /linshibin/blog/item/4a02ad3acd7c692fb8998f01.html”ResponSespectrumProgram（精確法求解）% 反應譜 精確法 程序 Begin With matlab6.5%clear% *讀入地震記錄*fid = fopen(CHI010.txt);Accelerate,count = fscanf(fid,%g); %count 讀入的記錄的量Accelerate=9.8*Accelerate; %單位統一為 m和stime=0:0.005:(count-1)*0.005; %單位 s% *精確法計算各反應*%初始化各儲存向量Displace=zeros(1,count); %相對位移Velocity=zeros(1,count); %相對速度AbsAcce=zeros(1,count); %絕對加速度% *A,B矩陣*DampA=0,0.05,0.1; %三個阻尼比TA=0.0:0.05:6; %TA=0.000001:0.02:6; %結構周期Dt=0.005; %地震記錄的步長%記錄計算得到的反應，MDis為某阻尼時最大相對位移，MVel為某阻尼%時最大相對速度，MAcc某阻尼時最大絕對加速度，用于畫圖MDis=zeros(3,length(TA);MVel=zeros(3,length(TA);MAcc=zeros(3,length(TA);j=1; %在下一個循環中控制不同的阻尼比for Damp=0,0.05,0.1 t=1; %在下一個循環中控制不同的結構自振周期 for T=0.0:0.05: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*Damp2-1)/(Frcy2*Dt); %計算公式化簡 d_3t=Damp/(Frcy3*Dt); B=zeros(2,2); B(1,1)=e_t*(d_f+Damp/Frcy)*s/DamFrcy+(2*d_3t+1/Frcy2)*c)-2*d_3t; B(1,2)=-e_t*(d_f*s/DamFrcy+2*d_3t*c)-1/Frcy2+2*d_3t; B(2,1)=e_t*(d_f+Damp/Frcy)*(c-Damp/sqrt(1-Damp2)*s)-(2*d_3t+1/Frcy2)*(DamFrcy*s+Damp*Frcy*c)+1/(Frcy2*Dt); B(2,2)=e_t*(1/(Frcy2*Dt)*c+s*Damp/(Frcy*DamFrcy*Dt)-1/(Frcy2*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)*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*Frcy*Velocity(i+1)-Frcy2*Displace(i+1); end MDis(j,t)=max(abs(Displace); MVel(j,t)=max(abs(Velocity); if T=0.0 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 allfigure %繪制地震記錄圖plot(time(:),Accelerate(:) title(PEER STRONG MOTION DATABASE RECORD-CHI010)xlabel(time(s)ylabel(acceleration(g)gridfigure %繪制位移反應譜 plot(TA,MDis(1,:),-.b,TA,MDis(2,:),-r,TA,MDis(3,:),:k)title(Displacement)xlabel(Tn(s)ylabel(Displacement(m)legend(=0,=0.05,=0.1)gridfigure %繪制速度反應譜 plot(TA,MVel(1,:),-.b,TA,MVel(2,:),-r,TA,MVel(3,:),:k)title(Velocity)xlabel(Tn(s)ylabel(velocity(m/s)legend(=0,=0.05,=0.1)gridfigure %繪制絕對加速度反應譜 plot(TA,MAcc(1,:),-.b,TA,MAcc(2,:),-r,TA,MAcc(3,:),:k)title(Absolute Acceleration)xlabel(Tn(s)ylabel(absolute acceleration(m/s2)legend(=0,=0.05,=0.1)gridfigure %繪制標準加速度反應譜 M=max(abs(Accelerate); %地震記錄最大值plot(TA,MAcc(1,:)/M,-.b,TA,MAcc(2,:)/M,-r,TA,MAcc(3,:)/M,:k)title(Normalized Absolu