大跨桥梁具有柔性大、刚度小等特点, 在风荷载作用下易发生风致振动, 熊耀清等^[1]对大跨悬索桥进行研究, 考虑其横风向及竖风向振动特征并提出抗风缆等抗风措施。对于悬索桥及斜拉桥而言, 由于拉索细长, 其在风荷载作用下容易发生各种形式的振动, 驰振是其常见的振动形式。驰振是指细长物体因气流作用发生自激的纯弯曲大幅振动, 呈发散、不稳定性。这种振动最初发现于覆冰输电线, 其振幅高达10倍输电线直径。而斜拉桥的主梁、桥塔、拉索均可能发生驰振现象。现在公认的拉索驰振现象主要分为三类:干索驰振、覆冰驰振和尾流驰振。拉索的振动会加剧索的疲劳损伤, 导致拉索的失效, 从而造成斜拉桥的破坏。

　　 近年来很多专家对覆冰导线的驰振进行了大量的研究, Den Hartog^[2]于1932年基于单自由度的驰振模型首次提出了横风向驰振理论, 认为驰振是一种不稳定的自激振动, 指出模型在风荷载作用下在横风向上产生运动的位移变量, 引起相对风向角发生变化, 由此产生的气动负阻尼是引起驰振的主要因素。Nigol^[3]于1981年提出模型发生扭转才是引起模型与来流风向的相对风向角变化的主要原因。Yu^[4]于1992年给出惯性耦合的失稳理论, 认为覆冰导线或覆冰拉索在风荷载作用下, 拉索出现失稳的情况大部分是水平、竖向和扭转这三个方向耦合运动的结果。

　　 温晓光^[5]利用CFD数值模拟及PⅣ风洞试验, 对三维斜拉索干索驰振机理进行研究。黄韬、李寿英等^[6,7]进行了静力天平测力风洞试验, 研究了新月形、D形和扇形覆冰斜拉索的驰振情况, 并利用CFD数值模拟技术对新月形、D形覆冰模型进行了二位绕流模拟。马文勇等^[8]通过高频天平测力风洞试验研究了8种扇形覆冰导线气动力特性, 讨论了扇形覆冰导线覆冰厚度、覆冰角度及来流紊流度对覆冰导线平均气动力及驰振稳定性的影响。

　　 国内外学者也对理想双圆柱的尾流驰振开展了较多的研究, 包括风洞试验和数值模拟两方面。H Sockel^[9]与胡建华^[10]等对并列刚性节段模型及三维全气弹模型进行风洞试验, 研究其不稳定区域。Tanaka^[11]指出尾流驰振可能发生在两索相距很近的情况也可能发生在两索相距很远的情况, 因此两索间距可分为近距失稳区、稳定区和远距失稳区。李永乐等^[12]利用风洞试验研究了斜拉桥并列拉索尾流驰振问题, 对比分析了拉索间距、来流风向角等对下流索振动特性的影响。马如进、唐浩俊^[13,14]等利用大涡数值模拟及CFD数值模拟方法对串列拉索进行研究, 计算了下游索的驰振不稳定区域。

　　 覆冰会改变拉索的气动外形, 致使串列双索中上游索后方的尾流形态发生变化, 而下游索正好位于上游索的尾流域中, 这种变化会对下游索的尾流驰振产生较大的影响。因此本文将利用CFD数值模拟三维覆冰串列双索流域, 研究覆冰串列双索的尾流驰振形态及驰振临界风速。

　　 气动力参数采用如下定义:

　　 $\begin{array}{l}C{}_{\text{L}}=2F{}_{\text{L}}/(\rho U{}^{2}LB)(1)\\ C{}_{\text{D}}=2F{}_{\text{D}}/(\rho U{}^{2}LB)(2)\end{array}$

　　 式中:C_L, C_D分别为升力系数和阻力系数;F_L为斜拉索模型的升力, 来流速度方向逆时针转动90°为升力的正方向;F_D为斜拉索模型的阻力, 沿来流方向为正;ρ为空气密度, 取1.225kg/m³;U为前方均匀来流风速;B为斜拉索模型截面特征长度;L为斜拉索模型长度。

　　 偏心覆冰拉索在风荷载激励下产生升力, 诱发驰振。驰振产生的条件是气动负阻尼的产生, 表现为平动升力系数出现负斜率, 即拉索会不断从外界吸收能量, 从而产生大幅振动^[15]。依据Den Hartog驰振理论, 当气动力阻尼小于零时结构发生失稳, 即驰振发生的条件为:

　　 $\frac{\partial C{}_{\text{L}}}{\partial \alpha }+C{}_{\text{D}}<0(3)$

