CN101539477B - 一种深水顶张式立管涡激振动与疲劳分析的方法 - Google Patents

Abstract

一种深水顶张式立管涡激振动与疲劳分析的方法，涉及深水立管设计领域，其具体步骤为：步骤1，获取流场数据；步骤2，将流场数据代入立管振动方程式6和7中；步骤3，采用有限元方法求解式6、7，计算结果包括位移、速度、加速度和应力时程；步骤4，根据所述计算结果采用雨流计数法统计出一段时间内，某个幅值的应力循环次数n_i，将n_i代入疲劳损伤计算公式8中，计算出疲劳损伤。本发明所述的深水顶张式立管涡激振动与疲劳分析的方法，采用实际的管中管模型，同时考虑顺流向振动和横流向振动，提高了应力计算的准确度。

Description

一种深水顶张式立管涡激振动与疲劳分析的方法 

技术领域

本发明涉及深水立管设计领域，具体的说是一种深水顶张式立管涡激振动与疲劳分析的方法。

背景技术

VIV是圆柱体在定常流作用下的一种奇特振动现象，具体地说VIV是涡旋脱落引起的一种圆柱体振动现象，是圆柱体尾流处的交替涡旋泄放引起的振动。对于刚性圆柱体和柔性不太大的圆柱体(一般的工程结构)，其垂直于流速方向(横流向)的振幅远远大于平行于流速方向(顺流向)的振幅，而振动频率为流速方向振动频率的二分之一。因此，在现有的VIV研究中，一般忽略顺流向振动，而只研究横流向振动。VIV的数学模型和涡激升力模型仅有横流向的模型，如式1所示的Iwan尾流振子模型：

$\frac{d^2{\stackrel{\‾}{y}}_n}{{\mathrm{dt}}^2}+2{\mathrm{\ζ}}_n^T{\mathrm{\ω}}_n\frac{d{\stackrel{\‾}{y}}_n}{\mathrm{dt}}+{\mathrm{\ω}}_n^2{\stackrel{\‾}{y}}_n=\left(\frac{{\mathrm{\α}}_4}{v_n}\right)\frac{d{\stackrel{\‾}{z}}_n}{\mathrm{dt}}---\left(1\right)$

式1中：y_n为结构横流向位移模态参数；z_n为尾流流体的模态参数。

再如式2所示的锁定区涡激升力模型：

$\frac{1}{2}{C}_L\mathrm{\ρ}{\mathrm{DU}}^2\sin {\mathrm{\ω}}_{s}t---\left(2\right)$

式2中：C_L为升力系数；ρ为流体密度；D为圆柱体直径；U为流速；ω_s为Strouhal频率， ${\mathrm{\ω}}_s=\frac{2\mathrm{\πStU}}{D},$ 其中：St为斯特罗哈数(Strouhal number)；t为时间。

VIV是深水立管设计的主要问题之一，由于我国的深水开发刚刚起步，尚没有深水立管的设计分析软件，目前均采用国外的同类型软件进行设计分析，如Orcaflex和Shear7等。Orcaflex是目前被海洋工程界普遍接受的深水立管分析软件之一，它的VIV分析模块采用了尾流振子模型，计算结果与工程实测结果有较大的误差，所以，目前主要用于立管的波浪响应分析，而VIV分析通常采用Shear7。Shear7是深水立管VIV分析的专业软件，它采用工程实测数据修正计算参数，因此，计算结果与工程实测结果比较吻合，普遍为工程界所接受。Shear7采用的是频域分析方法。即：将立管的时域运动方程式3利用坐标变换式4转化为式5的模态(振型)坐标方程，式3如下：

$m\frac{d^{2}y(z,t)}{{\mathrm{dt}}^2}+c\frac{\mathrm{dy}(z,t)}{\mathrm{dt}}---\left(3\right)$

$+\mathrm{EI}\frac{d^{4}y(z,t)}{{\mathrm{dt}}^4}-T\frac{d^{2}y(z,t)}{{\mathrm{dt}}^2}=\frac{1}{2}{C}_L\mathrm{\ρ}{\mathrm{DU}}^2\sin {\mathrm{\ω}}_{s}t$

