大跨悬索桥时域颤振有限元精细化分析

2023-12-23吴长青张志田王林凯

振动与冲击 2023年24期

吴长青, 张志田, 王林凯, 罗华

(1. 湖南理工学院土木建筑工程学院,湖南岳阳 414006; 2. 海南大学土木建筑工程学院,海口 570228;3. 安徽省交通规划设计研究总院股份有限公司公路交通节能环保技术交通运输行业研发中心,合肥 230000)

JTG/T 3360-01-2018《公路桥梁抗风设计规范》是以颤振临界风速评估桥梁的颤振稳定性能,因此准确计算颤振临界风速是桥梁颤振分析的关键。颤振分析方法主要有风洞试验法与数值法,两类方法各有优势且相辅相成。试验法可以为数值法提供必要的气动参数且可检验数值模拟的可靠性;而数值法的时间与经济成本相对较低,在全桥颤振分析中尤为突出。数值法包括频域法与时域法,其中频域法研究较早,较为成熟,已被广泛应用于桥梁颤振临界问题的研究[1-3];时域法起步较晚但发展迅速,在颤振稳定及非线性颤振问题的研究中均得到了广泛应用[4-8]。

在桥梁颤振问题的研究中,一般不考虑平均风效应对颤振性能的影响,然而风对桥梁结构产生的各类效应是同时存在的,它们之间相互作用与影响。随着跨度的增加,桥梁结构的平均风效应越来越显著,对颤振性能的影响不容忽视。近年来已有学者研究了平均风效应引起的附加攻角效应与结构刚度变化等因素对颤振临界风速的影响。欧阳克俭等[9]通过节段模型风洞试验得到了-3°～+4°(步长为1°)攻角内的颤振导数,在颤振分析中通过插值算法自动计入附加攻角对颤振导数的影响,研究表明附加攻角效应降低了洞庭湖二桥的颤振临界风速。朱乐东等[10]采用节段模型风洞试验研究了平均风效应引起的附加攻角对桥梁颤振性能的影响,结构表明附加风攻角效应对节段模型系统的阻尼比与气动导数产生了显著的影响,进而提高了象山港大桥的颤振临界风速。熊龙等[11]研究了平均风效应对大跨悬索桥的颤振性能与动力特性变化的影响,结果表明平均风效应在一定程度上降低了桥梁的颤振临界风速。然而,上述研究均是先独立计算桥梁的平均风效应与气弹效应后再叠加,并未实现两者的一体化分析,无法考虑平均风效应引起的结构位形、应力、刚度等变化对颤振及后颤振性能的影响。

基于Scanlan自激力模型,通过阶跃函数法建立了主梁断面的自激力时域模型,采用伪稳态气动力扣除法成功建立了气动自激力与平均风荷载一体化的气动力时域模型。基于该模型,在ANSYS平台上实现了某大跨悬索桥节段及全桥模型的时域颤振分析,评估了平均风效应、几何非线性以及时间步长对桥梁颤振及后颤振性能的影响。

1 一体化气动力时域模型

1.1 自激气动力模型

Scanlan[12]采用颤振导数建立了钝体桥梁断面的气动自激力模型,该模型将正弦运动引起的自激升力及升力矩表示为如下形式

(1)

(2)

Scanlan自激力模型是一种时频混合模型,可借助阶跃函数将Scanlan模型转换为纯时域模型(记为IF型自激力模型)[13],表达式如下

Lse(x,s)=Lseα(x,s)+Lseh(x,s)

(3)

其中,

(4)

其中,

式中:x为结构断面的位置;σ为时间积分变量;s为无量纲时间,s=Ut/B,t为时间;Lse(x,s)与Mse(x,s)分别为断面运动引起的主梁单位长度上的总升力与升力矩时程;Lseh(x,s)与Lseα(x,s)分别为竖向运动与扭转运动引起的升力时程;Mseh(x,s)与Mseα(x,s)分别为竖向运动与扭转运动引起的升力矩时程;C′L,C′M分别为升力系数CL与升力矩系数CM对攻角的导数;α′为扭转位移α对s的一阶导数;h″为竖向位移h对s的二阶导数;φmn为阶跃函数,下标“m”为升力L或升力矩M,下标“n”为竖向运动h或扭转运动α,如φLh为竖向运动引起的与升力相关的阶跃函数。所有阶跃函数的表达式可统一为如下形式

(5)

式中:ai与di为待识别参数;r为阶跃函数所含的指数项个数,若r=3,则对应有6个待识别的参数。参数识别后即确定主梁断面的时域自激力模型,参数识别的方法可参考文献[14]。

