低雷诺数下多根圆柱横向受迫振动的数值模拟研究

吴 浩， 孙大鹏， 吕 林

(大连理工大学 海岸和近海工程国家重点实验室， 辽宁 大连 116023)

低雷诺数下多根圆柱横向受迫振动的数值模拟研究

吴 浩， 孙大鹏， 吕 林

(大连理工大学 海岸和近海工程国家重点实验室， 辽宁 大连 116023)

采用三步有限元方法直接求解ALE描述的N-S方程组，结合ALE动网格技术，建立了2D横向受迫振动的数值模型。利用该模型计算了1个主圆柱体与4个附属小圆柱体的横向受迫振动问题，Re=185，振幅比Ae/D=0.2，频率比fe/fs=0.86～1.04，间隔0.01。通过分析总流体力及其频谱随频率比的变化，发现总流体力出现了突变，且流体力的振动模态可以分为4种不同的类型，其中包括基本锁定状态。最后分析了升力突变前后流体-结构能量传递关系、结构位移和升力的相位关系以及涡脱落模式的变化。

多根圆柱；横向受迫振动；三步有限元方法

0 引言

海洋立管的涡激振动问题是近年来学术界高度关注的热点，立管涡激振动的简化模型即流体力学中经典的圆柱绕流问题。圆柱体在稳定流场中的运动可分为自激振动和受迫振动两种类型，围绕着这两种运动类型，学术界开展了广泛而深入的研究，对研究成果比较全面的回顾可见Bearman[1]、Williamson[2]和Sarpkaya[3]的综述文献。自激振动的研究表明，随着缩减速度的变化，圆柱体的振幅、振动频率和涡脱落模态都迅速发生变化，因此，要确定这些影响因素之间的相互联系是非常困难的。而在受迫振动的研究中，圆柱体的振幅、振动频率可以精确的控制。

对于在流场中发生受迫振动的圆柱体，若振动频率fe接近圆柱体固定不动时的涡脱落频率fs，则圆柱体所受流体力和涡脱落模态都会发生显著变化，这是受迫振动最主要的水动力特性。Gu等[4]通过粒径追踪和粒子图像测速方法将受迫振动圆柱的近尾迹可视化，并获得了流线图和速度场分布。对于给定的雷诺数和振幅，随着fe/fs的增加，从圆柱一侧最初形成的尾涡逐渐靠近圆柱底部，最终达到极限位置，然后尾涡的生成位置和符号都发生突变，这种现象在低雷诺数和高雷诺数时都会发生。其后，Lu和Dalton[5]、Guilmineau和Queutey[6]分别采用不同的数值计算方法对Gu的实验结果进行了验证，给出了圆柱体上的流体力系数。Blackburn和Henderson[7]通过数值模拟方法研究了单圆柱受迫振动问题，在计算中固定雷诺数和振幅，关注fe/fs=1附近振动频率的变化对水动力特性和尾迹结构的影响，通过分析流体力系数和流体-结构能量传递关系，将受迫振动时的流动状态划分为4种类型：周期状态(Periodic)、准周期状态(Quasi-periodic)、弱混沌状态(Weaklychaotic)和混沌状态(Chaotic)，并发现在周期状态下存在不对称尾迹结构和非零的平均升力系数。

国内方面，潘志远[8]总结了国内外圆柱体自激振动和受迫振动的主要研究成果，并讨论了两种运动类型的区别和联系。樊娟娟等[9]在N-S方程和k-ω湍流模型的基础上，结合Fluent动网格技术，研究了较高雷诺数下大振幅比圆柱横向受迫振动问题。赵静等[10]在N-S方程和k-ω湍流模型的基础上，利用流线迎风有限元方法结合ALE动网格技术，研究了亚临界雷诺数下圆柱体的受迫振动问题，分析了流体力与振动频率的联系并给出了锁定区间。

