CN102507082A - 一种深水立管的时域涡激升力确定方法 - Google Patents

Abstract

本发明涉及海洋深水立管的研究方法，具体涉及一种深水立管的时域涡激升力确定方法。该方法将由深水立管横流向振动速度和加速度产生的粘性阻力和附加质量力引入了涡激升力的计算，并采用迭代方法计算时域的涡激升力，从而建立了一个完整的流固耦合时域涡激升力计算方法。由于本发明增加了深水立管振动速度和加速度产生的粘性阻力和附加质量力，比现有方法更符合圆柱体横流向涡激振动的受力状态。

Description

一种深水立管的时域涡激升力确定方法

技术领域

本发明涉及海洋深水立管的研究方法，具体涉及一种深水立管的时域涡激升力确定方法。

背景技术

涡激升力是圆柱体受到的垂直于流体流动方向的一种交变流体荷载，是由圆柱体尾流处的涡旋泄放引起的。对于弹性圆柱体而言，涡激升力将使圆柱体产生垂直于流体流动方向的振动。而对于刚性圆柱体，涡激升力并不引起圆柱体振动。这就是说，不论圆柱体是否振动，涡激升力始终是存在的。

涡激升力的交变频率等于涡旋泄放的频率，而涡旋泄放的频率取决于流体流经圆柱体的速度大小，涡激升力的大小取决于流体流经圆柱体的速度。因此，对于刚性圆柱体来说，涡激升力的大小和频率仅仅与流体的流动速度有关。而对于弹性圆柱体，由于圆柱体本身的振动，不仅使得流体流经圆柱体的速度发生变化，而且使得流体在圆柱体的扰动下产生粘性阻力和附加质量力。因此，对于弹性圆柱体的涡激升力计算不能仅仅考虑流场的流速一个因素(传统涡激升力计算方法的缺点)，而应考虑圆柱体振动的速度和加速度。

深水立管或海底管线属于弹性圆柱体，其涡激振动为弯曲振动，因此，沿圆柱体轴向的不同位置，其振动速度和加速度是不同的。如果采用相同的涡激升力计算，则误差是可想而知的。

总体来说，现有技术存在如下缺陷：

(1)没有考虑深水立管(弹性圆柱体)横流向振动的加速度效应，忽略了深水立管周围流场压力梯度产生的附加质量力。因此，计算得到的深水立管振动频率偏低；

(2)没有考虑深水立管横流向振动的速度效应，忽略了深水立管与流体相对速度变化引起的粘性阻力，因此，计算结果偏小。

发明内容

本发明的目的在于针对现有技术的缺陷，提供一种考虑深水立管横流向振动加速度引起的附加质量力和振动速度引起的粘性阻力的深水立管的时域涡激升力确定方法，使计算结果更符合涡激振动的自然规律。

本发明的技术方案如下：一种深水立管的时域涡激升力确定方法，建立的涡激升力时域模型如下：

$\mathrm{FL}=\frac{1}{2}{C}_L\mathrm{\ρD}{(U-\stackrel{.}{u})}^2\cos {\mathrm{\ω}}_s^{\′}t-\frac{1}{2}{C}_D\mathrm{\ρD}\stackrel{.}{v}\left|\stackrel{.}{v}\right|-\frac{{\mathrm{\π}}^2}{4}\mathrm{\ρD}\stackrel{..}{v}$

式中：FL--涡激升力；

C_L--升力系数；

C_D--拖曳力系数；

ρ--流体密度；

D--深水立管直径；

U--流速；

--深水立管顺流向振动速度；

--深水立管横流向振动速度；

ω′_s--涡旋泄放频率，

其中，St为斯特罗哈数，为圆柱体顺流向振动速度；

t--时间；

--深水立管横流向振动加速度；

针对深水立管顺流向振动速度、横流向振动速度和加速度，采用迭代方法对上式进行计算，得到深水立管的时域涡激升力。

