宽高比为4的矩形断面涡激振动响应数值模拟

刘志文，周 帅，陈政清

（湖南大学 风工程试验研究中心，长沙 410082）

涡激振动是由于结构尾流区旋涡脱落所引起的一种振动现象，具有强迫、自激、限幅等特点，是一种典型的非线性振动现象。实际工程结构如桥梁结构、天线、索结构、热交换管及钻井平台立管等都有可能发生涡激共振现象，若设计不当则会在气流作用下产生大幅涡激振动，从而引起结构疲劳或影响结构的使用性能，因此必须予以重视。

钝体绕流和涡激振动在理论研究和工程实际中都有着很重的意义，现有的研究成果大多是通过试验得到的［1，2］，钝体结构涡激振动的数值模拟研究相对较少。Nomura［3，4］采用基于任意拉格朗日 － 欧拉法计算了雷诺数分别为Re=100、Re=140时圆柱和H型柱体的涡激振动响应，该方法假定柱体和流体交界面网格点的运动速度和柱体表面的运动速度相同，在远离柱体的边界上假定网格点固定。该方法最大的缺点是柱体的运动容易引起计算网格的畸变。曹丰产等［5］采用动网格法分析圆柱和流体的耦合作用，采用二阶投影法和多重网格法求解N－S方程，柱体振动方程采用Newmark－β法求解，实现了圆柱涡激振动响应的数值模拟。Sun等［6］采用LES模型和分块迭代方法对雷诺数为Re=200时圆柱涡激振动响应进行了数值模拟。李广望［7］采用任意拉格朗日－欧拉法求解圆柱的涡激振动，结合分块耦合法提高圆柱涡激振动计算效率，但仍对大振幅条件下的网格畸变和缠绕等不能很好解决。方平治、顾明［8］针对二维方柱在高雷诺数条件下的涡激振动进行了数值模拟，计算得到了方柱的涡激振动“锁定”和“拍”现象。徐枫、欧进萍［9］基于动网格技术和滑移网格技术，对二维方柱涡激振动响应进行了数值模拟。Pan等［10］采用k－ωSST湍流模型对低质量－阻尼比参数的二维圆柱涡激振动响应进行了数值模拟，计算得到的流体力、圆柱振动响应以及圆柱尾流涡结构均与试验结果吻合良好。值得注意的是，气流在流经圆柱、方柱时一般会发生“分离”现象，而很少出现“再附”现象。在实际桥梁结构工程中，气流流经主梁断面后一般会同时出现“分离”和“再附”现象。对出现“分离”和“再附”现象的结构进行涡激振动数值模拟对工程实际而言具有重要的意义。

为此本文采用“刚性运动区域+动网格区域+静止网格区域”的思路建立网格［11］，来解决结构振动过程中可能出现的网格畸变与负体积的问题，并将求解结构振动响应的Newmark－β通过UDF嵌入到Fluent中，对宽高比为4的矩形断面结构涡激振动响应进行了数值模拟，该方法区别于文献［9］中的滑移网格技术。

1 计算模型

1.1 控制方程

粘性不可压缩流体的控制方程有质量守恒方程（连续性方程）、动量守恒方程（N－S方程）和能量守恒方程。在进行涡激振动响应计算时不考虑能量方程，在直角坐标系下，基于雷诺平均的连续性方程和N－S方程分别为：

式中，i，j=1，2；ρ为空气密度，即ρ=1.225 kg/m3，μ为动力粘性系数，即μ=1.789 4 ×10－5kg/m·s。

柱体振动控制方程为：

柱体周围流场采用Fluent6.2中的k－ωSST湍流模型进行计算。具体求解设置为：采用SIMPLEC（Semi implicit method for pressure linked equation consistent）算法求解动量方程中速度分量和压力的耦合问题；采用PRESTO（Pressure staggering option）求解速度分量；对流项采用QUICK（Quadratic upwind interpolation for convection kinetics）求解。在进行柱体涡激振动计算时，首先针对静止柱体绕流进行计算，计算物理时间步长为0.005 s，计算残差控制在 5 ×10－4，计算进行到充分绕流为止；然后将柱体突然释放，进行柱体涡激振动响应计算。进行涡激振动响应计算的每个时间步内，先求解流体控制方程，得到速度场、压力场以及作用在柱体上的升力FL（t），通过嵌入Fluent中的UDF程序来提取作用于柱体上的升力，并将该升力代入柱体振动方程（3）的右端，采用Newmark－β法求解柱体振动响应，并将柱体速度赋予柱体周围随柱体一起做刚体运动的随动流体，然后将柱体的速度通过Fluent中的动网格宏DEFINE_CG_MOTION进行传递，使柱体周围的动网格获得速度，从而更新动网格位置，待网格迭代收敛后，整个流场更新完成从而进行下一个时间步的计算，直到计算结束为止。涡激振动计算时间步设置为0.005 s，计算残差控制在2 ×10－3。

