0 引　言

立管作为将海底油田井口与油气生产勘探平台连接起来的细长结构物，被认为是海上平台系统的重要组成部分。而在一定流速的来流作用下，立管将在横流及顺流方向上同时出现涡激振动（vortex-induced vibration, VIV），引起结构疲劳损伤。

目前，涡激振动的主要研究方法分为实验方法、计算流体力学法（CFD）以及经验模型法。Trim^[1]记录并分析了大尺度立管在不同均匀流下的横流、顺流振动实验数据；Chaplin^[2]设计了阶梯流作用下立管的一系列涡激振动实验；Baarholm^[3]对实尺度立管模型在不同剪切流下立管横流和顺流向振动幅值及振动频率进行了分析；Song^[4-5]研究了均匀流下立管模型所受到流体激励力在顺流和横流方向的分布规律。数值模拟方面，Wang^[6]则采用CFD方法预报了均匀及线性剪切流下立管双向涡激振动并与实验进行了对比。而基于尾流振子模型方法的立管双向涡激振动响应预报工作也陆续展开^[7-8]。但对于常用于工业界的流体力分解模型，相关研究主要针对横流的振动^[9]，或在仅考虑横流响应对顺流向激励力放大作用进行顺流向振动预报^[10-12]。而近年的实验表明横流向与顺流向振动间存在明显的耦合效应^[13]。因此，有必要建立能够考虑横流-顺流向耦合效应的立管涡激振动预报模型，以更为准确地进行立管动力响应分析。

本文针对顶张式立管，基于双自由度受迫振动试验数据，采用有限元法和能量平衡方程建立了考虑双向耦合效应的涡激振动频域预报模型，通过与实验结果的对比验证了其合理性，并讨论附加质量系数变化效应对预报结果的影响规律。

1 数学模型 1.1 立管振动方程

将顶张式立管简化为受轴向张力作用的简支梁，坐标原点位于立管悬挂点，z轴向下为正，x轴为顺流方向，y轴为横流方向。将简支梁在空间上进行有限元离散，得到单元刚度矩阵 ${\boldsymbol K} = {{\boldsymbol K}_M} + {{\boldsymbol K}_T}$ ，其中：

$ {{\boldsymbol K}_M}{\text{ = }}\frac{{EI}}{{{l^3}}}\left[ \begin{aligned} {12}&{6l}&{ - 12}&{6l} \\ {6l}&{4{l^2}}&{ - 6l}&{2{l^2}} \\ { - 12}&{ - 6l}&{12}&{ - 6l} \\ {6l}&{2{l^2}}&{ - 6l}&{4{l^2}} \end{aligned} \right] ,$

(1)

$ {{\boldsymbol K}_T}{\text{ = }}\frac{T}{{30l}}\left[ \begin{aligned} {36}&{3l}&{ - 36}&{3l} \\ {3l}&{4{l^2}}&{ - 3l}&{{l^2}} \\ { - 36}&{ - 3l}&{36}&{3l} \\ {3l}&{{l^2}}&{3l}&{4{l^2}} \end{aligned} \right] 。$

(2)

梁单元集中质量矩阵可写为：

$ {\boldsymbol M}{\text{ = }}\frac{{({m_0} + {m_a})}}{{420}}\left[ \begin{aligned} {156}&{22l}&{54}&{ - 13l} \\ {22l}&{4{l^2}}&{13l}&{ - 3{l^2}} \\ {54}&{13l}&{156}&{ - 22l} \\ { - 13l}&{ - 3{l^2}}&{ - 22l}&{4{l^2}} \end{aligned} \right]。$

(3)

式中： $ {{\boldsymbol K}_M} $ 与 $ {{\boldsymbol K}_T} $ 分别立管单元抗弯刚度和轴向拉压刚度矩阵；EI为弯曲刚度；T为轴向张力；l为单元长度； $ {m_0} $ 为立管单位长度质量； $ {m_a}{\text{ = }}\dfrac{1}{2}{C_a}\rho {\text{π}} {D^2} $ ，为单位长度附加质量， $ {C_a} $ 为附加质量系数， $ \,\rho $ 为流体密度， $ D $ 为立管直径。