实际海洋工程中的立管周围都安装了多根附属管线，如节流压井管线、BOP液控管线和水龙带等。从局部看，当发生涡激振动时，主管与附属管线作为一个整体进行运动，附属管线将对主管周围流场的水动力特性造成明显影响。基于这一工程背景，Wu等[11]开展水池拖曳实验，研究了长细比为1 750的柔性立管上安装4根附属管线的涡激振动问题，通过改变间距和附属管线沿管长的覆盖率，详细分析了附属管线对立管位移振幅、主控频率和振动模态等结构响应特性的影响。利用三步有限元方法直接求解ALE描述的N-S方程组，结合ALE动网格技术，建立了2D横向受迫振动数值模型。首先求解经典的单圆柱横向受迫振动问题，对数值模型进行验证，然后进一步将该数值模型应用于计算1个主圆柱体与4个附属小圆柱体的横向受迫振动问题，在计算中给定雷诺数和振幅，关注fe/fs对圆柱体的总流体力系数和涡脱落模态的影响。

1 数值模型

对于密度为常数、忽略质量力、层流态的二维不可压缩N-S方程组简化形式，选取圆柱体直径D(对于多根圆柱体，指主圆柱体直径)、来流流速U和时间尺度D/U作为无量纲化尺度参数，并采用ALE方法处理圆柱体运动边界，可以得到ALE描述的无量纲形式N-S方程组：

式中：u′=u/U、x′=x/D、t′=tU/D、p′=p/ρU2、Re=UD/υ，υ为运动学粘性系数；u′m=um/U为无量纲ALE网格运动速度分量；i,j=1, 2为坐标分量。为方便计算，将U和D均取为常数1。

采用高阶时间格式和Galerkin有限元空间离散结合的三步有限元方法求解ALE描述的无量纲形式N-S方程组，该方法的具体离散求解过程、网格依赖性和计算稳定性分析可参考吕林[12]和唐国强[13]的工作。

在受迫振动的计算中，将圆柱体周围的有限范围指定为ALE网格更新区域，在该区域内将网格看作弹性模量为常数的弹性体，这样由于边界条件改变而引起的网格坐标变化可以用Laplace方程表示。在每一时间步，预先给定结构的运动位置，即网格的运动边界条件，然后在ALE区域内求解Laplace方程自动生成新的网格坐标，从而得到网格运动速度u′m，ALE区域外的网格则不用更新，有效地减小了计算量。测试表明，通过选择合适的ALE区域大小，网格在运动中不会发生畸变，且保持了良好的正交性。

2 无量纲参数

2.1 流体力系数

通过求解N-S方程得到流场数值解，即可根据圆柱表面的压力分布以及速度梯度计算圆柱受力，圆柱单位长度上的正向力Fx和横向力Fy分别为：

2.2 流体-结构能量传递系数

结构在发生横向受迫振动时的瞬时位移为y(t)，其无量纲形式为y′(t)=y(t)/D，在每一个振动周期内，流体向结构传递的机械能可以写成如下的无量纲形式：

式中：T为结构受迫振动的周期。

E被称为流体-结构能量传递系数，E为正值则流体对结构做功，E为负值则结构对流体做功。结构位移和升力的基本谐波之间的相位角φ也常用于衡量流体-结构能量传递，当E为正值时，0°<φ<180°。但是相位角仅仅反映了流体力的基本谐振动与结构之间的能量传递信息，当流场处于非周期状态时，升力的基本谐波可能是不存在的，因此在结果分析中采用能量传递系数更加准确直观。

Blackburn和Henderson[7]指出，单根圆柱体受迫振动时，尾涡生成和圆柱体运动之间的相位关系变化伴随着E符号的变化。在单根圆柱体受迫振动的基本锁定区，尾迹处于周期性的涡脱落状态，升力为稳定的周期振动，此时结构位移和升力的相平面图为一闭合环路，且E的符号与闭合环路的方向有关，对于正的能量传递，环路方向为顺时针，反之对于负的能量传递，环路方向为逆时针。

3 计算结果