进一步，如上所述的深水立管的时域涡激升力确定方法，其中，所述的采用迭代方法进行计算的具体步骤如下：

1)给定深水立管顺流向振动速度、横流向振动速度和加速度及计算时间的初值：

t_j＝t₀＝0， ${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_0^{\left(0\right)}=0,$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_0^{\left(0\right)}=0,$ ${\stackrel{.}{v}}_j^{\left(i\right)}={\stackrel{.}{v}}_0^{\left(0\right)}=0,$ ${\stackrel{..}{v}}_j^{\left(i\right)}={\stackrel{..}{v}}_0^{\left(0\right)}=0$

式中：j--时间步数，计算开始时j＝0；

i--迭代次数，每个时间步开始时i＝0；

2)计算给定流速下的约化速度：

$V_r=\frac{U}{f_{n}D}$

式中：V_r--约化速度；

      U--流体流速；

      f_n--深水立管的固有频率；

      D--深水立管直径；

3)当V_r＜5或V_r＞7时，将t_j，代入下式计算第j时间步内第i次迭代的脉动拖曳力：

${\mathrm{FD}}_j^{\left(i\right)}=\frac{1}{2}{C}_D\mathrm{\ρD}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^2\cos {\mathrm{\ω}}_s^{\′}{t}_0+$

$\frac{1}{2}{C}_D\mathrm{\ρD}(U-{\stackrel{.}{u}}_j^{\left(i\right)})|U-{\stackrel{.}{u}}_j^{\left(i\right)}|-\frac{{\mathrm{\π}}^2}{4}\mathrm{\ρD}{\stackrel{..}{u}}_j^{\left(i\right)}$

当5≤V_r≤7时，将t_j，

代入下式计算第j时间步内第i次迭代的脉动拖曳力：

${\mathrm{FD}}_j^{\left(i\right)}=\frac{1}{2}{C}_D\mathrm{\ρD}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^2\cos {\mathrm{\ω}}_s^{\′}{t}_0+$

$\frac{1}{2}{C}_D\mathrm{\ρD}(U-{\stackrel{.}{u}}_j^{\left(i\right)})|U-{\stackrel{.}{u}}_j^{\left(i\right)}|-\frac{{\mathrm{\π}}^2}{4}\mathrm{\ρD}{\stackrel{..}{u}}_j^{\left(i\right)}$

在以上两个公式中，

FD--脉动拖曳力；

C_D--拖曳力系数；

ρ--流体密度；

D--深水立管直径；

U--流体流速；

--深水立管顺流向振动速度；

ω′_s--涡旋泄放频率，

其中，St为斯特罗哈数；

t--时间；

ü--深水立管顺流向振动加速度；

4)将步骤3)中计算得到的脉动拖曳力

代入深水立管的振动方程式：

${m\stackrel{..}{u}}_j^{(i+1)}+c{\stackrel{.}{u}}_j^{(i+1)}+k{u}_j^{(i+1)}={\mathrm{FD}}_j^{\left(i\right)}$

式中：m--深水立管的质量；

      c--深水立管的阻尼系数；

      k--深水立管的弯曲刚度；

--第j时间步内第i+1次迭代的顺流向振动位移；

--第j时间步内第i+1次迭代的顺流向振动速度；

--第j时间步内第i+1次迭代的顺流向振动加速度；

计算第j时间步内第i+1次迭代的深水立管顺流向涡激振动的速度和加速度

5)如果

ε为预先设定的计算精度，则继续进行迭代计算，令：

${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_j^{(i+1)},$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_j^{(i+1)}$

然后，重复步骤3)～5)的计算；

如果 $\max \left\{\right|{\stackrel{..}{u}}_j^{(i+1)}-{\stackrel{..}{u}}_j^{\left(i\right)}|,|{\stackrel{.}{u}}_j^{(i+1)}-{\stackrel{.}{u}}_j^{\left(i\right)}\left|\right\}\≤\mathrm{\ϵ},$ 则令：

