一直以来, 深水立管涡激振动的横向振动(垂直于立管轴向的振动)得到了国内外学者的广泛深入的研究^[1]。目前国内外学者对参激振动和横向振动的耦合振动也有不少研究成果^[2-8], 但涉及考虑立管大变形的参激振动和横向振动的耦合振动研究较少。笔者在以往提出的考虑大变形深水立管涡激振动模型^[9-15]的基础上, 考虑深水立管顶端浮体的垂荡运动对立管涡激振动的影响, 提出考虑大变形的深水立管参激振动与涡激振动的耦合振动模型(以下统称耦合模型), 并在此模型的基础上研究考虑大变形(结构挠度与结构直径比例大于1.0即认为发生了大变形)时参数激扰对立管涡激振动响应的影响。

1 数学模型

图 1为深海立管坐标系的选取及弯曲示意图。由考虑大变形的深水立管的涡激振动模型^[9]得顺流向运动方程和横向运动方程。

顺流向运动方程为

$ \begin{array}{l} \bar m\frac{{{\partial ^2}x\left( {z,t} \right)}}{{\partial {t^2}}} + {E_{\rm{I}}}\frac{{{\partial ^4}x\left( {z,t} \right)}}{{\partial {z^4}}} - \left( {T - {E_{\rm{I}}}{\kappa ^2}} \right)\frac{{{\partial ^2}x\left( {z,t} \right)}}{{\partial {z^2}}} + \\ c\frac{{\partial x\left( {z,t} \right)}}{{\partial t}} = \frac{1}{2}{{C'}_{\rm{D}}}\rho D{\left( {U - \dot x} \right)^2}\cos \left( {4{\rm{ \mathsf{ π} }}{{f'}_{\rm{s}}}t} \right). \end{array} $

(1)

横向运动方程为

