考虑轴向张力时变效应的圆柱体涡激振动响应特性研究

2019-05-27袁昱超薛鸿祥唐文勇

振动与冲击 2019年9期

袁昱超，薛鸿祥，唐文勇

(1.上海交通大学海洋工程国家重点实验室，上海 200240;2.高新船舶与深海开发装备协同创新中心，上海 200240)

立管作为油气生产系统的重要组成部分，起到将油气资源由海底输送至顶端平台的作用，典型的平台-顶张式立管-海床系统可简化为图1所示布置形式。影响立管动力响应的结构刚度可分为由立管固有属性决定的弯曲刚度和由轴向张力提供的附加刚度两种成分。在复杂波浪环境中，浮式平台易发生垂荡运动，诱发张紧器压缩或拉伸，从而导致立管顶端张力随时间波动。因此，考虑张力时变效应的细长圆柱体涡激振动相较张力恒定条件更接近于实际海洋环境。

近年来，Franzini等[1-2]开展了张力简谐变化的小尺度立管模型试验研究。Karniadakis等[3-4]基于二维切片理论借助CFD方法研究了结构弯曲刚度可变条件下圆柱体涡激振动问题。Park等[5]认为由时变张力引发的参数激励会改变结构原有响应特征。Da Silveira等[6]发现张力时变条件下立管动力响应存在模态阶跃。王东耀等[7]、Wu等[8]和Chen等[9]均指出平台垂荡引起的张力波动可能导致立管激发更高阶模态及更大幅度动力响应。唐友刚等[10]研究了时变张力对于剪切流工况下TTR立管涡激振动响应的影响效应并得出结构响应存在对应0.5倍参激频率的亚谐振成分。

图1 平台-顶张式立管-海床系统示意图Fig.1 Sketch of platform-TTR-seabed system

流体力分解模型是广泛用于预报圆柱体涡激振动的一类方法，详见文献[11-14]。以往流体力分解模型并未考虑时变张力对结构动力响应的影响。本文基于Wang等提出的时域流体力分解模型，在每一分析步之初更新结构刚度矩阵以计及张力时变效应。采用Franzini等公布的小尺度模型试验结果验证本文方法的有效性。结合另一38 m细长圆柱体模型，研究了28组工况下结构涡激振动响应特性，讨论了时变张力的幅值和频率对响应特性的影响规律。

1 涡激振动时域模型

本文选取Cartesian坐标系，其中x轴顺应来流速度方向，y轴垂直于流速方向，z轴为圆柱体轴向。

1.1 涡激振动水动力载荷模型

基于流体力分解模型，涡激振动横流向流体力载荷Fhydro,y可由式(1)表示，主要包含三个载荷成分，即激励力Fv，阻尼力Fd和惯性力Fm。

(1)

式中：Cv为激励力系数，A为响应位移幅值，D为直径，f为响应频率，V为遭遇流速，t为时间，ρf为流体密度，cf为水动力阻尼系数，Ca为附加质量系数，本文设定Ca为1.0。

图2为激励力系数Cv云图，源于圆柱体受迫振荡试验数据[15]。该受迫振荡试验测得的St数约为0.193。

图2 涡激振动激励力系数云图Fig.2 VIV excitation force coefficient contour

1.2 锁定判定准则

本文选取无因次频率带宽[0.125，0.25]作为锁定区间，该锁定区间需根据式(2)对真实环境及试验条件下St数的差异进行修正。

(fD/VSt)test=(fD/VSt)actual

(2)

若响应频率落在涡激振动激发带宽内，锁定发生，本文假定结构将被锁定到涡激振动激励中心对应的无因次频率0.17，并根据式(2)进行修正，作为涡激振动主导激励频率。

当Cv为负时，水动力阻尼力作用于圆柱体，水动力阻尼系数可由式(3)计算。

(3)

当响应幅值和频率超出试验数据范围时，采用Venugopal提出的阻尼模型[16]，其中的经验系数Chf，Csw及Clf推荐选取0.18，0.25及0.2。

1.3 结构动力响应求解