${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{\·}{u}}_j^{\left(n\right)},$ ${\stackrel{.}{v}}_j^{\left(i\right)}={\stackrel{\·}{v}}_j^{\left(0\right)},$ ${\stackrel{..}{v}}_j^{\left(i\right)}={\stackrel{..}{v}}_j^{\left(0\right)}$

n为第j时间步内计算脉动拖曳力的总迭代次数；

对于第1个时间步即j＝0，

否则，

其中，m为第j-1时间步内计算涡激升力的总迭代次数；

6)将

代入下式计算第j时间步内第i次迭代的涡激升力：

${\mathrm{FL}}_j^{\left(i\right)}=\frac{1}{2}{C}_L\mathrm{\ρD}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^2\cos {\mathrm{\ω}}_s^{\′}{t}_j-\frac{1}{2}{C}_D\mathrm{\ρD}{\stackrel{.}{v}}_j^{\left(i\right)}\left|{\stackrel{.}{v}}_j^{\left(i\right)}\right|-\frac{{\mathrm{\π}}^2}{4}\mathrm{\ρD}{\stackrel{..}{v}}_j^{\left(i\right)}$

7)将步骤6)中计算得到的涡激升力代入深水立管的横流向涡激振动方程式：

${m\stackrel{..}{v}}_j^{(i+1)}+c{\stackrel{.}{v}}_j^{(i+1)}+k{v}_j^{(i+1)}={\mathrm{FL}}_j^{\left(i\right)}$

式中：m--深水立管的质量；

      c--深水立管的阻尼系数；

      k--深水立管的弯曲刚度；

--第j时间步内第i+1次迭代的横流向振动位移；

--第j时间步内第i+1次迭代的横流向振动速度；

--第j时间步内第i+1次迭代的横流向振动加速度；

计算第j时间步内第i+1次迭代的深水立管横流向涡激振动速度

和加速度

8)如果

ε为预先设定的计算精度，则继续进行迭代计算，令：

${\stackrel{.}{v}}_j^{\left(i\right)}={\stackrel{.}{v}}_j^{(i+1)},$ ${\stackrel{..}{v}}_j^{\left(i\right)}={\stackrel{..}{v}}_j^{(i+1)}$

然后，重复步骤6)～8)的计算；

如果

则开始下一个时间步的计算，令：

t_j＝t_j+1＝t_j+Δt， ${\stackrel{.}{v}}_j^{\left(i\right)}={\stackrel{.}{v}}_{j+1}^{\left(0\right)}={\stackrel{.}{v}}_j^{\left(m\right)},$ ${\stackrel{..}{v}}_j^{\left(i\right)}={\stackrel{..}{v}}_{j+1}^{\left(0\right)}={\stackrel{..}{v}}_j^{\left(m\right)}$

m为第j时间步内计算涡激升力的总迭代次数；

重复步骤3)～8)的计算，直至计算时长满足需要。

本发明的有益效果如下：本发明将由深水立管横流向振动速度和加速度产生的粘性阻力和附加质量力引入了涡激升力的计算，并采用迭代方法计算时域的涡激升力，从而建立了一个完整的流固耦合时域涡激升力计算方法。与现有的涡激升力计算方法相比，本发明具有下述优点：

(1)增加了深水立管振动速度和加速度产生的粘性阻力和附加质量力，比现有方法更符合圆柱体横流向涡激振动的受力状态；

(2)考虑了流固耦合作用，从而可以进行非线性的横流向涡激振动分析；

(3)可以计算深水立管轴向不同位置的涡激升力时程，因而可以得到涡激升力沿深水立管(弹性圆柱体)轴向的变化规律。

附图说明

图1为本发明的方法流程图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细描述。