引入简支梁顶端和底端的约束条件，按照单元节点坐标组装得到总刚度矩阵 $ \left[ {{{\boldsymbol K}_M}{\text{ + }}{{\boldsymbol K}_T}} \right] $ 和质量矩阵 $ \left[ {\boldsymbol M} \right] $ 后，可得振动方程：

$ \left[ {{{\boldsymbol K}_M}{\text{ + }}{{\boldsymbol K}_T}} \right]\varphi - {\omega ^2}\left[{\boldsymbol M} \right]\varphi = 0 {\text{。}}$

(4)

式中： $ \omega $ 为圆频率， $f{\text{ = }}{{2\pi } \mathord{\left/ {\vphantom {{2{{\text{π}}} } \omega }} \right. } \omega }$ 即为固有频率； $ \varphi $ 为模态振型。

通过求解矩阵方程即可得到简支梁不同模态的频率和振型。

1.2 涡激振动附加质量系数模型

研究表明当立管发生涡激振动时，其顺流向响应频率往往为横流向的2倍^[13]。Zheng^[14]基于此开展了一系列圆柱体横流-顺流双自由度耦合受迫振动实验，得到了不同横流向振幅 ${{{A_y}} \mathord{\left/ {\vphantom {{{A_y}} D}} \right. } D}$ ，顺流向振幅 ${{{A_x}} \mathord{\left/ {\vphantom {{{A_x}} D}} \right. } D}$ ，顺流-横流间位移相位角 $ \theta $ 及无因次频率 $\hat f$ 情况下附加质量系数 ${C_{a,CF}}$ 、 ${C_{a,IL}}$ 的数据，其中 $\hat f{\text{ = }}\dfrac{{f \cdot D}}{U}$ ，f为振动频率，U为来流流速。而考虑到附加质量系数主要受无因次频率的影响，本文将文献[14]中附加质量系数简化为二维 $\hat f \sim {C_a}$ 曲线，分别如图1和图2所示。本文暂未考虑相位角的影响，并将其固定为0。

而考虑到不同工况中附加质量将随流场分布发生变化，立管响应频率及振型可通过迭代计算得到。

在横流方向上，计算步骤为：

1）假设 ${C_{a,CF}}{\text{ = 1}}$ ，求解立管固有频率 ${f_{CF,0}}$ ；

2）根据流场流速分布，计算立管不同位置无因次频率 ${\hat f_{CF,k}}{\text{(}}z{\text{) = }}\dfrac{{{f_{CF,k}} \cdot D}}{{U(z)}}$ ( $k \in N$ )；

3）根据横流方向附加质量系数简化曲线得到附加质量分布，更新振动方程，计算得到新的振动频率 ${f_{CF,k + 1}}$ ；

4）重复步骤2和步骤3直至结果收敛，得到横流向实际振动频率 ${f_{CF}}$ 和振型 ${\varphi _{CF}}(z)$ 。

由于基于顺流方向附加质量系数计算得到的振动频率往往并不满足总是横流响应频率的2倍，常出现 $2{f_{i,CF}}$ 落于 ${f_{j,IL}}$ ， ${f_{j{\text{ + }}1,IL}}$ 两者之间情况。因此需要对顺流附加质量进行调整。在本文模型中通过减小顺流向附加质量的方式增大固有频率，以满足 ${f_{j,IL}} = 2{f_{i,CF}}$ ，其过程如下：

1）针对第i阶横流向振动模态，假设其相对应顺流振动频率满足 ${f_{IL,0}} = 2{f_{i,CF}}$ ，并根据流场分布计算不同位置无因次频率 ${\hat f_{IL,0}}{\text{(}}z{\text{) = }}\dfrac{{2{f_{i,CF}} \cdot D}}{{U(z)}}$ ；

2）根据附加质量系数分布曲线得到顺流向附加质量分布，建立振动方程，计算第i阶横流向振动模态对应的顺流向模态及频率；

3）对比 ${f_{j,IL}}$ 和 ${f_{i,CF}}$ ，调整顺流向附加质量，直至满足 ${f_{j,IL}} = 2{f_{i,CF}}$ ，并得到实际顺流向振型 ${\varphi _{j,}}_{IL}(z)$ 。

1.3 涡激振动流体力模型 1.3.1 激励力模型