1.2 计算区域与网格划分

考虑到工程实际中，气流流经桥梁主梁断面时一般会发生“分离”和“再附”现象，这区别于气流流经方柱和圆柱时的流动现象；同时为验证本文方法的精度，选择文献［12］中宽高比为4的矩形断面进行涡激振动响应数值模拟。矩形断面宽为B=0.3 m，矩形断面高为D=0.075 m，将矩形断面竖向涡激振动系统简化为沿竖向振动的弹簧－质量－阻尼系统，如图1所示。涡激振动计算参数取自文献［12］，自振频率fn=3.905 Hz，单位长度矩形柱体质量为1.362 kg/m，阻尼比为0.177%，对应的折算速度Vred=V∞/fnD（其中V∞为来流风速（m/s），fn为结构振动频率（Hz），D为矩形断面高（m））分别为 4、6、8、9、10、12、14 及 16，对应的雷诺数（以矩形断面高度D为参考尺寸）范围为6 015～24 060。

计算区域为：上游边界距矩形断面为10B，上下侧边界距矩形断面中心距离为10B，下游边界距矩形断面为20B。边界条件设置如下：上游为速度入口，即Velocity-inlet边界条件，湍流强度为Iu=0.5%，湍流粘性比为 10%；出口边界采用Outflow边界条件，即完全发展的出流边界条件；上下两侧采用对称边界条件，即Symmetry边界条件；矩形断面四周采用无滑移边界条件，即No-Slip Wall边界条件。

图1 结构涡激振动模型Fig.1 Vortex-induced vibration model of structure

图2 矩形断面网格分块划分及边界条件Fig.2 Grid blocks and boundary conditions of the rectangular section

为了解决动网格模拟中网格的畸变、缠绕等问题，采用“刚性运动区域+动网格区域+静止网格区域”的方法建立矩形断面涡激振动网格划分区域。在刚性运动区域，采用结构化网格对结构周围的区域网格进行加密，这部分网格在结构振动过程中随结构一起运动，不需要重新划分网格。动网格区域采用三角形网格进行划分，在结构振动过程中，该部分网格进行重新划分，由于该部分网格远离结构，故网格间距较大，网格不容易出现负体积，从而有效解决结构大幅涡激振动可能引起的网格畸变和负体积现象。为了提高计算效率，最外侧网格采用静止的结构化网格进行划分。在划分网格时，首先对矩形断面附近建立边界层网格，矩形断面附近第一层网格中心至矩形断面表面的距离为δ=0.001B，然后分块进行网格划分，整个计算区域网格总数为75 848。图3分别给出了矩形断面周围流场区域网格总体图及局部示意图。

图3 矩形断面网格划分Fig.3 Computational grid of the rectangular section

2 数值模拟结果

2.1 静止矩形断面绕流结果

限于篇幅，图4仅给出了折算风速为Vred=V∞/fnH=9.0时静止矩形断面升力系数时程曲线及升力系数FFT变换幅值频谱曲线。从图4中可以看出气流流经该矩形断面时所产生的升力系数出现明显的振幅不变的振荡现象，表明气流流经该断面时所产生的漩涡已达到稳定状态。表1给出了宽高比为4的矩形断面阻力系数、升力系数根方差以及斯托罗哈数计算结果。从表1中可以看出，本文阻力系数、斯托罗哈数结果与文献［14］试验结果比较接近，但升力系数根方差比文献［13］的结果要偏小。

图4 静止矩形断面绕流数值模拟结果Fig.4 Numerical simulation results of flow past the stationary rectangular section

表1 宽高比为4的矩形断面绕流数值模拟与试验结果Tab.1 Numerical and experimental results of flow past the fixed rectangular section with aspect ratio 4

2.2 弹性支承矩形断面涡激振动模拟结果

根据文献［12］的宽高比为4的矩形断面涡激振动振幅与速度图（即A－V图），分别对折算风速为4、6、8、9、10、12、14 及 16 时进行了涡激振动响应计算。限于篇幅，图5仅给出折算风速Vred=V/fnD分别为4、10及16时，矩形断面升力系数CL、涡激振动响应时程曲线。从图5中可以看出，不同折算风速条件下，弹性支承矩形断面的涡激振动响应振幅是不同的。