传统的涡激振动理论认为，深水立管的顺流向涡激振动非常小，可以忽略不计。因此，深水立管的涡激振动就是垂直流体流动方向的振动被人们所接受。“涡激振动”仅仅意味着深水立管的垂直流向振动，这可以从目前的多款商业软件的计算方法得到证明。Orcaflex、VIVA和Shear7是目前海洋工程界公认的深水立管涡激振动分析软件，被广泛应用于深水立管的涡激振动与疲劳分析。但是，这些软件只提供横向涡激振动分析结果。当然现有的其它分析方法中也没有顺流向涡激振动的分析方法，唯一的计算方法是本发明的申请人已获得的前一项发明专利中提出的顺流向涡激振动分析方法“一种深水顶张式立管涡激振动与疲劳分析的方法”ZL200910136583.X。由于本发明提出的涡激升力计算方法考虑了流固耦合效应，计算时需要用到深水立管的顺流向振动速度。因此，本发明提出的涡激升力计算方法包括深水立管顺流向振动的计算。

本发明同时考虑深水立管横流向振动速度和加速度效应、并考虑了涡旋泄放锁定区和非锁定区圆柱体顺流向振动的不同性质，提出了基于下述模型的圆柱体时域涡激升力计算方法：

$\mathrm{FL}=\frac{1}{2}{C}_L\mathrm{\ρD}{(U-\stackrel{.}{u})}^2\cos {\mathrm{\ω}}_s^{\′}t-\frac{1}{2}{C}_D\mathrm{\ρD}\stackrel{.}{v}\left|\stackrel{.}{v}\right|-\frac{{\mathrm{\π}}^2}{4}\mathrm{\ρD}\stackrel{..}{v}---\left(1\right)$

式中：FL--涡激升力；

      C_L--升力系数；

      C_D--拖曳力系数；

      ρ--流体密度；

      D--深水立管直径；

      U--流速；

--深水立管顺流向振动速度；

--深水立管横流向振动速度；

ω′_s--涡旋泄放频率，

其中，St为斯特罗哈数，

为圆柱体顺流向振动速度；

t--时间；

--深水立管横流向振动加速度。

上式中，第一项为涡旋泄放引起的升力，第二项为深水立管振动速度引起的粘性阻力，第三项为深水立管振动加速度引起的附加质量力。

由于公式(1)包含圆柱体的顺流向振动速度、横流向振动速度和加速度，因此，必须采用迭代方法计算，如图1所示，具体计算步骤如下：

1)给定深水立管顺流向振动速度、横流向振动速度和加速度及计算时间的初值(公知技术)：

t_j＝t₀＝0， ${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_0^{\left(0\right)}=0,$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_0^{\left(0\right)}=0,$ ${\stackrel{.}{v}}_j^{\left(i\right)}={\stackrel{.}{v}}_0^{\left(0\right)}=0,$ ${\stackrel{..}{v}}_j^{\left(i\right)}={\stackrel{..}{v}}_0^{\left(0\right)}=0$

式中：j--时间步数，计算开始时j＝0；

i--迭代次数，每个时间步开始时i＝0；

2)计算给定流速(流速是设计或分析给定的条件，为已知值)下的约化速度(公知技术)：

$V_r=\frac{U}{f_{n}D}$

式中：V_r--约化速度；

U--流体流速；

f_n--深水立管的固有频率；

D--深水立管直径；

3)当V_r＜5或V_r＞7时，将t_j，

代入下式计算第j时间步内第i次迭代的脉动拖曳力：

${\mathrm{FD}}_j^{\left(i\right)}=\frac{1}{2}{C}_D\mathrm{\ρD}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^2\cos {\mathrm{\ω}}_s^{\′}{t}_0+$

$\frac{1}{2}{C}_D\mathrm{\ρD}(U-{\stackrel{.}{u}}_j^{\left(i\right)})|U-{\stackrel{.}{u}}_j^{\left(i\right)}|-\frac{{\mathrm{\π}}^2}{4}\mathrm{\ρD}{\stackrel{..}{u}}_j^{\left(i\right)}$

