细长输流管内外流耦合振动特性研究

鲍 健，陈正寿

(1.浙江海洋大学 船舶与海运学院，浙江 舟山 316022；2.中国海洋大学 工程学院，山东 青岛 266100)

随着海洋资源的开采逐渐走向更深的海域，大长径比的深海立管/管线应用越来越广泛。深海输流管线主要运输高温高压油气，沿水深的不均匀流动容易激发细长弹性管的多模态振动响应，另外管内流体高速输运也会改变弹性管结构的动力特性[1-2]，内外流耦合作用对管线整体运动响应及疲劳损伤影响显著。内外流与细长柔性管线之间发生流固耦合的机理是不同的，外流以一定的速度流过管线会诱发尾涡脱落，进而引起弹性管体涡激振动(VIV)。当尾涡脱落频率与管体结构固有频率接近甚至相等时，“锁定”现象发生，从而增强管线振动并加速管线的疲劳破坏[3]。当内流沿管线轴向运动时，附加质量力、离心力和科氏力诱发管体流致振动(FIV)。这种双流固耦合工况(外流—管体，内流—管体)会导致管线动力学特性相较单流固耦合工况更加复杂。

在柔性立管/管线涡激振动方面，国内外学者已经做了相当多探索性的研究工作。Chaplin等[4]通过阶梯流下柔性圆柱体振动响应试验，获得了沿轴向和横流向(CF)的一系列激发波数。Trim等[5]对均匀剪切流动下的VIV响应进行了试验研究，发现顺流向的疲劳损伤程度与横流向相同。Wang和Xiao[6]基于大涡模拟湍流模型对竖直立管在均匀和线性剪切流下的涡激振动过程进行数值研究，得到的结果与试验数据吻合良好。当轴向内流沿非线性几何形状的弯曲管体运动时，额外的离心力和科氏力使管线强非线性内外流耦合响应过程更加复杂，考虑内流的弹性管涡激振动问题在近十几年来才得到了明显的关注。Guo和Lou[7]通过弹性管内流的模型试验，发现随着内流流速的增加，弹性管的振动幅度将增大，而振动频率将减小。Duan等[8]通过数值研究发现，行波在输流管线的横流向和顺流向(IL)VIV响应中都很明显。Yang等[9]通过参数化的数值研究发现，立管的最大位移和应力可以根据内部流速而增加或减小，并且临界内流流速会导致不同外流速度下立管的主振模态阶数增加。Chen和Kim[10]通过对弹性管在上行、下行内流和剪切外流下振动响应的数值仿真，发现内流的输送方向在决定大跨度管体的动力响应方面起着不同的作用。柳博瀚等[11]基于计算流体力学(CFD)方法开展管体内流对弹性管振动响应影响的研究，结果表明较高的内流流速能够显著影响管体的振动特性。

目前关于内外流耦合下弹性管涡激振动的研究主要关注内流流速的变化对立管振动频率、振幅，主振模态等特性的影响，如要全面评估弹性管在实际工作环境的动力响应和疲劳损伤，则需对内外流耦合诱发的多模态共振，行、驻波等特性开展更深入研究。文中首先介绍采用的数值方法，然后描述数值计算模型，并验证数值方法的可靠性，最后基于各工况的数值结果，分析内外流耦合下细长弹性管的动力学特性，并阐述管内涡结构的演变规律。

1 数值方法

1.1 流场控制方程

采用多物理场耦合分析软件STAR-CCM+[12]，研究内外流耦合下柔性管的强流固耦合特性。大涡模拟(LES)及其动态的Smagorinsky SGS模型[13]被用于求解流体域的水动力。过滤处理后的Navier-Stokes方程为：

(1)

(2)

采用有限体积法对流场进行离散化求解，压力速度耦合采用SIMPLE算法，动量方程的离散采用有界中心差分格式。

1.2 结构动力学控制方程

将柔性管体简化为两端简支的欧拉—伯努利梁模型[14]，可表示为：

