基于ANSYS平台的桥梁时域颤振阶跃函数算法

卿前志，张志田，肖 玮，朱明坤

(1.湖南大学 风工程试验研究中心，湖南长沙410082;2.上海市政工程设计研究总院(集团)有限公司，上海200092)

结构的颤振稳定性能是超大跨度桥梁设计中至关重要的控制因素之一。目前桥梁颤振性能的评价方法主要有3种，即经典理论法、直接风洞试验测试法和试验加理论法［1－2］。经典理论法是以Theodorsen机翼颤振理论为基础的，主要针对平板机翼的古典耦合颤振问题，后经 Bleich、Klöppel和 Thiele、Selberg和Van der Put等的努力，这一理论被应用于悬索桥颤振的近似计算;直接试验法一般通过节段模型或全桥气弹模型直接在风洞中测试桥梁的颤振性能;试验加理论法则通过强迫振动或自由振动装置测试断面的颤振导数［3］，再将颤振导数应用到结构颤振理论中分析结构的气动稳定性能。

迄今为止，颤振理论分析包含频域和时域两大类方法。频域法又可分为多模态法和全模态法，其通过将颤振运动方程转化为特征值问题，求解方程的特征值确定颤振临界风速。时域法则首先需获得时域化的气动自激力，再在动力有限元分析中进行求解［2］。频域法属于线性分析方法，与之相比，时域分析能方便地考虑各类非线性因素的影响，并能反映结构在颤振后的振幅演变特性，有利于基于过程性能的桥梁气动性能研究而日益受到重视。对于存在流动分离的钝体桥梁断面，以颤振导数表示的自激力是一种时频混合形式，因而通常不能直接用于桥梁时域颤振分析，需要将其转化为时域化的表达式。目前，时域分析气动自激力表达的方法主要有两种:①通过脉冲响应函数结合有理函数进行表达［4］;②沿用机翼理论的方法，采用阶跃函数结合桥梁断面气动导数进行表达［5－6］。采用有理函数法对断面自激力进行表达存在的主要缺陷有两个方面:①与断面的准定常风荷载特性(平均风荷载特性)不相容;②其极限特性不具备明确的物理意义，只适合于结构响应均值为零的情况［7］。鉴于此，笔者采用阶跃函数模拟自激力进行时域颤振分析。

自20世纪70年代以来，随着计算机技术的飞速发展，国内外相继出现了许多大型通用有限元分析软件。ANSYS软件因其强大的前后处理和计算分析功能而在桥梁工程中有着广泛的应用。大跨度桥梁的气弹问题通常不能在一般商业软件中进行分析，但由于ANSYS软件具有良好的二次开发接口功能而被应用于桥梁结构的频域颤振分析中［8］。但是，很少有利用ANSYS二次开发功能进行桥梁时域颤振分析的报道。文献［9］虽然提出了一种在ANSYS中实现颤振时程分析的方法，但是，其必须同时对频率和风速进行搜索，因此实质上仍然是一种频域方法。为了充分利用ANSYS强大的各类非线性分析功能进行大跨度桥梁的非线性气动稳定性研究，笔者提供了在ANSYS中实现颤振时程分析的另一方法，即在瞬态动力学分析过程中嵌入气动自激力的递推算法来考虑其运动状态历史记忆特性，得到结构的颤振临界风速。本文的方法也适合其他风振问题如颤抖振时域分析中的气动自激力模拟，为利用大型通用软件进行大跨桥梁结构抗风分析提供参考。

1 气动自激力阶跃函数表达

结构体系的运动方程为:

式中:［M］，［C］，［K］分别为结构质量矩阵、阻尼矩阵和刚度矩阵;{Fa}为结构上受到的荷载列向量。

忽略断面侧向振动的影响，则作用在断面单位长度上气动自激力可表示如下:

式中:ρ为空气密度;U为来流风速;K=Bω/U为折算频率;B为断面宽度;(i=1～4)为颤振导数，一般通过风洞试验的方法进行识别得到;˙h，˙α分别表示对实时间t的导数。

从式(2)、式(3)中可以看出，气动自激力是折算频率K及风速U的函数。

为了得到时域化自激力表达式，可采用阶跃函数对断面气动自激力进行表达。假设断面在某一时刻相对来流的攻角姿态突然产生一个阶跃位移α0，那么，形成的气动升力时程可表示如下:

式中:C'L为断面升力系数对攻角的导数，C'L=dCL/dα(CL为升力系数，是风攻角α的函数);φ(s)称为阶跃升力函数(s=Ut/B为无量纲时间变量)。

在桥梁风致振动问题研究中，可用如下更加灵活的形式进行表达［6］:

阶跃函数所描述的阶跃攻角变化引起的瞬态气动力变化形式确定后，对于任意形式的小幅振动时程，所引起的气动升力可根据线性叠加的原理表达如下:

式中:α'(σ)=dα/ds。

当采用阶跃函数来表示桥梁断面阶跃位移引起的气动力时，通常将竖向运动以及扭转运动对气动力的贡献分开处理［10－11］。于是，单位长度梁体所受的气动升力与扭矩分别有:

式(7)、式(8)即为采用阶跃函数表示的桥梁断面气动自激力表达式，考虑了扭转位移与竖向速度两种运动状态的影响，也有文献考虑了侧向运动对自激力的贡献，即将式(7)与式(8)中积分项数扩充至3项，分别增加了 φLp与 φMp的贡献［6］。

对比式(2)、式(3)与式(7)、式(8)桥梁断面颤振力表达式后，根据两种表达式频谱特性一致特点，可通过积分变换的方法求得颤振导数与阶跃函数中待定参数的关系，其详细过程可见文献［7，12］。

求出桥梁断面升力以及升力矩对应各运动状态的阶跃函数表达式后，采用动力有限元该法进行三维颤振分析过程中，可采用式(7)与式(8)进行自激力求解。由于两式中有关于运动状态过程的时间积分，因此须化解成为递推的计算式，否则自激力的计算将成为一个十分耗时的过程［7］。将自激升力与升力矩分解成与竖向以及扭转运动相关的两部分:

式中:Lseα(s)，Mseα(s)分别表示扭转运动引起的自激升力与升力矩时程;Lseh(s)，Mseh(s)分别表示竖向运动引起的自激升力与升力矩时程。

以式(10)最后一个等式右边的第1项，即扭转运动引起的自激升力为例，假设s时刻(s=Ut/B)的值已知，通过变换可得式(12):

式中:

从式(12)可知，对于任意时刻的自激力，其求解式中的第1项只与当前运动状态(扭转位移或者竖向速度)有关，而对于后面若干项积分式，则与从结构运动开始到当前的整个运动时间历程有关。因此，要根据 s时刻的 Lseα(s)递推出 s+Δs时刻的Lseα(s+Δs)，只需要解决式(13)的递推计算即可，该式可通过以下方法完成递推:

同理，可以得出 Lseh(s)，Msea(s)，Mseh(s)的递推表达式。

得到断面的气动自激力表达式后，即可方便地在ANSYS中实现颤振时程分析。值得一提的是ANSYS中节点速度和加速度的求取，参照ANSYS帮助文件［13］，对于采用Newmark法的时程分析，节点加速度和速度可表示为如下差分格式:

根据式(15)、式(16)，可通过前一时刻的速度、加速度及当前时刻的位移求得当前时刻的速度和加速度大小，逐步提高风速，便可得到结构的颤振临界风速。

图1给出了在ANSYS有限元分析软件中实现颤振时程分析的计算流程。

图1 ANSYS颤振时程分析流程Fig.1 Flow chart of time history flutter analysis

利用ANSYS软件进行桥梁颤振时域分析的基本步骤如下:

2 数值算例

2.1具有理想平板断面的简支梁

为验证本文模型和求解方法的正确性，首先采用文献［8］中具有理想平板断面的简支板梁进行验证。模型参数为:简支板梁长度L=300 m，宽40 m，约束两端扭转自由度。平板断面的竖向和横向刚度分别为 EIz=2.1× 106MPa·m4，EIy=1.8×107MPa·m4，扭转刚度GIx=4.1× 105MPa·m4。每延米的质量m=20 000 kg/m，质量惯矩 Im=4.5 × 106kg·m2/m，空气密度ρ=1.248 kg/m3。建模时质量矩阵采用集中质量矩阵形式。主梁选用BEAM4单元模拟，质量惯矩采用MASS 21单元模拟，整个模型具有30个桥面梁单元和29个扭转质量单元。为便于和文献［8］计算值进行对比，结构阻尼设为0。

首先对理论平板气动导数进行阶跃函数参数拟合，表1列出了阶跃函数参数拟合值。

表1 平板气动导数阶跃函数参数拟合值Table 1 Thin foil fitted parameters of indicial functions

图2为通过阶跃函数拟合参数反算得到的颤振导数值和理论颤振导数值的对比。

图2 理想平板颤振导数拟合值与理论值的比较Fig.2 Comparison between fitted flutter derivatives and theoretic results of thin foil

从图2可以看出，反算的颤振导数值和理论值之间差别很小，计算结果具有较好的精度。得到阶跃函数参数后，给结构施加初始激励做动力响应分析，逐步提高风速，便可得到结构的颤振临界风速。