1.2 平均风荷载

作用在主梁上的平均风荷载通常采用主梁断面的静风三分力系数来描述,其表达式如下

(6)

(7)

(8)

1.3 一体化的气动力模型

结构风致静力变形会对桥梁的颤振稳定性能产生不可忽略的影响,因此在桥梁颤振分析中需要建立平均风与气弹效应一体化的时域气动力模型。第一种建立方法是直接将平均风荷载与气动自激力叠加,即将式(3)与式(7)、式(4)与(8)对应叠加,表达式如下

(9)

(10)

这种方法简单,认为平均风效应与气弹效应满足叠加原理。已有研究表明,IF型自激力模型的极限值代表了结构平均响应引起的气动力,恰好平均风荷载模型中也包含了这部分气动力[15]。因此,平均风荷载与自激力的简单叠加会重复考虑结构平均响应所引起的那部分气动力。鉴于此,形成了气动力一体化的第二种方法,即采用伪稳态气动力扣除法先从IF型自激力模型中减去被重复计算的气动力部分,然后再叠加平均风荷载,具体表达式如下

(11)

(12)

为验证伪稳态气动力扣除法的有效性,以某大跨悬索桥为例,采用3种工况考察主梁中点的阶跃响应。3种荷载工况具体描述如表1所示。其中工况2与工况3的区别在于后者扣除了伪稳态自激力而前者未扣除。需要指出的是,气弹效应不会对结构最终的静力平衡位置或者平均结构响应产生任何影响。3种工况对应的主梁中点的响应时程,如图1所示。工况1与和工况3对应的响应逐渐收敛于同一极限值(真实的静平衡位置),而工况2最终收敛于一个错误的平衡位置,由此证明了伪稳态自激力扣除法可以有效地解决部分气动荷载被重复考虑的问题。

图1 平均风荷载的阶跃响应Fig.1 Step response to mean wind loads

表1 荷载模式描述Tab.1 Load pattern descriptions

2 数值分析

2.1 工程概况与气动参数

算例为一座主跨1 660 m的大跨度悬索桥,主梁采用闭口流线型钢箱梁断面,其宽度为49.7 m,高度为4.0 m,如图2所示。利用ANSYS软件建立“鱼骨式”的全桥有限元模型,如图3所示。其中主梁、刚臂、桥塔与桥墩均采用空间梁单元beam4模拟,主缆与吊杆采用link10单元模拟。基于风洞试验的方法得到了主梁断面的气动参数。基于遗传优化方法得到阶跃函数各参数值,如表2所示。然后根据参数可反算得到颤振导数的拟合值。三分力系数,如图4所示;颤振导数的试验值与拟合值,如图5所示。由图4、图5可知,两者吻合得较好,拟合精度满足工程要求。

图2 主梁的标准断面图(mm)Fig.2 Standard configuration of the bridge girder(mm)

图3 全桥有限元模型Fig.3 Finite element model of the whole bridge

图4 主梁断面的三分力系数Fig.4 Three-component coefficient of main-girder section

图5 主梁断面的颤振导数试验值与拟合值Fig.5 Tested and fitted values of flutter derivatives

表2 阶跃函数参数值Tab.2 Parameters value of indicial functions

2.2 节段模型颤振分析

基于ANSYS软件建立桥梁节段有限元模型,如图6所示。主梁节段采用beam4单元模拟,弹簧采用link8单元模拟,质量惯性矩采用mass21单元模拟,模型具有竖向与扭转双自由度。节段模型时域颤振分析的主要计算参数如表3所示。

图6 桥梁节段有限元模型Fig.6 Finite element model of the bridge section model

表3 节段模型时域颤振分析的主要计算参数Tab.3 Main calculation parameters of time-domain flutter analysis of the section model

在时域计算中,时间步长Δt的设置会对颤振临界风速的求解造成一定的误差,当Δt趋于零时,这种误差将会消失,然而在有限元计算中,Δt总会有一定的尺度。采用不同时间步长计算得到的颤振临界风速与其相对风洞试验结果的误差,如图7所示。相对误差定义式如下

图7 不同时间步长对应的临界风速以及相对误差Fig.7 Critical wind speeds of flutter and relative errors versus time step

(13)

式中:Ucr为采用某一时间步长计算得到的颤振临界风速;U0为节段模型风洞试验所得的颤振临界风速。

由图7可知,随着Δt的增大,求解的临界风速显著增大。当Δt≤0.035 s时,计算所得的临界风速已具有较小的误差。鉴此, Δt取0.025 s,相对误差为0.25%,已具有足够的精度。