(3)

式中：EI表示管体弯曲刚度；T为立管顶张力，m为管体系统的质量，c为结构阻尼，Fi为立管所受的水动力载荷分量。通过有限元法将式(3)离散化，再采用纽马克—贝塔法求解。

(4)

通过流场域和固体域之间的数据映射与交换实现了双向流固耦合，其中径向基函数方法[15]被用于流场网格边界和内部节点的移动插值。

2 数值仿真模型

2.1 计算域及网格划分

选取Chen等[16]开展的剪切流下玻璃纤维增强塑料(FRP)管涡激振动试验作为计算参考模型，并展开进一步的数值模拟研究。FRP弹性管模型主要参数如表1所示。

表1 FRP管模型的主要参数

图1为该FRP管数值模型示意，沿管体A端(最小剪切流速，Vmin=0)至B端(最大剪切流速，Vmax)施加线性剪切外流，弹性管模型A端固定，B端铰接，并且在B端施加一个恒定预张力，其中管内均匀流自A端运输至B端。如图2所示，外流场计算域网格划分采用多面体网格和定向网格技术[12]，管体近壁面流体网格划分较密，保证了y+<1，计算域远场的网格较为稀疏，以适应网格的大变形。外流沿x轴正向流动，流动入口距立管横截面中心15D，流动出口距立管横截面中心50D，上下边界分别距立管横截面中心15D(D为FRP管外径)。边界条件设置为：外流入口为速度入口边界，来流方向与边界垂直，外流出口设置为参考压力为0的压力出口边界，与立管外表面重合的流体域表面设置为无滑移壁面边界，前后两个侧面设置为对称边界，上下两个面设置为速度入口边界，速度大小方向与入口边界一致。在管体结构的有限元网格划分方面，立管横截面为环形结构，沿周长均匀划分40个网格节点，沿径向多层分布，并沿轴向以适当的层间距均匀层状拉伸。

图1 FRP管数值模型示意

图2 外流场网格划分

2.2 数值方法可靠性验证

通过傅里叶功率谱识别主频的方法，对采用不同度量单位的模型试验和数值模拟数据进行对比验证。图3分别给出了2种剪切流速下的FRP管跨中处振动位移功率密度谱对比结果。图3中，线1和线2分别代表模型试验中测得的横流向和顺流向应变时间历程功率谱密度，其对应的纵坐标为功率谱密度(模型试验)，横坐标为频率单位Hz。线3和线4分别表示在数值仿真中监测到的横流向和顺流向振动位移时间历程功率谱密度分布，其对应的纵坐标为功率谱密度(数值仿真)。横坐标为频率，单位Hz。

图3 试验值和数值模拟结果的频谱对比

从以上2个频谱对比分析可以发现，外流流速对确定主振模态的频率起重要作用。当最大剪切流速度较低时(Vmax=0.3 m/s)，在数值模拟和试验中观察到的CF振动过程功率谱密度函数几乎是单峰分布。当剪切流速增加到Vmax=0.9 m/s时，一些从振模态会被激发。可以得出结论，在数值模拟中观察到的主频分布与试验结果相符。由此证明了文中采用的流固耦合数值方法的可靠性。

其次，使用大涡模拟(LES)及其动态的Smagorinsky SGS模型，进行了雷诺数Re=3 900的单圆柱绕流数值模拟。采用端木玉等[18]的计算模型，圆柱直径D=0.01 m。计算域长度(x轴的方向)为15D，上游断面距离圆柱中心为5D，下游断面距离圆柱中心为15D，方向与来流方向一致；计算域宽度(z轴方向)为πD，方向沿着圆柱方向；计算域高度(y轴方向)为10D，上下面距离圆柱为5D，方向平行于圆柱横截面。采用多面体非结构化网格对计算域进行划分，为了满足大涡模拟对近壁面网格密度的要求，保证了近壁面处y+<1。从表2可知，文中计算的结果与前人计算结果比较接近，验证了数值模型的有效性。

