钢悬链线立管成本低、对浮体运动及高温高压环境的适应性强，被认为是最适合应用于深水油气开发的海洋立管；但其在顶部浮体和海洋环境的作用下会产生大幅度和大角度的运动，具有较强的几何非线性^[1]。传统的求解方法如线性弹簧法、悬链线方程法、准静力法等不能较好处理立管在大变形情况下的非线性动力学问题。细长柔性杆模型用弹性杆的轴线位置来描述杆的位置, 可以较好地解决立管的大变形问题^[2]。Nordgren^[3]用有限差分近似和半显式积分给出了有限振幅运动问题的细长柔性杆模型分析方法；Garrett^[4]用有限单元法求解了不可伸缩的细长柔性杆模型问题，并进行了静力分析和动力分析；Paulling等^[5]考虑杆的轴向伸长，进一步发展了柔性杆模型。以上研究都没有考虑内流的作用，目前尚未见到用细长柔性杆模型模拟海洋立管时考虑内流影响的研究报道。

近年来，内流对立管动力特性的影响已经引起国内外学者的重视。Wu等^[6]采用奇异摄动技术研究内流对顶张力立管动力特性的影响，发现内流会减弱顶部张力的作用，但当顶部张力较大时，内流的影响很小，当顶部张力较小时，内流的影响变得较为明显。Paidoussis^[7]研究了输流管道的动力学和稳定性问题，发现内部流体高速定常流动时，立管会发生弯曲和振颤。郭海燕等^[8-9]利用有限单元法研究顶张力立管在内流、海流和波浪作用下的侧向变形和动力特性，发现内流流速的增加会增大立管的变形和应力，同时会降低立管的固有频率，但当内流流速较小时，其影响较小。目前对于考虑内流作用的海洋立管动力特性的研究多集中于顶张力立管这种刚性立管，对钢悬链线立管这种大变形的柔性立管考虑较少。Vandiver等^[10]将钢悬链线立管等效成顶张力立管，顶张力立管的长度即钢悬链线立管的弧长，将钢悬链线立管的轴向张力线性化，并沿管长施加，然后用有限单元法及简化方法推导固有频率和振型，这种等效转换的方法比较繁琐，且存在一定的误差。Orcina公司^[11]编制的商业软件OrcaFlex可以分析钢悬链线立管在内流作用下的动力特性，但该软件对立管模型的建立采用了集中质量法。本文考虑内流的作用, 利用柔性杆模型建立钢悬链线立管的动力学模型，采用有限元法分析内流对钢悬链线立管动力特性的影响。

1 计算模型 1.1 柔性杆模型

柔性杆模型如图 1所示，杆的位形由杆轴线位置表示，即空间曲线r(s, t)。

考虑内流作用的微元体平衡方程可由动量守恒获得：

$ F'+q=\left( \rho +{{m}_{f}} \right)\ddot{r}。$

(1)

式中：F为轴线上的合力；q为立管单位长度施加的外力；ρ和m_f分别为立管和内流单位长度的质量。

忽略转动惯量和剪切变形的影响，由动量矩守恒得：

$ M'+r'\times F+m=0。$

(2)

式中：M为轴线上的合力矩；m为单位长度的外力矩。通常立管、锚线等细长结构物可以忽略扭矩和分布外力矩，则M=r′×EIr″, m=0。EI为抗弯刚度。

如果柔性杆可伸长且伸长量为小量，则变形条件：

$ \frac{1}{2}\left( r'\cdot r'-1 \right)\approx \frac{\lambda }{EA}。$

(3)

式中：λ=F·r′+(EIr″)′·r′；EA为抗拉刚度。则杆的运动方程为：

