TMD系统最优参数的设计方法

SchoolofCiviSchoolofCivi戲土木工程与力学学院NGINEERING&MECHANICS结构动力学小论文班级 土木卓越1201班学号U2O1210323姓名陈祥磊指导老师 叶昆2015.01.05

TMD系统最优参数的设计方法摘要：调谐质量阻尼器TMD由质块，弹簧与阻尼系统组成。即由将其振动频率调整至主结构频率附近，改变结构共振特性，以达到减震作用。将调谐质量阻尼器(TMD)装入结构的目的是减少在外力作用下基本结构构件的消能要求值。在该情况下，这种减小是通过将结构振动的一些能量传递给以最简单的形式固定或连接在主要结构的辅助质量一弹簧一阻尼筒系统构成的TMD来完成的。现在的建筑结构在地震作用下容易产生过大的反应进而发生破坏，因此TMD等减震结构显得非常重要，要将TMD应用于实际结构中，鉴于结构的空间都是有限的，所以TMD不能过大，即TMD的质量相对于结构而言应该很小。本文中选择mTMD二0.05xM，即TMD的质量为主体结构的5%。其次，TMD应该能够发挥明显的减震作用，因此我们需要对TMD的参数进行设计选择。本文对结构基底在受地震激励下的TMD参数设计进行了研究，并且用真实的地震波通过MATLAB编程的方法实现TMD的作用以搜索到最优的TMD参数。关键词：TMD阻尼比频率比参数优化一、TMD减震理论简介下图所示为两自由度体系的结构图，通过这个结构来研究TMD结构的减震机理。列出两个质点的平衡方程如下：(x—x)+k(x—x)=—mx2 2 1 2 2 1 2g■(x—x)+kx—k(x—x)=—mxmx+cxmx+cx—c写成矩阵形式即为：m0__xm0__x"c+c—cxk+k—kx11+1221+12210mc—x—c——c— ccc—x—c——k— ckc—x—c—整个结构的阻尼矩阵：C=^M+PK，要求出«、m1-2」0，通过结构的第一二主振频率求得： a=空#2 0=垄Wi+q ®+®1212由于直接用各质点相对与地面的位移值难以直接反应结构在地震下的层间位移,所以，将位移量进行变换，将各层间位移量作为基本未知量，即令x=卩X—x=卩11212再列出两个质点的平衡方程如下：TOC\o"1-5"\h\z■ ■■ ■■ I ・ -Jf ・・mx|Li +pi +c j-i +k |Li=—mx2 1 2 2 2 2 2 2gm+mmc0■k0■_i_m+m1221+11+11———12mmi0ci0kim1- 22」1-2」—2」1-2」—2」1-2」L2」写成矩阵形式为:对于多自由度的结构而言,此时的质量矩阵、刚度矩阵将会发生改变1121211mpi+m(pi+i)+ci+ki=—Cm+m上写成矩阵形式为:对于多自由度的结构而言,此时的质量矩阵、刚度矩阵将会发生改变1121211_ - 11 1 2g+•••+作m?+•••+^5m?+•••+^5m?+•••+作…咗…m5i2+00 …c2…・・・•0"0「叮卩2+00k2…0"…0i2=—作+•••+^5m^+•••+^5m5m5…m5-_is-00…c5-_卩5-00…k5-_i5--m5 -xgg其中的质量矩阵不再是对角矩阵，而是满秩矩阵，其表达式如下:m“+m•••+mm2+13Nm+m•••+mm+23N2m+•••+m3Nm一N编写程序形成M矩阵时，M矩阵符合下列表达式:mc…3•+mN••••••mNm…3•+mN••••••mN•••••••••mN•••••••••mNmNmN•••mN-mmk；(j?「M=ijM=

匚ij编程时即可形成满秩的质量矩阵。经过变换后，C矩阵是一个对角矩阵，原来的C矩阵为_C1+C2-C2000-C2C2+C3-C3000-C3C3+C40000C4C4+C5-C5000-C5C5经过变换后，C矩阵为_C10000_0C200000C300000C400000C5-阻尼矩阵的变化通过下面的程序实现：fori=l:ND讦i==NDC(i)=CMatrix1(i,i);elseC(i)=CMatrix1(i,i)+CMatrix1(i,i+1);endend刚度矩阵也是相同的变化。经过这样的变化之后便于我们研究结构在地震力作用下的各层层间位移。虽然TMD系统可以用于减震，但是应该选择合适的频率比f和阻尼比t，否则难以发挥则用，甚至起到反作用。对于理想的无阻尼的系统而言，可以求解出f和£的解析式，进而求出最优值。其最优值如下：f「F/(i+卩)；If石实际建筑工程中的系统都是有阻尼的，而且安装了TMD的系统激振频率一旦偏离TMD系统的固有频率，主结构的振幅将急剧增大。所以，研究有阻尼系统的TMD是很有实际意义的，此时求解的出的TMD参数才能真正应用于实际结构中。考虑了主结构的阻尼的TMD系统，将无法导出最优频率比和最优阻尼