表2 单圆柱模型验证

2.3 网格依赖性

网格密度是影响数值结果准确性的关键因素，因此在做进一步的数值预报前首先进行了圆柱绕流问题的流体网格依赖性验证。这里所用网格为多面体非结构化网格，并对柱体近壁区网格局部加密。对于来流速度U=1.5 m/s下刚性圆柱(D=34.8 mm，L=πD)绕流问题，采用了大涡模拟湍流模型，并考虑了沿柱体内壁周长方向节点的数量Nc，无量纲参数Δr/D(Δr是第一层网格高度)和无量纲参数Δz/D(Δz是轴向网格层间距)3种影响网格密度的关键参数，计算并比较了不同网格密度下(Mesh-1，Mesh-2，Mesh-3，Mesh-4)的平均阻力系数，均方根阻力系数和均方根升力系数，如表3所示。结果表明：Mesh-3和Mesh-4的均方根升力系数的误差值减小到5.4%，相应的平均阻力系数和均方根阻力系数随网格密度的增加几乎没有变化。可以得出结论：Mesh-3(其特征为Nc=160，Δr/D=0.000 5，Δz/D=0.005)适用于给定管体模型流体域的数值模拟。因此，在接下来的数值模拟工作中，采用Mesh-3中使用的外流场网格密度。使用STAR-CCM+[12]，在32个处理器(intel(R)Xeon(R)CPU E5-2620 v3 @ 2.40GHz)上，每一种工况耗时超过15×24 h。

表3 网格依赖性测试结果

3 计算结果及分析

3.1 振动强度的时间—空间分布

随着管体振型及振动轨迹发生改变，沿管体轴线的振动能量也发生变化，通过对比在不同内流流速工况下管体振动位移在时间空间上分布云图来分析内流对管体展向模态干涉和波形转变的影响。外流Vmax=1.5 m/s不变，内流分别为0.5 m/s,1.0 m/s,1.5 m/s,2.0 m/s的4个工况下管体顺流向与横流向的位移—时程空间分布如图4～7所示，横坐标为标准化时间，纵坐标为标准化的沿管体轴向各监测点的位置。其中，深色区域表示该瞬时振动位置偏离初始位置较大。

图4 内流Vin=0.5 m/s时，无量纲振动位移时程空间分布

随着内流流速由0.5 m/s增长到2.0 m/s，管体横流向振动表现为6阶主振模态，顺流向振动表现为9阶主振模态。这表明了在一定的内流流速范围内，外流流速对柔性管横流向和顺流向振型起决定性作用。另外，沿管跨方向振动位移的大小不但随时间变化，而且与内流流速密切相关。当内流流速从0.5 m/s增加到1.0 m/s，振动强度提升，但是当内流流速由1.0 m/s增加到2.0 m/s，振动强度反而有所下降。这些都归因于流体—固体和模态—模态间的连续能量传递。一个有趣的现象是，沿管跨方向局部的位移极值现象很明显，表现为深色区域。如图5(a)所示，在t/T=0.35～0.55的时间段中，存在一个位移极大值区域，在t/T=0.55～0.75的时间段中存在一个位移极小值区域。这意味着虽然这个工况下管体横流向振动表现为6阶振型，但是1阶模态的振动强度很显著。

图5 内流Vin=1.0 m/s 时，无量纲振动位移时程空间分布

另外在顺流向上，位移极值区域沿管长交替出现，如图4(b)至图7(b)。这表明了2阶甚至3阶模态开始参与振动能量的转移，多模态共存及跃迁现象显著。这种内流—管体和外流—管体的双重流固耦合系统使管体的振动响应更加复杂。

图7 内流Vin=2.0 m/s 时，无量纲振动位移时程空间分布