对于流体对立管的激励力，同样采用Zheng^[14]实验得到的横流向和顺流向涡激振动激励力系数。与单自由度强迫振动实验不同，Zheng所进行的实验可以反映出横流与顺流振动间存在的耦合效应，并在激励力系数的分布上有所体现。以0°位移相位角情况为例，典型激励力系数分布云图如图3和图4所示。可以发现，随顺流振幅增大，横流向激励力系数激励中心所对应无因次频率减小；而顺流向激励力系数将随横流振幅呈现增大趋势。参考文献[13]，将立管横流向无因次频率介于0.125和0.25之间的区域作为能量输入区。计算中根据无因次响应频率及幅值进行插值得到相应的横流及顺流向激励力系数。

1.3.2 阻尼模型

当立管结构响应超出能量输入区范围时，认为流体作用力将转变为水动力阻尼。本文使用由 Venugopal^[15]针对涡激振动提出的经验阻尼模型。

Venugopal阻尼模型可分为高流速区和低流速区，其中高流速区水动力阻尼系数为：

$ {r_{hh}} = {{\rho {U^2}(z){C_{vh}}} \mathord{\left/ {\vphantom {{\rho {U^2}(z){C_{vh}}} \omega }} \right. } \omega } ,$

(5)

低流速区水动力阻尼系数为：

$ {r_{hl}} = {r_{sw}}(z) + \rho DU(z){C_{vl}}。$

(6)

式中： $ {r_{sw}}(z){\text{ = }}\dfrac{{\omega {\text{ π }} \rho {D^2}}}{2}\left[ {2\sqrt {\frac{2}{{{{\omega {D^2}} \mathord{\left/ {\vphantom {{\omega {D^2}} \nu }} \right. } \nu }}}} {\text{ + }}{C_{sw}}{{\left( {\dfrac{A}{D}} \right)}^2}} \right] $ ，为静水阻尼系数； $ \nu $ 为流体运动粘性系数； $ {C_{sw}} $ 为静水阻尼经验系数； $ {C_{vh}} $ 和 $ {C_{vl}} $ 分别为高低流速区阻尼经验系数。

1.4 涡激振动响应幅值计算

参考SHEAR7方法，通过预估能量选取所振动激发模态^[7]。针对横流向振动模态i及对应顺流向振动模态j，计算预估能量：

$ {\varPi _{i\& j}} = {\varPi _{i,CF}} + {\varPi _{j,IL}} = \frac{{{{\left| {{Q_{i,CF}}} \right|}^2}}}{{2{R_{i,CF}}}} + \frac{{{{\left| {{Q_{j,IL}}} \right|}^2}}}{{2{R_{j,IL}}}}。$

(7)

其中：

$ Q{\text{ = }}\int_{Li} {\frac{1}{2}\rho {C_{V,\max }}(\hat f)D{U^2}(z)\left| {\varphi (z)} \right|{\text{dz}}}, $

(8)

$ R{\text{ = }}\int_L {{r_s}(z){\varphi ^2}(z)\omega {\rm{d}}z} + \int_{L - Li} {{r_h}(z){\varphi ^2}(z)\omega {\text{dz}}}{\text{。}} $

(9)

式中： $ {C_{V,\max }}(\hat f) $ 为无因次频率 $ \hat f $ 对应激励力系数最大值；Li为立管能量输入区长度。

选取激发振动模态时，根据最大预估能量值 ${\varPi _{\max }}$ ，计算各横流-顺流振动模态组合预估能量比 $\dfrac{{{\varPi _{i\& j}}}}{{{\varPi _{\max }}}}$ ，以0.7为能量阈值，选取横流及顺流向涡激振动激发模态^[7]。

当处于涡激振动状态时，在一段时间内立管的能量输入及输出可分别写为 $\left\langle {{\varPi ^{in}}} \right\rangle = \displaystyle\int_0^{nT} {P(t)V(t){\text{dt}}} $ 和 $\left\langle {{\varPi ^{out}}} \right\rangle = \displaystyle\int_0^{nT} {R(t){V^2}(t){\text{dt}}} $ 。

其中，模态激励力为：

