袁昱超， 薛鸿祥， 唐文勇

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

涡激振动(VIV)一直以来都是海洋细长结构物设计过程中面临的一个关键性问题，尤其对于顶张式立管的动力响应分析和疲劳寿命评估有着举足轻重的影响.由于其极强的非线性，涡激振动在流固耦合领域仍然是一个极具挑战性的研究课题.作为油气生产系统的重要组成部分，顶张式立管是连接海上浮式平台和海底井口的桥梁，扮演着将油气资源由海底输送至顶端平台的媒介.影响立管动力响应的结构刚度可分为由立管材料固有物理属性决定的弯曲刚度和由立管轴向有效张力提供的附加横向刚度2种成分，通常情况下，会在立管顶端施加一定的预紧力以提供必要的附加横向刚度.在复杂的海洋环境中，波浪与海流的联合作用易激发海上浮式平台的垂荡运动，平台的升沉带动张紧器的压缩或拉伸，从而导致立管轴向张力呈现随时间波动的变化特征.这不仅使得立管动力响应分析时需要引入动边界条件，更为本质的是结构刚度将随时间周期性变化.因此，计及顶端平台垂荡运动(即轴向张力时变效应)的顶张式立管涡激振动响应相较轴向张力恒定条件更接近于真实的海洋环境，也势必更为复杂.

试验研究和数值模拟是目前对涡激振动加深了解和认识的2类主要途径.近年来，随着相关试验技术的不断进步，Franzini等[1]针对小尺度立管模型计及顶端垂荡运动的涡激振动问题开展了相关试验研究.时域数值模拟同样是研究此类问题时一种可供选择的有效方法，且不易受到模型尺寸的限制.时域涡激振动数值模型根据其依托的核心算法可分为计算流体力学、尾流振子模型和流体力分解模型3类.Josefsson等[2]基于二维切片理论借助计算流体力学方法研究了结构弯曲刚度可变条件下圆柱体涡激振动响应.现有关于立管轴向张力时变(有时也被称为参数激励)条件下的涡激振动研究主要基于尾流振子模型[3-7].然而，已有的数值模拟研究鲜有与已公布的模型试验结果进行对比，其模拟时变张力工况的有效性及预报精度仍有待进一步论证.

流体力分解模型同样是被学术界和工程界广泛用于预报顶张式立管涡激振动的一类较为成熟的方法，Lie[8]早在1995年就已开展了相关研究工作.在近十年间，Sidarta 等[9]和Wang等[10]相继提出了各自的流体力分解模型用于开展立管在各种形式来流下涡激振动时域响应的相关研究.但上述流体力分解模型大多没有计及顶端平台垂荡运动对立管结构动力响应的影响，而是简单地认为轴向张力在整个涡激振动过程中保持恒定.本文推荐了一套可供选择的流体力分解模型，能够用于预报张力时变条件下顶张式立管涡激振动时域响应.涡激振动流体力载荷基于受迫振荡试验数据，结构刚度矩阵根据实时轴向张力沿管长方向的分布在每一分析步中进行更新.Franzini等[1]公布的模型试验被用于验证推荐数值模型在恒定及时变张力条件下的有效性.数值模拟与试验观测对比研究中发现，时变张力与涡激振动联合作用下立管结构响应呈现亚谐振成分、非对称振型及Mathieu型共振等特性，本文针对以上问题作进一步展开讨论，并给出了合理解释.

1 涡激振动时域模型

典型的平台-顶张式立管-海床系统可简化为如图1所示布置形式.图中：Cartesian坐标体系用以描述推荐涡激振动时域模型；x轴顺应来流方向；y轴垂直于流速方向；z轴为立管轴向即竖直向上方向.

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

1.1 涡激振动流体力载荷模型

基于涡激振动流体力分解模型，流体力载荷

Fhydro,y=Fv+Fd+Fm=

(1)

式中：Fv为激励力；Fd为阻尼力；Fm为惯性力；Cv为激励力系数；A为响应位移幅值；D为立管直径；f为响应频率；V为遭遇流速；t为时间；ρf为流体密度；cf为水动力阻尼系数；Ca为附加质量系数，本文假定Ca为经验值 1.0.

