低雷诺数下圆柱顺流向和横流向涡激振动响应分析

韦昱呈，鲁 丽

（西南交通大学力学与工程学院，成都 610031）

引 言

在工程中广泛存在由流体诱发的振动问题，如机翼的颤振、蒸汽发生器传热管的流致振动、海洋立管的涡激振动等。目前学术界内被认可的流致振动机理有：周期性漩涡脱落、流体弹性不稳定、声共振以及湍流抖振［1］。在反应堆结构中，圆柱体结构为发生流致振动的主要部件，例如在横流作用下的蒸汽发生器传热棒、燃料棒等。当流体介质以一定速度流经圆柱结构时，会在圆柱尾部两侧交替产生周期性涡脱落，从而使圆柱受到周期性变化的升力，诱发圆柱结构的振动；而圆柱结构的振动又会反作用于流体，改变流体的流动状态，从而改变尾涡模态。这种流体与结构间由于周期性漩涡脱落产生的相互作用被称为涡激振动［2］（Vortex-Induced Vibration，VIV）。

对于二维弹性支撑圆柱的实验研究较多。Feng［3］最早在风洞试验中观察到圆柱结构的涡激振动现象。Williamson等［4-5］在一系列经典水洞实验中发现，当质量阻尼参数（m*ζ）较大时，圆柱涡激振动响应存在两个分支，即初始分支（initialbranch）和下端分支（lowerbranch）；当质量阻尼参数较小时，响应存在3 个分支，即初始分支、上端分支（upperbranch）以及下端分支；各响应分支间存在跳跃现象。吕振等［6］采用肋板蒙皮结构降低圆柱质量比并进行风洞试验，得到了圆柱涡激振动的两个响应分支。目前对圆柱涡激振动的数值模拟较主流的方法是将圆管等三维细长柔性结构简化为一系列质点模型，利用静力学等效方法获得各质点的弹簧支撑参数，则每个质点可视为二维的质量-弹簧-阻尼结构，然后通过建立每个质点的二维涡激振动模型计算其涡激振动响应。该方法近似认为各个质点的涡激振动响应是横流作用下三维细长柔性圆管的涡激振动响应［7-9］。

由于圆柱在顺流方向的振动幅值远小于其在横流方向的振幅，因此在以往大多数研究中，都以圆柱在横流方向上的振动为主，忽略了顺流方向对横流方向振动的影响［10-14］。Jauvtis 等［15］认为，当圆柱的质量比m*<6时，必须考虑顺流向振动对圆柱涡激振动响应的影响。Martins 等［16］采用Williamson 的实验参数对单自由度及两自由度圆柱进行数值模拟，得出两自由度圆柱最大横向振幅为1.5D（单自由度为1.1D，其中D 为圆柱直径）。黄智勇等［17］分别建立单自由度与两自由度圆柱涡激振动数值模型，结果表明，当m*<3.5时，考虑了顺流方向振动的涡激振动模型比单自由度涡激振动模型在横流方向产生了更大的振幅。

基于实验模型和数值模型，较多学者对高雷诺数下单圆柱的单自由度与两自由度涡激振动展开研究，并取得了丰富的成果，而对低雷诺数下单圆柱的单自由度、两自由度涡激振动对比研究较少。基于此，本文采用CFD 商业软件Fluent 求解器以及动网格技术，嵌套自编UDF 程序实现双向流固耦合，分别建立单自由度及两自由度二维弹性支撑刚性圆柱的涡激振动模型，对雷诺数Re=150时的圆柱涡激振动响应进行数值计算分析，并讨论顺流方向的振动对圆柱涡激振动响应的影响。

1 数值模型的建立

通过求解二维非定常纳维-斯托克斯方程（N-S 方程）来获取单圆柱涡激振动的流场流动状态，不可压缩流体的连续性方程可表示为：

不可压缩粘性流体非定常流动的N-S 方程有如下形式：

其中，ρ 为不可压缩流体密度，ui为在i 方向上的瞬时速度分量，p 为压力，ν=μ/ρ 为运动粘度系数，uj为动力粘度系数，t、xi分别为时间、笛卡尔坐标系。

使用CFD 商业软件Fluent 求解器求解流体控制方程，获得流场流动状态，在低雷诺数（Re=150）下进行数值模拟，采用Laminar 层流模型，采用的边界条件为：流场左端为速度入口（Velocity-Inlet），右端为压力出口（Pressure-Outlet），上下采用对称边界（Symmetry），圆柱表面采用无滑移壁面条件（No-SlipWall）。选取的计算流域大小为20D×40D，尾涡轨迹区域长度为30D，圆柱前端以及上下边界距离圆柱均为10D，流体域模型如图1所示。

