深海立管涡激振动双稳态现象

李江涛， 吴志强， 王远岑， 张祥云

（天津大学a. 机械工程学院； b. 天津市非线性动力学与混沌控制重点实验室， 天津300350）

0 引 言

多稳态现象是非线性系统的重要动力学行为之一。 当系统中存在多稳态时，不同初始条件或扰动，可能导致系统行为趋于不同吸引子，最终振动特征往往有重要区别。 多稳态导致的这种振动行为的不确定性，会对工程结构的设计和运行带来严重危害。

管道涡激振动中结构动力与尾流水动力的相互作用非常复杂， 而且二者的耦合作用具有很强的非线性特征[1]。 其响应是否存在多稳态现象，若存在其有何特点，是一个值得关注的理论和工程问题。

已有小尺度实验研究表明， 流速增大或减小时管道振动存在迟滞现象， 在迟滞区有多解共存现象，即涡激振动多稳态行为确实存在。 Feng（1968）[2]首先在弹性支撑柱体风洞实验中发现折合速度约为6 时管道出现迟滞现象，其响应曲线同时包含初始分支（initial branch）和下端分支（lower branch）。Brika、Labeville 等[3]随后也在风洞中进行了相关的实验，选用了阻尼比与Feng（1968）[2]不同的空心管道，得到了更大的“锁频”区间和两倍于Feng 实验结果的滞后区域。 后来，Williamson[4-6]对中等雷诺数的低质量比弹性支撑圆柱体进行了一系列的实验，发现质量阻尼比m*ζ 不同时，涡激振动系统会出现两种不同的响应。 大m*ζ 时，系统的响应存在两个分支:初始激励分支（initial excitation branch）和下端分支（lower branch），最大振幅出现在初始激励分支上；小m*ζ 时，系统响应存在三个分支:初始分支（initial branch）、上端分支（upper branch）和下端分支（lower branch），最大振幅出现在上端分支上。系统响应在初始分支与上端分支、上端分支与下端分支之间都发生了跳跃现象，前者属于迟滞，后者属于间歇切换。

尽管涡激振动的多稳态现象已被实验证明，理论分析的研究却较少有相关的报道，焦点也多在预测模型上。 关于结构涡激振动响应的预测模型有多种，其中尾流振子模型因为其具有较为明确、合理的物理意义以及较好的计算精度，被工程界广泛采用。 该模型通过相互独立的方式分别建立圆柱振子运动方程和流体振子运动方程，然后利用它们共同预报流体-弹性系统的动力响应。Hartlen 和Currie[7]在Bishop 和Hasson[8]等人研究的基础上，首次提出了尾流动力可以使用Van der Pol 方程作为升力系数的控制方程，并且可以与结构的振动方程联立求解。后来Iwan[9]和Blevins[10]通过研究，给出了用于二维流场弹性支撑刚性圆柱的尾流振子模型。Facchinetti 和Langre[11-12]对前人的工作进行了总结，并且对比了位移、速度、加速度耦合尾流振子模型，认为加速度耦合模型可以更好地模拟结构对流体的作用，但是在预测锁频区域宽度时，仅适用于小质量比情况，随着质量比增大，简缩速度范围估计偏于保守。陈伟民和郑仲欣[1]在前人研究的基础上，结合并利用速度耦合模型和加速度耦合模型的优点，提出来一种新的非线性耦合模型，且与实验结果对照较为准确。 关于涡激振动中的多稳态现象，还有人做了一些相关的研究,陈威霖和及春宁等[13]研究了单圆柱涡激振动中的振幅不连续和相位切换现象，结果表明当阻流比小于0.05 后，锁定区间向更高的折合流速偏移，上端分支和下端分支之间的迟滞环宽度增大近4 倍，并解释了圆柱升力和位移之间相位差跳跃现象。

本文利用Van der Pol 尾流振子模型描述流体的作用， 选用非光滑平方非线性项描述流体对结构的阻尼作用，采用加速度耦合模型，建立深海立管的涡激振动动力学模型。应用Poincaré 映射方法，研究流速变化引起立管涡激振动行为的分岔现象，从而深入分析流速变化对涡激振动的影响。

1 深海立管涡激振动模型