式3中：m为立管单位长度的质量；y(z，t)为立管的横流向弯曲位移，其中，z为立管轴向坐标；c为结构阻尼系数；EI和T分别为立管截面抗弯刚度和有效张力。

式4如下：

$y(z,t)=\underset{r}{\mathrm{\Σ}}{Y}_r\left(z\right)q_r\left(t\right)---\left(4\right)$

式4中：Y_r(z)为立管第r阶振型；q_r(t)为立管第r阶振型的坐标。

式5如下：

$M_r{\stackrel{\·\·}{q}}_r\left(t\right)+C_r{\stackrel{\·}{q}}_r\left(t\right)+K_r{q}_r\left(t\right)=P_r\left(t\right),(r=\mathrm{1,2},...,n)---\left(5\right)$

式5中：M_r为立管第r阶模态质量， $M_r={\∫}_0^L{Y}_r^2\left(z\right)\mathrm{mdz};$ C_r为立管第r阶模态阻尼， $C_r={\∫}_0^L{Y}_r^2\left(z\right)(c+c_a)\mathrm{dz};$ K_r为立管第r阶模态刚度， $K_r={\∫}_0^L\mathrm{EI}{\left(\frac{d^2{Y}_r\left(z\right)}{{\mathrm{dz}}^2}\right)}^2\mathrm{dz}+{\∫}_0^{L}T{\left(\frac{{\mathrm{dY}}_r\left(z\right)}{\mathrm{dz}}\right)}^2\mathrm{dz};$ P_r(t)为第r阶模态力， $P_r\left(t\right)={\∫}_0^L\frac{1}{2}{C}_L\mathrm{\ρ}{\mathrm{DU}}^2{Y}_r\left(z\right)\mathrm{dz}\sin {\mathrm{\ω}}_{s}t.$

然后采用随机振动的谱分析方法求解式5的方程式，即采用付立叶变换将时域的力函数转换为频域的谱函数，从而求得立管的响应谱。因此，Shear7的计算结果是以响应的平均值和均方差表示的，不能给出响应时程，因而不能给出立管响应的最大值。

上述两个软件的方法是目前比较有代表性的两种VIV分析方法，它们都没有考虑立管的顺流向振动问题。完全基于传统的圆柱体VIV理论和分析方法，而深水立管由于柔性较大，其VIV为高阶(＞10阶)多模态(5～20个模态)振动(一般圆柱体的VIV为1～3阶模态振动)，发明人的研究表明，对于长细比较大的深水立管，其顺流向振动幅度已达到横流向的50％(深水立管的VIV顺流向振幅与横流向振幅之比大于50％)，由于其振动频率高于横流向一倍，因此，顺流向振动引起的疲劳损伤与横流向相当。此种情况下，深水立管VIV分析应考虑顺流向振动问题是显而易见的，但现有技术给出的分析方案明显无法满足这一需要。

深水顶张式立管是由两个或三个钢管组成的管中管结构，如图1所示，由内至外依次包括油管(图1.a)、内套管(图1.b)、外套管(图1.c)。目前的分析软件均采用弯曲刚度等效的单层管模型，使得应力计算与实际结构不符。按弯曲刚度等效得到较大的等效管内径，因此，轴向应力偏大。由于外层套管仅承受外压，因此，按等效管计算的环向应力和径向应力偏小。

综上所述，现有的深水顶张式立管的VIV分析存在以下缺点：

1、没有考虑顺流向振动响应，因此，位移响应和应力响应均较小，计算结果偏于不安全。

2、尾流振子模型计算的横流向响应误差较大，主要是由于尾流振子模型将尾流处的流体看作是与固体具有相同性质(具有固定形状)的物体做相同模态的振动，并与结构相互作用。

3、涡激升力模型没有考虑流固耦合效应，因此，只适用于锁定区的VIV分析。

4、涡激升力的频率采用Strouhal频率，而在锁定区，涡旋泄放频率不符合Strouhal频率计算式。因此，涡激升力用于锁定区的VIV分析，其频率与实际的涡旋泄放频率不同，造成计算结果的不准确。