图6给出了折算风速Vred=V/fnD分别为4、10及16时矩形断面升力系数CL的FFT变化幅值频谱曲线。从图6中可以看出，当折算风速为Vred=V/fnD=4时，矩形断面升力系数CL的FFT变化幅值频谱对应的卓越频率为fs=1.532 6 Hz，低于矩形断面结构振动频率，未发生锁定现象；当折算风速为Vred=V/fnD=10时，矩形断面升力系数CL的FFT变化幅值频谱对应的卓越频率为fs=3.759 1 Hz，与矩形断面的结构振动频率接近，即出现“锁定”现象；随着折算风速的增加，矩形断面的漩涡脱落频率逐渐增加，“锁定”现象消失。

图7 涡激振动幅值及漩涡脱落频率比随折算风速变化Fig.7 VIV amplitude and ratio of shedding frequency to natural frequency versus reduced wind velocity

图7分别给出了矩形断面涡激振动振幅随折算风速的变化曲线、漩涡脱落频率与结构振动频率之比随折算风速的变化曲线。从图7（a）中可以看出，宽高比为4的矩形断面涡激振动响应数值模拟结果（如锁定区间、最大振幅）均与文献［10］中的风洞试验结果吻合较好。由图7（b）可知，当折算风速为Vred=V/fnD=8～12时，矩形断面漩涡脱落频率与结构振动频率之比接近1，即发生了“锁定”现象。

3 结论

本文针对宽高比为4的矩形断面的涡激振动进行了数值模拟，并与已有试验结果进行了比较，主要结论如下：

（1）采用“刚性运动区域+动网格区域+静止网格区域”的思路建立网格，可以有效解决结构涡激振动响应计算中可能产生的网格畸变、负体积等问题；

（2）通过对宽高比为4的矩形断面的涡激振动数值模拟可知，采用本文所建立的方法进行具有“分离”、“再附”现象的钝体断面的涡激振动是可行的。

需要说明的是，本文仅仅针对形状较为简单的钝体断面进行了涡激振动数值模拟，对于断面形状较为复杂的钝体断面（如考虑防撞护栏、检修车轨道等附属设施的桥梁主梁断面），则需要进一步的研究与试验验证。

［1］鲜 荣，廖海黎，李明水.大跨主梁沿跨向涡振Scanlan非线性模型应用［J］.振动与冲击，2009，28（4）：54－58.

［2］欧阳克俭，陈政清，韩 艳，等.桥面中央开口悬索桥涡激共振与制涡试验研究［J］.振动与冲击，2009，28（7）：199－202.

［3］Nomura T. Finite elementanalysis ofvortex-induced vibrations of bluff cylinders［J］.Journal of Wind Engineering and Industrial Aerodynamics，1993，46，47：587 －594.

［4］ Nomura T.A numerical study on vortex-induced oscillations of bluff cylinders［J］.Journal of Wind Engineering and Industrial Aerodynamics，1993，50：75－83.

［5］曹丰产，项海帆.圆柱非定常绕流及涡激振动的数值计算［J］.水动力学研究与进展（A 辑），2001，16（1）：111－118.

［6］李广望，任安禄，陈文曲.ALE方法求解圆柱的涡激振动［J］.空气动力学报，2004，22（3）：283－288.

［7］ Sun D，Owen J S，Wright N G，et al.Fluid-structure interaction of prismatic line-like structures using LES and Block-Iterative coupling［C］//Aerodynamics Laboratory National Research Council Canada.5th International Colloquium on Bluff Body Aerodynamics and Applications，Ottawa，Canada：2004.133－136.

［8］方平治，顾 明.高雷诺数条件下二维方柱涡激振动的数值模拟［J］.同济大学学报（自然科学版），2008，36（2）：161－165.

［9］徐 枫，欧进萍.方柱非定常绕流与涡激振动的数值模拟［J］.东南大学学报（自然科学版），2005，35（S1）：35－39.

［10］ Pan Z Y，Cui W C，Miao Q M.Numerical simulation of vortex-induced vibration of a circular cylinder at low massdamping using RANS code［J］.Journal of Fluids and structures，2007，23：23 －37.

［11］ Huang L，Liao H L，Wang B，et al.Numerical simulation for aerodynamic derivatives of bridge deck［J］.Simulation Modelling Practice and Theory，2009，17：719－729.

［12］ Masaru M，Tomomi Y，Hitoshi T，et al.Vortex-induced vibration and its effect on torsional flutter instability in the case of B/D=4 rectangular cylinder［J］.Journal of wind Engineering and Industrial Aerodynamics，2008，96：971－983.

［13］ Sun D，Owen J S，Wright N G.Application of thek-ωturbulence model for a wind-induced vibration study of 2D bluff bodies［J］.Journal of Wind Engineering and Industrial Aerodynamics，2009，7，77 －87.

［14］ Nakaguchi H，Hashimoto K，Muto S，An experimental study on aerodynamic drag of rectangular cylinders［J］.J.Japan Soc.Aeronaut.Space Sci，1968，16：1 －5（in Japanese）.