　　 式中:$\frac{\partial C{}_{\text{L}}}{\partial \alpha }+C{}_{\text{D}}$为驰振力系数;α为偏心覆冰斜拉索迎风向角。

　　 当气动阻尼与结构阻尼之和为零时为驰振发生的临界状态, 可得驰振临界风速U_cr理论公式^[16]:

　　 $U{}_{\text{c}\text{r}}=-\frac{4m\omega \xi }{\rho B}\frac{1}{C{}_{{\text{L}}^{\prime }}+C{}_{\text{D}}}(4)$

　　 式中:m为振子质量;ω为结构振动圆频率;ξ为阻尼比;$C{}_{{\text{L}}^{\prime }}=\frac{\partial C{}_{\text{L}}}{\partial \alpha }$。

　　 由于拉索的一阶模态频率最小, 相应驰振临界风速最小, 因此, 由一阶模态频率计算拉索的驰振临界风速。拉索驰振临界风速U_{cr, 1}可由下式得到^[6]:

　　 $U{}_{\text{c}\text{r},1}=-\frac{\text{π}m{}_1\omega {}_1\xi {}_1}{\rho B}\frac{1}{C{}_{{\text{L}}^{\prime }}+C{}_{\text{D}}}(5)$

　　 式中:m₁为单位长度斜拉索质量;ω₁为斜拉索一阶模态圆频率, ω₁=2πf₁, 其中f₁取斜拉索一阶模态频率;ξ₁为斜拉索一阶模态的阻尼比。

　　 覆冰斜拉索选取扇形典型冰型, 斜拉索直径为120mm, 覆冰厚度为40mm, 拉索模型分为直索和斜拉索, 拉索长度为600mm, 斜拉索倾斜角度为70°, 其模型外形尺寸见图1。

　　 计算流域采用矩形区域, 计算区域的大小为4.5m×3.0m, 拉索的中心位于原点处, 距离上游流体入口为1.5m, 距离下游流体出口为3.0m, 距离左右流域壁面均为1.5m, 见图2, 图中X轴正方向为速度来流方向, Y轴正方向为升力方向, Z轴为拉索的长度方向。

　　 三维拉索绕流的网格采用结构化网格, 拉索周围边界层网格采用外O形网格, 对边界层网格进行加密, 边界层径向网格增长系数为1.05, 网格划分见图3。

　　 图1 三维覆冰拉索模型示意图

　　  

　　 图2 计算区域尺寸及坐标

　　  

　　 图3 覆冰拉索网格划分

　　 边界条件的定义:上游流体入口定义为速度入口 (Velocity inlet) , 下游流体出口定义为压力出口 (Pressure outlet) , 上下壁面定义为对称边界 (Symmetry) , 其他定义为自由流动边界 (Free slip wall) 。进口速度为12m/s, 湍流强度 (Turbulent Intensity) 设置为3.8%, 湍流粘性率 (Turbulent Viscosity Ratio) 设置为10。　  

　　 求解器 (Solver) 选用基于压力法的求解器 (Pressure Based) , 三维空间 (3D) , 采用非稳态计算方法, 时间步长取0.001s。湍流模型采用SST k-ω (Shear Stress Transport k-ω) 模型, 松弛因子采用默认, 动量、湍动能、比耗散采用二阶迎风格式, 选择SIMPLEC算法。

　　 模拟计算时取5°风向角递增, 风向角示意图见图4, 为了使模拟更精确, 针对每一工况进行建模、划分网格。

　　 图4 风向角示意图

　　 图5～8给出了0°风向角和90°风向角下直索和斜拉索的速度、压力云图。计算结果表明, 无论是直索还是斜拉索, 拉索的迎风面和背风面速度都较小, 在拉索的背风面出现速度的最小值, 拉索的两侧风速较大, 最大值达到17m/s, 直索和斜拉索都出现明显的尾流涡交替现象, 但斜拉索尾流涡的速度大小及发散形式都与直索有着细微差别, 说明拉索的倾斜角度对其自身的驰振有一定影响。从图6、图8可以看出, 在拉索的迎风面出现最大正压, 在拉索的背风面出现最大负压, 0°风向角下直索出现了明显的尾流涡脱落现象, 但0°风向角下斜拉索没有出现尾流涡脱落, 同时拉索四周压力大小及分布形式都有着明显差别。对于不同的风向角, 速度和压力云图都存在较大差别, 这说明, 无论是拉索的倾斜角度还是风向角都对拉索的驰振影响很大。

