TMD系統(tǒng)最優(yōu)參數(shù)的設(shè)計(jì)方法

1、 結(jié)構(gòu)動(dòng)力學(xué)小論文 班 級 土木卓越1201班 學(xué) 號 u201210323 姓 名 陳祥磊 指導(dǎo)老師 葉 昆 2015.01.05tmd系統(tǒng)最優(yōu)參數(shù)的設(shè)計(jì)方法摘要：調(diào)諧質(zhì)量阻尼器tmd由質(zhì)塊，彈簧與阻尼系統(tǒng)組成。即由將其振動(dòng)頻率調(diào)整至主結(jié)構(gòu)頻率附近，改變結(jié)構(gòu)共振特性，以達(dá)到減震作用。將調(diào)諧質(zhì)量阻尼器(tmd)裝入結(jié)構(gòu)的目的是減少在外力作用下基本結(jié)構(gòu)構(gòu)件的消能要求值。在該情況下，這種減小是通過將結(jié)構(gòu)振動(dòng)的一些能量傳遞給以最簡單的形式固定或連接在主要結(jié)構(gòu)的輔助質(zhì)量彈簧阻尼筒系統(tǒng)構(gòu)成的tmd來完成的?，F(xiàn)在的建筑結(jié)構(gòu)在地震作用下容易產(chǎn)生過大的反應(yīng)進(jìn)而發(fā)生破壞，因此tmd等減震結(jié)構(gòu)顯得非常重要，要將t

2、md應(yīng)用于實(shí)際結(jié)構(gòu)中，鑒于結(jié)構(gòu)的空間都是有限的，所以tmd不能過大，即tmd的質(zhì)量相對于結(jié)構(gòu)而言應(yīng)該很小。本文中選擇，即tmd的質(zhì)量為主體結(jié)構(gòu)的5%。其次，tmd應(yīng)該能夠發(fā)揮明顯的減震作用，因此我們需要對tmd的參數(shù)進(jìn)行設(shè)計(jì)選擇。本文對結(jié)構(gòu)基底在受地震激勵(lì)下的tmd參數(shù)設(shè)計(jì)進(jìn)行了研究，并且用真實(shí)的地震波通過matlab編程的方法實(shí)現(xiàn)tmd的作用以搜索到最優(yōu)的tmd參數(shù)。關(guān)鍵詞：tmd 阻尼比 頻率比 參數(shù)優(yōu)化1、 tmd減震理論簡介 下圖所示為兩自由度體系的結(jié)構(gòu)圖，通過這個(gè)結(jié)構(gòu)來研究tmd結(jié)構(gòu)的減震機(jī)理。列出兩個(gè)質(zhì)點(diǎn)的平衡方程如下： 寫成矩陣形式即為：整個(gè)結(jié)構(gòu)的阻尼矩陣：，要求出、，通過結(jié)構(gòu)的

3、第一二主振頻率求得： 由于直接用各質(zhì)點(diǎn)相對與地面的位移值難以直接反應(yīng)結(jié)構(gòu)在地震下的層間位移，所以，將位移量進(jìn)行變換，將各層間位移量作為基本未知量，即令 再列出兩個(gè)質(zhì)點(diǎn)的平衡方程如下： 寫成矩陣形式為：對于多自由度的結(jié)構(gòu)而言，此時(shí)的質(zhì)量矩陣、剛度矩陣將會(huì)發(fā)生改變其中的質(zhì)量矩陣不再是對角矩陣，而是滿秩矩陣，其表達(dá)式如下：編寫程序形成m矩陣時(shí)，m矩陣符合下列表達(dá)式： 編程時(shí)即可形成滿秩的質(zhì)量矩陣。經(jīng)過變換后，矩陣是一個(gè)對角矩陣，原來的矩陣為經(jīng)過變換后，矩陣為阻尼矩陣的變化通過下面的程序?qū)崿F(xiàn)：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剛度矩陣也是相同的變化。經(jīng)過這樣的變化之后便于我們研究結(jié)構(gòu)在地震力作用下的各層層間位移。雖然tmd系統(tǒng)可以用于減震，但是應(yīng)該選擇合適的頻率比和阻尼比，否則難以發(fā)揮則用，甚至起到反作用。對于理想的無阻尼的系統(tǒng)而言，可以求解出和的解析式，進(jìn)而求出最優(yōu)值。其最優(yōu)值如下： ； 實(shí)際建筑工程中的系統(tǒng)都是有阻尼的，而且安裝了tmd的系統(tǒng)激振頻率一旦偏離tmd系統(tǒng)的固有頻率，主結(jié)構(gòu)的振幅將急劇增大。所以，研究有阻尼系統(tǒng)的tmd是很有實(shí)際意義的，此時(shí)求解的出的tmd參數(shù)才能真正應(yīng)用于實(shí)際結(jié)構(gòu)中?？紤]了主結(jié)構(gòu)的阻尼的tmd系統(tǒng)，將無法導(dǎo)出最優(yōu)頻率比和