比的解析式。虽然可以通过非线性规划的方法得到最优频率比和最优阻尼比的近似解，但是过于繁杂，需要一种程序化的方法来简化求解过程。结构参数计算结构为下图所示的多自由度体系结构，研究此结构在地震动作用下的位移、速度以及加速度等的变化过程。结构计算参数如下:1、2、mi-m2mN二m ；ki二k3、结构参数中N=5；k=1.0。MtmdCtmdK.5M4C4KdM21、2、mi-m2mN二m ；ki二k3、结构参数中N=5；k=1.0。MtmdCtmdK.5M4C4KdM2CaMiCiKi结构计算简图4、其中每一层的阻尼比为E=0.05，TMD的质量为结构总质量的5%，即Mtmd=5%x5x1000xIO3=250x103kg

5、根据上面的结构图和计算参数，求得在无阻尼时TMD的最优阻尼C和刚度匸 i丿严J'd+pr0.1336f= +P/2/（I+卩）匸 i丿严J'd+pr0.1336CTMD二2mTMD^TMD^TMD二=819856= 》pM=376653571取K〜3.767x107；C=8.20x10 6 8 10 12 14 16 18 20 22Driftof4thStory 6 8 10 12 14 16 18 20 22Driftof4thStory在求解有阻尼结构的TMD最优频率比和阻尼比时，以求得的无阻尼情况下的最优阻尼C和刚度K作为初值，来求有阻尼系统的最优值。-ualuaoe-dsn0 5 10 15 20 25 30 35 40Drift-ualuaoe-dsn0 5 10 15 20 25 30 35 40Driftof4thStory放大之后细节图为:x10'3caLuaoe-dsocaLuaoe-dso第四层楼的层间位移图0 5 10 15 20 25 30 35 40Driftof5thStory放大之后细节图为:o-3-ualuaoe-dsQ6 8 10 12 14 16 18 20 22 24Drifto-3-ualuaoe-dsQ6 8 10 12 14 16 18 20 22 24Driftof5thStory-2第五层楼的层间位移图O.-015101520DriftofTMD253035 40TMD的相对位移图O-cO.-015101520DriftofTMD253035 40TMD的相对位移图O-caluaoe-dsQ由程序可得TMD的位移幅值为max(yTMD)=0.062802，其与层间最大位移的比值k=yTMDk=yTMD0.06280.0154=4.1再观察TMD的位移图可知，TMD的位移远大于结构本身各层的层间位移值,从能量角度而言，TMD结构通过吸收地震的能量从而达到减小结构地震反映的目的，所以自身的位移值会很大。实际结构中通常将TMD装置安装在结构的顶部,由于结构的承载能力和空间有限，所以其质量不能过大，一般都是结构总质量的5%以内。结构加上TMD前后的层间位移幅值对比表位移模式位移幅值withoutTMDwithTMD比值Y10.0154087280.0131248130.851777861Y20.0139448910.0118609580.850559356Y30.0115306660.0098419230.853543284Y40.0082223740.0070791670.860963836Y50.0042738240.0037749890.883281401将上述表格的数据制成图，从图上可以清楚地看到加上TMD之后的效果是很明显的。蓝色线条表示的是未加TMD的层间位移，红色线为加了TMD后的层间位移。由图表可知，对于结构的层间位移，加了TMD时的位移幅值大概都是未加TMD时位移幅值的85%,可见隔震效果比较明显，能够明显的减小各层的位移，对于提高结构的抗震性能非常有用。此时的频率比和阻尼比还不是待求的最优值，表明在最优条件下，TMD的作用会更加明显。我们还可以从表中看到,

