谢　炼　张文明　宗周红

(东南大学土木工程学院，南京　210096)



基于ANSYS的桥梁结构颤振时域重启动迭代法

谢炼张文明宗周红

(东南大学土木工程学院，南京210096)

摘要:为了简化桥梁结构颤振时域分析、提高计算效率，以结构自激力脉冲响应函数表达式为基础，基于ANSYS软件的重启动技术，提出了桥梁结构颤振时域计算的重启动迭代方法．分别对具有理想平板截面的简支梁和半漂浮体系斜拉桥灌河大桥进行颤振时域分析来验证该方法的有效性．结果表明:重启动迭代计算方法所获得的结构颤振临界响应与理论解吻合良好，且计算结果对于时间步长的选取没有依赖性，该方法能够真实反应所有参与模态的振动响应，可应用于实桥的颤振时域分析中．

关键词:桥梁结构;颤振时域分析;自激力脉冲响应函数;重启动迭代

引用本文:谢炼，张文明，宗周红．基于ANSYS的桥梁结构颤振时域重启动迭代法［J］．东南大学学报(自然科学版)，2016，46(1) : 133－139．DOI: 10．3969/j．issn．1001－0505．2016．01．022．

大跨度索支撑桥梁结构的颤振研究采用频域分析法和时域分析法．其中，频域分析法计算简便，但尚属于结构的线性分析范畴;时域分析法可以考虑结构的非线性特征，反映结构自激振动的时程变化趋势，故能更精确地得到桥梁结构的振动形态及幅度变化［1］．时域求解中需要利用结构自激力的时域表达式．目前，结构自激力主要存在2种时域表达形式:①阶跃函数结合桥梁断面气动导数的表达形式，可以直接应用于时域计算，但在考虑气动耦合时，气动力阶跃函数往往难以确定［2－4］;②脉冲响应函数与有理函数相结合的表达方式，根据气动导数与折减风速的关系，经一定的参数拟合，直接由桥梁断面非定常气动力传递函数计算结构的自激作用力［5－6］．

目前，诸多学者对桥梁结构颤抖振时域分析进行了研究，并基于有限元编程进行了实际分析计算［7－12］．在ANSYS平台中采用时域方法计算桥梁自激振动响应时，大多采用Matrix27单元在结构上施加气动力的方式．然而，采用Matrix27单元时需要输入结构自振频率，即要对桥梁结构的自振频率进行事先估计，这可能会导致频率估算误差的产生，且计算过程需要人工干预;由于输入频率的存在，采用Matrix27单元的计算方法只能相对准确地反应跟踪模态下的结构自激振动响应，不能准确反应结构多阶模态的耦合作用，利用此方法计算结构抖振时，会导致非常大的自激振动响应计算误差［9－11］．基于自激力的脉冲函数表达式，可将结构自激力以集中力的方式直接施加在结构特定节点上，无需事先估算结构频率便可直接计算结构时程响应信号，进而经Fourier变换得到结构振动频率，从而真实地反应所有参与模态的贡献，计算结果更为可信．

本文根据自激作用力时域表达式，采用ANSYS有限元软件对桥梁结构自激振动响应进行时域分析，提出了一种重启动迭代法．该方法直接将结构自激力以集中力的方式施加于结构的特定节点上，通过每个时间步长内的迭代计算出结构每一时刻的真实振动响应，进而得到桥梁结构自激振动时程响应，识别出结构颤振临界风速和颤振频率．

1　气动自激力时域表达

1．1自激力脉冲函数表达式

结构振动时自激力可以表示为［4－6，13］

式中，Lse(t)，Dse(t)，Mse(t)分别为结构的自激升力、阻力和升力矩; h，p，α分别表示与结构的竖向运动、侧向运动和扭转运动相应的运动状态量; B为桥面宽度; CLh，CLp，CLα，CDh，CDp，CDα，CMh，CMp，CMα为由颤振导数与折减风速关系拟合得到的参数，下文用向量C统一表示; t为时间; Fdr(C，x，t)为结构运动状态与自激力之间的传递函数时域形式，且

式中，ρ为空气密度; U为设定风速;τ为时间积分变量; x为结构的位移; C = { C1，C2，C3，d3，C4，d4}T为待定系数，由相应的颤振导数与折减风速的关系拟合得到，其中d3＞0，d4＞0．式(1)计算自激力时考虑了竖向、侧向和扭转运动，本文不考虑侧向运动的影响．

1．2自激力递推格式

