結(jié)構(gòu)動力學(xué)大作業(yè)

1、 結(jié)構(gòu)動力學(xué)大作業(yè) 班 級 土木卓越1201班 學(xué) 號 U201210323 姓 名 陳祥磊 指導(dǎo)老師 葉 昆 2014.12.30 結(jié)構(gòu)動力學(xué)大作業(yè) SDOF體系在任意荷載作用下的動力響應(yīng)一、結(jié)構(gòu)參數(shù) 計算結(jié)構(gòu)為右圖所示的1、 2、 3、結(jié)構(gòu)參數(shù)中；。2、 確定各階頻率和振型 多自由度體系自由振動時的運動方程為 . 寫成矩陣形式即為 假設(shè)此方程的解答為，帶入到運動方程中得到振動方程 此方程要有非零解必須滿足頻率方程，可解得各階主頻率再根據(jù) 可求出結(jié)構(gòu)的主振型。在主振型中，通常將最后一個位移值設(shè)定為1，只要在程序中加入下列語句：MDOF.YMatrix(:,i)=MDOF.YMatrix(:,

2、i)/MDOF.YMatrix(MDOF.ND,i)運行程序之后得到如下結(jié)果：1、各階頻率和周期W112.7290261T10.493610843W237.15584832T20.169103535W358.57252468T30.107271888W475.24400343T40.083504133W585.81966052T50.0732138222、 各階陣型Y1Y2Y3Y4Y50.284629677-0.8308300261.30972147-1.6825070661.9189859470.546200349-1.0881559210.372785591.397877389-3.228