当5≤V_r≤7时，将t_j，

代入下式计算第j时间步内第i次迭代的脉动拖曳力：

${\mathrm{FD}}_j^{\left(i\right)}=\frac{1}{2}{C}_D\mathrm{\ρD}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^2\cos {\mathrm{\ω}}_s^{\′}{t}_0+$

$\frac{1}{2}{C}_D\mathrm{\ρD}(U-{\stackrel{.}{u}}_j^{\left(i\right)})|U-{\stackrel{.}{u}}_j^{\left(i\right)}|-\frac{{\mathrm{\π}}^2}{4}\mathrm{\ρD}{\stackrel{..}{u}}_j^{\left(i\right)}$

在以上两个公式中，

FD--脉动拖曳力；

C_D--拖曳力系数；

ρ--流体密度；

D--深水立管直径；

U--流体流速；

--深水立管顺流向振动速度；

ω′_s--涡旋泄放频率，

其中，St为斯特罗哈数；

t--时间；

ü--深水立管顺流向振动加速度；

4)将步骤3)中计算得到的脉动拖曳力

代入深水立管的振动方程式：

${m\stackrel{..}{u}}_j^{(i+1)}+c{\stackrel{.}{u}}_j^{(i+1)}+k{u}_j^{(i+1)}={\mathrm{FD}}_j^{\left(i\right)}$

式中：m--深水立管的质量；

      c--深水立管的阻尼系数；

      k--深水立管的弯曲刚度；

--第j时间步内第i+1次迭代的顺流向振动位移；

--第j时间步内第i+1次迭代的顺流向振动速度；

--第j时间步内第i+1次迭代的顺流向振动加速度；

计算第j时间步内第i+1次迭代的深水立管顺流向涡激振动的速度和加速度

(公知技术)；

5)如果ε为预先设定的计算精度(根据需要确定，如ε＝1×10^-5)，则继续进行迭代计算，令：

${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_j^{(i+1)},$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_j^{(i+1)}$

然后，重复步骤3)～5)的计算；

如果 $\max \left\{\right|{\stackrel{..}{u}}_j^{(i+1)}-{\stackrel{..}{u}}_j^{\left(i\right)}|,|{\stackrel{.}{u}}_j^{(i+1)}-{\stackrel{.}{u}}_j^{\left(i\right)}\left|\right\}\≤\mathrm{\ϵ},$ 则计算第j时间步内第i次迭代的涡激升力，令：

${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_j^{\left(n\right)},$ ${\stackrel{.}{v}}_j^{\left(i\right)}={\stackrel{\·}{v}}_j^{\left(0\right)},$ ${\stackrel{..}{v}}_j^{\left(i\right)}={\stackrel{..}{v}}_j^{\left(0\right)}$

n为第j时间步内计算脉动拖曳力的总迭代次数；

对于第1个时间步，即j＝0，

否则，

其中，m为第j-1时间步内计算涡激升力的总迭代次数；

6)将

代入下式计算第j时间步内第i次迭代的涡激升力：

${\mathrm{FL}}_j^{\left(i\right)}=\frac{1}{2}{C}_L\mathrm{\ρD}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^2\cos {\mathrm{\ω}}_s^{\′}{t}_j-\frac{1}{2}{C}_D\mathrm{\ρD}{\stackrel{.}{v}}_j^{\left(i\right)}\left|{\stackrel{.}{v}}_j^{\left(i\right)}\right|-\frac{{\mathrm{\π}}^2}{4}\mathrm{\ρD}{\stackrel{..}{v}}_j^{\left(i\right)}$

上式中各参数的含义与公式(1)中对应参数的含义相同。

7)将步骤6)中计算得到的涡激升力

代入深水立管的横流向涡激振动方程式：