3.1 数值模型的验证

图1 单圆柱横向受迫振动的流体力系数，Re=185，Ae/D=0.2

图2 多根圆柱体计算模型

CDCLrmsSt该文1．3180．4470．193Lu和Dalton[5]1．3100．4220．195Guilmineau和Queutey[6]1．2870．4430．195

3.2 多根圆柱体横向受迫振动

多根圆柱体数值计算模型如图2所示，主圆柱的无量纲直径D=1，以主圆柱中心为原点，建立笛卡尔直角坐标系，x轴方向为顺流方向(IL)，y轴方向为横流方向(CF)。分别对应IL和CF方向，对称放置4个附属小圆柱，其无量纲直径d=0.25，主圆柱外表面到附属小圆柱外表面的无量纲距离为0.375。

如图3所示，计算域范围为60D×50D，主圆柱中心距离上、下、左边界25D，距离右边界35D，在主圆柱周围8D×8D范围内为ALE网格运动区域。整个计算域采用结构化网格进行剖分，总单元数为37 470，总节点数为37 840，主圆柱表面网格尺寸约为0.02D，附属小圆柱表面网格尺寸约为0.01D，圆柱体附近的计算网格剖分如图4所示。计算时间步长Δt′为无量纲时间0.001。

图3 计算域示意图 图4 圆柱体附近的网格剖分

图5 多根固定圆柱体绕流的总力系数时程曲线和FFT频谱，Re=185,Ae=0.2

图6 多圆柱横向受迫振动的总流体力系数，Re=185，Ae/D=0.2

在多根圆柱体横向受迫振动的计算中，主圆柱和4个附属小圆柱保持相对位置不变，整体沿横流方向做简谐运动，运动方程为y=Aesin(2πfet)。在计算中限定Re=185，Ae/D=0.2，则流体力和涡量场是唯一变量频率比fe/fs的函数，因此在fe/fs=0.86～1.04范围内，将fe/fs的计算间隔精细到0.01，从而能够准确地捕获流体力和涡脱落模式发生转变的时机。

选取总流体力系数CD和CL的振动稳定段进行傅里叶变换，则得到了总流体力系数频谱在频率比fe/fs上的分布情况，如图7所示。图7中f′ -fe/fs平面的虚线从左至右依次为f′=1feD/U,2feD/U,3feD/U,4feD/U，这意味着在这些虚线上，无量纲频率f′依次为结构无量纲振动频率feD/U的1～4倍。从图7中可以看出，在各频率比上，CD的主控无量纲频率为2feD/U，CL的主控无量纲频率为feD/U。这与单根圆柱体横向受迫振动是相同的，说明了多根圆柱体在受迫振动时，总流体力仍然受结构运动控制。与CDrms和CLrms相对应，CD主控频率的幅值在计算范围内随fe/fs的增加而增加。CL主控频率的幅值在fe/fs<0.94时随fe/fs的增加而减小，在fe/fs>0.94时随fe/fs的增加而增加。

CL频谱的频率组成在以下4个频率比区间上具有明显不同的特征：在fe/fs=0.86～0.89区间，除主控无量纲频率feD/U以外，在3倍频上也有单峰出现，但幅值非常小，同时在超低频段和1、2倍频附近都有宽频带小幅值分布，可见该区间的CL振动以feD/U的奇数倍频谐振为主，并具有较弱的随机性；在fe/fs=0.90～0.93区间，仅feD/U的1、3倍频上有单峰出现，且3倍频的幅值非常小，但不存在宽频带小幅值分布，可见该区间的CL振动是feD/U的周期性奇数倍频谐振；在fe/fs=0.94～0.99区间，与前两个区间相比，频率组成发生了明显的变化，主控无量纲频率feD/U的幅值突然减小，同时在多个非倍频位置也出现峰值，并且在整个频段上都有宽频带小幅值分布，可见该区间的CL振动具有明显的随机性；在fe/fs=1.00～1.04区间，feD/U的1、3倍频的幅值显著增加，而非倍频位置的幅值变化不大，整个频段的宽频特性突然消失了，可见该区间的CL振动是feD/U的奇数倍频谐振和非倍频谐振结合的形式，并具有一定的周期性。通过对CL频谱的分析可知，从低到高4个频率比区间上的CL振动分别对应Blackburn和Henderson[7]描述的弱混沌状态、周期状态、混沌状态和准周期状态，将这4个区依次标记为I、II、III和IV区。CD频谱随fe/fs的变化与CL频谱很相似，也可以根据频率组成划分为4个频率比区间，不同的是CD频谱除主控无量纲频率2feD/U以外，在4倍频上有单峰出现且幅值非常小，即CD振动中存在偶数倍频谐振。

