TMD小论文-cxl解读

结构动力学小论文
班级土木卓越1201班
学号 U201210323
姓名陈祥磊
指导老师叶昆
2015.01.05
TMD 系统最优参数的设计方法
摘要：
调谐质量阻尼器TMD 由质块，弹簧与阻尼系统组成。

即由将其振动频率调整至主结构频率附近，改变结构共振特性，以达到减震作用。

将调谐质量阻尼器(TMD)装入结构的目的是减少在外力作用下基本结构构件的消能要求值。

在该情况下，这种减小是通过将结构振动的一些能量传递给以最简单的形式固定或连接在主要结构的辅助质量—弹簧—阻尼筒系统构成的TMD 来完成的。

现在的建筑结构在地震作用下容易产生过大的反应进而发生破坏，因此TMD 等减震结构显得非常重要，要将TMD 应用于实际结构中，鉴于结构的空间都是有限的，所以TMD 不能过大，即TMD 的质量相对于结构而言应该很小。

本文中选择M m TMD ⨯=05.0，即TMD 的质量为主体结构的5%。

其次，TMD 应该能够发挥明显的减震作用，因此我们需要对TMD 的参数进行设计选择。

本文对结构基底在受地震激励下的TMD 参数设计进行了研究，并且用真实的地震波通过MATLAB 编程的方法实现TMD 的作用以搜索到最优的TMD 参数。

关键词：TMD 阻尼比 频率比 参数优化 一、TMD 减震理论简介
下图所示为两自由度体系的结构图，通过这个结构来研究TMD 结构的减震机理。

列出两个质点的平衡方程如下：
()()g x
m x x k x x c x 212212222m -=-+-+ ()()g x m x x k x k x x c x c x 1122111221111m -=--+--+ 写成矩阵形式即为：
g x m m x x k k k k k x x
c c c c c x x m m ⎥⎦
⎤
⎢⎣⎡-=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡--++⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡--++⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡2121222212122221212100 整个结构的阻尼矩阵：K M C βα+=，要求出α、β，通过结构的第一二主振频率求得： 2
1212ωωωξωα+=
212ωωξ
β+=
由于直接用各质点相对与地面的位移值难以直接反应结构在地震下的层间位移，
所以，将位移量进行变换，将各层间位移量作为基本未知量，即令
11μ=x 212x μ=-x
再列出两个质点的平衡方程如下：
()g x
m k c 22222212m -=+++μμμμ ()()g x m m k c 21111121211m m +-=++++μμμμμ 写成矩阵形式为：
g x m m m k k c c m m m m m ⎥⎦⎤⎢⎣⎡+-=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣
⎡+2212121212121222210000μμμμμμ
对于多自由度的结构而言，此时的质量矩阵、刚度矩阵将会发生改变
g x
m m m m m k k k c c c m m m m m
m m m m m m m m
⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡++++-=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢
⎢⎢⎢⎣
⎡++++++++5525152152152152152155555252
55
251000000000000μμμμμμμμμ
其中的质量矩阵不再是对角矩阵，而是满秩矩阵，其表达式如下：
⎥⎥⎥⎥
⎥⎥
⎦
⎤⎢⎢
⎢
⎢⎢⎢⎣
⎡⋅
⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅
⋅⋅+⋅⋅⋅+⋅⋅⋅⋅⋅⋅+⋅⋅⋅++⋅⋅⋅+⋅⋅⋅⋅⋅⋅+⋅⋅⋅++⋅⋅⋅+N N
N
N N N N N N
N N N N
m m m m m m m m m m m m m
m m m m m m m m m 33232
3231编写程序形成M 矩阵时，M 矩阵符合下列表达式： ()i j m M N j
k ij ≥=∑; ()i j m M N
i k ij ≤=∑;
编程时即可形成满秩的质量矩阵。

经过变换后，C 矩阵是一个对角矩阵，原来的
C 矩阵为
⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢
⎢⎢⎢⎣⎡--++--+--+5555
44433332222
1000
0000000000C C C C C C C C C C C C C C C C 经过变换后，C 矩阵为
⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡5432100000000000000000000C C C C C 阻尼矩阵的变化通过下面的程序实现： for i=1:ND
if i==ND
C(i)=CMatrix1(i,i); else
C(i)=CMatrix1(i,i)+CMatrix1(i,i+1); end end
刚度矩阵也是相同的变化。

经过这样的变化之后便于我们研究结构在地震力作用下的各层层间位移。

虽然TMD 系统可以用于减震，但是应该选择合适的频率比
f 和阻尼比ξ，
否则难以发挥则用，甚至起到反作用。

对于理想的无阻尼的系统而言，可以求解出
f 和ξ的解析式，进而求出最优值。