${m\stackrel{..}{v}}_j^{(i+1)}+c{\stackrel{.}{v}}_j^{(i+1)}+k{v}_j^{(i+1)}={\mathrm{FL}}_j^{\left(i\right)}$

式中：m--深水立管的质量；

      c--深水立管的阻尼系数；

      k--深水立管的弯曲刚度；

--第j时间步内第i+1次迭代的横流向振动位移；

--第j时间步内第i+1次迭代的横流向振动速度；

--第j时间步内第i+1次迭代的横流向振动加速度；

计算第j时间步内第i+1次迭代的深水立管横流向涡激振动速度和加速度

(公知技术)；

8)如果

ε为预先设定的计算精度(根据需要确定，如ε＝1×10^-5)，则继续进行迭代计算，令；

${\stackrel{.}{v}}_j^{\left(i\right)}={\stackrel{.}{v}}_j^{(i+1)},$ ${\stackrel{..}{v}}_j^{\left(i\right)}={\stackrel{..}{v}}_j^{(i+1)}$

然后，重复步骤6)～8)的计算；

如果则开始下一个时间步的计算，令：

t_j＝t_j+1＝t_j+Δt， ${\stackrel{..}{v}}_j^{\left(i\right)}={\stackrel{..}{v}}_{j+1}^{\left(0\right)}={\stackrel{..}{v}}_j^{\left(m\right)},$ ${\stackrel{..}{v}}_j^{\left(i\right)}={\stackrel{..}{v}}_{j+1}^{\left(0\right)}={\stackrel{..}{v}}_j^{\left(m\right)}$

m为第j时间步内计算涡激升力的总迭代次数；Δt一般取0.02秒；

重复步骤3)～8)的计算，直至计算时长满足需要(时长根据需要设定，如100秒)。

本发明相对于现有技术具有如下特点：

1、考虑了深水立管振动的惯性耦合效应

涡激振动是深水立管在流体定常流动作用下的一种振动形式，包括顺流向振动和横流向振动，其中的横流向振动是由涡激升力引起的。由于深水立管振动加速度对流场产生的扰动作用(流体的加速度为零)，深水立管周围的流场将产生压力梯度，该压力梯度对圆柱体的作用被称为附加质量力。因此，深水立管(弹性圆柱体)的涡激升力应在刚性圆柱体涡激升力的基础上增加由圆柱体的振动加速度引起的附加质量力。

2、考虑了深水立管振动的粘性耦合效应

深水立管的横流向振动速度对流场的扰动作用使深水立管与流体之间产生了大小和方向随时间变化的相对速度(流体本身在圆柱体振动方向的速度为零)，从而导致深水立管受到周期性变化的流体粘性阻力作用。因此，深水立管(弹性圆柱体)的涡激升力应在刚性圆柱体涡激升力的基础上增加由圆柱体与流体之间的相对速度引起的交变粘性阻力。

3、解决了深水立管弯曲振动的展向流固耦合效应

对于弹性圆柱体，如深水立管和海底管线，其涡激振动为弯曲振动，因此，沿圆柱体轴向的不同位置，其振动速度和加速度是不同的。如果采用相同的涡激升力计算，则误差是可想而知的。而本发明提出的涡激升力计算方法可以解决上述问题，使计算结果更符合涡激振动的本质特征。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若对本发明的这些修改和变型属于本发明权利要求及其同等技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (2)

1.一种深水立管的时域涡激升力确定方法，其特征在于：建立的涡激升力时域模型如下：

$\mathrm{FL}=\frac{1}{2}{C}_L\mathrm{\ρD}{(U-\stackrel{.}{u})}^2\cos {\mathrm{\ω}}_s^{\′}t-\frac{1}{2}{C}_D\mathrm{\ρD}\stackrel{.}{v}\left|\stackrel{.}{v}\right|-\frac{{\mathrm{\π}}^2}{4}\mathrm{\ρD}\stackrel{..}{v}$