$ \begin{split} & P(t) =\int_0^L {P(z,t)\left| {{\varphi _{CF/IL}}(z)} \right|{\text{dz}}} =\\ &\; \int_0^L {\frac{1}{2}\rho D{U^2}(z){C_V}[\hat f,Y(z),X(z)]\sin (\omega t)\left| {{\varphi _{CF/IL}}(z)} \right|{\text{dz}}} 。\end{split} $

(10)

式中： $ Y(z){\text{ = }}\dfrac{{{A_y}}}{D}\left| {{\varphi _{CF}}(z)} \right| $ 为立管各位置横流向振幅； $ X(z) = \dfrac{{{A_x}}}{D}\left| {{\varphi _{IL}}(z)} \right| $ 为立管各位置顺流向振幅。

立管模态速度为：

$ {V_i}(t) = {A_i}{\omega _i}\sin ({\omega _i}t), $

(11)

模态阻尼为：

$ R(t){\text{ = }}\int_0^L {\left[ {{r_s} + {r_h}(z)} \right]{\varphi ^2}(z){V_i}(t){\text{dz}}} {\text{。}}$

(12)

针对涡激振动响应中激发的横流向模态i及顺流向模态j建立能量平衡方程，可得模态振幅：

$ \begin{gathered} \frac{{{A_{x,i}}}}{D} = \hfill \\ \frac{{\dfrac{1}{2}\int_{Li} {\rho D{U^2}(z){C_{V,CF}}[{{\hat f}_{i,CF}},{Y_i}(z),{X_j}(z)]\left| {{\varphi _{i,CF}}(z)} \right|{\text{dz}}} }}{{\displaystyle\int_L {{r_s}\varphi _{i,CF}^2(z){\omega _{i,CF}}(t){\text{dz}}} + \displaystyle\int_{L - Li} {{r_h}(z)\varphi _{i,CF}^2(z){\omega _{i,CF}}(t){\text{dz}}} }} , \end{gathered} $

(13)

$ \begin{gathered} \frac{{{A_{y,j}}}}{D} = \hfill \\ \frac{{\dfrac{1}{2}\int_{Li} {\rho D{U^2}(z){C_{V,IL}}[{{\hat f}_{j,IL}},{Y_i}(z),{X_j}(z)]\left| {{\varphi _{j,IL}}(z)} \right|{\text{dz}}} }}{{\displaystyle\int_L {{r_s}\varphi _{j,IL}^2(z){\omega _{j,IL}}(t){\text{dz}}} + \displaystyle\int_{L - Li} {{r_h}(z)\varphi _{j,IL}^2(z){\omega _{j,IL}}(t){\text{dz}}} }} {\text{。}} \end{gathered} $

(14)

同时，对于柔性立管多阶模态同时被激发这一现象，本文通过模态叠加的方式进行计算处理。

综上，顶张式立管涡激振动响应幅值计算流程如图5所示。

2 涡激振动模型与试验结果对比验证

选取均匀流、阶梯流、剪切流3种流场工况进行计算，并与实验结果进行对比。为分析附加质量系数的影响，同时计算了附加质量系数固定为1情况下的涡激振动响应。

2.1 均匀流工况对比分析

采用Trim 等^[1]进行的大尺度立管模型均匀流试验。该试验立管模型长38 m，整体置于拖曳水池中，并在两端施加预张力。模型主要参数如表1所示。

图6与图7对预报所得涡激振动响应与试验数据进行了全面的对比。其中 ${C_a}{\text{ = 1}}$ 表示附加质量系数统一取值为1， ${C_a}(\hat f)$ 表示附加质量与无因次频率有关，附加质量沿立管管长取值不同。

在振动频率方面，与本文假设相同，试验得到顺流向主导振动频率为横流向的2倍。当来流流速低于1.0 m/s时，预报与试验结果间较为一致，而随着流速增大其间会产生一定的误差，但总体而言模型的预报误差处在可接受的范围内。在振幅方面， ${C_a}{\text{ = 1}}$ 与 ${C_a}(\hat f)$ 两种频域预报结果处于同一量值，横流向均方根位移平均值约为0.4D，与试验结果较为接近；而顺流向约为0.1D，小于实验数据。其原因在于，对于顺流向流体激励力而言，其存在2种能量输入区间^[11]：第一能量输入区对应折合速度约为1～2.5，而第二能量输入区对应折合速度约为2.5～4。由于强迫振动试验数据的限制，本文中的频域预报模型均未考虑第一能量输入区的激励力，导致预报结果较小。