　　 图9是三维FLUENT模拟计算得到的斜拉索0°风向角和90°风向角时阻力系数和升力系数时程曲线。在0°风向角时, 随着稳定涡街和旋涡脱落的形成, 阻力系数达到稳定, 约为0.85, 而升力系数在0°风向角附近呈稳定的周期变化, 随着时间的推移, 升力系数的周期性变化幅值在减小;在90°风向角时, 随着稳定涡街的形成, 阻力系数在0.8附近稳定, 升力系数在0.37附近呈周期性变化。将升力系数时程曲线进行傅里叶变换, 即可得到旋涡脱落频率和Strouhal数, 为涡激共振研究提供参考。

　　 通过对直索、斜拉索0°～180°风向角的阻力系数和升力系数取平均, 得到各风向角的平均阻力系数与平均升力系数, 见图10。直索的平均阻力系数曲线呈两端高中间低的形状, 在130°风向角和170°风向角左右出现最大值, 均约2.08, 在130°～170°风向角之间数值相对较低, 约0.9, 这主要是因为扇形覆冰拉索在两边风向角的迎风面积大, 在中间风向角的迎风面积小。直索的平均升力系数则呈现两端大中间小的形状, 升力系数除在0°风向角附近有波动外, 在0°～40°风向角之间趋于平缓, 在50°风向角附近上升到0.5左右, 在150°风向角达到最大的负峰值, 约-1.5, 在0°风向角和180°时均接近于0。斜拉索的平均阻力系数与平均升力系数呈现和直索相似的形状, 但平均阻力系数相对直索而言, 变化比较平缓, 说明拉索的倾斜角度会改变驰振气动力系数的大小, 但不会影响气动力改变的大致趋势。

　　 将图10 (a) 和图10 (b) 的气动力系数和黄韬^[6]的风洞试验值 (图11) 做对比, 发现三维模拟与试验数据曲线变化规律一致, 但三维模拟对应的平均阻力系数的最大及最小峰值都大于试验数据对应的峰值, 三维模拟获得的平均升力系数最大正峰值小于试验数值, 而获得的最大负峰值大于试验数值。三维模拟曲线比风洞试验曲线平缓。

　　 图9 风速12m/s时典型风向角的阻力及升力系数时程曲线

　　  

　　 图10 风速12m/s时三维模拟直索和斜拉索覆冰气动力系数

　　  

　　 图11 风速12m/s时风洞试验扇形覆冰斜拉索气动力系数

　　  

　　 图12 风速12m/s时三维模拟直索和斜拉索覆冰驰振力系数

　　 为了判断覆冰拉索发生驰振的可能性, 运用式 (3) 可得到驰振力系数, 但对于一个特定风向角而言, ∂C_L/∂α可能取决于C_L-α曲线的左侧曲线斜率或者右侧曲线斜率^[15]。本文将计算左驰振力系数和右驰振力系数, 当左、右驰振力系数同时小于0时, 才认为拉索处于不稳定状态。从图12可以看出, 本文所研究的扇形覆冰直索及斜拉索模型, 左右驰振力系数均大于0, 故拉索处于稳定状态。

　　 覆冰双索模型及流域尺寸保持和单索一致, 双索的尾流驰振可分为近距失稳和远距失稳, 本文研究覆冰双索的近距失稳, 两索间距为4倍的拉索直径, 为480mm, 具体模型见图13。

　　 三维覆冰双索绕流的网格采用结构化网格, 网格划分及数值模拟方法采用和单索相同的方式, 见图14。

　　 图15给出了同一时刻0°风向角下Z=100, 200, 300, 400, 500, 600mm时监控截面的压力云图。从图中可以看出, 双索没有出现明显的漩涡脱落, 对于上游索而言, 在拉索的迎风面产生最大正压, 在拉索的背风面出现最大负压, 对于拉索不同展向高度处, 压力的分布形式略有差别, 随着拉索高度变化, 相应截面的压力云图是存在相位差的, 说明了空间拉索绕流场的三维流动特性。而对于下游索而言, 由于下游索处于上游索的尾流域中, 所以, 下游索的迎风面出现较大的负压, 背风面压力为0或为负压。

　　 图16和图17为双索典型风向角的速度和压力云图。在0°风向角没有看到明显的尾流涡交替现象。而在90°风向角由于下游索已不再处于上游索的尾流域中, 所以两拉索的速度及压力分布和单索相似。上游索的速度、压力大小分布形式都和单索相似, 而下游索由于上游索的阻挡, 在该索的迎风面出现负压, 随着拉索展向高度的变化, 背风面压力为0或为负压。

　　 图18是三维FLUENT模拟计算得到的0°风向角和90°风向角的阻力系数和升力系数时程曲线。在0°风向角时, 随着周期性的旋涡脱落, 阻力系数在0附近呈周期性变化, 升力系数也在0附近呈稳定的周期变化, 且随着时间推移, 升力系数的变化幅值在减小;在90°风向角时, 随着周期性的旋涡脱落, 阻力系数在0.9附近呈周期变化, 升力系数在0.5附近呈周期性变化。