式中：FL--涡激升力；

      C_L--升力系数；

      C_D--拖曳力系数；

      ρ--流体密度；

      D--深水立管直径；

      U--流速；

--深水立管顺流向振动速度；

--深水立管横流向振动速度；

ω′_s--涡旋泄放频率，

其中，St为斯特罗哈数，

为圆柱体顺流向振动速度；

t--时间；

--深水立管横流向振动加速度；

针对深水立管顺流向振动速度、横流向振动速度和加速度，采用迭代方法对上式进行计算，得到深水立管的时域涡激升力。

2.如权利要求1所述的深水立管的时域涡激升力确定方法，其特征在于：所述的采用迭代方法进行计算的具体步骤如下：

1)给定深水立管顺流向振动速度、横流向振动速度和加速度及计算时间的初值：

t_j＝t₀＝0， ${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_0^{\left(0\right)}=0,$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_0^{\left(0\right)}=0,$ ${\stackrel{.}{v}}_j^{\left(i\right)}={\stackrel{.}{v}}_0^{\left(0\right)}=0,$ ${\stackrel{..}{v}}_j^{\left(i\right)}={\stackrel{..}{v}}_0^{\left(0\right)}=0$

式中：j--时间步数，计算开始时j＝0；

      i--迭代次数，每个时间步开始时i＝0；

2)计算给定流速下的约化速度：

$V_r=\frac{U}{f_{n}D}$

式中：V_r--约化速度；

      U--流体流速；

      f_n--深水立管的固有频率；

      D--深水立管直径；

3)当V_r＜5或V_r＞7时，将t_j，代入下式计算第j时间步内第i次迭代的脉动拖曳力：

${\mathrm{FD}}_j^{\left(i\right)}=\frac{1}{2}{C}_D\mathrm{\ρD}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^2\cos {\mathrm{\ω}}_s^{\′}{t}_0+$

$\frac{1}{2}{C}_D\mathrm{\ρD}(U-{\stackrel{.}{u}}_j^{\left(i\right)})|U-{\stackrel{.}{u}}_j^{\left(i\right)}|-\frac{{\mathrm{\π}}^2}{4}\mathrm{\ρD}{\stackrel{..}{u}}_j^{\left(i\right)}$

当5≤V_r≤7时，将t_j，

代入下式计算第j时间步内第i次迭代的脉动拖曳力：

${\mathrm{FD}}_j^{\left(i\right)}=\frac{1}{2}{C}_D\mathrm{\ρD}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^2\cos {\mathrm{\ω}}_s^{\′}{t}_0+$

$\frac{1}{2}{C}_D\mathrm{\ρD}(U-{\stackrel{.}{u}}_j^{\left(i\right)})|U-{\stackrel{.}{u}}_j^{\left(i\right)}|-\frac{{\mathrm{\π}}^2}{4}\mathrm{\ρD}{\stackrel{..}{u}}_j^{\left(i\right)}$

在以上两个公式中，

FD--脉动拖曳力；

C_D--拖曳力系数；

ρ--流体密度；

D--深水立管直径；

U--流体流速；

--深水立管顺流向振动速度；

ω′_s--涡旋泄放频率，

其中，St为斯特罗哈数；

t--时间；

ü--深水立管顺流向振动加速度；

4)将步骤3)中计算得到的脉动拖曳力

代入深水立管的振动方程式：

${m\stackrel{..}{u}}_j^{(i+1)}+c{\stackrel{.}{u}}_j^{(i+1)}+k{u}_j^{(i+1)}={\mathrm{FD}}_j^{\left(i\right)}$

式中：m--深水立管的质量；

      c--深水立管的阻尼系数；

k--深水立管的弯曲刚度；

--第j时间步内第i+1次迭代的顺流向振动位移；

--第j时间步内第i+1次迭代的顺流向振动速度；