由式(1)可知，涡激振动流体力载荷系数Cv和cf的确定是流体力分解模型的核心内容，其在每一分析步内均同时受控于响应幅值和频率.图2所示为本文选取的激励力系数云图[11].图中：粗实线标注出对应Cv=0的能量输入及输出边界；当Cv<0时，激励力将遵循下文规定法则转化为水动力阻尼；当超出图2所示的试验数据范围时，同样认为激励力将依据下文提及的经验阻尼模型转化为水动力阻尼.上述受迫振荡试验工况测得的斯坦顿数(St)约为 0.193.

图2 涡激振动激励力系数库[11]Fig.2 VIV excitation force coefficient database [11]

计算流体力载荷时，立管的振动幅值和频率是最为关键的决定因素，因此，在每一分析步中，需提取各节点位移(y)、速度(v)和对应时刻(t)求解振动幅值和频率分别为

A=|yb-ya|/2,f=1/[2(tb-ta)]

(2)

式中：(ya,yb)和(ta,tb)分别为相邻两个瞬时速度为零(即v=0)点a和b对应的位移和时刻.

1.2 锁定判定准则

锁定(Lock-in)是涡激振动问题中尤其需要引起重视的响应特性，本文参照文献[12]中双自由度圆柱体受迫振荡试验研究成果，选取无因次频率带宽[0.125，0.25]作为锁定区间，上述锁定区间需要考虑真实环境及试验条件下St的差异进行修正.修正方法为

(fr/St)test=(fr/St)actual

(3)

式中：fr为无因次频率且有fr=fD/V.

若无因次频率落在涡激振动激发带宽内，锁定将发生.本文假定立管单元的无因次频率将被锁定到涡激振动激励中心对应的 0.17 (见图2)以确保激励力系数取到峰值，此处的 0.17 同样需要根据式(3)进行修正，此时对应的频率则作为涡激振动主导频率用于根据图2计算激励力系数.

水动力阻尼同样是涡激振动数值模拟领域较为重要的一个课题，当激励力系数为正时，激励力与立管速度同向，整个系统内能量由流体传入立管，而当激励力系数为负时，水动力阻尼力作用于立管，能量也将由立管传回流体.参考涡激振动频域方法中的处理手段，本文假定流体传入立管的等效能量在一个振动周期内传回流体，根据能量守恒原理，水动力阻尼系数计算公式可表示为

(4)

若无因次频率位于激发带宽外，则可根据图2求得Cv<0的激励力系数，进而基于式(4)得到对应的水动力阻尼系数.当立管响应幅值和频率超出图2所示试验数据范围，采用Venugopal提出的涡激振动经验阻尼模型[13]模拟水动力阻尼力，其中的经验系数Chf、Csw和Clf推荐选取 0.18、0.25 和 0.2.

1.3 结构动力响应求解

顶张式立管涡激振动响应通常遵循小位移假定，立管模型可认为是满足Euler-Bernoulli梁理论的抗弯弹性结构物，其在横流向的运动控制微分方程

(5)

式中：m为立管单位长度质量；cs为结构阻尼系数且有cs=2mωξ，ξ为结构阻尼比，ω为圆频率；E为弹性模量；I为惯性矩；T0为恒定轴向有效张力；ΔT、PT和φ分别为时变张力的幅值、周期和初始相位.

式(5)左侧的第3和4项分别对应由立管固有属性EI决定的弯曲刚度和由轴向有效张力提供的附加横向刚度.为了更好地说明本文的推荐方法，图3给出涡激振动时域分析中每一分析步的完整流程图.

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

本文推荐的涡激振动时域数值模型基于直接积分动力分析，采用Hilber-Hughes-Taylor (HHT)方法[14]求解式(5)中的结构运动控制方程.HHT法作为一种隐式数值求解方法已较为成功地应用于经典动力学问题的仿真模拟，由于其具有更好的收敛性、计算精度和效率，也被认为是一种改进的Newmark-β法.