1.1 立管振动模型

立管顶部通过张紧器与海上平台相连接，张紧器可以简化为一个刚度为k 的弹簧。 其作用主要有两点:一是可以为管道提供一个较大的静张力从而避免管道产生较大的弯曲；二是可以减小由于平台的垂直运动产生的纵向应力对立管的影响[14]；立管底端简化为铰支连接，模型示意如图1 所示。

图1 立管模型示意图Fig.1 Diagram of the riser model

基于Euler 梁的弯曲振动理论，立管横向振动的基本方程可表示为

式中:y （z, t ）为立管的横向位移， 是关于垂向坐标轴z 和时间t 的连续函数；EI 为立管的抗刚度；mp=πρp（D2-d2）/4 为单位长度管道质量，mf=πρfd2/4 为内部流体质量，ma=CaπρwD/4 为管道附加质量，其中ρp、 ρf、 ρw分别为立管材料密度、内部流体密度和海水密度；D、d 分别为立管的内外直径；Ca为附加质量系数，一般取1； f （z, t ）为立管单位长度上的作用力；T （z, t ）为立管的有效张力，其表达式为

1.2 尾流振子模型

尾流振子模型即不考虑具体的流场结构，而是将流体和其中的振荡物体视为一个整体系统，把尾流部分看成一个非线性振子，尾流振子的振动引起结构的振动；反过来，结构的振动又对尾流有一个反馈的作用。 本文采用Van der Pol 方程来表示旋涡脱落的振动特性，同时用惯性耦合来表达其与结构运动间的耦合作用。 具体形式如下[11-12]:

式中:η 为流体对结构的瞬时升力系数CL和结构静态横向升力系数CL0之比，η=2CL/CL0；Ωf为漩涡脱落的圆频率，Ωf=2πStU/D，其中U 为水流流速，St为Strouhal 数，本文取St=0.2；A、ε 为经试验测定的耦合系数，分别取12 和0.3。

流体对立管的作用力包含两部分: f= fL+ fL′。 其中涡激升力为

阻尼力为

式中CD为流体阻力系数，取1.2。

1.3 流固耦合模型

引入如下无量纲量将方程（1）和（3）进行无量纲化:

式中，M=mp+mf+ma并对其进行伽辽金离散，假设解有如下形式:

其中离散系统中各系数具体表达式如下所示:

表1 计算模型参数[15]Tab.1 Calculation parameter of model

2 流速变化对系统稳定性影响

因此，分别计算了有阻尼、无阻尼两种情况下系统特征值随流速变化的情况见图2。 无阻尼情况下，特征值随流速变化，纯虚特征值逐渐接近至相等，然后变成两对虚部相等的复重根，再变为两对纯虚根。 U∈（0.025-0.085 m/s）间属于耦合颤振，为平衡点不稳定，而其他区间平衡点是中心稳定的；有阻尼时，U 约为0.05 m/s 时，涡激频率与结构频率大小交换（临界值可从虚部颜色变化看出）。 由于始终有一对复根实部大于0，平衡点始终是不稳定的，即立管处于涡激振动状态。

图2 特征值随流速变化图Fig.2 Variation of eigenvalues with flow velocity

3 流速变化对系统行为的影响

复杂非线性系统动力学理论分析通常比较困难， 因而数值方法成为通常采用的手段， 其中Poincaré 映射方法得到较广泛的应用。 根据Poincaré 截面上点的分布情况，即可对系统响应类型进行判断。 当系统响应达到稳态以后，若Poincaré 截面上的点数目有限时，对应的响应是周期的；若这些点形成闭环，则对应的运动是概周期的；若这些点在有限区域形成一定形态的吸引子，则对应的运动可能是混沌的。

3.1 流速变化导致的分岔及双稳态现象

Poincaré 映射方法用于分析参数引起的系统响应行为变化时，可得Poincaré 截面上点随参数的变化图，称为Poincaré 分岔图。 如果采用参数延拓法，同时获得参数递增和递减的两个分岔图，则可从是否存在滞后/跳跃现象来判断是否存在多稳态现象。 用Poincaré 方法分析系统分岔时，Poincaré 截面选择有重要影响。 为得到振动幅度随参数的变化图，本文选用截面定义为