由式(1)可知，自激力时程计算涉及函数积分，因此，在实际计算每一个离散时间点对应的自激力时均需要进行一次积分计算．为了加快计算速度、节省计算内存，采用递推公式代替式(2)中的积分项．以扭转传递函数FMα为例，在ti时刻，对于某一结构点，令

则ti +1= ti+Δt时刻有

其余函数的递推表达式类似于FMα．

2　ANSYS动力平衡方程求解

由自激力递推式(6)，利用ANSYS软件计算结构自激振动响应时，需要在每一个时间步内施加非定常的气动力．ti +1时刻的力与该时刻的运动状态(位移、速度、加速度)相关，只有求得该时刻的运动状态才能准确求出该时刻的自激力，而该时刻的运动状态又依赖于该时刻的自激力，故ti +1时刻的力和位移只能由一个耦合隐式方程同时求解，即需要进行迭代．具体步骤为:①通过迭代求解出某一时刻对应的结构运动状态，进而求解出结构的自激力;②将结构自激力以集中力的方式施加于结构上，继续迭代求解下一时刻的结构运动状态;③按预订时间步长增加时间，直至求得所需要的时程响应．其中，在迭代求解每一步结构准确振动状态时，采用了ANSYS重启动技术，因此本文称这种方法为重启动迭代法．采用重启动技术的优点在于，在ti～ti +1时间内迭代时，时间起点均为ti时刻，每次迭代可以回到相同的时间起点，而不会造成时间的消耗，可以较为精确地考虑时间积分效应．

Escalator system length Ls and split section length A1, A2(including manual input split section length A1, A2) are calculated from SAP system.

令ti时刻结构的位移、速度、加速度、自激力分别为，F(ti)，ti +1时刻结构的位移、速度、加速度、自激力分别为x(ti+Δt)，F(ti+Δt)，则ti和ti +1时刻的运动方程分别为

式中，M，Q，K分别表示质量矩阵、阻尼矩阵和刚度矩阵．