5、由于立管的外套管应力最大，因此，立管设计主要是外套管应力校核。而现有软件均采用弯曲刚度等效的单层管模型，造成轴向的内壁弯曲应力计算值偏小(内径偏差)、张力应力计算值偏大(截面积偏差)，环向应力计算值偏小(壁厚偏差)，径向应力计算值偏大(受力情况与实际不符)。

6、不能计算油管和内套管应力，因而无法校核油管和内套管的强度和疲劳寿命。

发明内容

针对现有技术中存在的缺陷，本发明的目的在于提供一种深水顶张式立管涡激振动与疲劳分析的方法，采用实际的管中管模型，同时考虑顺流向振动和横流向振动，提高了应力计算的准确度。

为达到以上目的，本发明采取的技术方案是：

一种深水顶张式立管涡激振动与疲劳分析的方法，其特征在于：其具体步骤为：

步骤1，获取流场数据；

步骤2，将流场数据代入立管振动方程式6和7

$(m+m_a)\frac{{\∂}^{2}x(z,t)}{{\∂t}^2}+(c+c_a)\frac{\∂x(z,t)}{\∂t}---\left(6\right)$

$+\mathrm{EI}\frac{{\∂}^{4}x(z,t)}{{\∂z}^4}-T\frac{{\∂}^{2}x(z,t)}{{\∂z}^2}=\frac{1}{2}{C}_D\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos {2\mathrm{\ω}}_s^{\′}t$

$(m+m_a)\frac{{\∂}^{2}y(z,t)}{{\∂t}^2}+(c+c_a)\frac{\∂y(z,t)}{\∂t}---\left(7\right)$

$+\mathrm{EI}\frac{{\∂}^{4}y(z,t)}{{\∂z}^4}-T\frac{{\∂}^{2}y(z,t)}{{\∂z}^2}=\frac{1}{2}{C}_L\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos {\mathrm{\ω}}_s^{\′}t$

中；

步骤3，采用有限元方法求解式6、7，计算结果包括位移、速度、加速度和应力时程；

步骤4，根据所述计算结果采用雨流计数法统计出一段时间内，某个幅值的应力循环次数n_i，将n_i代入疲劳损伤计算公式8

$D=\underset{i}{\mathrm{\Σ}}\frac{n_i}{N_i}---\left(8\right)$

中，计算出疲劳损伤。

本发明所述的深水顶张式立管涡激振动与疲劳分析的方法，采用实际的管中管模型，同时考虑顺流向振动和横流向振动，提高了应力计算的准确度。

附图说明

本发明有如下附图：

图1.a油管实际受力状态示意图

图1.b内套管实际受力状态示意图

图1.c外套管实际受力状态示意图

图2等效管模型受力状态示意图

图3雨流计数法示意图

具体实施方式

以下结合附图对本发明作进一步详细说明。

本发明的目的在于：

1、在深水立管的VIV分析中，考虑顺流向振动，从而使分析模型和结果与实际工程问题更吻合。

2、改进现有涡激升力模型，考虑流固耦合效应对涡激升力的影响，改善VIV分析精度。

3、改进现有涡旋泄放频率的计算模型，使其不仅适用于锁定区，也适用于非锁定区。

4、采用时域分析方法，直观地分析立管VIV的最大振幅及其应力和疲劳损伤。

5、采用立管的实际结构(管中管)计算应力，使得应力的计算符合结构的实际受力状态，则应力状态更符合实际工程结构实际的应力状态。可以计算油管和内套管的应力，从而可以对油管和内套管进行应力校核。

根据上述目的，本发明公开了一种深水顶张式立管涡激振动与疲劳分析的方法，其具体步骤为：

步骤1，获取流场数据；

步骤2，将流场数据代入立管振动方程式6和7

$(m+m_a)\frac{{\∂}^{2}x(z,t)}{{\∂t}^2}+(c+c_a)\frac{\∂x(z,t)}{\∂t}---\left(6\right)$