2 数值模型试验验证

本文选取由文献[1]中公布的一组小尺度立管模型(未拉伸长度L0为2.552 m)计及顶端垂荡运动的涡激振动试验结果用以验证推荐数值模型预报恒定及时变张力工况的有效性.

2.1 试验描述

试验立管模型采用竖直放置形式，部分浸没于水中，模型顶部由伺服电机装置施加垂荡幅值At/L0=1%的简谐运动.沿立管跨长分布43个测量点，测量点位移通过水下相机借助光学追踪技术直接获取.顶端布置测压元件用以测量轴向张力.图4给出试验装置布置示意图.表征立管模型物理特性的相关参数如表1所示.

结合上述信息，可估算出该试验模拟的轴向张力时变工况下ΔT/T0约为 0.51.除一个恒定张力工况外，文献[1]中仅公布了2个垂荡频率(即ωT/ωCT=2,3)下试验观测结果.其中：ωT为轴向张力变化圆频率；ωCT为恒定张力条件下由涡激振动诱发的立管响应圆频率.3个工况对应的雷诺数(Re)介于[700, 6 100]，本文选取实际St近似为 0.21.

图4 模型试验装置布置图[1]Fig.4 Experimental facilities arrangement of the model test[1]

模型参数数值未拉伸长度/m2.552拉伸长度/m2.602浸没长度/m2.257外径/m0.022弯曲刚度/(N·m2)0.056质量比3.48浸没质量/kg2.05定常预张力/N40

2.2 固有频率对比

自由衰减试验表明立管模型n阶固有频率fN,n为 0.84nHz[1].本文基于FEM方法进行模态分析得到模型固有频率预报值，与实测值对比情况如图5所示.计算所得第1、2、3阶固有频率分别为 0.85、1.70 和 2.57 Hz，与试验数据吻合较好.而从第4阶开始，计算值略大于实测值.事实上，对于一端受拉的简支梁，在立管固有弯曲刚度和轴向张力提供的附加横向刚度共同作用下，结构各阶固有频率理论上应呈现预报结果所示的增幅逐阶递增规律，而自由衰减试验揭示的严格线性关系更像是仅考虑附加横向刚度所得结果.由于本文涉及的3个工况下，主激发模态最高为3阶，因此固有频率对比所示误差对下文结构动力响应研究影响甚微.

图5 立管模型固有频率Fig.5 Natural frequencies of the riser model

2.3 响应位移时历与幅值谱对比

3个计算工况下监测点(z/L0=0.43)处立管结构响应的预报结果与试验观测对比情况如图6所示.图6(a)表征恒定张力工况，图6(b)和(c)分别表征ωT/ωCT=2,3的时变张力工况.其中：VR,1为约化速度，可通过VR,1=V/(fN,1D)计算得到.图6(a)表明，均匀流条件下恒定张力工况中立管涡激振动响应展现出明显的单一模态激发特征，位移时历严格服从正弦变化规律，预报和实测响应位移幅值均大致稳定在 0.6D，预报值略大于实测值.预报所得主导频率同样与试验数据吻合，均接近于结构一阶固有频率，由此可知，约化速度VR,1=5.63 时，涡激振动稳定地激发该立管模型一阶固有模态.

图6 监测点处响应位移时历曲线及对应幅值频响谱Fig.6 Time histories of response displacement and corresponding amplitude spectra at monitor point

由图6(b)可知，变张力使得预报与实测结果中的响应位移幅值放大约1倍.试验观测中捕捉到位移振幅剧烈调制的现象且调制过程无明显的规律性或周期性，而对应的预报结果中虽然也存在一定的振幅调制现象，但就整个时间历程而言位移响应表现得更为稳定.尽管如此，预报所得位移幅值频响谱仍与实测结果保持较高的一致性，能量响应集中在ω/ωCT=1处并伴随着少量成分出现在ω/ωCT=2, 3, 4, 5等处，两者之间细微的区别在于实测结果中其他频率成分弱于预报结果所示，且如图中虚圆圈标注所示，预报结果中ω/ωCT=3,5的能量响应较之ω/ωCT=2,4更强.