图3、图4分别为风速137 m/s及139 m/s简支梁跨中点位移响应曲线。

图3 跨中位移时程(U=137 m/s)Fig.3 Time history displacements of mid-span section(U=137 m/s)

图4 跨中位移时程(U=139 m/s)Fig.4 Time history displacements of mid-span section(U=139 m/s)

由图3可知风速为137 m/s时跨中点竖向及扭转运动具有等幅特性。从图4可知当风速为139 m/s时，结构振动位移出现明显的发散现象。得到结构的颤振临界风速137 m/s，对颤振临界状态位移响应时程做频谱分析得到结构的颤振频率为f=0.386 4 Hz，与文献［8］中提供的颤振临界风速精确解U=136.3 m/s及颤振频率f=0.391 4 Hz非常接近。

图5为给出了风速137 m/s及139 m/s时相轨迹。对比图5(a)、(b)的相轨迹，可知风速达到139 m/s时表征结构运动状态的极限环存在明显发散特性。

图5 相轨迹Fig.5 Phase diagrams

3 结论

桥梁时域颤振分析在非线性特性模拟、后颤振形态以及结构风振全过程再现等方面具有频域方法不可代替的优势，因此时域颤振分析在将来的大跨桥梁抗风研究中将会得到越来越广泛的应用。

笔者讨论了阶跃函数对结构断面自激力进行拟合的方法，在此基础上，利用有限元分析软件ANSYS的二次开发语言(APDL)对结构进行时域颤振瞬态动力有限元分析，实现了具有记忆特性的阶跃函数自激力递推算法与高效、可靠的大型商用软件的结合。

算例表明，基于拟合的气动自激力，在ANSYS三维有限元分析软件中能方便地实现颤振自激力的施加，计算得到的颤振临界风速和相关文献中的报道结果吻合良好。

研究表明，利用ANSYS平台的二次开发功能进行桥梁气动稳定性能的时域分析是可行的，这一方法可为桥梁颤振分析或颤抖振理论分析提高效率，降低自主开发大型软件过程中的算法错误或疏漏所带来的风险。

［1］陈政清.桥梁风工程［M］.北京:人民交通出版社，2005.

［2］项海帆.现代桥梁抗风理论与实践［M］.北京:人民交通出版社，2005.

［3］Scanlan R H，Tomko J J.Airfoil and bridge deck flutter derivatives［J］.Journal of Engineering Mechanics，ASCE，1971(97):1717－1737.

［4］Scanlan R H.Motion－related body force functions in two－dimensional low－speed flow［J］.Journal of Fluids and Structures，2000，14:49－63.

［5］Theodorsen T.NACA Report 496:General Theory of Aerodynamic Instability and the Mechanism of Flutter［R］.VA，USA:Advisory Committee for Aeronautics，1934.

［6］Caracoglia L，Jones N P.Time domain vs.frequency domain characterization of aeroelastic forces for bridge deck sections［J］.Journal of Wind Engineering and Industrial Aerodynamics，2003，91:371－402.

［7］张志田，陈政清，李春光.桥梁气动自激力时域表达式的瞬态与极限特性［J］.工程力学，2011，28(2):75－85.Zhang Zhitian，Chen Zhengqing，Li Chunguang.Limiting and transient characteristics of time－domain expressions for bridge selfexcited aerodynamic force［J］.Engineering Mechanics，2011，28(2):75－85.

［8］Hua X G，Chen Z Q，Ni Y Q，et al.Flutter analysis of long－span bridges using ANSYS［J］.Wind and Structures，2007，10:61 －82.

［9］华旭刚，陈政清，祝志文.一种在ANSYS中实现颤振时程分析的方法［J］.中国公路学报，2008，15(4):32 －34.Hua Xugang，Chen Zhengqing，Zhu Zhiwen.An approach of timehistory analysis of flutter in ANSYS［J］.China Journal of Highway and Transport，2008，15(4):32 －34.

［10］Costa C，Borri C.Application of indicial functions in bridge deck aero elasticity［J］.Journal of Wind Engineering and Industrial Aerodynamics，2006，94(9):859 －881.

［11］Salvatori L，Borri C.Frequency and time－domain methods for the numerical modeling of full－bridge aero elasticity［J］.Computers and Structures，2007，85:675 －687.

［12］Zhang Z T，Chen Z Q，Cai Y Y，et al.Indicial functions for bridge aero－elastic forces and time－domain flutter analysis［J］.Journal of Bridge Engineering，10.1061/(ASCE)BE.1943 －5592.0000176(Sep.6.2010).

［13］ANSYS Inc..ANSYS Theory Reference Documentation for ANSYS 11.0［M］.PA:ANSYS Inc.，2007.