图1 单圆柱二维涡激振动流场模型

为确保网格质量，流体域网格全部采用结构化网格进行划分。距离圆柱较远区域的网格划分采用指数形式分布，以减少整体网格数量，保证计算效率；在近壁面区域以及圆柱尾流区域对网格进行加密，靠近圆柱表面附近为边界层网格（Y+<1），以确保计算精度要求，整个计算域的网格数量为11 548。流体域网格模型如图2所示。

图2 弹性支撑圆柱流场计算网格模型

二维弹性支撑刚性圆柱结构如图3所示。

图3 二维弹性支撑刚性圆柱结构

根据牛顿第二定律，二维弹性支撑刚性圆柱的运动控制方程可以表示为：

其中，m 为圆柱的质量，c为结构阻尼系数，k为结构刚度系数，FD，FL分别为圆柱受到的阻力以及升力。式（3）还可以表示为：

其中，ω0= k m 为圆柱的固有圆频率系统阻尼比。

作用在圆柱上的阻力及升力可通过流场计算得到，阻力、升力与阻力系数CD、升力系数CL之间符合以下关系：

其中，U∞为流场来流速度，L 为圆柱的长度，对于二维平面，取圆柱长度为1 m；D为圆柱直径。

使用Runge-Kutta 法对式（4）进行离散求解，采用C语言将离散后的结构运动方程写入用户自定义函数UDF（User-Defined Function）中。具体的双向流固耦合算法实现步骤为：在一个时间步长内，先采用Fluent 求解器对流场进行求解，根据边界条件获取流场和二维圆柱表面的压力、速度等信息，利用DEFINE宏提取作用在圆柱表面的阻力及升力；然后代入圆柱的结构运动方程中，通过求解二维圆柱的运动方程，得到当前时间步下圆柱运动的位移和速度；再利用动网格技术，使用当前时间步下圆柱的位移和速度对网格进行更新，然后进入下一个时间步的求解。具体耦合算法流程如图4所示。

图4 双向流固耦合计算流程

2 响应计算及分析

2.1 固定圆柱绕流验证

为确保流场模型计算的准确性，分别在雷诺数Re=100，150，200 情况下进行固定圆柱绕流计算，并将各雷诺数下计算得到的升力系数最大幅值CmaxL同文献［18-19］进行对比，见表1。

表1 不同雷诺数下CmaxL 对比

斯特劳哈尔数St数（Strouhal number）表达式为：

通过对升力系数时程响应数据进行频谱分析，得到尾流涡脱频率fs，再通过式（6）计算不同雷诺数对应的斯托劳哈尔数，同文献［18-19］进行对比，见表2。从表1与表2中可看出，本文使用的CFD 求解模型与前人的计算结果吻合良好，确保了流场计算模型的有效性。

表2 不同雷诺数下St数对比

2.2 单自由度的响应计算

选取的圆柱涡激振动计算参数为：质量比m*=5，阻尼比ζ =0.0056，圆柱直径D=0.0554 m，来流速度U∞=0.5 m/s，雷诺数Re=150。无量纲约化速度为：

由式（7）可见，可通过改变圆柱固有频率fn的大小来实现Ur的变化，Ur的取值范围为2~10。

图5 所示为不同约化速度下横向无量纲位移y/D、升力系数以及阻力系数的响应时程曲线。从图5 中位移比与升力系数之间的关系可以看出，当约化速度较小时，涡激振动位移与升力系数的相位基本相同；随着约化速度的增加，位移与升力系数之间出现反相位。这种圆柱涡激振动位移与升力系数从同相位变为反相位的现象称为“相位开关”。另外，从图5 中升力系数与阻力系数之间的关系可以看出，升力系数的振动周期是阻力系数的两倍，这是由于在圆柱尾部两侧形成周期性涡脱落时，每发生两次涡脱落会使顺流方向的阻力产生一个周期变化，而每发生四次涡脱落才会使横流方向的升力产生一个周期变化。

图5 单自由度横向位移、升力系数、阻力系数时程曲线

图6 所示为单自由度涡激振动响应计算结果。图6（a）所示为圆柱涡激振动的横向无量纲振幅Ay/D 随约化速度Ur的变化曲线，从图中可以看出，随着Ur的变化，圆柱涡激振动横向振幅呈现出3 个不同区域，即初始分支（2≤Ur≤4）、下端分支（4.5≤Ur≤8）以及非同步区域（8.5≤Ur≤10），且下端分支的振幅远大于其余两个区域，表明此时响应进入“锁定”（Lock-In）区域。同时还可以观察到，在下端分支上横向振幅随着Ur的增大而减小，该变化规律同文献［20］中的变化规律一致。

