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

1、 结构动力学小论文 班 级 土木卓越1201班 学 号 u201210323 姓 名 陈祥磊 指导老师 叶 昆 2015.01.05tmd系统最优参数的设计方法摘要：调谐质量阻尼器tmd由质块，弹簧与阻尼系统组成。即由将其振动频率调整至主结构频率附近，改变结构共振特性，以达到减震作用。将调谐质量阻尼器(tmd)装入结构的目的是减少在外力作用下基本结构构件的消能要求值。在该情况下，这种减小是通过将结构振动的一些能量传递给以最简单的形式固定或连接在主要结构的辅助质量弹簧阻尼筒系统构成的tmd来完成的。现在的建筑结构在地震作用下容易产生过大的反应进而发生破坏，因此tmd等减震结构显得非常重要，要将t

2、md应用于实际结构中，鉴于结构的空间都是有限的，所以tmd不能过大，即tmd的质量相对于结构而言应该很小。本文中选择，即tmd的质量为主体结构的5%。其次，tmd应该能够发挥明显的减震作用，因此我们需要对tmd的参数进行设计选择。本文对结构基底在受地震激励下的tmd参数设计进行了研究，并且用真实的地震波通过matlab编程的方法实现tmd的作用以搜索到最优的tmd参数。关键词：tmd 阻尼比 频率比 参数优化1、 tmd减震理论简介 下图所示为两自由度体系的结构图，通过这个结构来研究tmd结构的减震机理。列出两个质点的平衡方程如下： 写成矩阵形式即为：整个结构的阻尼矩阵：，要求出、，通过结构的

3、第一二主振频率求得： 由于直接用各质点相对与地面的位移值难以直接反应结构在地震下的层间位移，所以，将位移量进行变换，将各层间位移量作为基本未知量，即令 再列出两个质点的平衡方程如下： 写成矩阵形式为：对于多自由度的结构而言，此时的质量矩阵、刚度矩阵将会发生改变其中的质量矩阵不再是对角矩阵，而是满秩矩阵，其表达式如下：编写程序形成m矩阵时，m矩阵符合下列表达式： 编程时即可形成满秩的质量矩阵。经过变换后，矩阵是一个对角矩阵，原来的矩阵为经过变换后，矩阵为阻尼矩阵的变化通过下面的程序实现：for i=1:nd if i=nd c(i)=cmatrix1(i,i); else c(i)=cmatri

4、x1(i,i)+cmatrix1(i,i+1); endend刚度矩阵也是相同的变化。经过这样的变化之后便于我们研究结构在地震力作用下的各层层间位移。虽然tmd系统可以用于减震，但是应该选择合适的频率比和阻尼比，否则难以发挥则用，甚至起到反作用。对于理想的无阻尼的系统而言，可以求解出和的解析式，进而求出最优值。其最优值如下： ； 实际建筑工程中的系统都是有阻尼的，而且安装了tmd的系统激振频率一旦偏离tmd系统的固有频率，主结构的振幅将急剧增大。所以，研究有阻尼系统的tmd是很有实际意义的，此时求解的出的tmd参数才能真正应用于实际结构中。考虑了主结构的阻尼的tmd系统，将无法导出最优频率比和

5、最优阻尼比的解析式。虽然可以通过非线性规划的方法得到最优频率比和最优阻尼比的近似解，但是过于繁杂，需要一种程序化的方法来简化求解过程。二、结构参数 计算结构为下图所示的多自由度体系结构，研究此结构在地震动作用下的位移、速度以及加速度等的变化过程。结构计算参数如下：1、 ；2、 ；3、结构参数中；。 结构计算简图 4、其中每一层的阻尼比为=0.05，tmd的质量为结构总质量的5%，即 5、根据上面的结构图和计算参数，求得在无阻尼时tmd的最优阻尼和刚度： =0.9642 ； =0.1336 =819856 =37665357取；。在求解有阻尼结构的tmd最优频率比和阻尼比时，以求得的无阻尼情况下

6、的最优阻尼和刚度作为初值，来求有阻尼系统的最优值。放大之后细节图为：第四层楼的层间位移图放大之后细节图为：第五层楼的层间位移图tmd的相对位移图由程序可得tmd的位移幅值为0.062802，其与层间最大位移的比值再观察tmd的位移图可知，tmd的位移远大于结构本身各层的层间位移值，从能量角度而言，tmd结构通过吸收地震的能量从而达到减小结构地震反映的目的，所以自身的位移值会很大。实际结构中通常将tmd装置安装在结构的顶部，由于结构的承载能力和空间有限，所以其质量不能过大，一般都是结构总质量的5%以内。 结构加上tmd前后的层间位移幅值对比表位移幅值 位移模式without tmdwith tm

