基于ANSYS的颤振频域分析方法

1．引 言
颤振是一种由于气动不稳定引起的自激发散振动。这一振动一旦发生，将导致结构彻底 破坏。1940 年，美国的旧 Tacoma 大桥因风振而坠落的事故震惊了桥梁工程界，桥梁气动稳 定问题的研究得到了普遍重视。尤其是现代桥梁日益向轻柔、大跨径方向发展，基本自振周 期已达到十秒以上，使桥梁气动稳定问题日益突出，确定桥梁颤振临界风速已成为大跨度桥 梁设计中的一个基本问题。
本文提出了一种利用 ANSYS 直接进行大跨桥梁颤振频域分析的有限元模型和方法。该 有限元模型采用 ANSYS 中自定义单元 Matrix27 来模拟桥面受到的自激气动力，而单元 Matrix27 的刚度或阻尼矩阵描述成风速和振动频率的函数。通过阻尼特征值分析，可确定系 统的各个复特征值，其中特征值的实部为阻尼而其虚部则为阻尼振动频率。因此颤振临界风 速的确定即转化为在某一个风速下系统的某个复特征值具有零或者接近与零的实部。最后通 过算例来说明本文提出的模型和方法的正确性和可靠性。本文提出的方法使得桥梁设计人员 和研究人员可以直接利用 ANSYS 来分析大跨桥梁的耦合颤振。
需要指出的是，在添加了 Matrix27 单元的系统模型后，系统的质量 M 跟原来的桥梁一 样而刚度已经变为 K − K ae 而不是原结构的 K 。因此，系统的瑞利阻尼为
4
C = αM + β (K − K ae )
(10)
所以，式（7）所描述的运动方程描述为
2
式中
X
e
和
X&
e
分别为单元
e
的节点位移和节点速度向量；K
e ae
和
C
e ae
分别为单元
e
的气动刚度
和气动阻尼矩阵。类似于一致质量矩阵和集中质量矩阵的，容易导出单元的一致气动力矩阵
和集中气动力矩阵[14]。当采用集中气动力矩阵时，气动刚度和气动阻尼矩阵具有如下形式：
K
e ae
=
⎡K ⎢
e ae1
⎣0
式（7）描述的为桥梁－风整合系统（简称系统）的参数化运动方程，其中风速和振动 频率为系统参数。其复模态特性可以根据复特征值分析得到。对于具有 n 自由度的系统，系
统 共 有 n 对 共 轭 特 征 值 ， λ j = σ j ± iω j 和 n 对 共 轭 特 征 向 量 ，
Φ j = p j ± iq j ( j = 1, 2, L, n) 。如果所有特征值的实部小于零，系统是动力稳定的；如果
Ce4
=
2C
e ae
(5b)
组集单元气动力矩阵后即可得到总体气动力矩阵
Fae = K ae X + Cae X&
(6)
式中 Cae 和 K ae 分别代表总体气动阻尼和气动刚度矩阵。
将式（6）代式（1）可得到系统的运动方程
MX&& + (C − Cae )X& + (K − K ae )X = 0
(7)
从 20 世纪 70 年代以来，随着计算技术的飞速发展，结构分析有了很大的突破，国外相 继出现了许多大型通用有限元分析程序，如 ANSYS[2]、ABAQUS、ADINA 和 MSC/NASTRAN 等，这些程序具有良好的界面、方便的前后处理和强大的计算分析功能以及开放的二次开发 系统。深受结构分析人员青睐，陈政清教授学术梯队从 1998 年引进了 ANSYS5.5 的版本， 在近几年，先后完成了包括洞庭湖三塔斜拉桥静、动力及稳定性分析，黑石铺湘江大桥中承
矩 M se 可表示为竖向位移 h ，水平位移 p 和扭转位移α 的函数[16][10]，如图 1 所示
( ) L se
=
1 ρU 2 2
2B
⎢⎣⎡KH1*
h&
U
+ KH * 2
Bα& U
+ K 2 H *α 3
+ K 2H* 4
h
B
+ KH * p&
5U
+
K 2H* 6
p⎤ ⎥
B⎦
(2a)
D se
自节点位移有关。
ANSYS 的 Matrix27 单元是一种功能很强的单元。如图 2 所示，该单元具有两个节点， 每个节点有 6 个自由度，其单元坐标系和总体坐标系平行；该单元没有固定的几何形状。跟 其它结构分析单元不同的是，它可以通过实常数的方式输入对称或不对称的质量、刚度或阻
尼矩阵[2]。为将上述自激力荷载在 ANSYS 中实现，本文采用了图 3 的计算图示（对某个桥
=
1 2
ρU 2 (2B)⎢⎣⎡KP1*
p&
U
+
KP2*
Bα& U
+
K 2 P3*α
+
K 2 P4*
p B
+
KP5*
h&
U
+
K 2 P6*
h⎤
B
⎥ ⎦
(2b)
( ) M se
=
1 ρU 2 2
2B 2
⎢⎣⎡ KA1*
h&
U
来自百度文库
+
KA2*
Bα& U
+
K 2 A3*α
+ K 2 A4*
h B
+
KA5* p&
，
它
们
是
无
量
纲
风
速
式（2）表示的是桥面单位长度上的受到的气动力，将这些分布荷载转化为作用于单元
两个节点的集中荷载。作用于单元 e 两节点的等效自激力可表示为：
Faee
=
K
e ae
X
e
+
C
e ae
X&
e
(3)
Lse
M se
U
D se
ϕ
α
Wind attack angle
B
图 1 桥梁断面的位移和气动力
当考虑结构阻尼时，式（7）需要适当修改。结构的阻尼通常以模态阻尼比的形式给出， 在进行阻尼复特征值分析时候，需要将模态阻尼比描述为瑞利阻尼形式
C = αM + βK
(8)
式中α 和 β 分别为瑞利阻尼系数。
当只考虑两阶模态 m 和 n 的模态阻尼比时候，瑞利阻尼系数和结构模态阻尼比的关系 可以描述为[6]：
1
式钢管混凝土拱桥动力特性及稳定性分析，益阳茅草街大桥中承式钢管混凝土拱桥静、动力 以及稳定性分析，秦沈客运专线预应力锚固区局部应力分析等较大项目的研究，并利用 ANSYS 的流体力学模块进行车辆－桥梁在横风中的动力特性研究、桥梁断面气动三分力计 算等，在 ANSYS 使用上积累了大量经验。此外通过对 ANSYS 二次开发，将颤振分析的多 模态参与单参数搜索颤振方法嵌入 ANSYS。
经过几十年的研究，桥梁的颤振分析方法，尤其在频域颤振分析方法，已日趋成熟。颤 振分析主要分为两类方法。第一类是基于物理坐标的颤振分析方法，其主要优点不需要预先 选取颤振参与模态，缺点在于计算量大；另一类是基于模态坐标的颤振分析方法，其主要优 点在于采用模态叠加法来描述结构的运动方程，计算量小，其缺点在于需要预先选取颤振参 与模态。Bleich[3]最早采用 Theodorsen 非平稳气动力表达式分析了旧 Tocoma 悬索桥的气动 稳定问题。Scanlan[15]~[16]和他的研究者最早提出了描述自激气动力的颤振导数的概念，并提 出了一套采用实测颤振导数来分析颤振的频域分析方法。谢霁明和项海帆[23]引入结构控制 中的状态空间概念，将自激气动力表达为有理函数形式，提出了桥梁颤振分析的状态空间法。 Agar[1]采用状态空间法，提出了一种双参数搜索的颤振分析方法。Miyata 和 Yamada[12]则直 接在物理坐标分析颤振，提出了不迭代的单参数搜索算法，但无法考虑结构阻尼的影响。 Namini[13]提出了一种双参数搜索颤振分析的 PK-F 法，该方法能够跟踪结构颤振前后的复特 征值的变化。程韶红[25]和张新军[24]等进一步发展了颤振分析的 PK-F 法，他们在颤振分析中 考虑了结构静风效应和风的攻角变化的影响。陈政清[4]将颤振特征运动用模态坐标表示，提 出了多模态参与单参数搜索法，在不考虑结构阻尼时不需迭代，并从能量角度描述了各个参 与模态的贡献。Tanaka[18]和 Kasuchi[11]同样将颤振转化为复特征值问题，采用双参数的行列 式搜索法求解颤振特征方程。Dung[8]等通过引入模态追踪技术，提高了基于物理坐标的颤振 分析方法的计算效率。Chen[5]等将气动力表达式用有理函数来描述，并在状态空间分析了桥 梁的颤振问题。葛耀君和 Tanaka[9]综述了比较了颤振分析的二维方法、多模态方法和全阶方 法，提出了能考虑结构阻尼的颤振全阶分析方法，需要对风速搜索和频率迭代的过程。丁泉 顺等[7]进一步发展了桥梁颤振分析的多模态和全模态方法，他们同样引入单参数搜索的算 法，能够跟踪结构颤振前后的复模态特性随风速的变化。
一种基于 ANSYS 的颤振频域分析方法
华旭刚 1，陈政清 2
（1．香港理工大学 香港； 2． 湖南大学风工程试验研究中心， 湖南长沙 410082）
摘 要：本文提出了一种在 ANSYS 直接进行大跨桥梁三维颤振频域分析的有限元模型和方法。该有 限元模型采用 ANSYS 中自定义单元 Matrix27 来模拟桥面受到的自激气动力，而单元 Matrix27 的刚 度或阻尼矩阵描述成风速和振动频率的函数。通过阻尼特征值分析，确定系统的各个复特征值，其中 复特征值的实部为阻尼而其虚部则为振动频率。因此颤振临界风速的确定即转化为在某一个风速下系 统的某个复特征值具有零或者接近与零的实部。最后通过算例来说明本文提出的模型和方法的正确性 和可靠性。 关键词：大跨桥梁；耦合颤振；复特征值分析；Matrix27；ANSYS
0⎤
K
e ae1
⎥ ⎦
,
C
e ae
=
⎡C ⎢
e ae1
⎣0
0⎤
C
e ae1
⎥ ⎦
(4a)
⎡0 0 0 0 0 0⎤
⎢⎢0 P6*
P4*
BP3* 0 0⎥⎥
Ke ae1
=
⎢0 a⎢⎢0
H
* 6
BA6*
H
* 4
BA4*
BH
* 3
B 2 A3*
0 0⎥ 0 0⎥⎥
(4b)
⎢0 0 0 ⎢
0 0 0⎥ ⎥
⎢⎣0 0 0 0 0 0⎥⎦
U
+ K 2 A6*
p⎤ ⎥
U⎦
(2c)
式中： U 为风速； ρ 为空气的质量密度； B 为桥面宽； K = Bω /U 为无量纲频率；ω 为
振
U~
动圆频率
= U /( fB)
；
H
* i
,
Ai* ,
(i = 1, 2, L, 6) 称 为 无 量 纲 气 动 导 数
或无量纲频率的函数，与桥梁截面的几何形状有关。
⎡0 0 0 0 0 0⎤
⎢⎢0 P5*
P1*
BP2* 0 0⎥⎥
Ce ae1
=
⎢0 b⎢⎢0
H
* 5
BA5*
H
* 1
BA1*
BH
* 2
B 2 A2*
0 0⎥ 0 0⎥⎥
(4c)
⎢0 0 0 ⎢
0 0 0⎥ ⎥
⎢⎣0 0 0 0 0 0⎥⎦
式中， a = ρU 2 K 2 Le / 2 ， b = ρUBKLe / 2 ， Le 为单元的长度。需要指出的是，式（4） 形式还依赖于总体坐标系的方向。由式（4）可知，单元 e 两个节点上受到的自激力只与各
2．基本理论
桥梁结构在均匀流中的运动方程可以描述为
MX&& + CX& + KX = Fse
(1)
式中 M, C , 和 K 分别为结构的总体质量、刚度和阻尼矩阵； X ， X& 和 X&& 分别代表节点
位移、速度和加速度； Fse 表示桥梁受到的自激气动力。
根据 Scanlan 自激气动力表达式，单位长度上受到的气动升力 Lse 、阻力 Dse 和气动扭
存在至少一个特征值的实部大于零，系统是动力不稳定的。颤振的临界状态为在某一无量纲
风速下，系统有且只有一个特征值的实部为零。 对于大跨桥梁，结构的有限元分析模型往往达到成千上万的自由度。计算结构所有复频
率将变得不切实际。但是桥梁的最低颤振临界风速往往发生在结构的低阶频率，所以一般只
需要跟踪结构的低阶复特征值随风速的变化。
( ) α = 2 ωmωn
ω
2 n
−
ω
2 m
ωnξm − ωmξn
(9a)
β
=
2
ω
ω
2 n
mωn −ω
2 m
⎜⎜⎝⎛
−
ξm ωn
+ ξn ωm
⎟⎟⎠⎞
(9b)
式中ωm ,ωn 分别为结构第 m 和 n 阶圆频率； ξm ,ξn 分别对第 m 和 n 阶模态阻尼比。当更
多的模态阻尼比需要等效时，瑞利阻尼系数可以通过最小二乘得到。
Z
Ze
i
Xe
K , C , or M
j
Ye
e
Y
X
图 2 Matrix27 单元示意图
3
e i
j
e1 e3
k
e2 e4
l
图 3 ANSYS 中颤振分析计算图示
从上述分析知，当桥面单元长度相等时，该对单元的刚度矩阵和阻尼矩阵分别为：
K e1
=
2K
e ae
Ce3
=
2C
e ae
(5a)
K e2
=
2K
e ae
面单元 e ）：在每个桥面主梁节点处添加 1 对 Matrix27 单元（包括一个刚度单元和一个阻尼 单元），该单元的一个节点为桥面节点 i 或 j ，另一个节点固定。一个单元 Matrix27 只能模
拟气动刚度或者气动阻尼，而不能同时模拟两者。例如，在节点 i 处，单元 e1 用于模拟气 动刚度而单元 e3 用于模拟节点 i 处受到的等效气动阻尼，单元 e1 和 e3 共用节点 i 和 k。