基于Euler-Bernoulli梁理论，横流涡激振动动力学方程可由式(4)表示。

(4)

式中：m为结构单位长度质量，cs为结构阻尼系数cs=4πmfξ，ξ为结构阻尼比，E为弹性模量，I为惯性矩，T0为张力恒定成分，ΔT，PT和φ分别为时变张力的幅值、周期和初始相位。

图3给出本文数值分析方法流程图。

图3 涡激振动时域分析流程图Fig.3 Flow-chart of VIV time domain analysis

2 数值模型的试验结果验证

Franzini等针对长度L0=2.552 m，直径D=0.022 m的立管模型开展了试验研究，本节通过对比该试验结果，验证本文数值方法的有效性。试验立管采用竖直放置形式，顶端施加40 N预张力，浸没水中长度为2.257 m。模型顶部施加垂荡幅值At/L0=1%的简谐运动激励。弯曲刚度为0.056 N·m2，质量比和浸没重量分别为3.48和7.88 N/m。自由衰减试验表明，该立管模型n阶固有频率fN,n为0.84nHz。本文采用有限元方法所得前三阶固有频率分别为0.85 Hz，1.70 Hz和2.57 Hz，与实测值吻合较好。

Franzini等给出了1个张力恒定工况和2个张力时变工况(ωT/ωCT=2和3)的试验结果，其中ωT为时变张力圆频率，ωCT为恒定张力时涡激振动响应圆频率。本节对以上3组工况的试验结果与数值结果进行对比，如图4和图5所示。

2.1 张力恒定工况

由图4可知，当张力恒定时，涡激振动呈现出明显的单一模态激发特征。预报所得位移幅值沿管长包络线与试验观测总体较为吻合，仅在立管上半部分略大于实测值。两者振幅峰值均出现在z/L=0.42附近，偏离立管中点。

监测点(z/L0=0.43)处位移时历服从正弦变化规律，预报和实测位移幅值均约为0.6D，幅值频响谱与试验结果吻合，该工况仅单一频率被激发，主导频率接近于结构一阶固有频率。由此可知，约化速度VR,1=5.63时，涡激振动稳定地激发该立管模型一阶固有模态。

图4 张力恒定工况结构响应对比Fig.4 Comparisons of structural response for constant tension case

2.2 张力时变工况

如图5所示，张力时变工况下，立管涡激振动呈现明显有别于张力恒定工况的响应特性，且受ωT/ωCT影响显著。在ωT/ωCT=2工况中，幅值包络线被显著放大(约1倍)。Franzini等推断这一放大现象是由于结构发生了由Mathieu不稳定引起的参数共振。实测结果不再是标准的一阶模态振型，尤其是在z/L位于[0.25,0.6]范围内，一些高阶模态明显参与立管涡激振动。预报结果同样显示出高阶模态被激发的响应特征，最大响应幅值略大于实测值。而对于ωT/ωCT=3工况，预报和实测响应幅值包络线吻合较好，均展现出明显的三阶模态振型。预报值仅在z/L介于[0,0.3]内略大于试验结果。预报与实测振幅峰值出现位置均偏离标准三阶模态振型对应位置，呈现明显的不对称性。

对于竖直放置的立管模型，由于结构自重存在，张力沿管长方向由顶端向底端线性递减，结构固有模态为满足类Bessel函数的非正交解集。因此，即使在均匀流条件下，结构整体振型在理论上也会呈现不对称性。针对张力沿轴向分布不均匀的工况，Da Silveira等借助自编程序及Orcaflex软件均得到与图5所示趋势一致的数值结果，即结构振型具有行波特性。