图6 单自由度涡激振动响应计算结果

图6 （b）给出了不同约化速度Ur所对应的升、阻力系数均方根值，从图中可以看出，在初始分支上随着Ur的增加迅速增加到0.9 附近；进入同步区域后，迅速减小至0.07 附近并趋于稳定；退出锁定区域后，跳跃至0.2 附近并随着Ur的增加而缓慢增加。略有不同，在初始分支前端（Ur=2~3.5），稳定保持在一个较小的数值附近；在初始分支末端（Ur=4），突然增加，当Ur=4.5 时，即响应刚进入下端分支时，增加至0.5附近，随后在锁定区域内随着Ur的增加逐渐减小；进入非同步区域后，趋于平稳。

2.3 两自由度的响应计算

在响应求解程序中同时考虑顺流方向和横流方向的自由度，实现两自由度圆柱涡激振动的响应计算。

单自由度与两自由度圆柱涡激振动响应计算结果对比如图7 所示。由图7（a）可以看出，两自由度涡激振动情况下横向振幅同样存在3 个区域，这个现象与单自由度一致，但是两自由度情况下的下端分支振幅大于单自由度，这表明考虑顺流向的振动后，横向振幅有所增加。由图7（b）可看出两自由度情况下升力系数及阻力系数的变化规律与单自由度一致，升力系数在进入下端分支前随Ur的增加而增加，进入下端分支后随Ur的增加而减小直至趋于稳定，进入非同步区域后跳跃至一个较大的值后随Ur的增加而缓慢增加；阻力系数在刚进入下端分支时增加至最大值，而后随Ur的增加而减小，最终趋于稳定。图7（c）所示为各个Ur对应的横向振动频率比fy/fn，从图中可以看出，单自由度和两自由度涡激振动响应的下端分支对应的频率比在1 附近，表明此时振动频率接近固有频率，发生共振，响应进入锁定区域；其余两个区域的频率比则保持在St=0.183 附近。图7（d）给出了不同约化速度下横向位移响应与升力系数间的相位角θ分布图，从图中可以看出，当响应处于初始分支时（2≤Ur≤4），横向位移与升力同相位；进入下端分支后，相位角在Ur=5 与Ur=5.5 之间发生跳跃，随着Ur的增加，相位角逐渐增加至165°；直至进入非同步区域（Ur=8.5~10）后，横向响应与升力系数完全反相位，即相位角为180°。

图7 单自由度与两自由度圆柱涡激振动响应计算结果

图8 所示为不同Ur下两自由度圆柱涡激振动的运动轨迹。图8（a）、图8（c）、图8（d）的约化速度分别位于初始分支、下端分支和非同步区域，从图中可以看出，圆柱在3 个不同响应区域内的运动轨迹呈现出“8”字形，这表明涡激振动具有一定的自限性。由图8（b）可看出，当Ur=4 时，圆柱的运动轨迹较紊乱。图9 所示为Ur=4 时单自由度和两自由度涡激振动各参数的响应时程曲线，图中可见当Ur=4 时，单自由度与两自由度的响应均呈现出“拍”现象。这表明当响应在即将离开初始分支进入下端分支时，圆柱的振动频率与尾流涡脱频率发生分离，圆柱在两个不同频率的相互作用下出现了拍振现象。

图8 两自由度涡激振动运动轨迹

图9 Ur=4时响应的“拍”现象

3 结论

利用CFD 商业软件Fluent 求解器，在定雷诺数（Re=150）下对单圆柱结构涡激振动进行数值模拟，并结合动网格技术以及UDF程序实现双向流固耦合，通过直接改变结构运动控制方程中的圆柱固有频率fn来模拟圆柱在不同约化速度下的涡激振动响应，得出结论如下：

（1）在较大的约化速度范围（Ur=2~10）对圆柱进行单自由度和两自由度流固耦合数值模拟，获得圆柱在不同约化速度下的涡激振动特性，并捕捉到响应的“初始分支”及“下端分支”，观察到“锁定”、“相位开关”和“拍”等现象，锁定区域对应的约化速度范围为Ur=4.5~8.0。

（2）圆柱结构在考虑了顺流方向上的振动后，会使横流方向的振幅有所增加，且在锁定区域内顺流方向的振幅也会增大，因此不能忽略其对横流方向振动造成的影响。

（3）圆柱涡激振动响应仅在即将离开初始分支进入下端分支时出现“拍”现象，其余各处圆柱涡激振动运动轨迹都呈稳定的“8”字形，表明圆柱涡激振动具有一定的自限性。