7、d比值y10.0154087280.0131248130.851777861y20.0139448910.0118609580.850559356y30.0115306660.0098419230.853543284y40.0082223740.0070791670.860963836y50.0042738240.0037749890.883281401将上述表格的数据制成图，从图上可以清楚地看到加上tmd之后的效果是很明显的。蓝色线条表示的是未加tmd的层间位移，红色线为加了tmd后的层间位移。 由图表可知，对于结构的层间位移，加了tmd时的位移幅值大概都是未加tmd时位移幅值的85%，可见

8、隔震效果比较明显，能够明显的减小各层的位移，对于提高结构的抗震性能非常有用。此时的频率比和阻尼比还不是待求的最优值，表明在最优条件下，tmd的作用会更加明显。我们还可以从表中看到，无论是否加了tmd，层间位移都是逐层递减的，底层的层间位移（也就是第一层的位移）最大，顶层的层间位移最小。分析时可以选择任意一个作为我们分析的目标值。 3、 优化算法 1、采用matlab中的优化函数 将结构本身看做一个质点，考虑tmd对这个系统的减震作用，结构分析如下： 加tmd后结构分析图上图是主结构-tmd系统模型。设质量为、，弹簧的刚度系数为、，阻尼器的粘性阻尼系数为、。设在地面的运动加速度作用下，在时刻t，

9、主结构和tmd系统相对于基底（轴）的位移分别为、，而其相对加速度分别为、，绝对加速度分别为，根据达朗贝尔原理可以建立“主结构tmd”系统的运动微分方程： 将结构的动力反应作为函数；而 实际上是一个关于、的隐函数，所以将求解结构层间位移的程序作为目标函数，其中的、值用待求的未知量表示。通过程序不断地搜索任意的值带入直到搜索到最小的位移，此时的目标函数可以选择是任意一层的层间位移，也可以选择层间位移之和，也就是最上面一层的对地位移，本文选择后面一种，即以最上一层的位移值作为目标函数，通过函数优化找出其最小值时tmd结构的、值。编写程序运行结果如下：设定的参数：x = 1.0e+07 *(3.766

10、525299961552 0.082213714221286)也就是；。也就是此种方法下求得的最优值。将此时的值带入到tmd程序中，算的各层的层间位移列表如下：结构加上tmd时优化前后的层间位移幅值对比表位移幅值位移模式without tmdwith tmdwith tmd（yh)比值y10.0154087280.0131248130.0083998320.545134673y20.0139448910.0118609580.0071615690.51356224y30.0115306660.0098419230.0053447690.463526501y40.0082223740.00707

11、91670.0037532540.456468411y50.0042738240.0037749890.0027321370.639272263 从图表上可以清楚地看到加上参数最优的tmd之后的效果是很明显的。黄色线条表示的是未加tmd的层间位移。由图表可知，对于结构的层间位移，加了优化的tmd时的位移幅值大概都是未加tmd时位移幅值的50%左右，对于此时的顶层位移比上原顶层位移为51.313814%，可见隔震效果比较明显，能够极大的减小各层的位移，对于提高结构的抗震性能非常有用。再看优化效果,对于加了tmd的第一种情况，通过参数优化使得各层间位移明显减小。减少程度均为30%多，可得优化是很明

12、显的。 2、采用遗传算法求解最优值由于 是一个关于tmd结构、的隐函数，无法采用直接求导方法求其极值。此外，该优化问题为不等式约束优化问题。对不等式约束优化问题，传统的方法采用罚函数法，此法常因边界条件及惩罚因子的设置不当，而无法收敛。另外一个解决算法就是遗传算法，遗传算法是基于生物遗传和进化机制、适合于复杂系统的自适应概率优化技术。通过遗传算法可以求得最优的、。由于目前所掌握的matlab知识有限，而且此方法特别复杂。所以还不能将此方法进行实际的编程运用，只是大概了解了一下思想。【附录】1、tuned_mass_damper_2dofclc;clear;close;%global ewave

13、file_name = e:matlab_codematlab_codeground_motions_libraryimpvallh-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

14、, 1 inf); %ewave.time=zeros(ewave.npts,1); for i=1:ewave.npts ewave.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);% parameters 1 - without tmdglobal evector1global mmatrix1global