--第j时间步内第i+1次迭代的顺流向振动加速度；

计算第j时间步内第i+1次迭代的深水立管顺流向涡激振动的速度

和加速度

5)如果

ε为预先设定的计算精度，则继续进行迭代计算，令：

${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{.}{u}}_j^{(i+1)},$ ${\stackrel{..}{u}}_j^{\left(i\right)}={\stackrel{..}{u}}_j^{(i+1)}$

然后，重复步骤3)～5)的计算；

如果 $\max \left\{\right|{\stackrel{..}{u}}_j^{(i+1)}-{\stackrel{..}{u}}_j^{\left(i\right)}|,|{\stackrel{.}{u}}_j^{(i+1)}-{\stackrel{.}{u}}_j^{\left(i\right)}\left|\right\}\≤\mathrm{\ϵ},$ 则令：

${\stackrel{.}{u}}_j^{\left(i\right)}={\stackrel{\·}{u}}_j^{\left(n\right)},$ ${\stackrel{.}{v}}_j^{\left(i\right)}={\stackrel{\·}{v}}_j^{\left(0\right)},$ ${\stackrel{..}{v}}_j^{\left(i\right)}={\stackrel{..}{v}}_j^{\left(0\right)}$

 n为第j时间步内计算脉动拖曳力的总迭代次数；

对于第1个时间步，即j＝0，

否则，

其中，m为第j-1时间步内计算涡激升力的总迭代次数；

6)将

代入下式计算第j时间步内第i次迭代的涡激升力：

${\mathrm{FL}}_j^{\left(i\right)}=\frac{1}{2}{C}_L\mathrm{\ρD}{(U-{\stackrel{.}{u}}_j^{\left(i\right)})}^2\cos {\mathrm{\ω}}_s^{\′}{t}_j-\frac{1}{2}{C}_D\mathrm{\ρD}{\stackrel{.}{v}}_j^{\left(i\right)}\left|{\stackrel{.}{v}}_j^{\left(i\right)}\right|-\frac{{\mathrm{\π}}^2}{4}\mathrm{\ρD}{\stackrel{..}{v}}_j^{\left(i\right)}$

7)将步骤6)中计算得到的涡激升力

代入深水立管的横流向涡激振动方程式：

${m\stackrel{..}{v}}_j^{(i+1)}+c{\stackrel{.}{v}}_j^{(i+1)}+k{v}_j^{(i+1)}={\mathrm{FL}}_j^{\left(i\right)}$

式中：m--深水立管的质量；

      c--深水立管的阻尼系数；

      k--深水立管的弯曲刚度；

--第j时间步内第i+1次迭代的横流向振动位移；

--第j时间步内第i+1次迭代的横流向振动速度；

--第j时间步内第i+1次迭代的横流向振动加速度；

计算第j时间步内第i+1次迭代的深水立管横流向涡激振动速度

和加速度

8)如果

ε为预先设定的计算精度，则继续进行迭代计算，令：

${\stackrel{.}{v}}_j^{\left(i\right)}={\stackrel{.}{v}}_j^{(i+1)},$ ${\stackrel{..}{v}}_j^{\left(i\right)}={\stackrel{..}{v}}_j^{(i+1)}$

然后，重复步骤6)～8)的计算；

如果

则开始下一个时间步的计算，令：

t_j＝t_j+1＝t_j+Δt， ${\stackrel{.}{v}}_j^{\left(i\right)}={\stackrel{.}{v}}_{j+1}^{\left(0\right)}={\stackrel{.}{v}}_j^{\left(m\right)},$ ${\stackrel{..}{v}}_j^{\left(i\right)}={\stackrel{..}{v}}_{j+1}^{\left(0\right)}={\stackrel{..}{v}}_j^{\left(m\right)}$

m为第j时间步内计算涡激升力的总迭代次数；

重复步骤3)～8)的计算，直至计算时长满足需要。

Images