由图6(c)可知，预报和实测位移幅值均在0.6D附近轻微波动，而对应的幅值频响谱中揭示出有别于前2个工况的新现象.预报和实测所得频响谱在ω/ωCT=1，3处均出现突出的能量响应峰值，且ω/ωCT=3处的能量峰值大于ω/ωCT=1处对应值，即立管动力响应不再主要展现一阶固有模态而变为由结构3阶固有模态主导.一些其他的频率成分同样出现在预报所得频响谱中，而实测频响谱似乎并未发现类似的亚谐振响应.Franzini等[15]指出，根据模型试验研究成果，在ωT/ωCT=3的纯顶端垂荡运动条件(即流速为0，无涡激振动响应)下发现了稳定的横向振动响应，但仅限于第3阶模态而对于第1、2阶模态并不适用.由此认为，在VR,1=5.63 这一约化速度下第3阶模态响应伴随着第1阶涡激振动同时发生，对于造成这一特殊现象的深层次原因，文献[15]中并未提及.上述事实也在一定程度上解释了为何ωT/ωCT=3的时变张力工况中体现的频响特征与前2个工况有所区别.

图7 响应位移幅值沿管长包络线Fig.7 Amplitude envelops of response displacement along the riser

2.4 响应幅值包络线对比

图7给出上述3个工况下预报和实测所得响应位移幅值沿管长包络线对比情况.图中：L为立管拉伸长度.对于恒定张力工况，预报所得响应位移幅值包络线与试验观测较为吻合，均为规则的1阶振型，仅在立管上半部分略大于实测值.两者振幅峰值均出现在z/L=0.42 附近，偏离立管中点.相比之下，在ωT/ωCT=2的时变张力工况中，幅值包络线被显著放大(1倍甚至更多).Franzini等[1]推断这一放大现象是由于结构发生了由Mathieu不稳定引起的参数共振.实测结果不再是标准的1阶模态振型，尤其是在z/L位于[0.25, 0.6]范围内，一些额外模态明显参与立管涡激振动.预报结果同样显示出高阶模态被激发的响应特征，在z/L=0.4 附近(即1阶振型极值出现位置)最大响应幅值略大于实测值，呈现1阶模态主导3阶模态参与的整体振型，与图6(b)中幅值频响谱所示频响特性对应.而对于ωT/ωCT=3的时变张力工况，预报和实测所得响应幅值包络线吻合较好，均展现出明显的3阶模态振型，预报值仅在z/L=[0, 0.3] 范围内略大于试验数据.预报与实测振幅峰值出现位置同样偏离标准3阶模态振型对应位置，且两者振幅包络线均呈现一定程度的非对称性.对比恒定张力工况，ωT/ωCT=3的时变张力对结构响应同样具有放大效应，只是不如ωT/ωCT=2条件下剧烈.

3 联合作用下结构响应特性研究

本节将就前文对比研究中发现的时变张力与涡激振动联合作用下结构动力响应展现的亚谐振响应、非对称振型及Mathieu型共振等特性展开讨论，探究其深层次原因.

3.1 亚谐振响应

由图6(b)中幅值频响谱所示，除涡激振动锁定频率外，ω/ωCT=3，5处同样存在着明显的能量响应，而这两处恰好对应ωCT+ωT和ωCT+2ωT.事实上，时变张力与涡激振动联合作用下ωCT±kωT的亚谐振频响成分的确可能被激发，现基于解析推导给出相应的理论解释.式(5)中的水动力阻尼主要起到限制立管振动幅值的作用，对于结构共振条件的推导影响较小，现忽略水动力阻尼项cf∂y/∂t并假定y(z,t)=φ(z)Y(t)对动力学方程进行解耦.其中：φ(z)为模态振型函数；Y(t)为幅值函数.对于两端简支梁可表示为sin(iπz/L)，式(5)可整理为

(6)

将Y写为级数展开形式

Y=Y0+εY1+ε2Y2…

(7)