图7 多圆柱横向受迫振动的总流体力系数频谱(虚线从左至右依次为f′=1, 2, 3, 4 feD/U)

图9给出了流体力的振动稳定段上能量传递系数的均方根Erms随fe/fs的变化。从图9中可以看出，只有在II区(周期状态)Erms接近于0，说明升力处于稳定的周期振动状态。III区(混沌状态)的Erms是最大的，说明升力振动的随机性很强。I区(弱混沌状态)和IV区(准周期状态)的Erms则介于II和III区之间。

图8 能量传递系数平均值随fe/fs的变化 图9 能量传递系数均方根随fe/fs的变化

图10给出了流体-结构能量传递关系发生突变前后，两个频率比fe/fs=0.93和0.94上的总升力系数CL时间过程线。从图10中可以看出，当fe/fs=0.93时，升力是稳定的周期振动，且如前所述，因为存在较弱的奇数阶谐振动，在每半个周期上时程线并不关于极值左右对称。频率比fe/fs增加了0.01，升力就变为随机振动。相应的，在这两个频率比上能量传递系数E时间过程线以及结构位移和升力的相平面图也具有完全不同的特征。如图11、图12所示，当fe/fs=0.93时，E为负值且几乎恒定不变，相平面图为一闭合环路，环路方向为逆时针与E的符号对应，而fe/fs=0.94时，E绝大多数时间为正值且波动很明显，相平面图不是闭合环路，其方向大体上为顺时针仍然与E的符号对应。因此，升力系数变化趋势和振动模态的突变与结构位移和升力的能量传递关系以及相位关系的突变有关。

图10 总升力系数时间过程线，fe/fs = 0.93, 0.94

图11 能量传递系数时间过程线 图12 横向位移和升力的相平面图

图13 一个运动周期内的瞬时涡量场，fe/fs = 0.93, 0.94 (实线：涡量为正；虚线：涡量为负)

图13给出了fe/fs=0.93和0.94时，一个运动周期内的无量纲位移y′和总升力系数CL时间过程线，并给出四个代表性时刻：y′=0(向上运动)、y′=0.2、y′=0(向下运动)和y′=-0.2的瞬时涡量场。从图13中可以看出，在流体-结构能量传递关系发生突变前后，两种频率比上的涡脱落模式差别很大。当fe/fs=0.93，即发生基本锁定时(图13(a)～图13(d))，涡脱落模式与单圆柱体卡门涡街非常相似，但是在近尾流处的涡结构更加复杂。(a)与(c)，(b)与(d)的涡量场近似于上下对称结构，由主圆柱和各附属小圆柱上形成的正负涡相互干涉，在近尾流区的上下两侧交替形成一个由3个正(负)涡围绕1个负(正)涡组成的涡团，然后该涡团逐渐变形融合，形成一个由多个正(负)涡组成的不规则涡团并脱落。在一个运动周期内，多根圆柱体的上下两侧对称的各脱落一个涡团，从而造成总升力的周期性振动，且振动主控频率等于结构振动频率。fe/fs=0.94时(图13(e)～图13(h))，涡脱落模式变得极不稳定，近尾流区上下两侧的涡脱落形态完全是不对称的，正负涡的干涉非常严重，涡团的结构和脱落频率在一个周期内不断地随机变化，无明显规律可循。结合之前的分析，fe/fs由0.93增加到0.94时，涡脱落模式由周期性的稳定脱落变为随机的不稳定脱落，能量传递关系由结构向流体传递变为流体向结构传递，正是这些水动力特性的突变造成了总流体力的突变。