15、cmatrix1global kmatrix1global mvecglobal kvec%nd = 5;%mvec = zeros(nd,1);kvec = zeros(nd,1);for i = 1:nd mvec(i) = 1000e3;endfor i = 1:nd kvec(i) = 1.0*2000e6;end%mmatrix1 = zeros(nd, nd);kmatrix1 = zeros(nd, nd);%for i = 1:nd mmatrix1(i,i) = mvec(i);end%for i = 1:nd if i = 1 kmatrix1(i,i ) = kvec(i

16、 ) + kvec(i+1); kmatrix1(i,i+1) = -kvec(i+1); else if i = nd kmatrix1(i,i-1) = -kvec(i); kmatrix1(i,i ) = kvec(i); else kmatrix1(i,i-1) = -kvec(i ); kmatrix1(i,i ) = kvec(i ) + kvec(i+1); kmatrix1(i,i+1) = -kvec(i+1); end endendevector1 = zeros(nd,1);%eigvecs, eigvals = eig(mmatrix1kmatrix1);%wvecto

17、r = sqrt(eigvals);tvector = 2*pi./sqrt(diag(eigvals);%beta = 2.0*0.05/(wvector(1,1) + wvector(2,2);alpha = beta*wvector(1,1)*wvector(2,2);%cmatrix1 = alpha*mmatrix1 + beta*kmatrix1;%for i=1:nd if i=nd c(i)=cmatrix1(i,i); else c(i)=cmatrix1(i,i)+cmatrix1(i,i+1); endend%cmatrix1 = zeros(nd,nd);mmatrix

18、1 = zeros(nd,nd);for i = 1:nd for j=1:nd if j=i for n=j:nd mmatrix1(i,j)=mmatrix1(i,j)+mvec(n); end else for n=i:nd mmatrix1(i,j)=mmatrix1(i,j)+mvec(n); end end endend %kmatrix1 = zeros(nd,nd);for i = 1:nd kmatrix1(i,i) = kvec(i);end%for i = 1:nd cmatrix1(i,i) = c(i);end%for i = 1:nd evector1(i) =mm

19、atrix1(i,i);end% parameters 2 tmd - with tmdglobal evector2global mmatrix2global cmatrix2global kmatrix2global mvec2%evector2 = zeros(nd+1,1);mmatrix2 = zeros(nd+1,nd+1);cmatrix2 = zeros(nd+1,nd+1);kmatrix2 = zeros(nd+1,nd+1);%for i = 1:nd+1 if i=nd+1 mvec2(i)=2.5e5; else mvec2(i) = mvec(i); endendf

20、or i = 1:nd+1 for j=1:nd+1 if j=i for n=j:nd+1 mmatrix2(i,j)=mmatrix2(i,j)+mvec2(n); end else for n=i:nd+1 mmatrix2(i,j)=mmatrix2(i,j)+mvec2(n); end end endend %for i = 1:nd+1 if i=nd+1 kmatrix2(i,i)=3.766e7; else kmatrix2(i,i) = kvec(i); endend%for i = 1:nd+1 if i=nd+1 cmatrix2(i,i)=819856; else cm

21、atrix2(i,i) =cmatrix1(i,i); endend%for i = 1:nd+1 evector2(i) =mmatrix2(i,i);end% execute the time-history analysis 01%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.ic, sol

22、ver1.opt);%subplot(3,1,1)plot(t1, v1(:,4), red ); grid on; hold on;xlabel(drift of 4th story);ylabel(displacement);%subplot(3,1,2)plot(t1, v1(:,5), blue); grid on; hold on;xlabel(drift of 5th story);ylabel(displacement);% subplot(3,1,3)% plot(t1, v1(:,6), blue); grid on; hold on;% xlabel(drift of 6t

23、h story);% ylabel(displacement);% execute the time-history analysis 02%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.ic, solver2.opt);%subplot(3,1,1)plot(t2, v2(:,4),blue,li

24、newidth,2); grid on; hold on;hleg1 = legend(without tmd,with tmd);%subplot(3,1,2)plot(t2, v2(:,5),yellow,linewidth,2); grid on; hold on;hleg2 = legend(without tmd,with tmd);%subplot(3,1,3)plot(t2, v2(:,6),blue,linewidth,2); grid on; hold on;xlabel(drift of tmd);ylabel(displacement);2、 youhua.m优化程序op

25、tions=optimset(options,tolfun,1e-10);x0=37665357 819856;x,fval=fminunc(fun1,x0,options);xfval3、function f=fun1(x)%global evector2global mmatrix2global cmatrix2global kmatrix2global cmatrix1global mvec2global ndglobal mvecglobal kvecglobal ewave%evector2 = zeros(nd+1,1);mmatrix2 = zeros(nd+1,nd+1);cmatrix2 = zeros(nd+1,nd+1);kmatrix2 = zeros(nd+1,nd+1);%for i = 1:nd+1 if i=nd+1 mvec2(i)=250000; else mvec2(i) = mvec(i); endendfor i = 1:nd+1 for j=1:nd+1 if j=i for n=j:nd+1 mmat