$+\mathrm{EI}\frac{{\∂}^{4}x(z,t)}{{\∂z}^4}-T\frac{{\∂}^{2}x(z,t)}{{\∂z}^2}=\frac{1}{2}{C}_D\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos {2\mathrm{\ω}}_s^{\′}t$

$(m+m_a)\frac{{\∂}^{2}y(z,t)}{{\∂t}^2}+(c+c_a)\frac{\∂y(z,t)}{\∂t}---\left(7\right)$

$+\mathrm{EI}\frac{{\∂}^{4}y(z,t)}{{\∂z}^4}-T\frac{{\∂}^{2}y(z,t)}{{\∂z}^2}=\frac{1}{2}{C}_L\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos {\mathrm{\ω}}_s^{\′}t$

中；

步骤3，采用有限元方法求解式6、7，计算结果包括位移、速度、加速度和应力时程；

步骤4，根据所述计算结果采用雨流计数法统计出一段时间内，某个幅值的应力循环次数n_i，将n_i代入疲劳损伤计算公式8

$D=\underset{i}{\mathrm{\Σ}}\frac{n_i}{N_i}---\left(8\right)$

中，计算出疲劳损伤。

分析深水顶张式立管涡激振动与疲劳需要流场数据(流场数据包括：流速及其沿水深的变化和不同流速发生概率，代入式6、7的右端项中字母U的位置和ω′_s计算式 ${\mathrm{\ω}}_s^{\′}=\frac{2\mathrm{\π}[U\left(z\right)-\stackrel{\·}{x}(z,t)]\mathrm{St}}{D}$ 中U的位置)，如流速及其沿水深的分布。流速沿水深的变化体现在公式6、7的右端项，即 $\frac{1}{2}{C}_L\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos {\mathrm{\ω}}_s^{\′}t$ 及其 ${\mathrm{\ω}}_s^{\′}=\frac{2\mathrm{\π}[U\left(z\right)-\stackrel{\·}{x}(z,t)]\mathrm{St}}{D}$ 中，流速U(z)中的z即为水深坐标。

不同流速发生的概率是不体现在公式中的，是计算疲劳损伤的依据。结构的疲劳分析是通过计算一段时间内(通常为一年)各种荷载引起的疲劳损伤来判断结构的疲劳强度，对于立管等海洋结构物，一般是计算一年时间内的疲劳损伤，其倒数就是疲劳寿命。而一年内，海流的流速是不断变化的，通过观测，人们可以得到不同流速发生的时段，统计后就得到了一年内，每个流速发生的总时段(小时数)，除以一年的时段(小时数)就得到了该流速的发生概率。这些数据是立管设计时的已知数据，输入后，程序在计算某个流速的疲劳损伤时，将计算结果乘以该流速的发生概率，就得到了一年内该流速引起的疲劳损伤，将所有流速引起的疲劳损伤累加起来，就是立管一年的疲劳损伤。

这些数据是油田开发时就已经通过海洋环境监测获得的，以便作为立管设计的依据，并提供给设计人员。

将流场数据代入涡激升力模型和立管振动方程中：

$(m+m_a)\frac{{\∂}^{2}x(z,t)}{{\∂t}^2}+(c+c_a)\frac{\∂x(z,t)}{\∂t}---\left(6\right)$

$+\mathrm{EI}\frac{{\∂}^{4}x(z,t)}{{\∂z}^4}-T\frac{{\∂}^{2}x(z,t)}{{\∂z}^2}=\frac{1}{2}{C}_D\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos {2\mathrm{\ω}}_s^{\′}t$

$(m+m_a)\frac{{\∂}^{2}y(z,t)}{{\∂t}^2}+(c+c_a)\frac{\∂y(z,t)}{\∂t}---\left(7\right)$

$+\mathrm{EI}\frac{{\∂}^{4}y(z,t)}{{\∂z}^4}-T\frac{{\∂}^{2}y(z,t)}{{\∂z}^2}=\frac{1}{2}{C}_L\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos {\mathrm{\ω}}_s^{\′}t$