5、最優(yōu)阻尼比的解析式。雖然可以通過非線性規(guī)劃的方法得到最優(yōu)頻率比和最優(yōu)阻尼比的近似解，但是過于繁雜，需要一種程序化的方法來簡化求解過程。二、結(jié)構(gòu)參數(shù) 計(jì)算結(jié)構(gòu)為下圖所示的多自由度體系結(jié)構(gòu)，研究此結(jié)構(gòu)在地震動(dòng)作用下的位移、速度以及加速度等的變化過程。結(jié)構(gòu)計(jì)算參數(shù)如下：1、 ；2、 ；3、結(jié)構(gòu)參數(shù)中；。 結(jié)構(gòu)計(jì)算簡圖 4、其中每一層的阻尼比為=0.05，tmd的質(zhì)量為結(jié)構(gòu)總質(zhì)量的5%，即 5、根據(jù)上面的結(jié)構(gòu)圖和計(jì)算參數(shù)，求得在無阻尼時(shí)tmd的最優(yōu)阻尼和剛度： =0.9642 ； =0.1336 =819856 =37665357??；。在求解有阻尼結(jié)構(gòu)的tmd最優(yōu)頻率比和阻尼比時(shí)，以求得的無阻尼情況下

6、的最優(yōu)阻尼和剛度作為初值，來求有阻尼系統(tǒng)的最優(yōu)值。放大之后細(xì)節(jié)圖為：第四層樓的層間位移圖放大之后細(xì)節(jié)圖為：第五層樓的層間位移圖tmd的相對位移圖由程序可得tmd的位移幅值為0.062802，其與層間最大位移的比值再觀察tmd的位移圖可知，tmd的位移遠(yuǎn)大于結(jié)構(gòu)本身各層的層間位移值，從能量角度而言，tmd結(jié)構(gòu)通過吸收地震的能量從而達(dá)到減小結(jié)構(gòu)地震反映的目的，所以自身的位移值會(huì)很大。實(shí)際結(jié)構(gòu)中通常將tmd裝置安裝在結(jié)構(gòu)的頂部，由于結(jié)構(gòu)的承載能力和空間有限，所以其質(zhì)量不能過大，一般都是結(jié)構(gòu)總質(zhì)量的5%以內(nèi)。 結(jié)構(gòu)加上tmd前后的層間位移幅值對比表位移幅值 位移模式without tmdwith tm

7、d比值y10.0154087280.0131248130.851777861y20.0139448910.0118609580.850559356y30.0115306660.0098419230.853543284y40.0082223740.0070791670.860963836y50.0042738240.0037749890.883281401將上述表格的數(shù)據(jù)制成圖，從圖上可以清楚地看到加上tmd之后的效果是很明顯的。藍(lán)色線條表示的是未加tmd的層間位移，紅色線為加了tmd后的層間位移。 由圖表可知，對于結(jié)構(gòu)的層間位移，加了tmd時(shí)的位移幅值大概都是未加tmd時(shí)位移幅值的85%，可見