其最优值如下：
)1/(2/1μμ++=f ； ）
（μμξ+=
183
实际建筑工程中的系统都是有阻尼的，而且安装了TMD 的系统激振频率一旦偏离TMD 系统的固有频率，主结构的振幅将急剧增大。

所以，研究有阻尼系统的TMD 是很有实际意义的，此时求解的出的TMD 参数才能真正应用于实际结构中。

考虑了主结构的阻尼的TMD 系统，将无法导出最优频率比和最优阻尼
比的解析式。

虽然可以通过非线性规划的方法得到最优频率比和最优阻尼比的近似解，但是过于繁杂，需要一种程序化的方法来简化求解过程。

二、结构参数
计算结构为下图所示的多自由度体系结构，研究此结构在地震动作用下的位移、速度以及加速度等的变化过程。

结构计算参数如下： 1、kg m 3101000⨯= ；m N K /1020006
⨯=
2、m m m m N =⋅⋅⋅⋅⋅⋅==21 ；k k k k N λ==⋅⋅⋅⋅⋅⋅==21
3、结构参数中5=N ；0.1=λ。

结构计算简图
4、其中每一层的阻尼比为ξ=0.05，TMD 的质量为结构总质量的5%，即
kg M TMD 33102501010005%5⨯=⨯⨯⨯=
5、根据上面的结构图和计算参数，求得在无阻尼时TMD 的最优阻尼C 和刚度
K ：
)1/(2/1μμ++=f =0.9642 ； ）
（μμ
ξ+=
183=0.1336
==TMD TMD TMD TMD m C ξω2=819856
()M f m K TMD TMD μωω22
1===37665357
取710767.3⨯≈K ；5
1020.8⨯=C 。

在求解有阻尼结构的TMD 最优频率比和阻尼比时，以求得的无阻尼情况下的最优阻尼C 和刚度K 作为初值，来求有阻尼系统的最优值。

放大之后细节图为：
第四层楼的层间位移图
放大之后细节图为：
第五层楼的层间位移图
TMD 的相对位移图
由程序可得TMD 的位移幅值为=)max (TMD y 0.062802，其与层间最大位移的比值
1.40154
.00628.01===y y k TMD
再观察TMD 的位移图可知，TMD 的位移远大于结构本身各层的层间位移值，从能量角度而言，TMD 结构通过吸收地震的能量从而达到减小结构地震反映的目的，所以自身的位移值会很大。

实际结构中通常将TMD 装置安装在结构的顶部，由于结构的承载能力和空间有限，所以其质量不能过大，一般都是结构总质量的5%以内。

结构加上TMD 前后的层间位移幅值对比表
位移幅值
位移模式
without TMD with TMD 比值 Y1 0.015408728 0.013124813 0.851777861 Y2 0.013944891 0.011860958 0.850559356 Y3 0.011530666 0.009841923 0.853543284 Y4 0.008222374 0.007079167 0.860963836 Y5
0.004273824
0.003774989
0.883281401
将上述表格的数据制成图，从图上可以清楚地看到加上TMD 之后的效果是很明显的。

蓝色线条表示的是未加TMD 的层间位移，红色线为加了TMD 后的层间位移。

由图表可知，对于结构的层间位移，加了TMD 时的位移幅值大概都是未加TMD 时位移幅值的85%，可见隔震效果比较明显，能够明显的减小各层的位移，对于提高结构的抗震性能非常有用。

此时的频率比和阻尼比还不是待求的最优值，表明在最优条件下，TMD 的作用会更加明显。

我们还可以从表中看到，
无论是否加了TMD，层间位移都是逐层递减的，底层的层间位移（也就是第一层的位移）最大，顶层的层间位移最小。

分析时可以选择任意一个作为我们分析的目标值。

三、优化算法
1、采用Matlab中的优化函数
将结构本身看做一个质点，考虑TMD对这个系统的减震作用，结构分析如下：
加TMD后结构分析图
上图是主结构---TMD 系统模型。

设质量为s m 、T m ，弹簧的刚度系数为s K 、
T K ，阻尼器的粘性阻尼系数为s C 、T C 。

设在地面的运动加速度)(t x
g 作用下，在时刻t ，主结构和TMD 系统相对于基底（x x -轴）的位移分别为)(t x s 、)(t x T ，
而其相对加速度分别为)(t x s 、)(t x g ，绝对加速度分别为)()(t x
t x g s +，根据达朗贝尔原理可以建立“主结构—TMD ”系统的运动微分方程：
()()()()()()()()0000s s s s s s T T s T T s T T T T T T T T T T x t x t x t x t m C C C k k k m x t x t x t x t m C C k k m ⎡⎤⎡⎤⎡⎤⎡⎤
+-+-⎡⎤⎡⎤⎡⎤⎡⎤++=-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥--⎣⎦⎣⎦⎣⎦⎣
⎦⎣⎦⎣⎦⎣⎦⎣⎦ 将结构的动力反应作为函数()C K Y ,；而()C K Y , 实际上是一个关于K 、
C 的隐函数，所以将求解结构层间位移的程序作为目标函数，其中的K 、C 值
用待求的未知量()C K x ,表示。

通过程序不断地搜索任意的()C K x ,0值带入直到搜索到最小的位移，此时的目标函数可以选择是任意一层的层间位移，也可以选择层间位移之和，也就是最上面一层的对地位移，本文选择后面一种，即以最上一层的位移值作为目标函数，通过函数优化找出其最小值时TMD 结构的K 、C 值。

编写程序运行结果如下：
设定的参数：x = 1.0e+07 *(3.766525299961552 0.082213714221286)
也就是37665253=K ；822137=C 。

也就是此种方法下求得的最优值。

将此时的值带入到TMD 程序中，算的各层的层间位移列表如下：
结构加上TMD时优化前后的层间位移幅值对比表
位移模式without TMD with TMD with TMD（yh) 比值
位移幅值
Y1 0.015408728 0.013124813 0.008399832 0.545134673 Y2 0.013944891 0.011860958 0.007161569 0.51356224
Y3 0.011530666 0.009841923 0.005344769 0.463526501 Y4 0.008222374 0.007079167 0.003753254 0.456468411 Y5 0.004273824 0.003774989 0.002732137 0.639272263
从图表上可以清楚地看到加上参数最优的TMD之后的效果是很明显的。