$ -\left( EIr'' \right)''+\left( \lambda r' \right)'+q=(\rho +{{m}_{f}})\ddot{r}。$

(4)

1.2 荷载分析

海洋立管处于复杂的海洋环境中，受到自身重力和周围流体的作用。

$ q=w+{{F}^{s}}+{{F}^{d}}。$

(5)

式中：w、F^s、F^d分别为立管单位长度受到的重力、水静力和水动力。

$ {{F}^{s}}=B+\left( Pr' \right)'。$

(6)

式中：B为立管单位长度受到的浮力；P为管内外压强差。

水动力利用Morison方程计算^[12]：

$ \begin{align} &{{F}^{d}}={{\rho }_{w}}\frac{\pi {{D}^{2}}}{4}\dot{V}+{{C}_{a}}{{\rho }_{w}}\frac{\pi {{D}^{2}}}{4}({{{\dot{V}}}^{n}}-{{{\ddot{r}}}^{n}})+ \\ &\frac{1}{2}{{C}_{D}}{{\rho }_{w}}D\left| {{V}^{n}}-{{{\dot{r}}}^{n}} \right|\left( {{V}^{n}}-{{{\dot{r}}}^{n}} \right)。\\ \end{align} $

(7)

式中：ρ_w为海水密度；D为立管外径；C_a为附加质量系数；C_D为拖曳力系数；V为水流速度。

钢悬链立管内部通过高温高压的油气时，气液两相段塞流的模拟非常复杂。为了简单起见，本文采用细长体活塞流模型来模拟内流作用，即假设立管内部流体为一具有无限柔性的细长状柱体，其断面上每一点都具有相同的流速v，这样内流对立管微元体的作用力为^[13]：

$ f\left( s, t \right)=-{{m}_{f}}\left( \frac{{{\partial }^{2}}r}{\partial {{t}^{2}}}+2v\frac{{{\partial }^{2}}r}{\partial s\partial t}+v\frac{^{2}{{\partial }^{2}}r}{\partial {{s}^{2}}} \right)。$

(8)

式中：第一项为惯性力；第二项为科氏力；第三项为内流保持与立管相同的曲率引起的作用力。

1.3 有限元离散

代入荷载作用力，将运动方程和变形条件写成张量形式:

$ \begin{align} & -\left( \rho \text{ }+\text{ }{{m}_{f~}} \right){{{\ddot{r}}}_{i~}}-{{C}_{A~}}\ddot{r}_{i}^{n}-2{{m}_{f~}}v\dot{r}{{'}_{i~}}-\left( EIr{{''}_{i~}} \right)''+ \\ & \left( \bar{\lambda }r{{'}_{i}} \right)'+{{{\bar{\omega }}}_{i}}+{{{\bar{F}}}^{d}}_{i}=0. \\ \end{align} $

(9)

$ \frac{1}{2}\left( r{{'}_{n}}r{{'}_{n}}-1 \right)-\frac{\bar{\lambda }-P+{{m}_{f}}{{v}^{2}}}{EA}=0。$

(10)

利用有限单元法将立管沿全长离散成N个单元，引入Hermite插值函数A_l(s)和P_m(s)，用单元两端节点处的变量U_il、λ_m来表示单元内部变量：

$ {{r}_{i}}\left( s, t \right)={{A}_{l}}\left( s \right){{U}_{il}}\left( t \right)。$

(11)

$ \bar{\lambda }\left( s, t \right)={{P}_{m}}\left( s \right){{{\bar{\lambda }}}_{m}}\left( t \right)。$

(12)

式中：A_l(s)、P_m(s)分别为二次、三次Hermite形函数；U_il(t)为单元两端节点处的位置向量和切线向量；λ_m(t)为单元两端节点和中点处的有效张力向量。

最终得到立管运动方程的矩阵形式：

$ \begin{align} &({{M}_{ijlk}}+M_{ijlk}^{a}){{{\ddot{U}}}_{jk}}+C_{ijlk}^{f}{{{\dot{U}}}_{jk}}+(K_{ijlk}^{1}+{{{\bar{\lambda }}}_{n}}K_{nijlk}^{2}){{U}_{jk}}-\\ &\ \ \ \ \ {{F}_{il}}=0 。\\ \end{align} $

(13)

$ {{A}_{mil}}{{U}_{ki}}{{U}_{kl}}-{{B}_{m}}+{{C}_{mn}}({{h}_{n}}-{{{\bar{\lambda }}}_{n}})=0。$

(14)

式中：C_ijlk^f为内流科氏力产生的阻尼矩阵；F_il为荷载向量；M_ijlk为立管与内流的质量矩阵，${M_{ijlk}} = \int_0^L {(\rho + {m_f}){A_l}{A_k}{\delta _{ij}}{\rm{d}}s} $；M_ijlk^a为附加质量矩阵, ${{M}_{_{ijlk}^{a}}}={{C}_{A~}}\left[\int_{0}^{L}{{{A}_{l~}}{{A}_{k~}}{{\delta }_{ij~}}\text{d}s}-\left(\int_{0}^{L}{{{A}_{l~}}{{A}_{k~}}A{{'}_{s~}}A{{'}_{t~}}\text{d}s} \right){{U}_{it~}}{{U}_{js}} \right]$，C_A为单位长度附加质量；K_ijlk¹为材料抗弯刚度引起的刚度矩阵，$K_{ijlk}^{1}=\int_{0}^{L}{EIA{{''}_{l}}A{{''}_{k}}{{\delta }_{ij}}\text{d}s}$；K_nijlk²为张力和曲率引起的几何刚度矩阵，$K_{nijlk}^{2}=\int_{0}^{L}{{{P}_{n}}A{{'}_{l}}A{{'}_{k}}{{\delta }_{ij}}\text{d}s}$。

1.4 模态分析

为描述模态振型与钢悬链线立管位形之间的相对位置关系，首先求解立管在重力、内流等荷载作用下的位形^[14]，即图 2、3、4中黑色实线表示的曲线。舍去钢悬链线立管运动方程中的惯性力项和阻尼力项，则

$ (K_{ijlk}^{1}+{{{\bar{\lambda }}}_{n}}K_{nijlk}^{2}){{U}_{jk}}-{{F}_{il}}=0。$

(15)

采用Newton-Raphson法求解非线性方程，将所有单元集成并施加边界条件，通过迭代计算，即可求解未知量U、λ，进而求得钢悬链线立管的位形。

忽略钢悬链线立管运动方程中的阻尼力项和荷载项，得：

$ ({{M}_{ijlk}}+M_{ijlk}^{a}){{{\ddot{U}}}_{jk}}+(K_{ijlk}^{1}+{{{\bar{\lambda }}}_{n}}K_{nijlk}^{2}){{U}_{jk}}=0。$

(16)

将方程的解U=φe^iωt代入运动方程，得到特征值问题：

$ (K_{ijlk}^{1}+{{{\bar{\lambda }}}_{n}}K_{nijlk}^{2})\left\{ \varphi \right\}={{\omega }^{2}}({{M}_{ijlk}}+{{M}_{_{ijlk}^{a}}})\left\{ \varphi \right\}。$

(17)

式中：ω为立管的固有频率；{φ}为相应的模态。

2 立管动力特性分析 2.1 模型验证

基于以上考虑内流作用的细长柔性杆模型，本文用Matlab编写了钢悬链线立管的计算分析程序SCRModal，并对钢悬链线输流立管进行了分析计算。算例参数如下；立管密度为7 850 kg/m³，弹性模量2.07×10¹¹ Pa，海水密度为1 025 kg/m³，内部流体密度为998 kg/m³, 立管内、外径分别为0.20和0.26 m，立管总长1 040 m，顶部偏移190 m，水深1 000 m。

为了验证模型的准确性，将本文计算得到的钢悬链线立管的前六阶固有频率与相同参数条件下商业软件OrcaFlex的计算结果进行对比，结果如表 1。

从以上比较结果可以看出，本文模型与OrcaFlex的计算结果吻合较好。

2.2 动力特性分析

应用上述算例，探讨内流流速、顶张力和管内顶端压强对钢悬链线立管动力特性的影响。为了更好地描述模态振型与钢悬链线立管位形之间的相对位置关系，对立管模态位移的增加量进行放大，放大倍数10 000，然后与上一迭代步的位移相加，得到放大的立管模态。

2.2.1 内流流速的影响

根据实际工程，内流流速分别取5、10、15 m/s，计算钢悬链线立管的前六阶固有频率列成表 2；将前三阶模态绘成图 2。

从表中可以看出，随着内流流速的增加，钢悬链线立管的固有频率降低，虽然实际海洋工程中立管内流的流速不会很高，但随着内流流速的增加，立管的固有频率可能降低到“锁振”频率范围，使得立管产生动力放大，所以立管内流流速对钢悬链线立管固有频率的影响应引起足够重视。

内流流速的改变不仅会影响立管的固有频率，还会影响立管的模态形状。从图中可以看出，随着内流流速的增加，模态有下沉趋势，高阶模态反转点向立管底端移动。随着模态数的增加，这种变化更加明显。

2.2.2 顶张力的影响

顶张力分别取1 547、1 687、1 828 kN，计算钢悬链线立管的前十五阶固有频率和前三阶模态并绘成图 3。

从图中可以看出，随着顶张力的增加，钢悬链线立管的固有频率增加。这是由于顶张力的增加引起立管各单元的有效张力增加，立管刚度随之增加，固有频率因此增加。随着顶张力的增加，高阶模态的反转点向立管顶部移动。

2.2.3 管内顶端压强的影响

管内顶端压强分别取0、20、50 MPa，计算钢悬链线立管的前六阶固有频率列成表 3；将前三阶模态绘成图 4。

从表 3中可以看出，随着管内顶端压强的增加，立管的固有频率有所降低，但是变化不大。从图 4中可以看出，立管的模态有下沉趋势，随着模态数的增加，这一现象更加明显。这是由于管内顶端压强的增大引起立管的轴向张力减小，由管内张力和曲率引起的几何刚度矩阵随之减小，越是高阶模态，这种影响越明显。深水环境下的钢悬链线立管能达到几千米长，且相比浅水环境，深水环境更容易激发钢悬链线立管的高阶模态，所以对处于深水环境的钢悬链线立管，应该更加重视管内顶端压强对其动力特性的影响。

3 结论

本文以钢悬链线立管为研究对象，考虑立管轴向变形和内流作用的影响，采用细长柔性杆模型建立管内流体荷载和管外海洋环境荷载共同作用下钢悬链线立管的运动方程，运用有限元法求解，从而计算分析立管的动力特性, 探讨其在不同内流流速、顶张力和管内顶端压强等条件下的变化规律。研究结果表明：

(1) 随着内流流速的增加，钢悬链线立管的固有频率降低，高阶模态的反转点向立管底端移动。

(2) 随着顶张力的增加，钢悬链线立管的固有频率增加，高阶模态的反转点向立管顶部移动。

(3) 随着管内顶端压强的增加，钢悬链线立管的固有频率有所降低。