上式中仅右端(涡激升力，即公式的右端项： $\frac{1}{2}{C}_L\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos {\mathrm{\ω}}_s^{\′}t)$ 的涡激升力与这些数据有关，即公式中的流速一项。流速沿水深的变化体现在公式6、7的右端项 $\frac{1}{2}{C}_D\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos {2\mathrm{\ω}}_s^{\′}t$ 和 $\frac{1}{2}{C}_L\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos {\mathrm{\ω}}_s^{\′}t$ 中，而流速的发生概率是在上述方程计算结束后，计算疲劳寿命时需要的数据，具体过程前面已经介绍，这里不再赘述。

计算过程是解微分方程，采用有限元数值解来求解顺流向和横流向两个动力方程。

计算结果包括位移、速度、加速度和应力时程。式6、7是以公知的计算公式为基础，本发明没有改变公式的基本结构，只是加入了一些修正参数，如m_a、C_a。下面以式7为例给出推导过程：

令：y(z，t)＝Y(z)q(t)                        (a)

其中：Y(z)为单元插值函数，q(t)为结点位移函数。代入式7得：

$(m+m_a)Y\left(z\right)\stackrel{\·\·}{q}\left(t\right)+(c+c_a)Y\left(z\right)\stackrel{\·}{q}\left(t\right)$

$+\mathrm{EI}\·Y^{\mathrm{IV}}\left(z\right)q\left(t\right)-T\·Y^{\′\′}\left(z\right)q\left(t\right)=\frac{1}{2}{C}_L\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos {\mathrm{\ω}}_s^{\′}t---\left(b\right)$

(b)式两端同时乘以Y(z)并沿单元积分：

$(m+m_a){\∫}_0^{l}Y\left(z\right)Y\left(z\right)\mathrm{dz}\stackrel{\·\·}{q}\left(t\right)+(c+c_a){\∫}_0^{l}Y\left(z\right)Y\left(z\right)\mathrm{dz}\stackrel{\·}{q}\left(t\right)$

$+[\mathrm{EI}{\∫}_0^{l}Y\left(z\right)Y^{\mathrm{IV}}\left(z\right)\mathrm{dz}-T{\∫}_0^{l}Y\left(z\right)Y^{\′\′}\left(z\right)\mathrm{dz}]q\left(t\right)---\left(c\right)$