$ \begin{array}{l} \bar m\frac{{{\partial ^2}y\left( {z,t} \right)}}{{\partial {t^2}}} + {E_{\rm{I}}}\frac{{{\partial ^4}y\left( {z,t} \right)}}{{\partial {z^4}}} - \left( {T - {E_{\rm{I}}}{\kappa ^2}} \right)\frac{{{\partial ^2}y\left( {z,t} \right)}}{{\partial {z^2}}} + \\ c\frac{{\partial y\left( {z,t} \right)}}{{\partial t}} = \frac{1}{2}{C_{\rm{L}}}\rho D{\left( {U - \dot x} \right)^2}\cos \left( {2{\rm{ \mathsf{ π} }}{{f'}_{\rm{s}}}t} \right). \end{array} $

(2)

式中, $ \bar m$为立管微元段平均质量; E_I为抗弯刚度; T为有效张力; t为时间变量; κ(z, t)为立管弯曲振动时的轴线曲率; c为结构阻尼; C′_D为漩涡释放引起的脉动拖曳力系数; C_L为横向涡激升力系数; ρ为海水密度; D为立管直径; U为均匀海流流速; ${\dot x} $为立管横向振动速度; f′_s为涡泄频率; x(z, t)为平面内的挠曲线(顺流向); y(z, t)为出平面的挠曲线(横向); z为水深坐标。

参激振动与涡激振动的耦合振动模型^[14]的顺流向运动方程为

$ \begin{array}{l} \bar m\frac{{{\partial ^2}x\left( {z,t} \right)}}{{\partial {t^2}}} + {E_{\rm{I}}}\frac{{{\partial ^4}x\left( {z,t} \right)}}{{\partial {z^4}}} - T\frac{{{\partial ^2}x\left( {z,t} \right)}}{{\partial {z^2}}} + c\frac{{\partial x\left( {z,t} \right)}}{{\partial t}} = \\ {G_{{\rm{As}}}}\kappa \frac{{\partial u\left( {z,t} \right)}}{{\partial z}} + \frac{1}{2}{{C'}_{\rm{D}}}\rho D{\left( {U - \dot x} \right)^2}\cos \left( {4{\rm{ \mathsf{ π} }}{{f'}_{\rm{s}}}t} \right). \end{array} $

(3)

横向运动方程为

$ \begin{array}{l} \bar m\frac{{{\partial ^2}y\left( {z,t} \right)}}{{\partial {t^2}}} + {E_{\rm{I}}}\frac{{{\partial ^4}y\left( {z,t} \right)}}{{\partial {z^4}}} - T\frac{{{\partial ^2}y\left( {z,t} \right)}}{{\partial {z^2}}} + c\frac{{\partial y\left( {z,t} \right)}}{{\partial t}} = \\ {G_{{\rm{As}}}}\kappa \frac{{\partial u\left( {z,t} \right)}}{{\partial z}} + \frac{1}{2}{C_{\rm{L}}}\rho D{\left( {U - \dot x} \right)^2}\cos \left( {2{\rm{ \mathsf{ π} }}{{f'}_{\rm{s}}}t} \right). \end{array} $

(4)

式中, G_As为立管的抗剪刚度。

模型假设立管的轴线不可伸长, 现将两模型通过几何变形数值叠加得到考虑大变形的深水立管参激与涡激振动耦合模型, 其顺流向运动方程为

$ \begin{array}{l} \bar m\frac{{{\partial ^2}x\left( {z,t} \right)}}{{\partial {t^2}}} + {E_{\rm{I}}}\frac{{{\partial ^4}x\left( {z,t} \right)}}{{\partial {z^4}}} - \left( {T - {E_{\rm{I}}}{\kappa ^2}} \right)\frac{{{\partial ^2}x\left( {z,t} \right)}}{{\partial {z^2}}} + \\ c\frac{{\partial x\left( {z,t} \right)}}{{\partial t}} = {G_{{\rm{As}}}}\kappa \frac{{\partial u\left( {z,t} \right)}}{{\partial z}} + \frac{1}{2}{{C'}_{\rm{D}}}\rho D\left( {U - } \right.\\ {\left. {\dot x} \right)^2}\cos \left( {4{\rm{ \mathsf{ π} }}{{f'}_{\rm{s}}}t} \right). \end{array} $

(5)

横向运动方程为

$ \begin{array}{l} \bar m\frac{{{\partial ^2}y\left( {z,t} \right)}}{{\partial {t^2}}} + {E_{\rm{I}}}\frac{{{\partial ^4}y\left( {z,t} \right)}}{{\partial {z^4}}} - \left( {T - {E_{\rm{I}}}{\kappa ^2}} \right)\frac{{{\partial ^2}y\left( {z,t} \right)}}{{\partial {z^2}}} + \\ c\frac{{\partial y\left( {z,t} \right)}}{{\partial t}} = {G_{{\rm{As}}}}\kappa \frac{{\partial u\left( {z,t} \right)}}{{\partial z}} + \frac{1}{2}{C_{\rm{L}}}\rho D\left( {U - } \right.\\ {\left. {\dot x} \right)^2}\cos \left( {2{\rm{ \mathsf{ π} }}{{f'}_{\rm{s}}}t} \right). \end{array} $

(6)

2 数值分析 2.1 数值求解

分析提出的模型特点发现, 参数激扰对结构影响与结构的挠度、支座位移有关, 而结构的挠度和支座位移是随时间不断变化的, 系统的非线性特征较明显, 所以对于运动方程的数值求解采用适用于非线性系统的Newmark-β法的增量形式进行求解。依据Newmark-β逐步积分法, δ和β是与精度和稳定性有关的参数。当δ>0.5时, 将产生算法阻尼, 从而使振幅人为衰减; 当δ < 0.5时, 产生负阻尼, 积分计算过程中振幅逐步增长, 通常取临界值δ=0.5, β=0.25×(0.5+δ)²=0.25, 此时的积分无条件稳定^[13]。

当支座位移是定常值时, 则问题回到梁的复杂弯曲振动。因此考虑立管顶端边界支座位移随时间变化的情况。

设

$ u\left( {z,t} \right) = {u_0}\lambda \left( t \right), $

(7)

则

$ \frac{{\partial u}}{{\partial z}} \approx \frac{{{u_0}}}{l}\lambda \left( t \right). $

(8)

式中, u(z, t)为立管顶端位移; u₀为立管顶端位移幅值; l为立管长度; λ(t)为一个已知时间函数。

λ(t)的取值根据参数激扰结构的运动规律得到, 对于深水立管结构而言, 主要考虑立管顶部浮体垂荡运动的影响。本文中主要研究在浮体的垂荡运动影响下立管的涡激振动。

此时, 立管的涡激振动模型式(5)和(6)化为

$ \begin{array}{l} \bar m\frac{{{\partial ^2}x\left( {z,t} \right)}}{{\partial {t^2}}} + {E_{\rm{I}}}\frac{{{\partial ^4}x\left( {z,t} \right)}}{{\partial {z^4}}} - \left( {T - {E_{\rm{I}}}{\kappa ^2}} \right)\frac{{{\partial ^2}x\left( {z,t} \right)}}{{\partial {z^2}}} + \\ c\frac{{\partial x\left( {z,t} \right)}}{{\partial t}} = {G_{{\rm{As}}}}\kappa \frac{{{u_0}}}{l}\lambda \left( t \right) + \frac{1}{2}{{C'}_{\rm{D}}}\rho D\left( {U - } \right.\\ {\left. {\dot x} \right)^2}\cos \left( {4{\rm{ \mathsf{ π} }}{{f'}_{\rm{s}}}t} \right), \end{array} $

(9)

$ \begin{array}{l} \bar m\frac{{{\partial ^2}y\left( {z,t} \right)}}{{\partial {t^2}}} + {E_{\rm{I}}}\frac{{{\partial ^4}y\left( {z,t} \right)}}{{\partial {z^4}}} - \left( {T - {E_{\rm{I}}}{\kappa ^2}} \right)\frac{{{\partial ^2}y\left( {z,t} \right)}}{{\partial {z^2}}} + \\ c\frac{{\partial y\left( {z,t} \right)}}{{\partial t}} = {G_{{\rm{As}}}}\kappa \frac{{{u_0}}}{l}\lambda \left( t \right) + \frac{1}{2}{C_{\rm{L}}}\rho D\left( {U - } \right.\\ {\left. {\dot x} \right)^2}\cos \left( {2{\rm{ \mathsf{ π} }}{{f'}_{\rm{s}}}t} \right). \end{array} $

(10)

分析式(9)和(10)可知, 在有横向扰动的条件下, 参数激扰将与横向激扰组合对立管涡激振动产生影响。由于参数激扰的大小与立管横向弯曲的挠度有关, 因此考虑参数激扰时, 系统的弯曲振动方程是一个隐式方程。式(9)和式(10)可分别表示为

$ \begin{array}{l} \bar m\frac{{{\partial ^2}x\left( {z,t} \right)}}{{\partial {t^2}}} + {E_{\rm{I}}}\frac{{{\partial ^4}x\left( {z,t} \right)}}{{\partial {z^4}}} - \left( {T - {E_{\rm{I}}}{\kappa ^2}} \right)\frac{{{\partial ^2}x\left( {z,t} \right)}}{{\partial {z^2}}} + \\ c\frac{{\partial x\left( {z,t} \right)}}{{\partial t}} = {G_{{\rm{As}}}}\kappa \frac{{\partial u\left( {z,t} \right)}}{{\partial z}}\frac{{{u_0}}}{l}\lambda \left( t \right) + \frac{1}{2}{{C'}_{\rm{D}}}\rho D\left( {U - } \right.\\ {\left. {\dot x} \right)^2}\cos \left( {4{\rm{ \mathsf{ π} }}{{f'}_{\rm{s}}}t} \right), \end{array} $

(11)

$ \begin{array}{l} \bar m\frac{{{\partial ^2}y\left( {z,t} \right)}}{{\partial {t^2}}} + {E_{\rm{I}}}\frac{{{\partial ^4}y\left( {z,t} \right)}}{{\partial {z^4}}} - \left( {T - {E_{\rm{I}}}{\kappa ^2}} \right)\frac{{{\partial ^2}y\left( {z,t} \right)}}{{\partial {z^2}}} + \\ c\frac{{\partial y\left( {z,t} \right)}}{{\partial t}} = {G_{{\rm{As}}}}\frac{{{\partial ^2}x\left( {z,t} \right)}}{{\partial {z^2}}}\frac{{{u_0}}}{l}\lambda \left( t \right) + \frac{1}{2}{C_{\rm{L}}}\rho D\left( {U - } \right.\\ {\left. {\dot x} \right)^2}\cos \left( {2{\rm{ \mathsf{ π} }}{{f'}_{\rm{s}}}t} \right). \end{array} $

(12)

将式(11)、(12)分离变量后得到的形式为

$ \begin{array}{l} \bar m\sum\limits_{n = 1}^\infty {{\phi _n}} \left( z \right){{\ddot q}_n}\left( t \right) + {E_{\rm{I}}}\sum\limits_{n = 1}^\infty {{{\phi ''''}_n}\left( z \right){q_n}\left( t \right)} - \left( {T - } \right.\\ \left. {{E_{\rm{I}}}{\kappa ^2}} \right)\sum\limits_{n = 1}^\infty {\phi _n^{\prime \prime }} \left( z \right){q_n}\left( t \right) = f_{\rm{D}}^\prime (z,t) + \\ {G_{{\rm{As}}}}\sum\limits_{n = 1}^\infty {\phi _n^{\prime \prime }} \left( z \right){q_n}\left( t \right)\frac{{{u_0}}}{l}\lambda \left( t \right), \end{array} $

(13)

$ \begin{array}{l} \bar m\sum\limits_{n = 1}^\infty {{\phi _n}} \left( z \right){{\ddot q}_n}\left( t \right) + {E_{\rm{I}}}\sum\limits_{n = 1}^\infty {{{\phi ''''}_n}\left( z \right){q_n}\left( t \right)} - \left( {T - } \right.\\ \left. {{E_{\rm{I}}}{\kappa ^2}} \right)\sum\limits_{n = 1}^\infty {\phi _n^{\prime \prime }} \left( z \right){q_n}\left( t \right) = {f_{\rm{L}}}(z,t) + \\ {G_{{\rm{As}}}}\sum\limits_{n = 1}^\infty {\phi _n^{\prime \prime }} \left( z \right){q_n}\left( t \right)\frac{{{u_0}}}{l}\lambda \left( t \right). \end{array} $

(14)

利用振型的正交性将式(13)、(14)简化为

$ {M_n}{{\ddot q}_n}\left( t \right) + {K_n}{q_n}\left( t \right) + {N_n}{q_n}\left( t \right) = {F_n}\left( t \right) - {K_{{\rm{G}}n}}{q_n}\left( t \right)\lambda \left( t \right). $

(15)

其中

$ {M_n} = \bar m\int_0^l {\phi _n^2} \left( z \right){\rm{d}}z, $

$ {K_n} = {E_{\rm{I}}}\int_0^l {\phi {''}_n^2} \left( z \right){\rm{d}}z, $

$ {N_n} = \left( {{E_{\rm{I}}}{\kappa ^2} - T} \right)\int_0^l {{{\phi ''}_n}} \left( z \right){\rm{d}}z, $

$ {K_{{\rm{G}}n}} = {G_{{\rm{As}}}}\frac{{{u_0}}}{l}\int_0^l {\phi {'}_n^2} \left( z \right){\rm{d}}z, $

$ {F_n}\left( t \right) = \int_0^l {{\phi _n}} \left( z \right){{f'}_{\rm{D}}}\left( z \right){\rm{d}}z/\int_0^l {{\phi _n}} \left( z \right){f_{\rm{L}}}\left( z \right){\rm{d}}z. $

式(15)右端的第二项中包括未知的广义坐标响应, 因此只要支座位移不是常数, 就必须采用数值方法求解。式(15)的非迭代格式可表示为

$ \begin{array}{l} {M_n}{{\ddot q}_n}\left( {t + \Delta t} \right) + {K_n}{q_n}\left( {t + \Delta t} \right) + {N_n}{q_n}\left( {t + \Delta t} \right) = {F_n}\left( {t + \Delta t} \right) - \\ {K_{{\rm{G}}n}}{q_n}\left( t \right)\lambda \left( {t + \Delta t} \right). \end{array} $

(16)

而迭代格式可表示为

$ \begin{array}{l} {M_n}\ddot q_n^{i + 1}\left( t \right) + {K_n}q_n^{i + 1}\left( t \right) + {N_n}q_n^{i + 1}\left( t \right) = {F_n}\left( t \right) - \\ {K_{{\rm{G}}n}}q_n^i\left( t \right)\lambda \left( t \right). \end{array} $

(17)

2.2 模型验证

以提出的考虑参数激扰的深水立管涡激振动模型(式(9)和(10))为计算原理, 编制相应的深水立管涡激振动响应计算程序TTRPD1.0。

将提出的深水立管涡激振动模型应用于深水顶张力立管的涡激振动响应分析, 并同立管涡激振动分析软件shear7的计算结果进行对比(图 2)。shear7只能计算立管横向振动响应的均方根位移, 所以程序TTRPD1.0只计算了横向振动的均方根位移。

3 算例分析 3.1 基础参数及边界和初始条件

基础参数:立管长度为1500 m, 外径为0.3239 m, 内径为0.285 m, 弹性模量为207 GPa, 剪切模量为79 kPa, 材料密度为7850 kg/m³, 立管的顶张力为5400 kN, 阻尼系数为0.05。环境荷载模拟海流流速分别为0.06、0.1、0.2、0.4、0.6、0.8和1.0 m/s。响应参数:拖曳系数为1.0, 附加质量系数为1.0, 升力系数为0.9。这里只计算深水立管在匀速海流作用下的涡激振动响应。深水顶张力立管在海底是通过深水插式连接器与井口连接, 其上为锥形的应力接头, 因此应力接头处设置固定端约束较为合理。本文中主要研究深水半潜式钻井平台的垂荡运动对深水立管涡激振动的影响, 将顶张力立管的顶端视为可以沿着立管轴向发生位移的简支端。

立管的响应是由顶部平台和海流等环境荷载等引起的强迫运动, 所以可以忽略运动初始条件对稳态运动的影响, 因此初始条件取为零, 即:初始位移为零(x₀=0, y₀=0), 初始速度为零($ {{\dot x}_0}$ =0, $ {{\dot y}_0}$ =0)和初始加速度也为零($ {{\ddot x}_0}$ =0, ${{\ddot y}_0} $ =0)。

3.2 动力响应

选取立管两向涡激振动发生最大位移的节点作为参考点研究其响应特征, 研究不同流速下, 参激激扰对立管响应的影响, 此时立管顶部浮体的垂荡运动幅值取为A=3.0 m, 垂荡频率取为ω=1.6 rad/s。图 3为不同均匀流荷载作用下立管耦合模型与不考虑参数激扰的非耦合模型的对比位移时程。

3.2.1 位移时程

(1) 横向涡激振动响应。由图 3可见, 流速为0.06 m/s时, 耦合振动模型与不考虑参激振动的非耦合模型在参考点处的时程曲线几乎重合。由式(9)知, 参数激扰与立管结构的曲率密切相关, 在低流速下, 立管结构的变形较小, 相应结构的局部曲率较小, 此时参数激扰对涡激振动几乎没有影响, 这一点与理论分析相对应。此时参数激扰对立管振动响应的影响可忽略。同时也证明了计算程序的正确性。

表 1为立管响应位移最大的点在两种模型下的振幅变化情况。由表 1和式(9)分析可知, 参数激扰对立管涡激振动的影响不仅与环境荷载有关, 而且与参数激扰的幅值和频率相关。立管顶部平台的垂荡运动对立管横向涡激振动响应振幅的影响较小, 当流速为0.4 m/s时, 振幅增加最大为5%。并且立管顶部平台的垂荡运动对立管涡激振动响应的影响并不是随着流速的增大而增大, 因此立管顶部平台的垂荡运动对立管横向涡激振动响应的影响是环境荷载的函数, 但两者并非线性关系。

(2) 顺流向涡激振动响应。图 4为不同均匀流荷载作用下立管涡激振动顺流向的响应时程。表 2为参数激扰对立管顺流向振幅的影响。从图 4中可以看出, 参数激扰对立管的涡激振动顺流向响应的影响同样是随着流速变化而变化。

分析立管顺流向振动位移时程图 4和表 2可知, 参数激扰对顺流向响应的影响较大。平台的垂荡运动对立管顺流向涡激振动响应幅值的影响较大。并且在不同流速下, 平台的垂荡运动的影响是立管环境荷载的函数, 但两者之间并不呈线性关系。

将立管顺流向振动响应的时程经傅里叶变化得到立管顺流向振动的频谱图, 如图 5所示。

分析立管在不同流速下顺流向振动的频谱图可知, 平台的垂荡运动影响了立管的响应频率, 并且在不同流速下, 其影响程度不同, 耦合模型的响应频率多于非耦合模型的响应频率。

3.2.2 位形分析

图 6为立管在不同流速下发生最大位移时的位形图。

分析立管响应的两向位形图可知, 参数激扰对顺流向响应的影响, 包括振幅和频率的影响都明显大于对横向振动响应的影响。考虑参数激扰的耦合振动模型也出现了多阶模态响应和位形不完全对称的现象。

3.3 参数激扰相关参数对响应的影响

对于平台垂荡运动引起的参数激扰, 主要由平台垂荡的幅值和平台的垂荡频率。接下来分别研究这两个参数对立管动力响应的影响。

图 7为在顶部平台不同垂荡幅值情况下立管各节点的最大位移。此时, 立管的顶部平台的垂荡频率为ω=1.60 rad/s。

由图 7可以看出, 当平台垂荡频率不变时, 立管的涡激振动响应的位移随着平台垂荡幅值的增大而增大。因此深水平台对立管产生的参数激扰对立管具有危害性, 其作用不可忽略。

图 8为不同垂荡频率下立管的动力响应(垂荡幅值A=3.0 m)。由图 8可以看出, 保持平台的垂荡幅值不变A=3.0 m, 立管结构的动力响应节点位移随着平台垂荡频率的增大而减小。这一现象与TomoFUJIWARA在国家海洋研究院所做的长28.5 m的柔性顶张力立管涡激振动试验的结论相同^[16]。

4 结论

(1) 立管固有频率随着顶端平台的运动而变化。平台的垂荡运动使立管涡激振动的振幅增大, 并且顺流向振幅增大的幅度大于横向振幅的幅度, 参数激扰对立管的涡激振动响应的频率也产生影响。

(2) 当平台的垂荡频率和环境荷载不变时, 立管的动力响应的位移随着浮体垂荡幅值的增大而增大。当平台的垂荡幅值和环境荷载不变时, 立管结构的动力响应位移随着浮体垂荡频率的增大而减小。立管的涡激振动响应是环境荷载、立管结构参数、顶部平台垂荡运动的函数。