无论是否加了TMD,层间位移都是逐层递减的，底层的层间位移（也就是第一层的位移）最大，顶层的层间位移最小。分析时可以选择任意一个作为我们分析的目标值。层间位移幅值表1 2 3 4 5WithoutTMD1 2 3 4 5WithoutTMD-TMD三、优化算法1＞采用Matlab中的优化函数将结构本身看做一个质点，考虑TMD对这个系统的减震作用，结构分析如下：加TMD后结构分析图上图是主结构一-TMD系统模型。设质量为ms、mT,弹簧的刚度系数为K、ST sKt，阻尼器的粘性阻尼系数为Cs、CT。设在地面的运动加速度xg(t)作用下,在时刻t，主结构和TMD系统相对于基底(x-x轴)的位移分别为xs(t)、XT(t),而其相对加速度分别为xs(t)、xg(t)，绝对加速度分别为xs(t)+xg(t)，根据sgsg达朗贝尔原理可以建立“主结构一TMD”系统的运动微分方程：m0x(t)[C+C-C一「x(t)]k+k-kx(t)[m0一x(t)]S0mxS(t)'-'tj+s T-CTTCT」xS(t)l**t 」+s T-k匚 TTkT」xS(t)l**t js0mT_lxS(t)l**t 」将结构的动力反应作为函数Y(K,C)；而Y(K,C)实际上是一个关于K、C的隐函数，所以将求解结构层间位移的程序作为目标函数，其中的K、C值用待求的未知量x(K,C)表示。通过程序不断地搜索任意的x0(K,C)值带入直到搜索到最小的位移，此时的目标函数可以选择是任意一层的层间位移，也可以选择层间位移之和，也就是最上面一层的对地位移，本文选择后面一种，即以最上一层的位移值作为目标函数，通过函数优化找出其最小值时TMD结构的K、C值。编写程序运行结果如下：CommandWindow1.0e+07*3.766525299961552 0.082213714221286设定的参数：x=1.0e+07*(3.766525299961552 0.082213714221286)也就是K=37665253；C=822137。也就是此种方法下求得的最优值。将此时的值带入到TMD程序中，算的各层的层间位移列表如下：

结构加上TMD时优化前后的层间位移幅值对比表位移模式位移幅值withoutTMDwithTMDwithTMD（yh）比值Y10.0154087280.0131248130.0083998320.545134673Y20.0139448910.0118609580.0071615690.51356224Y30.0115306660.0098419230.0053447690.463526501Y40.0082223740.0070791670.0037532540.456468411Y50.0042738240.0037749890.0027321370.639272263泡 层间位移幅值表从图表上可以清楚地看到加上参数最优的TMD之后的效果是很明显的。黄色线条表示的是未加TMD的层间位移。由图表可知，对于结构的层间位移，加了优化的TMD时的位移幅值大概都是未加TMD时位移幅值的50%左右，对于此时的顶层位移比上原顶层位移为51.313814%，可见隔震效果比较明显，能够极大的减小各层的位移，对于提高结构的抗震性能非常有用。再看优化效果，对于加了TMD的第一种情况，通过参数优化使得各层间位移明显减小。减少程度均为30%多，可得优化是很明显的。2、采用遗传算法求解最优值由于Y(K,C)是一个关于TMD结构K、C的隐函数，无法采用直接求导方法求其极值。此外，该优化问题为不等式约束优化问题。对不等式约束优化问题，传统的方法采用罚函数法，此法常因边界条件及惩罚因子的设置不当，而无法收敛。另外一个解决算法就是遗传算法，遗传算法是基于生物遗传和进化机制、适合于复杂系统的自适应概率优化技术。通过遗传算法可以求得最优的K、C。由于目前所掌握的Matlab知识有限，而且此方法特别复杂。所以还不能将此方法进行实际的编程运用，只是大概了解了一下思想。【附录】1、Tuned_Mass_Damper_2DOFclc;clear;close;%%globalEWaveFile_Name ='E:\Matlab_Code\Matlab_Code\Ground_Motions_Library\IMPVALL\H-E01140.AT2';%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);%DATA=transpose(fscanf(fid,'%g%g',[1inf]));%EWave.Time=zeros(EWave.NPTs,1);fori=1:EWave.NPTsEWave.Time(i)=(i-1)*EWave.DT;endEWave.Acel=DATA(:,1);EWave.AMax=max(abs(EWave.Acel));%EWave.Acel=EWave.Acel*0.3*9.81/EWave.AMax;%EWave.DT=EWave.Time(2)-EWave.Time(1);%%Parameters1——WithoutTMDglobalEVectorlglobalMMatrixlglobalCMatrixlglobalKMatrixlglobalMVecglobalKVec%ND=5;%MVec=zeros(ND,1);KVec=zeros(ND,1);fori=1:NDMVec(i)=1000E3;endfori=1:NDKVec(i)=1.0*2000E6;end%%MMatrix1=zeros(ND,ND);KMatrix1=zeros(ND,ND);%fori=1:NDMMatrix1(i,i)=MVec(i);end%fori=1:ND讦i==1KMatrix1(i,i)=KVec(i)+KVec(i+1);KMatrix1(i,i+1)=-KVec(i+1);else讦i==NDKMatrix1(i,i-1)=-KVec(i);KMatrix1(i,i)=KVec(i);elseKMatrix1(i,i-1)=-KVec(i);KMatrix1(i,i)=KVec(i)+KVec(i+1);KMatrix1(i,i+1)=-KVec(i+1);endendendEVector1=zeros(ND,1);%[EIGVecs,EIGVals]=eig(MMatrix1\KMatrix1);%WVector=sqrt(EIGVals);TVector=2*pi./sqrt(diag(EIGVals));%BETA=2.0*0.05/(WVector(l,l)+WVector(2,2));ALPHA=BETA*WVector(1,1)*WVector(2,2);%CMatrix1=ALPHA*MMatrix1+BETA*KMatrix1;%fori=1:ND讦i==NDC(i)=CMatrix1(i,i);elseC(i)=CMatrix1(i,i)+CMatrix1(i,i+1);endend%CMatrix1=zeros(ND,ND);MMatrix1=zeros(ND,ND);fori=1:NDforj=1:ND讦j>=iforn=j:NDMMatrix1(i,j)=MMatrix1(i,j)+MVec(n);endelseforn=i:NDMMatrix1(i,j)=MMatrix1(i,j)+MVec(n);endendendend%KMatrix1=zeros(ND,ND);fori=1:NDKMatrix1(i,i)=KVec(i);end%fori=1:NDCMatrix1(i,i)=C(i);end%fori=1:NDEVector1(i)=MMatrix1(i,i);end%%Parameters2TMD——WithTMDglobalEVector2globalMMatrix2globalCMatrix2globalKMatrix2globalMVec2%EVector2=zeros(ND+l,l);MMatrix2=zeros(ND+1,ND+1);CMatrix2=zeros(ND+1,ND+1);KMatrix2=zeros(ND+1,ND+1);%fori=1:ND+1讦i==ND+1MVec2(i)=2.5E5;elseMVec2(i)=MVec(i);endendfori=1:ND+1forj=1:ND+1讦j>=iforn=j:ND+1MMatrix2(i,j)=MMatrix2(i,j)+MVec2(n);endelseforn=i:ND+1MMatrix2(i,j)=MMatrix2(i,j)+MVec2(n);endendendend%fori=1:ND+1讦i==ND+1KMatrix2(i,i)=3.766E7;elseKMatrix2(i,i)=KVec(i);endend%fori=1:ND+1讦i==ND+1CMatrix2(i,i)=819856;