假定ti时刻的运动状态和自激力均为精确解，则式(7)成立．计算ti +1时刻的运动状态时，首先假定F1(ti+Δt) = F(ti)，代入式(8)可以计算得到ti +1时刻的运动状态．但由于实际情况下F1(ti+Δt)与 F(ti)并不相等，故通过运动状态x1(ti+Δt)计算出的F2(ti+Δt)并不等于F1(ti+Δt)，需将F2(ti+Δt)代入式(8)计算，依次循环，直至前后2次计算出的运动状态xj(ti+满足一定要求时，便可求得ti +1时刻结构的真实运动状态和自激力．采用欧几里得范数控制迭代误差，计算式如下:

式中，xj，k为结构中第k个节点第j次位移迭代的结果; N为结构自激力振动计算时需要施加自激力的总节点数;εk1，εk2，εk3为允许的欧几里得范数误差限．竖向和扭转运动状态须同时计入式(9)～(11)．

计算的具体步骤如下:

①假定ti +1时刻F1(ti+Δt) = F(ti) ;

②根据式(8)计算ti +1时刻的xj(ti+Δt)，并根据自激力计算式(1)计算Fj +1(ti+Δt) ;

③利用求得的结构自激力Fj +1(ti+Δt)重新计算ti +1时刻的运动状态

⑤得到ti +1时刻的运动状态和自激力，按时间步长Δt增加时间，重复步骤①～步骤④，直至得到所需总时间为止．

3　算例验证

3．1理想平板截面简支梁模型

为了考察所提方法的精确性、稳定性及可靠性，采用文献［14］中具有理想平板截面的简支梁模型进行验证．设定x轴沿桥纵向，z轴沿桥横向，y轴为竖向．该简支梁跨长L = 300 m，宽度B = 40 m;材料弹性模量E = 210 GPa，泊松比ν= 0．3，竖向抗弯惯性矩Iz= 10 m4，横向抗弯惯性矩Iy= 85．714 m4，抗扭惯性矩It=5．076 2 m4;质量m = 20 000 kg，质量惯性矩Im=4．5×106kg·m2，空气密度ρ=1．225 kg /m．

将简支梁等分为30个单元建立有限元模型，质量集中在节点上．模型边界条件为:左端约束ROTX，UX，UY，UZ;右端约束ROTX，UY，UZ．采用LANB方法进行动力特性分析，前10阶模态计算结果见表1．

3．2简支梁颤振时域计算结果

进行零度攻角下简支梁的颤振响应分析，不考虑侧向运动引起的自激力，与竖向运动、扭转运动相对应的8个颤振导数随折减风速的变化曲线如图1所示．图中，U为来流风速; f为结构频率; H1，H2分别为竖向运动、扭转运动产生的气动阻尼关于自激升力的颤振导数; H3，H4分别为扭转运动、竖向运动产生的气动惯性和气动刚度关于自激升力的颤振导数; A1，A2分别为竖向运动、扭转运动产生的气动阻尼关于自激升力矩的颤振导数; A3，A4分别为扭转运动、竖向运动产生的气动惯性和气动刚度关于自激升力矩的颤振导数．

令时间步长为0．05 s，利用ANSYS软件计算得到的简支梁模型的自激振动响应时程结果见图2．由图可知，简支梁颤振临界风速为137．8 m /s，对颤振临界状态下的时域数据进行Fourier变换，得到相应的颤振频率为0．38 Hz．有限元计算结果与理论值的对比见表2．由表可知，计算值与理论值吻合良好．

表1　简支梁前10阶模态计算结果

图1　理想平板颤振导数

图2　简支梁结构自激振动响应时程曲线

表2　简支梁颤振计算结果

为研究时间步长的变化对计算结果的影响，分别计算了不同时间步长下该简支梁的颤振临界风速值，并与理论值进行比较，结果见表3．由表可知，时间步长的变化不会对计算结果产生较大影响，即不论时间步长为多少，对应每个时间点的运动状态与其自激力总是相协调对应的．因此，利用所提的重启动迭代方法时，结构自激振动响应的有限元计算结果对时间步长的选取没有依赖性．

表3　不同计算时间步长下颤振临界风速的比较

4　灌河大桥颤振响应时域分析

4．1灌河大桥模型

灌河大桥位于江苏省沿海高速(G15)连盐段，处于连云港市和盐城市交界处，于2006年11月2日建成通车．主桥采用双塔双索面半飘浮五跨连续组合梁斜拉桥，跨径组成为32．9 + 115．4 + 340 + 115．4 +32．9 m．主梁采用工字梁，钢纵梁、钢横梁、小纵梁通过节点板及高强螺栓连接形成空间钢架，钢架上架设预制桥面板，现浇膨胀混凝土湿接缝，与钢梁上的抗剪栓钉形成整体，组成组合梁体系．主梁全高3．08～3．41 m，两工字钢纵梁梁肋间距34 m．索塔采用C50型混凝土，上塔柱高42．0 m，中塔柱高61．8 m，下塔柱高15．829 m，索塔在桥面以上高度为96．548 m，塔底中心横向间距为35．6 m．塔柱采用空心箱形截面，单箱单室，塔柱顺桥向尺寸为7．0 m，上、中塔柱横桥向尺寸为4．0 m．斜拉索采用OVM250系列环氧涂层钢绞线拉索，钢绞线标准强度为1．86 GPa，弹性模量为195 GPa，全桥共104根斜拉索．桥垮布置及主梁标准断面如图3所示．

采用鱼骨梁模型在ANSYS软件中建立灌河大桥三维有限元模型，主梁和桥塔采用beam4单元，斜拉索采用link10单元，结果如图4所示．参考横隔梁位置将主梁划分为168个单元，主梁质量集中于节点处．有限元计算得到的结构自振频率与环境振动试验实测值［15］对比见表4．

图3　灌河大桥桥垮布置及主梁截面(单位: cm)

图4　灌河大桥三维有限元模型

表4　灌河大桥自振频率

4．2灌河大桥颤振响应分析

进行零度攻角下灌河大桥的颤振响应分析，不考虑主梁侧向运动引起的自激力，与竖向运动、扭转运动相对应的8个颤振导数随折减风速的变化曲线如图5所示．

图6　灌河大桥自激振动响应时程曲线

令阻尼比为1%，时间步长为0．10 s，利用本文方法计算灌河大桥自激振动响应，结果见图6．由图可知，灌河大桥的颤振临界风速约为172 m/s，与《灌河大桥建设文集》［16］给出的规范估计值170 m /s接近．对临界状态下结构时域振动信号进行Fourier变换处理，得到灌河大桥颤振频率约为0．58 Hz．

为研究时间步长的变化对计算结果的影响，计算了时间步长为0．05 s时灌河大桥颤振临界风速及颤振频率．表5列出了2种时间步长下的计算结果．由表可知，2种时间步长下有限元计算结果非常接近，说明计算结果对于时间步长的选取没有依赖性．

表5　灌河大桥颤振计算结果

5　结论

1)基于桥梁结构自激力的脉冲函数形式，给出了自激力计算的递推表达式，并利用ANSYS软件的重启动技术，提出了一种结构颤振时域分析的重启动迭代法，可用于求解结构颤振临界风速和颤振频率．