4 结 论

[ 1 ] Bearman P W. Vortex shedding from oscillating bluff bodies [J]. Annual Review of Fluid Mechanics, 1984, 16:195-222.

[ 2 ] Williamson C H K, Govardhan R. Vortex-induced vibrations [J]. Annual Review of Fluid Mechanics, 2004. 36:413-55.

[ 3 ] Sarpkaya T. A critical review of the intrinsic nature of vortex-induced vibrations [J]. Journal of Fluids and Structures, 2004, 19:389-447.

[ 4 ] Gu W, Chyu C, Rockwell D. Timing of vortex formation from an oscillating cylinder [J]. Physics of Fluids, 1994, 6:3677-3682.

[ 5 ] Lu X Y, Dalton C. Calculation of the timing of vortex formation from an oscillating cylinder [J]. Journal of Fluids and Structures, 1996, 10:527-541.

[ 6 ] Guilmineau E, Queutey P. A numerical simulation of vortex shedding from an oscillating circular cylinder [J]. Journal of Fluids and Structures, 2002, 16(6):773-794.

[ 7 ] Blackburn H M, Henderson R D. A study of two-dimensional flow past an oscillating cylinder [J]. Journal of Fluid Mechanics, 1999, 38(5):255-286.

[ 8 ] 潘志远. 海洋立管涡激振动机理与预报方法研究 [D]. 上海: 上海交通大学, 2006.

[ 9 ] 樊娟娟, 唐友刚, 张若瑜, 等. 高雷诺数下圆柱绕流与大振幅比受迫振动的数值模拟 [J]. 水动力学研究与进展, 2012, 27(1): 24-32.

[10] 赵静, 吕林, 董国海, 等. 亚临界雷诺数下圆柱受迫振动的数值研究[J]. 计算力学学报, 2012, 29(1): 74-79.

[11] Wu H, Sun D P, Lu L, et al. Experimental investigation on the suppression of vortex-induced vibration of long flexible riser by multiple control rods [J]. Journal of Fluids and Structures, 2012, 30:115-132.

[12] 吕林. 海洋工程中小尺度物体的相关水动力数值计算 [D]. 大连: 大连理工大学, 2006.

[13] 唐国强. 立管涡激振动数值模拟方法及物理模型实验 [D]. 大连: 大连理工大学, 2011.

Numerical Simulation of Lateral Forced Vibration of Multiple Circular Cylinders at Low Reynolds Number

(State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Liaoning Dalian 116023, China)

Three steps finite element method is used to directly solve the N-S equations described by ALE method, and combining with ALE dynamic grid technique, a numerical model for 2D lateral forced vibration is developed. The lateral forced vibration of one main circular cylinder with additional four smaller circular cylinders is calculated, at Re=185, amplitude ratioAe/D=0.2andfrequencyratiofe/fs=0.86～1.04withintervalof0.01.Byanalyzingthetotalfluidforcesandtheirspectrumsasafunctionoffe/fs,anabruptchangeontotalliftforceisfound,andtheoscillatingmodeoffluidforcecanbedividedintofourtypesincludingtheprimarylock-inregime.Thenthefluid-structureenergytransfer,thephasebetweenstructuralmovementandliftforceandthevortexsheddingmodeareanalyzedtorealizethenatureonaccountforabruptchangeoffluidforce.

multiple circular cylinders; forced vibration; three steps finite element method

2014-07-01

国家自然科学基金项目(项目编号：51409041，51279027)。

吴 浩(1980-)，男，工程师。

1001-4500(2015)02-0048-09

P

A