内流不仅可以改变管体振动强度，还能影响管体展向的驻波行波特性[20]。弹性FRP管的涡激振动过程表现为混合驻波行波主导，并且在选取时间范围内，驻波行波交替出现。驻波模式比行波模式更占优势，后者在管跨局部位置间歇性出现，尤其是在横流向的振动中。这一发现表明由于管体长度有限且激发的振动模态不够高，大部分振动能量仍集中在驻波区域。如图4所示，在外流Vmax=1.5 m/s，内流0.5 m/s的工况下，管体横流向和顺流向都表现出驻波主导的混合振动模式，表现为在管体的两端主要由驻波主导，且局部的行波组分在管中段间歇性出现。在t/T=0.65～0.85的时间段内，管体振动位移峰值位置不断波动，并且振型从6阶转化到了5阶。这种相邻模态间转变的出现体现出多模态间能量的传递。图5为在外流Vmax=1.5 m/s，内流1.0 m/s的工况下管体顺流向与横流向的位移时程空间分布。相比于内流0.5 m/s的工况，横流向振动表现出更多的行波组分和更加难以辨识的位移峰谷值位置，这说明此时内流流速的增加使得横流向振动变得更加强烈。但是在顺流向振动上却发现驻波主导模式更加显著。图6为在外流Vmax=1.5 m/s，内流1.5 m/s的工况下管体顺流向与横流向的位移时程空间分布。由图6可知，此工况下管体的振动情况与内流为0.5 m/s的工况相似，横流向和顺流向振动都表现出驻波主导，局部行波间歇性出现的模式。在t/T=0.50～0.60的时间段内，横流向振型由6阶转化成了5阶，多模态转变现象间歇性出现。图7为在外流Vmax=1.5 m/s，内流2.0 m/s的工况下管体顺流向与横流向的位移时程空间分布。由图7可知，此工况下管体横流向和顺流向的振动情况均为驻波主导，虽然存在较少的行波成分，但不足以导致模态阶数的转变。总的来说，管体的内流流速增加可以明显改变管体的振型和振幅，并且使其振动响应特性更加复杂化。

图6 内流Vin=1.5 m/s 时，无量纲振动位移时程空间分布

3.2 顺流向的平均偏移及最大变形处的振动轨迹

外流作用下，管体受到拖曳力作用，在顺流向上会产生比较大的静态变形，同时又因为交替变更的升力作用和瞬变的拖曳力分量，产生比较显著的局部振动变形。柔性管体结构的非线性振动必将影响管内流的运输路径。管内流路径的变化，又必然导致沿管体的科氏力、离心力、附加质量力发生变化，使管体的非线性振动响应更加复杂化。为了更全面立体地揭示大变形输流管体的振动过程，这里将研究内流流速的变化对于弹性管最大偏移位置及相应运动轨迹的影响。图8为外流为1.0 m/s，内流分别为0.0 m/s,0.5 m/s，1.0 m/s，1.5 m/s和2.0 m/s的5个工况下管体顺流向的平均偏移，由图可知，管体顺流向最大平均偏移量在5D～6D之间，最大偏移位置为y/L=0.375。内流流速的增大对管体顺流向最大平均偏移量有一定的影响。图9选择性地展示了外流为1.0 m/s，内流分别为0.5 m/s和1.5 m/s时最大偏移位置y/L=0.375处的管体运动轨迹，其中横坐标为管体顺流向上各点的相对位置，纵坐标为管体横流向上各点的相对位置。由于在开始一段时间内两端约束的柔性管体受到水流的冲击会产生一个较大的顺流向偏移，图中数据只截取管体产生稳定静态偏移后的运动轨迹。由图9可知，每个周期的“8”字形轨迹发生了变形，且相邻周期的运动轨迹可能会产生一个较大的位置跳跃，体现了大长径比管体在较大外流作用下强烈振动的过程。由于各个周期内的“8”字形轨迹变形程度不一，且振动位置跳跃较大，所以多个运动周期叠加之后的轨迹规律较不清晰。这也是细长输流管在内外流共同作用下多模态交涉的体现。