就监测点(z/L0=0.43)处立管局部结构响应而言，ωT/ωCT=2的时变张力激励使得预报与实测位移幅值相较张力恒定工况放大约一倍。试验观测存在位移振幅剧烈调制的现象且调制过程无明显的规律性。对应的预报结果中虽然也存在一定的振幅调制，但就整个时间历程而言位移响应表现得更为稳定。预报所得位移幅值频响谱与实测结果保持较高一致性，能量响应集中在ω/ωCT=1处并伴随着少量成分出现在ω/ωCT=2,3,4,5等处。实测结果中其他频率成分弱于预报结果所示。而对于ωT/ωCT=3工况，预报和实测响应位移时历同样呈现明显的振幅调制特性，幅值在0.6D附近轻微波动。预报和实测频响谱在ω/ωCT=1和3处均出现突出的能量响应峰值，且ω/ωCT=3处的能量峰值大于ω/ωCT=1处。ω/ωCT=2,4,5,6的频率成分同时出现在预报和实测频响谱中，只是实测频谱中响应能量相对微弱。针对该试验立管，Franzini等[17]在纯时变张力(ωT/ωCT=3)激励的试验工况中发现稳定的三阶横向振动响应，而无第一、二阶响应成分，即这一约化速度下第三阶模态响应伴随着第一阶涡激振动同时发生。这也在一定程度上解释了为何该工况体现的频响特征与下文大尺度圆柱体模型设计工况有所区别。

(a)At/L0=1%,ωT/ωCT=2,VR,1=5.8

(b)At/L0=1%,ωT/ωCT=3,VR,1=5.63

通过与模型试验对比显示，本文采用的数值预报方法可较有效地用于时变张力工况下圆柱体涡激振动响应分析。

3 时变轴向张力条件下涡激振动响应特性研究

与大多数真实海洋环境下细长结构物相比，小尺度立管具有一定的局限性，最显著的一点在于涡激振动仅激发结构一阶模态。为深入研究时变张力下涡激振动响应特性，有必要挑选另一尺度较大的圆柱体模型，使得在低流速条件下也可激发结构高阶模态。本文旨在研究时变张力这单一要素对涡激振动的影响，故选取一水平放置圆柱体模型作为研究对象。

Trim等[18]公布了Norwegian Deepwater Programme (NDP)采用的38 m长圆柱体模型，模型横向放置在试验水池中，表1给出模型的主要参数。

3.1 计算工况

38 m细长圆柱体模型在不同端部预紧力下前几阶固有频率分布如图6所示，可以发现预紧力对固有频率具有显著的影响，其本质原因在于预紧力直接影响结构刚度。

涡激振动发生时，结构将被锁定到距离主导频率最近的结构固有频率。而主导频率仅与St数、圆柱体直径和来流速度有关。因此，即便在不同预紧力下涡激振动主导频率均保持恒定，但激发模态阶数却未必不变。考虑张力时变效应时，可根据变化幅值的相对大小将工况分为两类。变化幅值较小时激发模态阶数保持不变，而变化幅值较大时激发模态则有可能随时间发生阶跃。

本文选取ΔT=500 N和4 000 N分别代表张力变幅较小和较大条件，14个时变张力频率涵盖高频区及低频区，时变张力初始相位取φ=0。设计的28组计算工况参数见表2，ΔT/T0和ωT/ωCT分别为张力无因次变化幅值和频率。流速均为0.4 m/s。

3.2 时变与恒定轴向张力结构响应对比

本节选取ΔT/T0=0.1，ωT/ωCT=1/20工况，对张力时变及恒定条件下涡激振动响应特性进行分析。采用连续小波变换(CWT)方法获取圆柱体中点处(z=19 m)响应位移时频分布云图，母小波选取Morlet小波。位移幅值频响谱根据快速傅里叶变换(FFT)得到。张力时变及恒定条件下圆柱体中点响应对比如图7所示，其中前三行分别为端部张力、响应位移和响应频率时历，第四行为位移幅值频响谱。

图7 时变与恒定张力下结构响应对比(z=19 m)Fig.7 Response comparison with time-varying and constant tension (z=19 m)

由图7可知，张力恒定时，涡激振动具有明显的单模态激发特征，响应幅值和频率在时间维度上保持不变。位移幅值频响谱显示，单一被激发的频率成分为2.14 Hz。根据文献[18]，与之对应的试验值为2.134 Hz，两者吻合度较高。