选取流速为0.4 m/s工况作为典型工况对立管响应进行具体的分析，频域预报得到横流向无量纲均方根位移及顺流向无量纲均方根位移沿管长分布如图8所示。结果表明，频域模型在横流及顺流方向上均预测得到的与实验一致的第三阶和第五阶模态。其中在横流向振幅上，2种预报模型的结果较为一致。而在顺流方向上，由于 ${C_a}{\text{ = 1}}$ 情况下顺流振频未经过调整，导致计算中无因次频率偏离激励中心对应的无因次频率，令计算得到结果与实验间偏差较大。

2.2 阶梯流工况对比分析

Chaplin等^[2]设计了一种特殊的试验装置，立管模型下半部分在拖曳水池中，为阶梯流流速段，高于水池水面部分置于与拖车同时拖动的装水容器中，为静水段，并在立管顶端施加张力。试验参数如表2所示。

选取3种流速工况进行计算分析，流速段流速分别为0.16 m/s，0.31 m/s，0.54 m/s，对应立管顶张力为405 N，457 N，598 N。响应幅值分布如图9所示。计算结果表明，在该类型阶梯流作用下，立管呈现单一模态主导的响应特性，同时2种附加质量处理方法下预报得到的振型、振幅沿立管分布基本一致。与实验值相比，横流向响应吻合度较好；而顺流向响应除在流速0.31 m/s工况下在0.6L至底端区域内预报结果相对较小，其余均能够较好地吻合试验数据。

2.3 剪切流工况对比分析

MARINTE公司于挪威Hanoytangen海域建立了实尺度顶张式立管模型涡激振动试验装置。立管模型的一端与拖车相连，另一端靠浮箱牵引。当拖车运动时，立管将受到线型剪切流作用。试验参数如表3所示。

针对顶部最大流速为0.54 m/s的剪切流工况，模型预报所得激发模态、振动频率、模态振幅及权重如表4，表5所示。由表可见，在剪切流下立管多个振动模态将被激发。而由于附加质量系数取值的不同导致频率预报存在差异，2种方法在横流及顺流向涡激振动激发模态上有明显不同。

通过模态叠加得到立管涡激振动响应幅值分布如图10所示。Baarholm等^[3]给出了该工况下实验得到的涡激振动响应振幅沿立管分布。

在横流方向上，试验数据表明立管涡激振动含有明显的第11阶模态成分；而 ${C_a}{\text{ = 1}}$ 情况下并未预报出该模态，并且所得到振动幅值分布未能体现靠近底部位置振幅相对较大的分布特点。在顺流方向上，由于附加质量对振型影响较小，2种预报所得结果均未出现底部振幅较大的特点。同时由于激发模态的不同，两者在分布存在一定差异。

由于立管曲率与其弯曲应力及疲劳损伤密切相关，本文同样针对剪切流工况下横流向和顺流向曲率进行分析，如图11所示。横流方向上，在立管0～0.3L范围内，2种附加质量情况下计算所得曲率与试验数据间存在一定差异，而其余部分均与试验结果较为一致，其中 ${C_a}{\text{ = 1}}$ 情况所得结果较为保守。顺流方向上，预报结果则均略小于实验值。同时需要关注的是，在振幅分布中顺流向约为横流向的1/4～1/5，但是在弯曲曲率上两者处在同一量级上，并且由于其振动频率为横流向的2倍，顺流向响应同样会对立管造成不可忽视的疲劳损伤，在实际工程设计中应予以关注。

3 结　语

本文基于双自由强迫振动试验数据，通过有限元方法及能量平衡理论，建立了顶张式立管的双向涡激振动频域预报模型，与不同流场试验观测进行了对比，并讨论了附加质量系数对预报结果的影响。研究结果表明：

1）本文建立的双向耦合涡激振动频域预报方法可较为准确地得到立管在横流向及顺流向涡激振动响应频率及位移分布，计及附加质量的变化效应能够更好地反映出立管在不同流场下的涡激振动响应特性。

2）涡激振动引起的立管顺流向响应幅值虽小于横流向，但由于激发的模态阶数高、振动频率大，顺流向涡激振动引起的疲劳损伤同样应引起足够关注。