黄色线条表示的是未加TMD的层间位移。

由图表可知，对于结构的层间位移，加了优化的TMD时的位移幅值大概都是未加TMD时位移幅值的50%左右，对于此时的顶层位移比上原顶层位移为51.313814%，可见隔震效果比较明显，能够极大的减小各层的位移，对于提高结构的抗震性能非常有用。

再看优化效果,对于加了TMD的第一种情况，通过参数优化使得各层间位移明显减小。

减少程度均为30%多，可得优化是很明显的。

2、采用遗传算法求解最优值
由于()C K Y , 是一个关于TMD 结构K 、C 的隐函数，无法采用直接求导方法求其极值。

此外，该优化问题为不等式约束优化问题。

对不等式约束优化问题，传统的方法采用罚函数法，此法常因边界条件及惩罚因子的设置不当，而无法收敛。

另外一个解决算法就是遗传算法，遗传算法是基于生物遗传和进化机制、适合于复杂系统的自适应概率优化技术。

通过遗传算法可以求得最优的K 、C 。

由于目前所掌握的Matlab 知识有限，而且此方法特别复杂。

所以还不能将此方法进行实际的编程运用，只是大概了解了一下思想。

【附录】
1、Tuned_Mass_Damper_2DOF
clc;
clear;
close;
%%
global EWave
File_Name = 'E:\Matlab_Code\Matlab_Code\Ground_Motions_Library\IMPV ALL\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', [1 inf]));
%
EWave.Time=zeros(EWave.NPTs,1);
for i=1:EWave.NPTs
EWave.Time(i)= (i-1)*EWave.DT;
end
EWave.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 TMD
global EVector1
global MMatrix1
global CMatrix1
global KMatrix1
global MVec
global KVec
%
ND = 5;
%
MVec = zeros(ND,1);
KVec = zeros(ND,1);
for i = 1:ND
MVec(i) = 1000E3;
end
for 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 ) + 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
end
end
EVector1 = zeros(ND,1);
%
[EIGVecs, EIGVals] = eig(MMatrix1\KMatrix1);
%
WVector = 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);
end
end
%
CMatrix1 = zeros(ND,ND);
MMatrix1 = 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
end
end
%
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) =MMatrix1(i,i);
end
%% Parameters 2 TMD ---- With TMD
global EVector2
global MMatrix2
global CMatrix2
global KMatrix2
global 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);
end
end
for 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
end
end
%
for i = 1:ND+1
if i==ND+1
KMatrix2(i,i)=3.766E7;
else
KMatrix2(i,i) = KVec(i);
end
end
%
for i = 1:ND+1
if i==ND+1
CMatrix2(i,i)=819856;
else
CMatrix2(i,i) =CMatrix1(i,i);
end
end
%
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, Solver1.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 6th 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','LineWidth',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优化程序
options=optimset(options,'tolfun',1e-10);
x0=[37665357 819856];
[x,fval]=fminunc(@fun1,x0,options);
x
fval
3、function f=fun1(x)
%%
global EVector2
global MMatrix2
global CMatrix2
global KMatrix2
global CMatrix1
global MVec2
global ND
global MVec
global KVec
global 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);
end
end
for 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
end
end
%
for i = 1:ND+1
if i==ND+1
KMatrix2(i,i)=x(1);
else
KMatrix2(i,i) = KVec(i);
end
end
%
for i = 1:ND+1
if i==ND+1
CMatrix2(i,i)=x(2);
else
CMatrix2(i,i) =CMatrix1(i,i);
end
end
%
for i = 1:ND+1
EVector2(i) =MMatrix2(i,i);
end
%% 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);
%
f=max(V2(:,1))+max(V2(:,2))+max(V2(:,3))+max(V2(:,4))+max(V2(:,5));
end
4、Aequilate
options=optimset(options,'tolfun',1e-10);
x0=[37665357 819856];
[x,fval]=fminunc(@fun1,x0,options);
x
Fval
（注：运行程序时，先运行TMD程序，再运行Aequilate程序，最后运行Youhua.m 优化程序）
【参考文献】
1、李春祥、刘艳霞、熊学玉 TMD系统最优参数的实用设计方法
2、强跃何建何运祥等改进遗传算法的结构调谐质量阻尼器风振控制的优化设计。