当张力时变时，涡激振动幅值随时间呈现规律性变化，其变化周期基本等于时变张力周期(即10 s)。伴随着张力的变化，能量将由外界传入涡激振动系统。根据能量守恒原理，响应位移趋向于变大以消散输入的能量，而输入能量大致与时变张力变化速度的平方成正比，故位移幅值变化应遵循张力变化速度绝对值的变化趋势。因此，在一个变化周期内，响应位移按照类似于“山峰”的形状调制两次，对应张力变化的两次加速和减速过程。如图中浅蓝色标记所示，最小幅值出现的时刻点并非严格对应于张力变化速度VT=0，而是略滞后于后者，即存在迟滞现象。

对于张力时变工况，响应频率同样展现出时变特性，随张力的增大而变高，反之亦然。由图6可知，端部预紧力位于4 500～5 500 N时，激发模态阶数始终保持第三阶，而结构固有频率与张力存在正相关关系，因此响应频率的变化与张力变化契合。根据位移幅值频响谱，时变张力下涡激振动为多频叠加响应，其中能量最强的峰频仍为2.14 Hz，对应恒定张力工况下涡激振动主导频率，同时激发了一些其他频率成分，激发区间为[1.8，2.4]Hz，且主导频率强度较张力恒定工况明显变弱。

3.3 时变张力幅值效应

张力时变效应对涡激振动响应的影响由变化幅值和频率共同决定，本节挑选ωT/ωCT=1和1/50代表时变张力频率较高和较低两种情况，研究时变张力幅值效应，圆柱体中点处结构响应对比如图8和图9所示，左右两列分别对应ΔT/T0=0.1和0.8。

图8 ωT/ωCT=1时结构响应对比(z=19 m)Fig.8 Response comparison with ωT/ωCT=1 (z=19 m)

图9 ωT/ωCT=1/50时结构响应对比(z=19 m)Fig.9 Response comparison with ωT/ωCT=1/50 (z=19 m)

张力高频变化条件下，ΔT/T0=0.1时，时频分析及快速傅里叶变换结果均显示结构动力响应特性与恒定张力工况大致相似，不同之处在于峰频2.14 Hz的强度变弱，且额外出现0.14 Hz被激发的现象；而ΔT/T0=0.8时，结构响应变得相对复杂，位移时历曲线呈现不稳定性，响应频率也随着时间发生明显的变化，能量最强的频率成分为1.99 Hz，该频率值非涡激振动主导频率2.14 Hz，位移幅值显著放大，超过恒定张力工况下对应值的两倍，达到0.055 m，这一响应放大现象将在下文中着重讨论。

张力低频变化条件下，随着ΔT/T0的增加，主激发频率(2.14 Hz)的强度变弱，频率激发带宽相应变宽。由于更多亚谐振频率成分参与涡激振动，ΔT/T0增大时响应位移和频率时历呈现出更为复杂的规律，结构响应具有明显的周期性，遵循的周期与时变张力周期一致。

本节涉及的四个计算工况对应的响应位移时空分布如图10所示，第一、二行分别对应ωT/ωCT=1和1/50，左右两列对应ΔT/T0=0.1及0.8。张力高频变化条件下，ΔT/T0较小时，结构响应几乎与恒定张力工况一致；而当ΔT/T0较大时，整体振动形态仍大致满足三阶振型，但响应位移时空分布不具有明显的规律。大幅值与小幅值振动交替出现，整个动力响应过程变得不再稳定。张力低频变化条件下，变化幅值效应呈现出新的特性。ΔT/T0较小时，结构在张力达到极值时被激发出有别于恒定张力工况的非标准三阶振型；而ΔT/T0较大时，结构响应出现直观的模态阶跃现象。张力高频变化条件下，由于周期太短，结构无法在不同模态间形成完整的转换过程，故模态阶跃并未出现。

图10 不同ΔT/T0下响应位移时空分布Fig.10 Temporal-spatial distributions of displacement with different ΔT/T0

3.4 时变张力频率效应

时变张力频率对圆柱体动力响应的影响同样显著。本节选取ΔT/T0=0.1和0.8两组工况分别讨论张力变幅较小和较大时变化频率对涡激振动响应的影响。结构中点处动力响应对比情况如图11和图12所示，由左至右分别对应ωT/ωCT=2，1/4和1/60。