2)具有理想平板截面的简支梁数值算例结果表明，根据重启动迭代法所获得的颤振临界风速及颤振频率与理论解吻合良好，且计算对于时间步长的选取没有依赖性．

3)灌河大桥颤振响应分析结果表明，颤振临界风速与《灌河大桥建设文集》给出的规范估计值接近，证明该方法可以较好地应用于实桥的颤振时域分析中．

4)重启动迭代法能够真实反应所有参与模态的振动响应，可以进一步应用于桥梁结构颤抖振时域分析的自激力计算中．

参考文献(References)

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

［2］Scanlan R H．Motion-related body force functions in two-dimensional low-speed flow［J］．Journal of Fluids and Structures，2000，14(1) : 49－63．

［3］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(2) : 371－402．

［4］丁泉顺．大跨度桥梁耦合颤抖振响应的精细化分析［D］．上海:同济大学土木工程学院，2001．

［5］Lin Y K，Li Q C．New stochastic theory for bridge stability in turbulent flow［J］．Journal of Engineering Mechanics，1993，119(1) : 113－127．

［6］Li Q C，Lin Y K．New stochastic theory for bridge stability in turbulent flowⅡ［J］．Journal of Engineering Mechanics，1995，121(1) : 102－116．

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

［8］Ahmad N．Finite element-based flutter analysis of cablesuspended bridges［J］．Journal of Structural Engineering，1992，118(6) : 1509－1526．

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

［10］Zhang W M，Ge Y J，Levitan M L，et al．Aerodynamic flutter analysis of a new suspension bridge with double main spans［J］．Wind＆Structures，2011，14 (3) : 187－208．

［11］曾宪武，韩大建．大跨桥梁风致抖振时域分析及在ANSYS中的实现［J］．桥梁建设，2004，34(1) : 9－12．DOI: 10．3969/j．issn．1003-4722．2004．01．003．Zeng Xianwu，Han Dajian．Time-domain analysis of wind-induced buffet on long-span bridges and implementation of analysis in ansys［J］．Bridge Construction，2004，34 (1) : 9－12．DOI: 10．3969/j．issn．1003-4722．2004．01．003．(in Chinese)

［12］郭增伟，葛耀君．桥梁自激力脉冲响应函数及颤振时域分析［J］．中国公路学报，2013，26(6) : 103－109．Guo Zengwei，Ge Yaojun．Impulse response functions of self-excited force and flutter analysis in time domain for bridge［J］．China Journal of Highway and Transport，2013，26(6) : 103－109．(in Chinese)

［13］Starossek U．Complex notation in flutter analysis［J］．Journal of Structural Engineering，1998，124(8) : 975 －977．

［14］张文明．多主跨悬索桥抗风性能及风致灾变全过程研究［D］．上海:同济大学土木工程学院，2011．

［15］樊星辰．基于健康监测的结合梁斜拉桥安全预警方法研究［D］．南京:东南大学土木工程学院，2014．

［16］田克平，韩大建．灌河大桥建设文集［M］．北京:人民交通出版社，2008: 52－57．

Restart iteration method for time-domain flutter analysis of bridge structures based on ANSYS

Xie Lian Zhang Wenming Zong Zhouhong
(School of Civil Engineering，Southeast University，Nanjing 210096，China)

Abstract:In order to simplify the time-domain flutter analysis of bridge structures and improve the computational efficiency，a restart iteration method is proposed based on the impulse function expression of self-excited force by using the restart technology of ANSYS．The effectiveness of the proposed method is proved by time-domain analysis of a simply supported beam with ideal flat section and a semi-floating cable-stayed bridge—Guanhe Bridge．The results show that the critical flutter response of bridge structures obtained by the proposed method agrees well with the theoretical solution．And the calculation results are independent of the time step．The restart iteration method can consider all vibration modes of structures，thus it can be applied to the flutter analysis of the bridge structure in the time domain．

Key words:bridge structure; time-domain flutter analysis; self-excited impulse response function; restart iteration

基金项目:国家自然科学基金资助项目(51178101，51208104，51378112)、江苏省自然科学基金资助项目(BK2012344)．

收稿日期:2015-06-11．

作者简介:谢炼(1988—)，男，硕士生;宗周红(联系人)，男，博士，教授，博士生导师，zongzh@ seu．edu．cn．

DOI:10．3969/j．issn．1001－0505．2016．01．022

中图分类号:U441．3

文献标志码:A

文章编号:1001－0505(2016) 01-0133-07