elseCMatrix2(i,i)=CMatrixl(i,i);endend%fori=1:ND+1EVector2(i)=MMatrix2(i,i);end%%Executethetime-historyanalysis01%Solver1.TD=[min(EWave.Time)max(EWave.Time)];Solver1.IC =zeros(10,1);Solver1.Opt =odeset('MaxStep',EWave.DT);%[T1,V1]=ode45(@Tuned_Mass_Damper_2DOF_ODE1,Solver1.TD,Solver1.Opt);%subplot(3,1,1)plot(T1,V1(:,4),'red');gridon;holdon;xlabel('Driftof4thStory');ylabel('Displacement');%subplot(3,1,2)plot(T1,V1(:,5),'blue');gridon;holdon;xlabel('Driftof5thStory');ylabel('Displacement');%%subplot(3,1,3)%plot(T1,V1(:,6),'blue');gridon;holdon;%xlabel('Driftof6thStory');%ylabel('Displacement');%%Executethetime-historyanalysis02%Solver2.TD=[min(EWave.Time)max(EWave.Time)];Solver2.IC =zeros(12,1);Solver2.Opt =odeset('MaxStep',EWave.DT);%[T2,V2]=ode45(@Tuned_Mass_Damper_2DOF_ODE2,Solver2.TD,Solver2.Opt);%subplot(3,1,1)plot(T2,V2(:,4),'blue','LineWidth',2);gridon;holdon;hleg1=legend('WithoutTMD','WithTMD');Solver1.IC,Solver2.IC,%Solver1.IC,Solver2.IC,subplot(3,l,2)plot(T2,V2(:,5),'yellow','LineWidth',2);gridon;holdon;hleg2=legend('WithoutTMD','WithTMD');%subplot(3,1,3)plot(T2,V2(:,6),'blue','LineWidth',2);gridon;holdon;xlabel('DriftofTMD');ylabel('Displacement');2、 Youhua.m优化程序options=optimset(options,'tolfun',1e-10);x0=[37665357819856];[x,fval]=fminunc(@fun1,x0,options);xfval3、 functionf=fun1(x)%%globalEVector2globalMMatrix2globalCMatrix2globalKMatrix2globalCMatrix1globalMVec2globalNDglobalMVecglobalKVecglobalEWave%%EVector2=ze