小ΔT/T0条件下，当ωT/ωCT=2(即PT=0.25 s)时，结构响应并不遵循0.25 s的变化周期，与恒定张力工况相比，振幅调制现象明显，位移显著变大，极值达到0.03 m；当ωT/ωCT=1/4和1/60时，动力响应的周期性逐渐显现，峰频2.14 Hz始终占据主导，其他频率成分同时被激发，这些被额外激发的亚谐振频率成分的具体数值严格遵循ωCT±kωT(k=1,2,3,……)的规律。

大ΔT/T0条件下，当ωT/ωCT=2时，响应位移幅值急剧增大，响应能量集中在1.99 Hz，而非一直保持的2.14 Hz，可认为此工况已诱发时变张力与涡激振动的结构联合共振，详见下文；当ωT/ωCT=1/60时，周期性重新出现，结构动力响应回归稳定状态；而对于ωT/ωCT=1/4工况，响应位移和频率均显得较为无序，位移幅值较小，各激发频率间的强弱关系也并不固定，可认为其处于高频共振区与低频稳定区之间的过渡区域。ΔT/T0较大时，同样存在着ωCT±kωT的频率成分可能被激发的情况。由图10可知，时变张力频率对涡激振动响应的另一影响效应体现在高频变化抑制模态阶跃的发生。

图11 ΔT/T0=0.1时结构响应对比(z=19 m)Fig.11 Response comparison with ΔT/T0=0.1 (z=19 m)

图12 ΔT/T0=0.8时结构响应对比(z=19 m)Fig.12 Response comparison with ΔT/T0=0.8 (z=19 m)

4 联合激励共振讨论

时变张力和涡激振动联合激励下，结构可能发生共振，响应位移幅值将显著增大。出于结构安全考虑，共振状态是海洋细长结构物在工程设计阶段希望极力避免的，因此本节将就这一可能出现的共振现象展开讨论。

图13给出各工况下最大位移分布，图中蓝色实线标注出恒定张力工况对应值。当ωT/ωCT大于1/2时(高频区)，最大响应位移普遍大于恒定张力工况，且极值出现在ωT/ωCT=2处，与文献[2]中试验结果一致。由此可推断出ωT=2ωCT是时变张力与涡激振动联合激励下结构最易发生共振的核心频率。而ωT=2ωCT正是触发结构Mathieu失稳的最危险条件，这种由于Mathieu失稳造成的结构共振伴随着涡激振动锁定现象的发生被进一步放大。因此，本文认为此类联合激励下发生的共振可定义为Mathieu型VIV共振。当ωT/ωCT处于高频区某一固定值时，响应位移的放大效应随着ΔT/T0的增大而越发明显，这一规律可以基于能量守恒原理给出相应解释。当ωT相同时，时变张力向涡激振动系统输入的能量与ΔT/T0正相关，故ΔT/T0较大时，响应位移趋向于增大更多以消除更多输入的能量。图中A/D>20的结果超出普遍认知，表明在极端工况条件下Mathieu型共振对结构响应的放大效应极为显著。文献[9]基于尾流振子模型研究此类极端工况时也得出与图13所示数量级相当的预报结果。

图13 不同工况下响应位移最大值Fig.13 Maximum displacements under different cases

随着ωT由共振中心2ωCT向两侧偏离，最大位移也将逐渐减小并趋向于张力恒定工况对应值。ωT/ωCT小于1/2时的结果在图13左侧的红色方框内被局部放大。由于ΔT/T0较大时，高频共振区与低频稳定区之间存在着结构响应较为混乱的过渡区域，故ωT/ωCT位于[0.05，0.25]时，ΔT/T0=0.8时的最大响应位移小于ΔT/T0=0.1对应值。如图8和图12所示，Mathieu型共振发生时，主导结构响应的激发频率不再为涡激振动锁定频率2.14 Hz，而是偏移至其附近的一个新值，约为2 Hz，因此可认为Mathieu型共振诱发了失谐现象。