本节利用Poincaré 映射方法分析水流流速在0.01-0.20 m/s 之间变化时立管振动行为的变化。 图3 给出了流速递增、递减两种情况下的Poincaré 分岔图计算结果，而竖直方向的箭头表示解发生跳跃的方向。（其中“实线箭头→” 表示正向，“虚线箭头 ”表示反向）。从图3 可以看出，增大流速和减小流速时算得的立管振动响应基本一致，立管的振动响应随着流速的增大呈现出先增大后减小的趋势，并且出现了涡激振动特有的“锁频”现象。当折合速度小于3 时，立管振动响应较小，此时涡激振动对管道的影响较小，尾流和管道之间的相互影响也较小；当折合速度处于4～11 之间时，管道出现较大的振动响应，振幅均大于0.6D，最大可达到1.2D。 此时管道振动频率接近漩涡脱落频率并出现了“锁频”现象；随着流速继续增大，立管振幅也随之减小。 在“锁频”区域两侧，出现了因解跳跃导致的滞后现象。 这与Feng （1968）[2]和Williamson（1996-1999）[4-6]等学者在实验中发现的现象类似，在滞后区间内有双稳态现象。

图3 增大流速和减小流速时立管振动响应对照图Fig.3 Comparison of riser vibration response with increasing flow velocity and decreasing flow velocity

特别需指出锁频解存在范围U∈（0.024-0.086 m/s）与前文得到的耦合颤振范围U∈（0.025-0.085 m/s）较为接近，也就是说可用线性分析得到的耦合颤振区来估计涡激共振导致的锁频区。这一结论与Langre（2006）[18]所得结论一致。

为了进一步说明这类行为的特点，又分别计算了典型流速下立管振动响应的时间历程、相图和频谱。 在锁频区左侧，当折合速度为3.02 时，增大和减小流速情况下得到的时间历程、相图、频谱如图4所示。 从两图中可以很明显地看出在增大流速时，其时间历程呈现出明显的“拍频”现象，属概周期运动，振动幅值最大值也仅为0.08D，频谱中含两种不同的频率成分；而在减小流速的情况下，其时间历程呈现出较为规整的正弦图像，响应振幅达到0.4D，而频谱中仅含一种频率成分，接近立管线性固有频率，属周期运动。说明当折合速度为2.77-3.02 区间内存在周期运动、概周期运动共存的双稳态现象。

图4 增大和减小流速时立管的时间历程图、相图和谱图（折合速度3.02）Fig.4 The history, phase diagram and spectrum of riser with increasing and decreasing flow velocity at a reduced speed of 3.02

在锁频区域右侧，当折合速度为11.06 时，增大和减小流速情况下的结果如图5 所示。在该折合速度下，响应的不同主要体现在振动幅值大小和频谱中主要频率的不同。 当增大流速时，系统的振动幅值可以达到0.93D 左右，频谱分析得到的主要频率为0.186 rad/s；而当减小流速时，系统的振动幅值只有0.17D 左右，频谱分析得到的主要频率为0.249 rad/s；显然，当折合速度为10.81-11.32 区间内双稳态现象是两类不同周期运动的共存。

图5 增大和减小流速时立管的时间历程图、相图和谱图（折合速度为11.06）Fig.5 The history, phase diagram and spectrum of riser with increasing and decreasing flow velocity at a reduced speed of 11.06

通过具体对各流速下立管振动响应的分析， 可以更清楚地看出立管在这些流速下的振动响应均产生了双稳态现象。 其差异既有振动形态的不同，也包含振动幅值、频率的差别。

3.2 响应频率和相位差变化规律

为了进一步研究立管响应频率在选取流速范围内变化规律， 本文分别计算绘制了立管响应频谱和尾流振子频谱随流速变化的三维图形，如图6-7 所示。

图6 1 000 m 立管振动响应三维频谱图Fig.6 Three-dimensional spectrum of 1 000 m riser’s vibration response

根据谱分析结果提取了主要频率成分随流速变化结果见图7， 为反应响应的本质图中还给出了有阻尼情况下特征值频率的结果。

图6-7 清晰地反映了其振动响应频率随流速变化的规律:

（1） 当流速较小时（0.010-0.022 m/s 之间），频谱中有两个主要频率成分: 较小的频率非常接近漩涡脱落频率， 随着流速增加而增加且幅值较小；较大的频率在0.15 rad/s 附近波动，基本和管道固有频率 （1 000 m 管道的固有频率为0.154 rad/s）相同，且随流速变化不明显，其与线性化矩阵的有阻尼特征值结果吻合较好。 也就是说在小流速时，立管的振动为概周期运动，其频率成分包括:涡激频率和立管频率，此时相当于漩涡脱落激励立管振动。 因流体阻尼作用小，流体与立管相互作用较弱，响应中保留了立管的频率成分；

图7 响应频率变化Fig.7 Variation of response frequency

（2） 当流速处在较大范围时（0.024-0.086 m/s 之间），此流速范围立管振动处在“锁频”区域，振幅较大，响应频率成分以基频为主，其频率值由U=0.024 m/s 时的0.148 rad/s 增大到了U=0.086 m/s 时的0.183 rad/s，略小于线性化矩阵的有阻尼特征值结果。 也就是说在中等流速下，立管的振动为周期运动，漩涡脱落频率与立管频率接近而发生“锁频”现象， 即发生共振；

（3） 随着流速的增大（0.088-0.1 m/s 之间），此流速范围管道振动脱离“锁频”区域，立管振幅幅值较小。 其频谱图像又发生了进一步的变化，主要频率从0.088 m/s 时的0.249 rad/s 跃变为0.1 m/s 时的0.302 rad/s，得到的频率接近于漩涡脱落频率Ωf=2πStU/D（Ωf=0.340-0.387 rad/s），处于线性化矩阵有阻尼特征值之间，接近较大的特征值。 也就是说在小幅周期振动区间，即在高流速下，由于系统的耦合性，立管振动响应的主要频率接近漩涡脱落频率，并且立管频率部分被衰减掉。

同时又提取了立管和尾流振子频谱分析中主频率对应的相位差， 并分析了其值随流速变化的规律，如图8 所示。 从图中可以看出，在选取流速范围内， 二者振动相位之间的关系是随着流速逐渐变化的。 当流速为0.03-0.048 m/s 之间时，二者振动相位是同相的；当流速为0.048-0.096 m/s 之间时，二者振动相位是相反的。当v=0.096 m/s 时，立管位移-尾流振子位移的相位差达到了180°，二者相位完全相反。

图8 立管和尾流振子相位差随流速变化图Fig.8 Phase difference of riser and wake oscillator with flow velocity

4 结 论

本文利用Van der Pol 尾流振子模型描述流体的作用，建立了深海立管的涡激振动动力学模型，利用Poincaré 映射方法，研究了流速变化引起立管振动行为的分岔现象，并具体分析了振动响应频率和相位的变化规律，主要结论如下:

（1） 无阻尼情况下，系统特征值随流速变化，纯虚特征值逐渐接近至相等，然后变成两对虚部相等的复重根，再变为两对纯虚根。 U∈（0.025-0.085 m/s）间属于耦合颤振平衡点不稳定，而其他区间平衡点是中心稳定的；有阻尼时，由于始终有一对复根实部大于0，平衡点始终是不稳定的，即立管处于涡激振动状态。

（2） 通过本文的计算和研究表明:1000m 长管道系统响应存在三个区域，小流速下概周期运动，中流速下大幅锁频周期运动，高流速下小幅周期运动。 在锁频区两侧，分别存在两种类型的双稳态现象:概周期运动与大幅周期运动共存，大幅与小幅周期运动共存。

（3） 通过对系统响应频率进行计算分析可知，系统处于小流速时，振动为概周期运动，频率成分包括:涡激频率和立管频率，相当于涡脱落激励立管振动；中等流速下，振动为大幅周期运动，涡脱落频率与立管频率接近而发生锁频现象，即发生共振；高流速下，振动为小幅周期运动，仅含涡脱落频率，立管频率部分被衰减掉。

（4） 在选取流速范围内，立管和尾流振子振动相位之间的关系是随着流速逐渐变化的，由小流速下的同相位变成了大流速下相反的相位。