3、707410.763521118-0.594351144-1.20361560.5211085583.5133370920.9189859470.309721468-0.7153703-1.830830026-2.6825070611111由主振型可以求得廣義質(zhì)量矩陣和廣義剛度矩陣，均為對角矩陣。程序如下所示：MDOF.MGMatrix = zeros(MDOF.ND,MDOF.ND);for i = 1:MDOF.ND MDOF.MGMatrix(i,i)= (transpose(MDOF.YMatrix(:,i)*MDOF.MMatrix*MDOF.YMatrix(:,i);endMDOF

4、.KGMatrix = zeros(MDOF.ND,MDOF.ND);for i = 1:MDOF.ND MDOF.KGMatrix(i,i)= (transpose(MDOF.YMatrix(:,i)*MDOF.KMatrix*MDOF.YMatrix(:,i);end三、使用能量法計算近似的一階頻率； 使用能量法求解一階頻率就是瑞利法，Rayleigh法主要用于求的近似解，原理是能量守恒定律：；最終求解的表達(dá)式為： 因為采用了近似解答，相當(dāng)于在精確解條件下添加了約束，提高了結(jié)構(gòu)的整體剛度，所以第一頻率會大于精確值。在實際計算時，假設(shè)的位移幅值函數(shù)決定最終計算的精確度，越接近實際的位移模式，

5、Rayleigh法的計算結(jié)果越精確。假設(shè)的位移幅值函數(shù)必須滿足以下幾點：必須滿足運動邊界條件，即幾何和自然邊界條件；所設(shè)位移幅值函數(shù)必須與實際振型形狀大致相近。通常可取結(jié)構(gòu)在耨個靜位移作用下的彈性曲線作為的近似表達(dá)式，此式應(yīng)變能可用相應(yīng)荷載所做的功來代替。 計算時采用列表計算和編程計算兩種方法互為佐證以證明結(jié)果的正確性，計算結(jié)果如下：1、列表計算法 層數(shù)質(zhì)量m MiXiMiXi2W11000000490000002000*E60.02450.024524500600.2521000000392000002000*E60.01960.0441441001944.81310000002940000

6、02000*E60.01470.0588588003457.4441000000196000002000*E60.00980.0686686004705.965100000098000002000*E60.00490.0735735005402.25匯總26950016110.7112.80 如上表所示得出的第一頻率的近似值為，與第一題中的結(jié)果比較，相對誤差為，即計算結(jié)果很接近精確解。 2、編程法 Rayleigh法計算過程簡單，易于通過程序?qū)崿F(xiàn)。變成后進行計算可得結(jié)果為 上面分別通過兩種不同的方法求得頻率值，現(xiàn)在考察計算結(jié)果與實際的吻合程度，一般建筑結(jié)構(gòu)的自振周期 則五層建筑結(jié)構(gòu)的第一主振型

7、的周期大概介于之間，計算結(jié)果中的符合一般經(jīng)驗。另外據(jù)經(jīng)驗可知一二階主振頻率之間有如下關(guān)系： 結(jié)算結(jié)果中的，接近1/3，進一步說明計算結(jié)果符合實際工況。4、 振型分解法 結(jié)構(gòu)有固定的自振頻率和周期，取決于結(jié)構(gòu)的參數(shù)和形式。結(jié)構(gòu)任意的自由振動模式均可以分解為主振型的組合，類似于高階方程組的基本解組，將結(jié)構(gòu)的位移通過正則坐標(biāo)表示為基本振型的組合 選定地震動數(shù)據(jù)（文件為Ground_MotionsIMPVALLH-E01140.at2）7807個，間隔時間為0.005s,將地震動幅值調(diào)整為0.3g，使用振型分解法計算結(jié)構(gòu)的地震響應(yīng)。地震動數(shù)據(jù)給定的地震波在各個時刻的加速度值。由給定的數(shù)據(jù)可以繪出地震動

8、的加速度和速度關(guān)于時間的曲線：加速度與時間關(guān)系圖：速度與時間關(guān)系圖：運行程序SDOF_Time_History_Analysis之后，得到結(jié)構(gòu)的位移與時間關(guān)系圖放大之后可得細(xì)節(jié)圖5、 矩陣迭代法 矩陣迭代法是采用逐步逼近的方法確定結(jié)構(gòu)的頻率和振型。體系做自由振動時各質(zhì)點的位移幅值為 這是一個迭代關(guān)系式，通過假定最初的帶入計算，可以求出，繼而不斷求得，當(dāng)與足夠接近時即可將作為第一振型。其運算過程如下：首先假定一個振型，通常假定為代入上式右邊，進行求解后即可得到和主振型的第一次近似值；然后以第一次近似值代入上式進行計算，得到和主振型的第二次近似值，如此下去，直到最后兩次的計算結(jié)果接近為止。當(dāng)一個振

9、型求得以后，利用振型的正交性，求出高階的頻率和振型。 當(dāng)結(jié)構(gòu)層數(shù)較多且計算精度要求較高時，宜用編程法來計算。思路一：由于在求解高階頻率和高階振型時，例如第K振型時，需要先利用主振型的正交性解出關(guān)于的表達(dá)式，以帶入到中以求出新的階迭代關(guān)系式；同理通過迭代即可求出第K階主振型。編程計算時，如果可以求出在考慮正交性后的新柔度矩陣，即可按照求解第一振型時的語言計算，此時只要設(shè)定一系列的循環(huán)語句即可。求解第二振型的程序如下：for i=1:4 for j=1:4 l=0; for n=5:MDOF.ND q=MDOF.XMatrix(j,1)* MDOF.DMatrix(i,n)/MDOF.XMatri

10、x(n,1);l=l+q; end MDOF.DMatrix(i,j)= MDOF.DMatrix(i,j)-l; endendMDOF.XMatrix(1:4,2)=1 -1 -1 -1;n2=0;while(max(abs(MDOF.XMatrix(1:4,2)-MDOF.SMatrix(1:4,2)10e-8)MDOF.SMatrix(1:4,2) = MDOF.XMatrix(1:4,2);MDOF.XMatrix(1:4,2) = MDOF.DMatrix(1:4,1:4)*MDOF.MMatrix(1:4,1:4)*MDOF.XMatrix(1:4,2);n2=n2+1;endMD

11、OF.WMatrix(2,1)=sqrt(MDOF.DMatrix(1:4,1:4)*MDOF.MMatrix(1:4,1:4)*MDOF.XMatrix(1:4,2)MDOF.XMatrix(1:4,2);MDOF.YMatrix(1:4,2) = (MDOF.WMatrix(2,1)(2*n2)*MDOF.XMatrix(1:4,2);其中二階時新的柔度矩陣的表達(dá)式：MDOF.DMatrix(i,j)= MDOF.DMatrix-MDOF.XMatrix(j,1)*MDOF.DMatrix(i,n)/MDOF.XMatrix(n,1) 這種方法可以求解出第二振型的頻率和主振型，但是在更高階

12、的振型求解中，由于關(guān)于的表達(dá)式是一個K階矩陣，求出的表達(dá)式很復(fù)雜，并且在代入的迭代關(guān)系式時，難以求出新柔度矩陣，因此難以應(yīng)用。思路二：W112.72902610T10.493610843W237.15584827T20.169103535W358.57252470T30.107271888W475.24400344T40.083504133W585.81966023T50.073213822Y1Y2Y3 Y4Y50.284636896-0.830799561.309716066-1.682287411.9178054020.546210873-1.088078280.3726322871.39

13、8035292-3.22719520.763529809-0.59425090-1.20365760.520635183.5123626270.9189898260.309781116-0.71524594-1.83059827-2.6822307011111上面迭代法求出的頻率與直接求杰出的頻率幾乎一模一樣，在限定的與差距足夠小時，求出的頻率可作為真值。6、 繪制反應(yīng)譜 研究單自由度體系在地震力作用下的地震反應(yīng)，繪制出不同周期的結(jié)構(gòu)在地震力作用下的地震反應(yīng)的最大加速度，即可得到一條地震動的反應(yīng)譜。單自由度體系的運動方程為 表示為：，由此方程可得慣性力表達(dá)式 絕對加速度：；得 若給定地震時地面

14、運動的加速度記錄和體系的阻尼比，可計算出單自由度體系的最大絕對加速度與體系自振周期的關(guān)系曲線，即為加速度反應(yīng)譜。7、 求解地震力 在上面求出了地震動作用下結(jié)構(gòu)的主振型以及頻率、周期值之后，可以求出振興參與系數(shù)： 再根據(jù)自振周期求出，由公式，即可求出N組地震作用，分別計算每組地震力作用下的內(nèi)力，按計算總效應(yīng)。Y1Y2Y3 Y4Y5F1205189.4527454445.2597472504.5893302947.955792603.59466F2393755.6763595196.7126134488.8284-251698.2579-155806.2023F3550422.1568325096

15、.6524-434225.0776-93829.48553169541.4416F4662496.7076-169410.6478-258082.1718329654.6118-129448.4711F5720899.7151-546977.4147360767.2325-180057.464248256.52569 F1207770.29997171.86 791617.74 552291.86 284043.83 【附錄】1、MDOF_Eigen_Analysis%clc;clear;close;%global MDOF%MDOF.ND = 5;%MDOF.MVec = zeros(MDO

16、F.ND,1);MDOF.KVec = zeros(MDOF.ND,1); for i = 1:5 MDOF.MVec(i) = 1000E3;end for i = 1:5 MDOF.KVec(i) = 1.0*2000E6;end%MDOF.MMatrix = zeros(MDOF.ND, MDOF.ND);MDOF.KMatrix = zeros(MDOF.ND, MDOF.ND);%for i = 1:MDOF.ND MDOF.MMatrix(i,i) = MDOF.MVec(i);end%for i = 1:MDOF.ND if i = 1 MDOF.KMatrix(i,i ) =

17、MDOF.KVec(i ) + MDOF.KVec(i+1); MDOF.KMatrix(i,i+1) = -MDOF.KVec(i+1); else if i = MDOF.ND MDOF.KMatrix(i,i-1) = -MDOF.KVec(i); MDOF.KMatrix(i,i ) = MDOF.KVec(i); else MDOF.KMatrix(i,i-1) = -MDOF.KVec(i ); MDOF.KMatrix(i,i ) = MDOF.KVec(i ) + MDOF.KVec(i+1); MDOF.KMatrix(i,i+1) = -MDOF.KVec(i+1); en

18、d endend%MDOF.Eigen_Vec, MDOF.Eigen_Val = eig(MDOF.MMatrixMDOF.KMatrix);%MDOF.WVec = sqrt(diag(MDOF.Eigen_Val);MDOF.TVec = 2*pi./MDOF.WVec;MDOF.YMatrix = zeros(MDOF.ND,MDOF.ND );for i=1:MDOF.NDMDOF.YMatrix(:,i)=null(MDOF.KMatrix-MDOF.WVec(i,1)*MDOF.WVec(i,1)*MDOF.MMatrix);endMDOF.MGMatrix = zeros(MD

19、OF.ND,MDOF.ND);for i = 1:MDOF.ND MDOF.MGMatrix(i,i) = (transpose(MDOF.YMatrix(:,i)*MDOF.MMatrix*MDOF.YMatrix(:,i);endMDOF.KGMatrix = zeros(MDOF.ND,MDOF.ND);for i = 1:MDOF.ND MDOF.KGMatrix(i,i) = (transpose(MDOF.YMatrix(:,i)*MDOF.KMatrix*MDOF.YMatrix(:,i);endMDOF.DMatrix = zeros(MDOF.ND,MDOF.ND);MDOF

20、.FGMatrix = zeros(MDOF.ND,MDOF.ND);%Test_01 = transpose(MDOF.Eigen_Vec)*MDOF.MMatrix*MDOF.Eigen_Vec;Test_02 = transpose(MDOF.Eigen_Vec)*MDOF.KMatrix*MDOF.Eigen_Vec;2、SDOF_Time_History_Analysis%clc;clear;close;%global EWave%File_Name=E:Matlab_CodeMatlab_CodeGround_Motions_LibraryIMPVALLH-E01140.AT2;%

21、路徑修改%fid = fopen(File_Name,r);% EWave.Str1 = fgetl (fid);EWave.Str2 = fgetl (fid);EWave.Str3 = fgetl (fid);EWave.NPTs=fscanf(fid, %i ,1);EWave.DT =fscanf(fid, %f ,1);EWave.Str4 = fgetl (fid);%EWave.Acel=zeros(EWave.NPTs,1); EWave.Time=zeros(EWave.NPTs,1); for i=1:EWave.NPTs EWave.Time(i)= (i-1)*EWav

22、e.DT; EWave.Acel(i)= fscanf(fid, %f ,1);end%EWave.AMax=max(abs(EWave.Acel);%EWave.Acel=EWave.Acel*0.3/EWave.AMax;% plot(EWave.Time,EWave.Acel);%global MDOF%MDOF.ND = 5;%MDOF.MVec = zeros(MDOF.ND,1);MDOF.KVec = zeros(MDOF.ND,1); for i = 1:5 MDOF.MVec(i) = 1000E3;end for i = 1:5 MDOF.KVec(i) = 1.0*200

23、0E6;end%MDOF.MMatrix = zeros(MDOF.ND, MDOF.ND);MDOF.KMatrix = zeros(MDOF.ND, MDOF.ND);%for i = 1:MDOF.ND MDOF.MMatrix(i,i) = MDOF.MVec(i);end%for i = 1:MDOF.ND if i = 1 MDOF.KMatrix(i,i ) = MDOF.KVec(i ) + MDOF.KVec(i+1); MDOF.KMatrix(i,i+1) = -MDOF.KVec(i+1); else if i = MDOF.ND MDOF.KMatrix(i,i-1)

24、 = -MDOF.KVec(i); MDOF.KMatrix(i,i ) = MDOF.KVec(i); else MDOF.KMatrix(i,i-1) = -MDOF.KVec(i ); MDOF.KMatrix(i,i ) = MDOF.KVec(i ) + MDOF.KVec(i+1); MDOF.KMatrix(i,i+1) = -MDOF.KVec(i+1); end endend%MDOF.Eigen_Vec, MDOF.Eigen_Val = eig(MDOF.MMatrixMDOF.KMatrix);%MDOF.GMMatrix= transpose(MDOF.Eigen_V

25、ec)*MDOF.MMatrix*MDOF.Eigen_Vec;MDOF.GKMatrix= transpose(MDOF.Eigen_Vec)*MDOF.KMatrix*MDOF.Eigen_Vec;MDOF.GEVector= transpose(MDOF.Eigen_Vec)*MDOF.MMatrix*ones(MDOF.ND,1);%global SDOF%Result.Time=EWave.Time;Result.Disp=zeros(length(Result.Time),MDOF.ND);Result.Vel=zeros(length(Result.Time),MDOF.ND);

26、for i=1:MDOF.NDSDOF.T=2*pirt(MDOF.Eigen_Val(i,i);SDOF.W=2*pi/SDOF.T;SDOF.D=0.00;%SDOF.AMatrix=zeros(2,2);%SDOF.AMatrix(1,2)=1;SDOF.AMatrix(2,1)=-2*SDOF.D*SDOF.W;SDOF.AMatrix(2,2)=-SDOF.W2;%SDOF.BVector=zeros(2,1);%SDOF.BVector(2)=-MDOF.GEVector(i)/MDOF.GMMatrix(i,i);% Execute the time-history analys

27、isSolver.TD=min(EWave.Time) max(EWave.Time);Solver.IC=zeros(2,1);Solver.Opt=odeset(MaxStep,2.0*EWave.DT);%T,V=ode45(SDOF_Time_History_Analysis_ODE,Solver.TD,Solver.IC,Solver.Opt);%for j=1:length(EWave.Time) Result.Disp(j,i)=interp1(T,V(:,1),EWave.Time(j); Result.Vel(j,i)=interp1(T,V(:,2),EWave.Time(

28、j);end%disp(i);% subplot(2,1,1)% plot(T, V(:,1),blue,LineWidth,2); grid on; hold on;% subplot(2,1,2)% plot(T, V(:,2),red,LineWidth,2); grid on; hold on;end% Result.Disp_ALL=transpose(MDOF.Eigen_Vec*transpose( Result.Disp); Result.Disp_1st=transpose(MDOF.Eigen_Vec(:,1)*transpose( Result.Disp(:,1); Re

29、sult.Disp_2nd=transpose(MDOF.Eigen_Vec(:,1:2)*transpose( Result.Disp(:,1:2); Result.Disp_3rd=transpose(MDOF.Eigen_Vec(:,1:3)*transpose( Result.Disp(:,1:3); % plot(Result.Time, Result.Disp_ALL(:,MDOF.ND),red); hold on; plot(Result.Time, Result.Disp_1st(:,MDOF.ND),blue); hold on;grid on;3、瑞利法程序ruilifa

30、%clc;%clc;clear;close;%global MDOF%MDOF.ND = 5;%MDOF.MVec = zeros(MDOF.ND,1);MDOF.KVec = zeros(MDOF.ND,1); for i = 1:5 MDOF.MVec(i) = 1000E3;end for i = 1:5 MDOF.KVec(i) = 2000E6;end%MDOF.X = zeros(MDOF.ND,1);MDOF.X1 = zeros(MDOF.ND,1);for i=1:MDOF.ND s=0; for n=i:MDOF.ND s=s+MDOF.MVec(n); end MDOF.

31、X1(i)=s*9.8/(MDOF.KVec(i);endfor i=1:MDOF.ND if i=1 MDOF.X(i)=MDOF.X1(i); else MDOF.X(i)=MDOF.X(i-1)+MDOF.X1(i); endendMDOF.WVec=sqrt(9.8*transpose(MDOF.X)*MDOF.MVec)/(transpose(MDOF.MVec)*MDOF.X.2);4、矩陣迭代法Diedaifa%clc;clear;close;%global MDOF%MDOF.ND = 5;%MDOF.MVec = zeros(MDOF.ND,1);MDOF.KVec = ze

32、ros(MDOF.ND,1); for i = 1:5 MDOF.MVec(i) = 1000E3;end for i = 1:5 MDOF.KVec(i) = 1.0*2000E6;end%MDOF.MMatrix = zeros(MDOF.ND, MDOF.ND);MDOF.KMatrix = zeros(MDOF.ND, MDOF.ND);%for i = 1:MDOF.ND MDOF.MMatrix(i,i) = MDOF.MVec(i);end%for i = 1:MDOF.ND if i = 1 MDOF.KMatrix(i,i ) = MDOF.KVec(i ) + MDOF.K

33、Vec(i+1); MDOF.KMatrix(i,i+1) = -MDOF.KVec(i+1); else if i = MDOF.ND MDOF.KMatrix(i,i-1) = -MDOF.KVec(i); MDOF.KMatrix(i,i ) = MDOF.KVec(i); else MDOF.KMatrix(i,i-1) = -MDOF.KVec(i ); MDOF.KMatrix(i,i ) = MDOF.KVec(i ) + MDOF.KVec(i+1); MDOF.KMatrix(i,i+1) = -MDOF.KVec(i+1); end endendMDOF.FMatrix =

34、 zeros(MDOF.ND, MDOF.ND);MDOF.FMatrix = inv(MDOF.KMatrix);MDOF.WMatrix = zeros(MDOF.ND,1);MDOF.DMatrix = zeros(MDOF.ND,MDOF.ND);MDOF.DMatrix = MDOF.FMatrix*MDOF.MMatrix;MDOF.YMatrix = zeros(MDOF.ND,MDOF.ND);MDOF.SMatrix = ones(MDOF.ND,MDOF.ND); MDOF.WMatrix(1,1) = sqrt(MDOF.DMatrix*MDOF.SMatrix(:,1)

35、MDOF.SMatrix(:,1); MDOF.YMatrix(:,1) = MDOF.WMatrix(1,1)2*MDOF.DMatrix*MDOF.SMatrix(:,1); while(max(abs(MDOF.YMatrix(:,1)-MDOF.SMatrix(:,1)0.0001) MDOF.SMatrix(:,1) = MDOF.YMatrix(:,1); MDOF.WMatrix(1,1) = sqrt(MDOF.DMatrix*MDOF.SMatrix(:,1) MDOF.SMatrix(:,1); MDOF.YMatrix(:,1) = MDOF.WMatrix(1,1)2*

36、MDOF.DMatrix*MDOF.SMatrix(:,1); MDOF.YMatrix(:,1) = MDOF.YMatrix(:,1)/MDOF.YMatrix(MDOF.ND,1); endfor i=2:5 MDOF.MGMatrix = zeros(MDOF.ND,MDOF.ND); MDOF.MGMatrix(i-1,i-1) = (transpose(MDOF.YMatrix(:,i-1)*MDOF.MMatrix*MDOF.YMatrix(:,i-1); MDOF.DMatrix = MDOF.DMatrix-1/(MDOF.WMatrix(i-1,1)2)*MDOF.YMat

37、rix(:,i-1)*transpose(MDOF.YMatrix(:,i-1)*MDOF.MMatrix/MDOF.MGMatrix(i-1,i-1); MDOF.WMatrix(i,1) = sqrt(MDOF.DMatrix*MDOF.SMatrix(:,i)MDOF.SMatrix(:,i); MDOF.YMatrix(:,i) = MDOF.WMatrix(i,1)2*MDOF.DMatrix*MDOF.SMatrix(:,i); while(max(abs(MDOF.YMatrix(:,i)-MDOF.SMatrix(:,i)0.0001) MDOF.SMatrix(:,i) =

38、MDOF.YMatrix(:,i); MDOF.WMatrix(i,1) = sqrt(MDOF.DMatrix*MDOF.SMatrix(:,i) MDOF.SMatrix(:,i); MDOF.YMatrix(:,i) = MDOF.WMatrix(i,1)2*MDOF.DMatrix*MDOF.SMatrix(:,i); MDOF.YMatrix(:,i) = MDOF.YMatrix(:,i)/MDOF.YMatrix(MDOF.ND,i); endend5、6、求解地震力程序Dizhenli%clc;clear;close;%global MDOF%MDOF.ND = 5;%MDOF

39、.MVec = zeros(MDOF.ND,1);MDOF.KVec = zeros(MDOF.ND,1); for i = 1:5 MDOF.MVec(i) = 1000E3;end for i = 1:5 MDOF.KVec(i) = 1.0*2000E6;end%MDOF.MMatrix = zeros(MDOF.ND, MDOF.ND);MDOF.KMatrix = zeros(MDOF.ND, MDOF.ND);%for i = 1:MDOF.ND MDOF.MMatrix(i,i) = MDOF.MVec(i);end%for i = 1:MDOF.ND if i = 1 MDOF.KMatrix(i,i ) = MDOF.KVec(i ) + MDOF.KVec(i+1); MDOF.KMatrix(i,i+1) = -MDOF.KVec(i+1); else if i = MDOF.ND MDOF.KMatrix(i,i-1) = -MDOF.KVec(i); MDO