风雨激振中斜拉索倾角对水线及拉索振动的影响
王剑 毕继红 何旭辉 关健 乔浩玥 邵倩
摘要：應用基于滑移理论的两自由度风雨激振模型，借助Comsol软件求解水膜运动方程中的风压力系数和风摩擦力系数，通过数值求解耦合的水膜运动方程和拉索振动方程，得到了不同倾角的斜拉索表面的上水线运动规律、拉索气动力变化规律及其振动响应；通过分析三者之间的内在联系探究拉索倾斜角度对风雨激振的影响。结果表明：斜拉索倾角对于起振风速区间的影响不大；拉索横风向及顺风向振幅均随倾角的增大而减小；上水线振荡频率在拉索倾角较小时接近拉索的自振频率，而当倾角较大时远大于拉索的自振频率；上水线在拉索自振频率附近的振荡幅度随拉索倾角的增大而逐渐减弱，导致气动力变化在拉索自振频率附近的能量越来越小，致使拉索振幅不断减小。
关键词：风雨激振；斜拉索；滑移理论；水线；倾角
中图分类号：TU312+.1；U443.38 文献标志码：A 文章编号：1004-4523（2018）01-0057-10
DOI：10.16385/j.cnki.issn.1004-4523.2018.01.007
引言
斜拉桥上的倾斜拉索随着桥梁跨度的不断增大而越来越长，刚度和阻尼越来越小，在特定风雨条件下极易发生大幅度的低频振动，即风雨激振现象。国内外诸多学者通过一系列的现场观测和风洞试验研究，发现斜拉索表面上水线的形成与振荡是风雨激振现象的重要影响因素。同时，研究人员在人工降雨风洞试验中发现拉索倾斜角度的改变对其振动有明显影响，即拉索振幅随着倾角的增大而显著减小。
由于风雨激振过程中水线厚度较小、形态复杂，加之降雨的影响，难以通过现场观测和风洞试验准确测量水线的运动状态，因此有必要采用数值模拟方法进行相应的研究。自2007年Lemaitre等首次应用滑移理论研究水平静止拉索表面上的水膜形态变化并模拟水线的形成以来，经过Taylor和Robertson、许林汕和葛耀君、毕继红和王剑等的不断完善，滑移理论逐渐发展成为数值模拟研究水线运动的重要方法。
现阶段已有的关于风雨激振中斜拉索倾角的研究还主要局限于拉索振幅，并未触及水线运动，而水线运动是引发风雨激振的重要因素。为此，本文应用已有的基于滑移理论的两自由度二维断面模型，分析不同倾斜角度的拉索表面水膜形态变化（水线运动）、拉索气动力及其振动之间的相互关系，研究斜拉索倾角对水线运动和拉索振动的影响，初步揭示风雨激振的产生机理。
1模型
受水平方向气流和重力共同作用的斜拉索，如图1（a）所示，拉索半径和倾角分别为R，α（0°≤α≤90°），风速为U0，风偏角为β（0°≤β≤90°）。
取图1（a）中的斜拉索A-A断面为研究对象，如图1（b）所示。忽略重力沿斜拉索轴向的分量，则作用在斜拉索断面内的重力分量为
gN=gcosα （1）
忽略轴向气流的影响，只考虑垂直于斜拉索的气流作用，则斜拉索A-A断面内的风速为
（2）
依据滑移理论，通过假设斜拉索表面有一层连续的水膜来模拟降雨，建立斜拉索风雨激振两自由度模型。耦合的水膜运动方程如下
从水膜运动方程（式（4））、拉索振动方程（式（5））及拉索气动力表达式（式（6））可以看出，水膜厚度在重力、表面张力和气流（Cp，Cf）的共同作用下会发生变化并形成水线，导致拉索升力和阻力发生变化（升力和阻力表达式中均有水膜厚度），进而导致拉索发生振动；而拉索振动又会反过来引起水膜形态（水膜运动方程中的拉索加速度）及风压力系数CP和风摩擦力系数Cf发生变化（考虑拉索在横风向及顺风向振动引起的相对速度的变化）；同时，水膜形态的变化也会导致Cp和Cf发生变化。
2数值求解
2.1风压力系数Cp与风摩擦力系数Cf
水膜运动方程中的风压力系数Cp和风摩擦力系数Cf是影响水膜运动的关键参数。水膜在重力、表面张力、风压力和风摩擦力作用下会发生形状改变，而水膜形状的改变必然会导致风压力系数Cp和风摩擦力系数Cf发生变化。由于风雨激振中的水线位置和形状千变万化，不可能通过试验来确定每一时刻的风压力系数Cp和风摩擦力系数Cf，为此，毕继红等在求解方程时应用计算流体力学方法，借助Fluent软件求得表面带有不同形状水膜的斜拉索的Cp和f。这一方法虽然可以求解耦合的水膜运动方程和拉索振动方程，但计算耗时极长，仅能获得前20s的运动状态；而顾明等的试验数据表明，从试验开始到风雨激振真正形成，需要数十秒甚至上百秒的时间。对此，本文采用有限元软件COMSOL计算随时间变化的风压力系数和风摩擦力系数，大幅提高了计算速度；同时由于COMSOL软件可以更好地与MATLAB软件结合，减少了先前计算中所需软件的数目，提高了计算稳定性。
本文采用COMSOL计算时仍采用与原先相同的假设，即假设水膜所受的气流作用与干燥拉索表面所受的气流作用相同。每一时间步内均根据上一步计算得到的水膜形状，将水膜考虑为固态，建立流场中的绕流障碍物，应用Spalart-Allmaras湍流模型，采用稳态计算此时刻的Cp和Cf。
设干燥拉索直径为d，计算区域及网格划分如图3所示。边界条件设置如下：
左侧进口处采用速度边界条件（u=UN，v=0）；右侧出口处采用压强边界条件（p=0）；上下边界处采用完全滑移边界条件；圆柱表面处采用无滑移边界条件（u=0，v=0）。
將应用COMSOL软件计算的结果与Achen-bach的实验数据、Fage和Flanker的实验数据、Celik和Shaffer的数值计算结果及Flu-ent软件的数值模拟结果进行对比。图4为Reg=105时的圆柱表面风压力系数Cp和风摩擦力系数Cf分布图。可以看出，应用COMSOL软件计算得到的Cp与Fage和Flanke测得的实验结果、Ce-lik和Shaffer的数值计算结果及Fluent软件的模拟结果相近，但在背风侧与Achenbach的实验数据有一定差别；Cf则与Achenbach的实验数据、Fage和Flanker的实验数据及Fluent软件的模拟结果接近，而与Celik和Shaffer的数值计算结果在一些位置有一定的偏差。综上所述，应用本文提出的数值模拟方法可以比较精确地获取风压力系数Cp和风摩擦力系数Cf。
2.2数值求解流程
水膜运动方程（式（4））与斜拉索振动方程（式（5a）和（5b））通过斜拉索的升力、阻力及横风向、顺风向加速度相耦合。采用差分法求解方程组，应用MATLAB软件进行数值计算并利用有限元软件COMSOL求解各个时间步不同水膜形态下的风压力系数Cp和风摩擦力系数Cf，基本流程如图5所示。
2.3基本参数及数值计算工况
3数值计算结果
由于风雨激振的真正形成需要数十秒甚至上百秒的时间，本文对上述各工况均进行了140s的数值计算，并选取100～140s的计算结果进行分析研究。
3.1拉索振动响应
图6显示了各个风速工况下斜拉索的横风向及顺风向振幅与Li等的风洞试验的对比。数值计算得到的横风向起振风速区间是6.76m/s<u0</u
3.3单位长度拉索气动力
从图15（a）所示拉索气动升力及阻力最大变化幅度与拉索倾角关系曲线可以看出，倾角α<30°时，升力及阻力变化幅度随α的增大而减小；而当α>30°时，气动力变化幅度却随α的增大而增大。这一变化趋势与上水线位置处的水膜最大厚度变化趋势（图10）相类似，说明上水线厚度的变化幅度对气动力变化幅度有重要影响。
对各倾角拉索的气动力变化时程进行频谱分析，发现各工况下的升力和阻力变化频率大致相同，均接近拉索的自振频率（f0=0.952Hz），如图16（a）所示。
图16（b）显示了拉索自振频率（f0=0.952Hz）附近气动力变化幅值随倾角α的变化关系。随着倾角α的增大，升力和阻力变化幅值不断减小，说明气动力在f0=0.952Hz频率的能量越来越小，从而导致拉索振动逐渐减弱。因此，决定拉索振动强弱的主要因素是气动力变化在拉索自振频率附近幅值，而非整体的最大幅值。
与前述上水线运动分析相对应，分别对图17和18所示的α=30°及α=55°时的升力和阻力时程进行频谱分析，如图19和20所示。两个工况下的升力和阻力的变化频率均接近拉索自振频率，但α=30°时的主频附近的变化幅值明显大于α=55°时，且能量更为集中，因此拉索振动幅度也更大。
4结论
本文基于滑移理论，应用两自由度风雨激振模型，分析拉索倾斜角度對风雨激振的影响，研究不同倾角拉索的气动力变化规律及其振动响应与上水线运动之间的内在联系，得到以下结论：
（1）与风洞试验结果相同，在特定风偏角工况下，拉索两个方向的振幅均随着倾角的增大而减小，但上水线位置处水膜最大厚度却先减小后增大，说明上水线处水膜最大厚度对拉索振动的影响有限。
（2）上水线振荡频率在拉索倾角较小时接近拉索的自振频率；而当倾角较大时远大于拉索的自振频率。上水线在拉索自振频率附近的振荡幅度随拉索倾角的增大而逐渐减弱，与拉索振幅的变化趋势相同。
（3）上水线与拉索之间的共振是导致拉索发生风雨激振的主要原因。随着拉索倾角的增大，上水线在拉索自振频率附近的振荡原来越弱，导致拉索自振频率附近的气动力变化的能量越来越小，致使拉索振动逐渐减弱。
摘要：應用基于滑移理论的两自由度风雨激振模型，借助Comsol软件求解水膜运动方程中的风压力系数和风摩擦力系数，通过数值求解耦合的水膜运动方程和拉索振动方程，得到了不同倾角的斜拉索表面的上水线运动规律、拉索气动力变化规律及其振动响应；通过分析三者之间的内在联系探究拉索倾斜角度对风雨激振的影响。结果表明：斜拉索倾角对于起振风速区间的影响不大；拉索横风向及顺风向振幅均随倾角的增大而减小；上水线振荡频率在拉索倾角较小时接近拉索的自振频率，而当倾角较大时远大于拉索的自振频率；上水线在拉索自振频率附近的振荡幅度随拉索倾角的增大而逐渐减弱，导致气动力变化在拉索自振频率附近的能量越来越小，致使拉索振幅不断减小。
关键词：风雨激振；斜拉索；滑移理论；水线；倾角
中图分类号：TU312+.1；U443.38 文献标志码：A 文章编号：1004-4523（2018）01-0057-10
DOI：10.16385/j.cnki.issn.1004-4523.2018.01.007
引言
斜拉桥上的倾斜拉索随着桥梁跨度的不断增大而越来越长，刚度和阻尼越来越小，在特定风雨条件下极易发生大幅度的低频振动，即风雨激振现象。国内外诸多学者通过一系列的现场观测和风洞试验研究，发现斜拉索表面上水线的形成与振荡是风雨激振现象的重要影响因素。同时，研究人员在人工降雨风洞试验中发现拉索倾斜角度的改变对其振动有明显影响，即拉索振幅随着倾角的增大而显著减小。
由于风雨激振过程中水线厚度较小、形态复杂，加之降雨的影响，难以通过现场观测和风洞试验准确测量水线的运动状态，因此有必要采用数值模拟方法进行相应的研究。自2007年Lemaitre等首次应用滑移理论研究水平静止拉索表面上的水膜形态变化并模拟水线的形成以来，经过Taylor和Robertson、许林汕和葛耀君、毕继红和王剑等的不断完善，滑移理论逐渐发展成为数值模拟研究水线运动的重要方法。
现阶段已有的关于风雨激振中斜拉索倾角的研究还主要局限于拉索振幅，并未触及水线运动，而水线运动是引发风雨激振的重要因素。为此，本文应用已有的基于滑移理论的两自由度二维断面模型，分析不同倾斜角度的拉索表面水膜形态变化（水线运动）、拉索气动力及其振动之间的相互关系，研究斜拉索倾角对水线运动和拉索振动的影响，初步揭示风雨激振的产生机理。
1模型
受水平方向气流和重力共同作用的斜拉索，如图1（a）所示，拉索半径和倾角分别为R，α（0°≤α≤90°），风速为U0，风偏角为β（0°≤β≤90°）。
取图1（a）中的斜拉索A-A断面为研究对象，如图1（b）所示。忽略重力沿斜拉索轴向的分量，则作用在斜拉索断面内的重力分量为
gN=gcosα （1）
忽略轴向气流的影响，只考虑垂直于斜拉索的气流作用，则斜拉索A-A断面内的风速为
（2）
依据滑移理论，通过假设斜拉索表面有一层连续的水膜来模拟降雨，建立斜拉索风雨激振两自由度模型。耦合的水膜运动方程如下
从水膜运动方程（式（4））、拉索振动方程（式（5））及拉索气动力表达式（式（6））可以看出，水膜厚度在重力、表面张力和气流（Cp，Cf）的共同作用下会发生变化并形成水线，导致拉索升力和阻力发生变化（升力和阻力表达式中均有水膜厚度），进而导致拉索发生振动；而拉索振动又会反过来引起水膜形态（水膜运动方程中的拉索加速度）及风压力系数CP和风摩擦力系数Cf发生变化（考虑拉索在横风向及顺风向振动引起的相对速度的变化）；同时，水膜形态的变化也会导致Cp和Cf发生变化。
2数值求解
2.1风压力系数Cp与风摩擦力系数Cf
水膜运动方程中的风压力系数Cp和风摩擦力系数Cf是影响水膜运动的关键参数。水膜在重力、表面张力、风压力和风摩擦力作用下会发生形状改变，而水膜形状的改变必然会导致风压力系数Cp和风摩擦力系数Cf发生变化。由于风雨激振中的水线位置和形状千变万化，不可能通过试验来确定每一时刻的风压力系数Cp和风摩擦力系数Cf，为此，毕继红等在求解方程时应用计算流体力学方法，借助Fluent软件求得表面带有不同形状水膜的斜拉索的Cp和f。这一方法虽然可以求解耦合的水膜运动方程和拉索振动方程，但计算耗时极长，仅能获得前20s的运动状态；而顾明等的试验数据表明，从试验开始到风雨激振真正形成，需要数十秒甚至上百秒的时间。对此，本文采用有限元软件COMSOL计算随时间变化的风压力系数和风摩擦力系数，大幅提高了计算速度；同时由于COMSOL软件可以更好地与MATLAB软件结合，减少了先前计算中所需软件的数目，提高了计算稳定性。
本文采用COMSOL计算时仍采用与原先相同的假设，即假设水膜所受的气流作用与干燥拉索表面所受的气流作用相同。每一时间步内均根据上一步计算得到的水膜形状，将水膜考虑为固态，建立流场中的绕流障碍物，应用Spalart-Allmaras湍流模型，采用稳态计算此时刻的Cp和Cf。
设干燥拉索直径为d，计算区域及网格划分如图3所示。边界条件设置如下：
左侧进口处采用速度边界条件（u=UN，v=0）；右侧出口处采用压强边界条件（p=0）；上下边界处采用完全滑移边界条件；圆柱表面处采用无滑移边界条件（u=0，v=0）。
將应用COMSOL软件计算的结果与Achen-bach的实验数据、Fage和Flanker的实验数据、Celik和Shaffer的数值计算结果及Flu-ent软件的数值模拟结果进行对比。图4为Reg=105时的圆柱表面风压力系数Cp和风摩擦力系数Cf分布图。可以看出，应用COMSOL软件计算得到的Cp与Fage和Flanke测得的实验结果、Ce-lik和Shaffer的数值计算结果及Fluent软件的模拟结果相近，但在背风侧与Achenbach的实验数据有一定差别；Cf则与Achenbach的实验数据、Fage和Flanker的实验数据及Fluent软件的模拟结果接近，而与Celik和Shaffer的数值计算结果在一些位置有一定的偏差。综上所述，应用本文提出的数值模拟方法可以比较精确地获取风压力系数Cp和风摩擦力系数Cf。
2.2数值求解流程
水膜运动方程（式（4））与斜拉索振动方程（式（5a）和（5b））通过斜拉索的升力、阻力及横风向、顺风向加速度相耦合。采用差分法求解方程组，应用MATLAB软件进行数值计算并利用有限元软件COMSOL求解各个时间步不同水膜形态下的风压力系数Cp和风摩擦力系数Cf，基本流程如图5所示。
2.3基本参数及数值计算工况
3数值计算结果
由于风雨激振的真正形成需要数十秒甚至上百秒的时间，本文对上述各工况均进行了140s的数值计算，并选取100～140s的计算结果进行分析研究。
3.1拉索振动响应
图6显示了各个风速工况下斜拉索的横风向及顺风向振幅与Li等的风洞试验的对比。数值计算得到的横风向起振风速区间是6.76m/s<u0</u
3.3单位长度拉索气动力
从图15（a）所示拉索气动升力及阻力最大变化幅度与拉索倾角关系曲线可以看出，倾角α<30°时，升力及阻力变化幅度随α的增大而减小；而当α>30°时，气动力变化幅度却随α的增大而增大。这一变化趋势与上水线位置处的水膜最大厚度变化趋势（图10）相类似，说明上水线厚度的变化幅度对气动力变化幅度有重要影响。
对各倾角拉索的气动力变化时程进行频谱分析，发现各工况下的升力和阻力变化频率大致相同，均接近拉索的自振频率（f0=0.952Hz），如图16（a）所示。
图16（b）显示了拉索自振频率（f0=0.952Hz）附近气动力变化幅值随倾角α的变化关系。随着倾角α的增大，升力和阻力变化幅值不断减小，说明气动力在f0=0.952Hz频率的能量越来越小，从而导致拉索振动逐渐减弱。因此，决定拉索振动强弱的主要因素是气动力变化在拉索自振频率附近幅值，而非整体的最大幅值。
与前述上水线运动分析相对应，分别对图17和18所示的α=30°及α=55°时的升力和阻力时程进行频谱分析，如图19和20所示。两个工况下的升力和阻力的变化频率均接近拉索自振频率，但α=30°时的主频附近的变化幅值明显大于α=55°时，且能量更为集中，因此拉索振动幅度也更大。
4结论
本文基于滑移理论，应用两自由度风雨激振模型，分析拉索倾斜角度對风雨激振的影响，研究不同倾角拉索的气动力变化规律及其振动响应与上水线运动之间的内在联系，得到以下结论：
（1）与风洞试验结果相同，在特定风偏角工况下，拉索两个方向的振幅均随着倾角的增大而减小，但上水线位置处水膜最大厚度却先减小后增大，说明上水线处水膜最大厚度对拉索振动的影响有限。
（2）上水线振荡频率在拉索倾角较小时接近拉索的自振频率；而当倾角较大时远大于拉索的自振频率。上水线在拉索自振频率附近的振荡幅度随拉索倾角的增大而逐渐减弱，与拉索振幅的变化趋势相同。
（3）上水线与拉索之间的共振是导致拉索发生风雨激振的主要原因。随着拉索倾角的增大，上水线在拉索自振频率附近的振荡原来越弱，导致拉索自振频率附近的气动力变化的能量越来越小，致使拉索振动逐渐减弱。