$=\frac{1}{2}{C}_L\mathrm{\ρ}D\cos {\mathrm{\ω}}_s^{\′}t{\∫}_0^{l}Y\left(z\right){[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\mathrm{dz}$

(c)式可表示为： $M\stackrel{\·\·}{q}\left(t\right)+C\stackrel{\·}{q}\left(t\right)+\mathrm{Kq}\left(t\right)=F\left(t\right)---\left(d\right)$

其中： $M=(m+m_a){\∫}_0^{l}Y\left(z\right)Y\left(z\right)\mathrm{dz}$

$C=(c+c_a){\∫}_0^{l}Y\left(z\right)Y\left(z\right)\mathrm{dz}$

$K=\mathrm{EI}{\∫}_0^{l}Y\left(z\right)Y^{\mathrm{IV}}\left(z\right)\mathrm{dz}-T{\∫}_0^{l}Y\left(z\right)Y^{\′\′}\left(z\right)\mathrm{dz}$

$F\left(t\right)=\frac{1}{2}{C}_L\mathrm{\ρ}D\cos {\mathrm{\ω}}_s^{\′}t{\∫}_0^{l}Y\left(z\right){[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\mathrm{dz}$

(d)式的求解可采用Newmark-β法，这是一个公知的方法。(d)式的计算结果为立管的位移、速度和加速度。将位移计算结果代入(a)式即可得到立管的位移响应。然后根据梁弯曲的应力计算公式：

$\mathrm{\σ}=E\frac{d^{2}y(z,t)}{d{z}^2}\frac{D}{2}---\left(e\right)$

计算出立管横截面的最大弯曲应力。

然后根据应力计算结果采用雨流计数法(雨流计数法是公知的方法，它需要应力时程(如图3所示)数据)，统计出一段时间内(图3为100秒)，某个幅值的应力循环次数n_i，这是雨流计数法的计算结果，再代入疲劳损伤计算公式(公知的公式)：

$D=\underset{i}{\mathrm{\Σ}}\frac{n_i}{N_i}---\left(8\right)$

计算出疲劳损伤。式8中：N_i为第i个应力循环幅值对应的材料疲劳破坏的循环次数，由材料的S-N曲线查得。材料的S-N曲线给出的是材料的疲劳应力与循环次数的关系，是传统的疲劳分析工具和方法，俗称为S-N曲线法。

本发明所述的深水顶张式立管涡激振动与疲劳分析的方法，采用了式6和式7所示的同时考虑顺流向振动和流固耦合的VIV分析模型：

$(m+m_a)\frac{{\∂}^{2}x(z,t)}{{\∂t}^2}+(c+c_a)\frac{\∂x(z,t)}{\∂t}---\left(6\right)$

$+\mathrm{EI}\frac{{\∂}^{4}x(z,t)}{{\∂z}^4}-T\frac{{\∂}^{2}x(z,t)}{{\∂z}^2}=\frac{1}{2}{C}_D\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos 2{\mathrm{\ω}}_s^{\′}t$

$(m+m_a)\frac{{\∂}^{2}y(z,t)}{{\∂t}^2}+(c+c_a)\frac{\∂y(z,t)}{\∂t}---\left(7\right)$

$+\mathrm{EI}\frac{{\∂}^{4}y(z,t)}{{\∂z}^4}-T\frac{{\∂}^{2}y(z,t)}{{\∂z}^2}=\frac{1}{2}{C}_L\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos {\mathrm{\ω}}_s^{\′}t$

式6、7中：x(z，t)和y(z，t)分别为顺流向和横流向振动响应函数；z为立管轴向坐标；t为时间；x(z，t)、y(z，t)、z和t四个参数为非输入参数，是程序计算的参数和坐标系统；m和c分别为立管单位长度的质量和结构阻尼系数，立管单位长度的质量是根据结构设计计算得到的，结构阻尼系数是经验系数，对于一定的材料，其值有一个取值范围，根据经验确定；m_a和c_a分别为附加质量系数和附加阻尼系数，附加质量系数和附加阻尼系数可由公知的Morison公式的相应参数计算得到；EI和T分别为立管截面抗弯刚度和有效张力，EI和T可由结构设计得到；C_D和C_L分别为拖曳力系数和升力系数，C_D和C_L为经验系数，可由相应的图表查得；ρ为流体密度，有表可查；U为流速，设计的基础数据，由油田的开发商提供；

为立管顺流向振动速度，由式6计算得到，因此，式6、7的求解顺序是先解6式，再解7式；ω′_s表示考虑流固耦合的涡旋泄放频率， ${\mathrm{\ω}}_s^{\′}=\frac{2\mathrm{\π}[U\left(z\right)-\stackrel{\·}{x}(z,t)]\mathrm{St}}{D},$ St为斯特罗哈数，D为立管外径。

由上述可见，式6、7所示的模型考虑了流固耦合效应对涡激升力和涡旋泄放频率的影响。

本发明给出的上述分析模型与现有分析模型的区别在于：

1、通过附加质量系数m_a和附加阻尼系数c_a来考虑流固耦合效应；

2、通过采用流体与立管的相对速度(式6和式7右端项的来考虑流固耦合对涡激升力和涡旋泄放频率的影响，这与现有的涡激升力模型：

不同，现有的涡激升力模型中采用的是流体的绝对速度U和Strouhal频率： ${\mathrm{\ω}}_s=\frac{2\mathrm{\π}\·U\·\mathrm{St}}{D};$

3、采用式6来计算顺流向的VIV。

本发明的应力计算与现有软件的计算方法有较大的区别，现有软件按等效管模型计算立管的应力，与立管的实际受力状态和应力状态不符。因为立管的最内层管是油管，内套管与油管之间的环形空间是气举线(采集天然气)，而外套管是隔水管，它与内套关之间的环形空间是常压空气。因此，油管和套管的受力是不同的，如图1所示，油管受内压p_O(原油压力)、外压p_G(天然气压力)作用，内套管仅受内压p_G(天然气压力)作用，而外套管仅受外压p_w(海水压力)作用。而等效管模型则采用式8的弯曲刚度等效：

EI_eq＝EI_tu+EI_in+EI_ou    (8)

式8中，EI_eq为等效管模型截面抗弯刚度；EI_tu为油管截面抗弯刚度；EI_in为内套管截面抗弯刚度；EI_ou为外套管截面抗弯刚度。等效管模型将三根管等效为一根管，计算应力时按照一根管同时受内压(原油压力)和外压(海水压力)作用计算，参见图2。应力计算不仅与外力有关，也与立管的壁厚有关。而等效管模型的壁厚和受力状态都与管中管不同，因此，计算出的应力与管中管不同，从而与实际工程不符。

本发明与现有软件计算方法不同之处在于：

1、现有软件没有系数m_a和c_a，本发明通过附加质量系数(式7中的m_a)和附加阻尼系数(式7中的c_a)考虑了横流向的流固耦合效应。流固耦合问题通常包含在流体荷载中，而VIV的横向振动是垂直于流动方向的振动，因此，其荷载中并不包括横向的流固耦合效应。本发明考虑到这个问题，在振动方程中增加了由于流固耦合而产生的附加质量和附加阻尼，从而解决了VIV横流向的流固耦合问题。

2、传统的VIV分析方法和软件仅采用式7计算横流向振动，而本发明增加了式6计算顺流向振动。

3、本发明的涡激升力模型通过采用流体与立管的相对速度考虑了流固耦合效应，即式6、7右端项的

4、本发明的涡激升力模型采用了考虑流固耦合的涡旋泄放频率，即采用流体与立管的相对速度计算式6、7右端项的涡旋泄放频率 ${\mathrm{\ω}}_s^{\′}=\frac{2\mathrm{\π}[U\left(z\right)-\stackrel{\·}{x}(z,t)]\mathrm{St}}{D},$ ω′_s表示涡旋泄放频率，它就是涡激升力的频率。结构振动响应的大小不仅与结构的固有频率有关，也与外荷载的频率(扰力频率)有关，因此，ω′_s对于立管的涡激振动响应是非常重要的参数，也是本发明的主要特点之一。原来的表达式为 ${\mathrm{\ω}}_s=\frac{2\mathrm{\πUSt}}{D},$ 称为斯特罗哈(Strouhal)频率。两个表达式的区别在于：本发明的表达式中包含了立管振动的速度

它是一个与水深和时间有关的参数，即式6中的方程变量。从而把原表达式的流体绝对速度U改为流体与结构的相对速度

即流固耦合效应。表达式中的St是斯特罗哈数(Strouhal number)，是一个经验系数，在10³＜Re＜10⁵范围(Re为雷诺数)，取0.2。表达式中的D为立管直径，是结构设计时确定的，在本发明中作为已知数输入。

5、本发明采用外套管的实际内外径计算弯曲应力，采用三根管的截面积之和计算张力引起的应力，而现有软件采用弯曲刚度等效的单层管模型计算弯曲应力和张力应力。张力引起的应力等于立管的有效张力除以立管的截面积，即： $\mathrm{\σ}=\frac{T}{A}.$ 由于立管的有效张力是施加在三根管上的，因此，上式的A应为三根管的截面积之和。

立管轴向应力是由立管的有效张力和弯矩以及立管的截面几何性质决定的，而按照弯曲刚度等效得到的等效管内径大于按截面积等效的等效管内径。因此，等效管的截面积小于三根管的实际截面积。从而使张力引起的应力计算结果偏大，而等效管内径小于外套管内径，因此，内壁弯曲应力的计算结果远小于外套管的内壁弯曲应力，用于校核外套管应力是偏于不安全的。

6、本发明采用外套管内、外径和实际受力状态(仅受外压)计算径向应力和环向应力，而现有软件则采用等效管的内外径和与实际不符的受力状态(同时承受内、外压)计算径向和环向应力。

7、本发明采用内套管内、外径和实际受力状态(仅受内压)计算内套管的应力，而现有软件不能计算内套管应力。

8、本发明采用油管内、外径和实际受力状态(受内、外压)计算油管的应力，而现有软件不能计算油管应力。

立管是由油管和套管组成的，根据套管是单层套管或双层套管，立管又分为单屏立管(单层套管)和双屏立管(双层套管)。国外现有软件只能计算外套管的应力而不能计算内套管和油管应力，因此，不能校核内套管和油管的强度与疲劳寿命。本发明增加了油管和内套管的应力计算。

本发明与现有VIV分析软件相比具有下述优点：

1、增加了VIV顺流向振动分析计算，比现有软件仅计算横流向VIV更加符合深水立管的运动和受力状态。

2、采用立管实际结构的几何参数和受力状态计算立管的各应力分量，比现有软件采用等效管的几何参数和受力状态(与实际结构不符)计算出的应力更准确。

3、本发明不仅可以计算外套管的强度和疲劳损伤，而且可以计算油管和内套管的强度和疲劳损伤，弥补了现有软件不能计算油管和内套管应力和疲劳损伤的不足。

4、在VIV运动方程中增加了附加质量系数和附加阻尼系数，考虑了流固耦合效应对立管VIV的影响，使VIV分析更加符合工程实际。

5、将涡激升力中的流体速度项改为流体与立管的相对速度，从而考虑了流固耦合效应对涡激升力的影响。

6、将Strouhal频率中的流体速度项改为流体与立管的相对速度，从而考虑了流固耦合效应对涡旋泄放频率的影响，使涡旋泄放频率与实验结果吻合得更好。

Claims (1)

1.一种深水顶张式立管涡激振动与疲劳分析的方法，其特征在于：其具体步骤为：

步骤1，获取流场数据，所述流场数据包括：流速及其沿水深的变化和不同流速发生概率；

步骤2，将流场数据代入立管振动方程式6

$(m+m_a)\frac{{\∂}^{2}x(z,t)}{\∂t^2}+(c+c_a)\frac{\∂x(z,t)}{\∂t}$

$+\mathrm{EI}\frac{{\∂}^{4}x(z,t)}{\∂z^4}-T\frac{{\∂}^{2}x(z,t)}{\∂z^2}=\frac{1}{2}{C}_D\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos {2\mathrm{\ω}}_s^{\′}t$

和7

$(m+m_a)\frac{{\∂}^{2}y(z,t)}{\∂t^2}+(c+c_a)\frac{\∂y(z,t)}{\∂t}$

$+\mathrm{EI}\frac{{\∂}^{4}y(z,t)}{\∂z^4}-T\frac{{\∂}^{2}y(z,t)}{\∂z^2}=\frac{1}{2}{C}_L\mathrm{\ρD}{[U\left(z\right)-\stackrel{\·}{x}(z,t)]}^2\cos {\mathrm{\ω}}_s^{\′}t$

中；

方程式6、7中，m和c分别为立管单位长度的质量和结构阻尼系数，m_a和c_a分别为附加质量系数和附加阻尼系数，x(z，t)和y(z，t)分别为顺流向和横流向振动响应函数，t为时间，z为立管轴向坐标，EI和T分别为立管截面抗弯刚度和有效张力，C_D和C_L分别为拖曳力系数和升力系数，D为立管外径，U(z)为流速，

为立管顺流向振动速度，表示考虑流固耦合的涡旋泄放频率，

步骤3，采用有限元方法求解式6、7，计算结果包括位移、速度、加速度和应力时程；

步骤4，根据所述计算结果采用雨流计数法统计出一段时间内，某个幅值的应力循环次数n_i，将n_i代入疲劳损伤计算公式8

$D=\underset{i}{\mathrm{\Σ}}\frac{n_i}{N_i}$

中，计算出疲劳损伤，

公式8中，N_i为第i个应力循环幅值对应的材料疲劳破坏的循环次数，n_i为某个幅值的应力循环次数。

Images