　　 通过对0°～90°风向角的阻力系数和升力系数取平均, 得到各个风向角下的平均阻力系数与平均升力系数, 其随风向角的变化规律如图19所示。平均阻力系数曲线呈两端高中间低的形状, 在0°～40°风向角之间呈波动状态, 在10°风向角和30°风向角时达到峰值, 约为0.1, 在35°～50°风向角之间急剧下降, 在50°风向角时达到最小值, 约为-0.2, 在50°～90°风向角之间上升。而对于升力系数而言, 呈现两端稍低中间高的山丘状, 在45°风向角达到最大值, 约为0.26, 在0°风向角和60°风向角出现最小值, 约为-0.05。

　　 将图19中三维模拟的串列平行双索在0°～90°风向角之间的气动力系数系数与三维模拟单索的气动力系数相比较, 可以看出串列下游索的气动力系数变化幅度较大, 和单索的气动力变化曲线有较大差异, 这说明下游索的尾流驰振和单索的驰振存在较大差异。

　　 为了判断覆冰拉索发生驰振的可能性, 运用式 (3) 计算得到双索的左驰振及右驰振系数。如图20所示, 当某风向角处左驰振系数及右驰振系数同时小于0时, 才认为拉索产生了驰振运动。拉索在50°, 55°, 60°, 65°, 70°风向角时, 左右驰振力系数都小于0, 发生驰振。

　　 根据图20计算得到的驰振力系数, 本文以某大跨斜拉桥为工程背景计算了拉索的驰振临界风速。表1给出了大跨斜拉桥部分斜拉索的结构参数^[17]。由阻尼比1%及表中各参数, 通过式 (5) 获得50°, 55°, 60°, 65°, 70°风向角下覆冰斜拉索在两侧曲线斜率下的驰振临界风速, 其数值见表2。

　　 从表2可以看出, 在同一风向角下, 左侧曲线斜率和右侧曲线斜率分别对应的斜拉索驰振临界风速差别较大, 在65°, 70°风向角处, 左侧曲线斜率和右侧曲线斜率分别对应的驰振临界风速均较小, 最小临界风速仅为0.738m/s, 易发生覆冰驰振。采用相同的数值模拟方式, 得到新月形覆冰拉索的驰振最小临界风速为7.9m/s, D形覆冰拉索驰振的最小临界风速为16.1m/s, 由于扇形覆冰拉索表面存在尖锐的棱角边, 这些棱角边往往是气流的分离点, 所以, 与新月形及D形覆冰拉索的气动力会有较大不同。马如进等^[13]通过对直径为180mm、间距为6.67倍的串列双索进行数值模拟, 得出斜拉索的最小临界风速为1.52m/s。尾流驰振下游索的临界驰振风速与两索间距及风向角有很大关系, 在本文的工况下, 又由于拉索覆冰的存在, 使得驰振临界风速较小, 所以, 对双索的尾流驰振, 应采取有效的措施。

　　 本文利用FLUENT软件对三维扇形覆冰单索及串列平行拉索进行数值模拟, 得到了覆冰拉索的气动力参数、驰振力系数以及驰振临界风速, 通过数值模拟结果分析得到如下结论:

　　 (1) 在本文分析的扇形覆冰单索风致振动中, 计算得到的单索的驰振力系数都大于0, 即单索不发生驰振。通过对比分析扇形覆水与新月形、D形覆冰拉索的最小驰振临界风速可得出, 斜拉索的覆冰形状及覆冰厚度对拉索的驰振性能影响较大。

　　 (2) 本文分析了平行串列扇形覆冰双索尾流驰振问题, 计算结果表明, 下游索在0°～40°风向角之间不发生尾流驰振, 在40°～70°风向角之间, 有发生驰振的可能, 在70°～90°风向角之间, 驰振力系数在0附近, 不发生驰振。

　　 (3) 对比单索和双索风致振动, 在两串列拉索相互干扰情形中, 上游索的响应情形和单索相似, 而下游索的响应和单索差异较大, 说明在双索相互干扰中, 对上游索影响较小, 而对下游索影响较大。

　　 (4) 本文索间距为4倍的拉索直径, 属于近距失稳。计算得出的某工程的覆冰斜拉索尾流驰振临界风速较小, 一方面和选取的斜拉索结构参数有关, 另一方面跟选取的索间距有关。长索由于其频率低, 尾流驰振风速较低, 下游索发生尾流驰振的可能性很高, 为减小拉索的不利振动, 应采取增加拉索系统阻尼的措施。