3种工况对应的颤振临界风速值,如表4所示。其中工况A对应的临界风速远大于风洞试验值,而工况B对应的临界风速与风洞试验值很接近。由此可知,忽略平均升力与升力矩的作用将会显著高估桥梁的颤振性能。

表4 节段模型颤振临界风速Tab.4 Critical wind speeds of flutter of the section model

2.3 全桥模型颤振分析

基于ANSYS软件对图3所示的有限元模型进行全桥时域颤振分析,主要计算参数如表5所示。其中竖弯、扭转阻尼比与节段模型时域颤振分析采取的一致,时间步长取0.025 s,风攻角取+3°。全桥时域颤振分析的4种工况描述,如表6所示。所有工况均考虑了自激升力与升力矩;除工况G1外,工况G2、工况G3与工况G4均考虑了几何非线性效应,在ANSYS求解中通过打开大变形开关来考虑几何非线性的影响;工况G3考虑了主梁的平均升力与平均升力矩;工况G4在工况G3的基础上还考虑主梁的平均阻力;所有工况均暂未考虑横向运动引起的气动自激力。

表5 全桥模型时域颤振分析的主要计算参数Tab.5 Parameters of time-domain flutter analysis of whole bridge model

表6 全桥模型颤振分析加载工况描述及颤振临界风速Tab.6 Load pattern descriptions and critical wind speeds of flutter analysis of the whole bridge

不考虑平均风效应时主梁中点的颤振扭转响应时程曲线,如图8所示。当风速U=95.0 m/s时,主梁的位移响应随时间衰减,见图8(a);当风速U=101.8 m/s时,主梁的位移响应首次出现等幅振动现象,见图8(b);当风速U=120.0 m/s时,主梁的扭转响应随时间不断增长,呈现出扭转发散的趋势,见图8(c);由此可确定工况G1的颤振临界风速值为101.8 m/s。工况G2对应的主梁中点位移响应曲线,见图8(d)。由图8可知,当考虑几何非线性效应时,桥梁的后颤振形态最终演变为限幅的软颤振(极限环振动)现象。

图8 不考虑平均风荷载时主梁中点的颤振响应时程Fig.8 Time-histories of self-excited torsional motion of the mid-span deck section with the mean wind effects excluded

工况G3对应的主梁中点颤振扭转响应时程曲线,如图9所示。由图9可知: 工况G3的颤振临界风速为77 m/s,见图9(b);当风速U=90 m/s时,主梁的扭转振幅最终呈现出小振幅的软颤振现象,最大的扭转振幅仅为0.71°,见图9(c)。

图9 考虑平均风效应时(工况G3)主梁中点的颤振响应时程Fig.9 Time-histories of self-excited torsional motion of the mid-span deck section with the mean lift and lift-moment included

工况G4对应的主梁中点颤振扭转响应时程曲线,如图10所示。当风速U=77 m/s,U=90 m/s与U=100 m/s时,主梁中点的位移响应均是先出现一段时间的软颤振随后逐渐衰减,且随着风速的增大,主梁的软颤振振幅增大,但之后振幅衰减的速度也越快。

图10 考虑平均风阻力时(工况G4)主梁中点的颤振响应时程Fig.10 Time-histories of self-excited torsional motion of the mid-span deck section with the mean wind drag included

全桥时域颤振分析各工况的颤振临界风速,见表6。由表6可知: 工况G1与工况G2的颤振临界风速几乎一致,主要原因是颤振临界状态下主梁处于微振状态,几何非线性十分微弱,不会影响桥梁的颤振临界风速;工况G3的颤振临界风速为77.0 m/s,比风洞试验结果小4 m/s,然而节段模型颤振分析中工况B的颤振临界风速与风洞试验值很接近,这可能是由于节段模型与全桥模型在结构刚度、阻尼与边界约束等方面的差异所造成的。工况G2的颤振临界风速为101.6 m/s,远大于工况G3的颤振临界风速,表明忽略平均风效应将大大高估桥梁的颤振稳定性能。需要指出的是,工况G3中U=90 m/s对应的主梁平均扭转角仅为0.092°,表明在所研究的风速范围内主梁的附加风攻角效应并不明显。工况G4在较高风速作用下振动响应最终呈现衰减趋势,其可能原因是在平均阻力的作用下,主梁产生了较大的侧向位移,结构体系的几何非线性效应与主缆系统扭转刚度增强[17],导致结构的振动逐渐衰减。综上所述,考虑平均风效应后,结构的应力状态、刚度特性以及变形形态等因素均发生了显著改变,正是这些因素的改变最终影响了悬索桥的颤振性能与响应特性,但具体的影响机制有待于进一步研究。