图8 外流Vmax=1.0 m/s时，不同内流速度Vin下顺流向的平均偏移

图9 y/L=0.375处管体运动轨迹

图10为外流为1.5 m/s，内流分别为0.0 m/s,0.5 m/s，1.0 m/s，1.5 m/s和2.0 m/s的5个工况下管体顺流向的平均偏移，与图8对比可知，在较大的外流作用下，管体顺流向最大平均偏移量高达13D～15D，且内流流速对最大平均偏移量的影响更加显著，但是最大偏移位置y/L=0.375几乎没有发生变化。图11选择性地展示了外流1.5 m/s，内流0.5 m/s和1.5 m/s时最大偏移位置y/L=0.375处的管体运动轨迹。由图11可知，管体最大偏移位置的运动轨迹在多个运动周期叠加下也变得混乱，反映了管体的振动能量在内外流共同作用下变化显著。这种显著的振动能量交换现象也与3.1小结中管体多模态共存和振动模态跃迁息息相关。

图10 外流Vmax=1.5 m/s时，不同内流速度下顺流向的平均偏移

图11 y/L=0.375处管体运动轨迹

3.3 二次流的涡结构演变

由于扭曲的运输路径弯曲和不断振动的管体，内流在柔性管体中运输，并随着管体不断振动的弯曲和旋转等因素及其相互交涉使管体内部产生复杂的二次流效应[21]。单独提取管体内流场数据，对管跨中截面流向涡的时程变化以及沿管长不同截面处的流向涡结构进行比较分析。从图12可知，在t=0.5～3.0 s的时间范围内，流向涡一开始附着在管壁上，随着管体的持续振荡，流向涡开始脱落，并在内部区域以涡对或单涡结构形式存在并运动。随着时间的推进，涡的强度、涡对数量以及涡对位置都在发生变换，这对柔性管体的涡激振动起到了显著的作用。最后可以看到大量离散的涡结构分布在各个位置，这些涡可能是由本身从管壁处脱落，也可能是管体的高频振动导致涡的分离。对于自由振动的柔性管体，横流向、顺流向以及轴向的弯曲和振动都会对内流场结构产生显著影响。

图12 管中跨截面内部涡结构的时间演变历程

4 结 语

基于CFD方法，对内外流作用下的柔性管体涡激振动流固耦合特性开展数值模拟研究。首先参考不含内流的涡激振动试验工况开展验证计算，数值结果与试验结果吻合良好，验证了数值方法的可靠性。为研究不同内流流速下大变形柔性管体的涡激振动响应，将不同工况计算得出的标准化振动位移的时—空分布、顺流向平均偏移、振动轨迹及内流场横流向涡的演化进行分析，得出了以下结论：

1)外流速度的增大导致更高的振动模态在管体横流向和顺流向响应中激发出来。随着内流流速的增加，主振模态阶数不变，但振动强度会发生明显改变，这表明了内流的运动带来的向心力、科氏力和附加质量力使整个管体系统的振动更加复杂。

2)输流管体的振动表现出多模态响应特性，除了主振模态，其他潜在的固有模态也被激发出来，并具备一定的振动能量。驻波主导了管体两个正交方向的振动，行波成分易于在顺流向上局部出现，这些都体现出管体涡激振动的模态转变和跃迁。内流流速对顺流向最大偏移位置有一定的影响，在多模态响应的影响下，该位置的运动轨迹有明显的变形和跳跃。

3)观察到内流场的二次流效应，其中流向涡沿着管壁不断生成并脱落，形成对涡或者单涡。随着持续性的局部振动变形，大量离散的涡被发现，这是涡激振动对管体内流场结构的影响。

文中主要研究了有限长弹性管模型在内外流作用下的振动响应问题，为一定长径比弹性管体的流固耦合分析提供了研究手段，分析结果可为较大长径比的弹性管涡激振动预报提供参考，但受限于当前数值计算能力，应用于工程实际中大长径比弹性管涡激振动计算模拟尚有一定难度。