结构动力学大作业

1、 结构动力学大作业 班 级 土木卓越1201班 学 号 U201210323 姓 名 陈祥磊 指导老师 叶 昆 2014.12.30 结构动力学大作业 SDOF体系在任意荷载作用下的动力响应一、结构参数 计算结构为右图所示的1、 2、 3、结构参数中；。2、 确定各阶频率和振型 多自由度体系自由振动时的运动方程为 . 写成矩阵形式即为 假设此方程的解答为，带入到运动方程中得到振动方程 此方程要有非零解必须满足频率方程，可解得各阶主频率再根据 可求出结构的主振型。在主振型中，通常将最后一个位移值设定为1，只要在程序中加入下列语句：MDOF.YMatrix(:,i)=MDOF.YMatrix(:,

2、i)/MDOF.YMatrix(MDOF.ND,i)运行程序之后得到如下结果：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由主振型可以求得广义质量矩阵和广义刚度矩阵，均为对角矩阵。程序如下所示：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法主要用于求的近似解，原理是能量守恒定律：；最终求解的表达式为： 因为采用了近似解答，相当于在精确解条件下添加了约束，提高了结构的整体刚度，所以第一频率会大于精确值。在实际计算时，假设的位移幅值函数决定最终计算的精确度，越接近实际的位移模式，

5、Rayleigh法的计算结果越精确。假设的位移幅值函数必须满足以下几点：必须满足运动边界条件，即几何和自然边界条件；所设位移幅值函数必须与实际振型形状大致相近。通常可取结构在耨个静位移作用下的弹性曲线作为的近似表达式，此式应变能可用相应荷载所做的功来代替。 计算时采用列表计算和编程计算两种方法互为佐证以证明结果的正确性，计算结果如下：1、列表计算法 层数质量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 如上表所示得出的第一频率的近似值为，与第一题中的结果比较，相对误差为，即计算结果很接近精确解。 2、编程法 Rayleigh法计算过程简单，易于通过程序实现。变成后进行计算可得结果为 上面分别通过两种不同的方法求得频率值，现在考察计算结果与实际的吻合程度，一般建筑结构的自振周期 则五层建筑结构的第一主振型

7、的周期大概介于之间，计算结果中的符合一般经验。另外据经验可知一二阶主振频率之间有如下关系： 结算结果中的，接近1/3，进一步说明计算结果符合实际工况。4、 振型分解法 结构有固定的自振频率和周期，取决于结构的参数和形式。结构任意的自由振动模式均可以分解为主振型的组合，类似于高阶方程组的基本解组，将结构的位移通过正则坐标表示为基本振型的组合 选定地震动数据（文件为Ground_MotionsIMPVALLH-E01140.at2）7807个，间隔时间为0.005s,将地震动幅值调整为0.3g，使用振型分解法计算结构的地震响应。地震动数据给定的地震波在各个时刻的加速度值。由给定的数据可以绘出地震动

8、的加速度和速度关于时间的曲线：加速度与时间关系图：速度与时间关系图：运行程序SDOF_Time_History_Analysis之后，得到结构的位移与时间关系图放大之后可得细节图5、 矩阵迭代法 矩阵迭代法是采用逐步逼近的方法确定结构的频率和振型。体系做自由振动时各质点的位移幅值为 这是一个迭代关系式，通过假定最初的带入计算，可以求出，继而不断求得，当与足够接近时即可将作为第一振型。其运算过程如下：首先假定一个振型，通常假定为代入上式右边，进行求解后即可得到和主振型的第一次近似值；然后以第一次近似值代入上式进行计算，得到和主振型的第二次近似值，如此下去，直到最后两次的计算结果接近为止。当一个振

9、型求得以后，利用振型的正交性，求出高阶的频率和振型。 当结构层数较多且计算精度要求较高时，宜用编程法来计算。思路一：由于在求解高阶频率和高阶振型时，例如第K振型时，需要先利用主振型的正交性解出关于的表达式，以带入到中以求出新的阶迭代关系式；同理通过迭代即可求出第K阶主振型。编程计算时，如果可以求出在考虑正交性后的新柔度矩阵，即可按照求解第一振型时的语言计算，此时只要设定一系列的循环语句即可。求解第二振型的程序如下：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);其中二阶时新的柔度矩阵的表达式：MDOF.DMatrix(i,j)= MDOF.DMatrix-MDOF.XMatrix(j,1)*MDOF.DMatrix(i,n)/MDOF.XMatrix(n,1) 这种方法可以求解出第二振型的频率和主振型，但是在更高阶

12、的振型求解中，由于关于的表达式是一个K阶矩阵，求出的表达式很复杂，并且在代入的迭代关系式时，难以求出新柔度矩阵，因此难以应用。思路二：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、 绘制反应谱 研究单自由度体系在地震力作用下的地震反应，绘制出不同周期的结构在地震力作用下的地震反应的最大加速度，即可得到一条地震动的反应谱。单自由度体系的运动方程为 表示为：，由此方程可得惯性力表达式 绝对加速度：；得 若给定地震时地面

14、运动的加速度记录和体系的阻尼比，可计算出单自由度体系的最大绝对加速度与体系自振周期的关系曲线，即为加速度反应谱。7、 求解地震力 在上面求出了地震动作用下结构的主振型以及频率、周期值之后，可以求出振兴参与系数： 再根据自振周期求出，由公式，即可求出N组地震作用，分别计算每组地震力作用下的内力，按计算总效应。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