将式(7)代入式(6)，可得到关于ε的各阶方程.出于简洁考虑，仅以第0阶和第1阶表达式为例进行推导：

通过求解式(8)和式(9)可得各自稳态解

Y0=A0cosωCTt+B0sinωCTt

(10)

(11)

式中：

关于ε的一阶方程的稳态解表明，额外的频率成分ωCT±ωT满足谐振条件.若求解关于ε的k阶方程，参照上述理论推导不难得出在涡激振动与时变张力联合作用下亚谐振频率成分ωCT±kωT的确可能被激发的结论.

3.2 非对称振型

由图7可知，3个工况(尤其是时变张力工况)下立管整体振型均呈现一定的非对称性，主要表现为振幅峰值出现位置偏离理论值且立管底端与顶端振型明显不对称.本文认为上述现象主要是由于该试验立管模型采用竖直放置形式，立管结构自重使得轴向张力定常成分沿管长(即垂向)分布不均匀，自上而下线性减小.因此，各分段对整体附加横向刚度矩阵的贡献也不相同.事实上，竖直放置立管的固有模态满足具有非正交解的Bessel类方程，模态振型峰值出现位置的确会向底端偏移.另外，T0(z) 沿管长不等，ΔT/T0自上而下逐渐变大，因此时变张力对立管底端的影响效应也会比顶端更为明显.Chen等[7]研究表明，随着ΔT/T0的增加，时变张力对结构响应的放大效应也越发显著.这也能够解释为何底端响应幅值大于顶端.文献[4]中采用尾流振子模型同样得出与本文预报结果趋势相同的非对称整体振型，一定程度上也佐证了本文数值模型的正确性.立管模型部分浸没水中导致的附加质量沿管长分布不均也会对结构整体振型产生一定影响，但由于该试验中浸没比达到87%，本文认为该影响效应与上述两点相比是次要的.

3.3 Mathieu型共振

图7(b)表明，ωT/ωCT=2的时变张力工况下的位移振幅包络线相较恒定张力工况放大1倍甚至更高，有理由认为此时的时变张力与已有的涡激振动协同作用下，立管结构发生强烈共振，出于结构安全性的考虑，共振状态是立管在工程设计阶段希望极力避免的，本文对于此类共振现象查阅相关文献给出内在机理解释.事实上，ωT=2ωCT的时变张力正是触发结构Mathieu失稳的最危险条件，文献[15]中试验研究表明仅ωT=2ωCT的时变张力即可导致结构发生参数共振.这种由于Mathieu不稳定造成的结构共振伴随着涡激振动锁定现象的发生被进一步放大，最终导致结构动力响应显著加剧.因此，本文认为此类时变张力与涡激振动联合作用下发生的共振现象可定义为Mathieu型共振.张力变化频率在ωT=2ωCT附近均可能诱发Mathieu型共振，故ωT/ωCT=3的时变张力同样对立管振动幅值具有一定的放大效应.

4 结论

计及顶端平台垂荡运动的顶张式立管涡激振动问题仍然是现阶段值得深入研究的前沿课题.本文基于涡激振动流体力分解模型提出一套可用于预报顶张式立管在海流及时变轴向张力联合作用下结构动力响应的时域数值方法.在对其恒定及时变张力条件下的有效性进行试验验证的基础上，详细讨论了对比研究中发现的若干响应特性，并得到以下结论：

(1) 基于本文数值方法的预报结果与试验观测吻合较好，推荐数值模型可有效预报轴向张力恒定及时变条件下顶张式立管涡激振动时域响应.

(2) 通过理论推导与数值模拟均可发现，当计及轴向张力时变效应时，除涡激振动锁定频率外，对应ωCT±kωT的亚谐振频率成分也可能被激发.

(3) 立管自重导致轴向张力定常成分沿管长分布不均匀，加之张力时变效应的影响，立管结构响应呈现非对称振型.

(4) 时变轴向张力与涡激振动联合作用下，结构会在张力变化频率ωT=2ωCT附近发生Mathieu型共振，响应位移幅值将明显放大.