8、隔震效果比較明顯，能夠明顯的減小各層的位移，對于提高結(jié)構(gòu)的抗震性能非常有用。此時(shí)的頻率比和阻尼比還不是待求的最優(yōu)值，表明在最優(yōu)條件下，tmd的作用會(huì)更加明顯。我們還可以從表中看到，無論是否加了tmd，層間位移都是逐層遞減的，底層的層間位移（也就是第一層的位移）最大，頂層的層間位移最小。分析時(shí)可以選擇任意一個(gè)作為我們分析的目標(biāo)值。 3、 優(yōu)化算法 1、采用matlab中的優(yōu)化函數(shù) 將結(jié)構(gòu)本身看做一個(gè)質(zhì)點(diǎn)，考慮tmd對這個(gè)系統(tǒng)的減震作用，結(jié)構(gòu)分析如下： 加tmd后結(jié)構(gòu)分析圖上圖是主結(jié)構(gòu)-tmd系統(tǒng)模型。設(shè)質(zhì)量為、，彈簧的剛度系數(shù)為、，阻尼器的粘性阻尼系數(shù)為、。設(shè)在地面的運(yùn)動(dòng)加速度作用下，在時(shí)刻t，

9、主結(jié)構(gòu)和tmd系統(tǒng)相對于基底（軸）的位移分別為、，而其相對加速度分別為、，絕對加速度分別為，根據(jù)達(dá)朗貝爾原理可以建立“主結(jié)構(gòu)tmd”系統(tǒng)的運(yùn)動(dòng)微分方程： 將結(jié)構(gòu)的動(dòng)力反應(yīng)作為函數(shù)；而 實(shí)際上是一個(gè)關(guān)于、的隱函數(shù)，所以將求解結(jié)構(gòu)層間位移的程序作為目標(biāo)函數(shù)，其中的、值用待求的未知量表示。通過程序不斷地搜索任意的值帶入直到搜索到最小的位移，此時(shí)的目標(biāo)函數(shù)可以選擇是任意一層的層間位移，也可以選擇層間位移之和，也就是最上面一層的對地位移，本文選擇后面一種，即以最上一層的位移值作為目標(biāo)函數(shù)，通過函數(shù)優(yōu)化找出其最小值時(shí)tmd結(jié)構(gòu)的、值。編寫程序運(yùn)行結(jié)果如下：設(shè)定的參數(shù)：x = 1.0e+07 *(3.766

10、525299961552 0.082213714221286)也就是；。也就是此種方法下求得的最優(yōu)值。將此時(shí)的值帶入到tmd程序中，算的各層的層間位移列表如下：結(jié)構(gòu)加上tmd時(shí)優(yōu)化前后的層間位移幅值對比表位移幅值位移模式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 從圖表上可以清楚地看到加上參數(shù)最優(yōu)的tmd之后的效果是很明顯的。黃色線條表示的是未加tmd的層間位移。由圖表可知，對于結(jié)構(gòu)的層間位移，加了優(yōu)化的tmd時(shí)的位移幅值大概都是未加tmd時(shí)位移幅值的50%左右，對于此時(shí)的頂層位移比上原頂層位移為51.313814%，可見隔震效果比較明顯，能夠極大的減小各層的位移，對于提高結(jié)構(gòu)的抗震性能非常有用。再看優(yōu)化效果,對于加了tmd的第一種情況，通過參數(shù)優(yōu)化使得各層間位移明顯減小。減少程度均為30%多，可得優(yōu)化是很明

12、顯的。 2、采用遺傳算法求解最優(yōu)值由于 是一個(gè)關(guān)于tmd結(jié)構(gòu)、的隱函數(shù)，無法采用直接求導(dǎo)方法求其極值。此外，該優(yōu)化問題為不等式約束優(yōu)化問題。對不等式約束優(yōu)化問題，傳統(tǒng)的方法采用罰函數(shù)法，此法常因邊界條件及懲罰因子的設(shè)置不當(dāng)，而無法收斂。另外一個(gè)解決算法就是遺傳算法，遺傳算法是基于生物遺傳和進(jìn)化機(jī)制、適合于復(fù)雜系統(tǒng)的自適應(yīng)概率優(yōu)化技術(shù)。通過遺傳算法可以求得最優(yōu)的、。由于目前所掌握的matlab知識有限，而且此方法特別復(fù)雜。所以還不能將此方法進(jìn)行實(shí)際的編程運(yùn)用，只是大概了解了一下思想。【附錄】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優(yōu)化程序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