CN108491615B - 一种三层顶张式立管动力响应的有限元分析方法 - Google Patents

Abstract

本发明公开了一种基于有限元分析的三层顶张紧式海洋立管的全耦合动力响应计算方法，属于深水油气开发工程管道技术分析领域。本发明采用有限元分析方法，分别建立三层管的有限元模型，将单元模型简化为欧拉梁单元；分别计算三个不同单元的质量矩阵与刚度矩阵，分别组合装配三层立管的质量矩阵与刚度矩阵；根据立管的质量矩阵与刚度矩阵，分别分析计算立管的频率与振型，进而选取扶正器的布置形式；根据立管所在海域环境，计算作用在外管上的水平拖曳力与水平惯性力；根据计算得出的油管与内管的横向位移响应，计算油管与内管各节点的应力、弯矩以及扶正器的横向作用力。本发明能正确分析出多层顶张紧式立管的横向受力与横向位移响应，具有概念明确、计算精度高等特点。

Description

一种三层顶张式立管动力响应的有限元分析方法

技术领域

本发明公开了一种基于有限元分析的三层顶张紧式海洋立管的全耦合动力响应计算方法，属于深水油气开发工程管道技术分析领域。

背景技术

深水顶张式立管是联系海底油井与水面设施的纽带与桥梁，立管为海底油井与水面设施之间提供油、气、水的输送，是深水油气开发的特有装备，尽管浅水油气开发装备中也有立管，但与深水立管不可同日而语，其结构形式与水动力性能几乎发生了质的变化。

深水顶张式立管由油管和套管组成，油管与套管最常用的方式为管中管形式，最外层管为外管，外管的作用为保护内管与油管免受海流载荷、隔绝海水腐蚀、支撑油管的作用，在深水采油的过程中，伴随着天然气的产生，内管与油管之间形成的通道用来输送天然气，油管的作用为连接水面采油树和海底井口、输送石油。

目前的顶张式立管设计的分析方法都是按照截面刚度等效的方法将油管、内管和套管等效为外径与外套管相同(保证水动力荷载相等)的单层管来计算的，但是基于拉压刚度(EA)等效与弯曲刚度(EI)等效得到的单层管截面参数是不同的，因此，采用等效管替代三层管耦合的计算结果与立管各功能管的实际受力状态是有差异的，采用等效管的方法不能模拟出各层管的实际受力状态、实际横向位移、各层管之间的实际耦合作用。

因此，提供一种正确解决各层立管之间耦合作用的方法对于深水油气开发具有重要的意义。

发明内容

本发明的主要目的是提供一种正确解决各层立管之间全耦合动力响应分析的方法，以解决现有技术不能很好模拟各层管之间全耦合作用的问题。

本发明的三层顶张式立管动力响应的有限元分析方法依次按以下步骤进行：

(1)首先根据下列公式1建立立管的动力特性方程表达式：

式中：n＝1,2,3分别表示油管、内管、外管；

[M_n]、[K_n]分别表示第n层管的质量矩阵与刚度矩阵；

ω_n为第n层管的固有频率；

为第n层管的振型向量；

(2)根据下列公式2建立立管的动力响应方程表达式：

式中：n＝1,2,3分别表示油管、内管、外管；

[M_n]、[C_n]、[K_n]分别表示第n层管的质量矩阵、阻尼矩阵与刚度矩阵；

{F}_t为n层立管所受外力的大小；

(3)建立立管的有限元模型，根据下列公式求解每根立管的单元质量矩阵与刚度矩阵，进而装配形成每根立管的质量矩阵[M_n]与刚度矩阵[K_n]：

其中，N₁N₂N₃N₄的表达式如下：

N₁(x)＝1-3γ²+2γ³ N₂(x)＝(γ-2γ²+γ³)l

N₃(x)＝3γ²-2γ³ N₄(x)＝(-γ²+γ³

其中，γ的表达式如下：

式中：m_nk_n为单元的质量矩阵与刚度矩阵；

l_n为第n层管划分的单元长度；

ρ_n、A_n分别为单元的密度与截面积，E_n、I_n分别为单元的弹性模量与截面惯性矩，当单元内的质量密度ρ_nA_n和刚度E_nI_n为常量时，单元质量矩阵与刚度矩阵可以直接表达出来；

N_i(x)为Hermite形函数，该函数可保证相邻单元在边界上横向位移与转角连续；

x_i为单元节点的起点位置，i表示单元编号；

(4)根据下式计算每层立管的阻尼矩阵：

[C_n]＝a₀[M_n]+a₁[K_n]

式中：δ为阻尼比；ω_ni、ω_nj分别为第n层管的i阶、j阶频率；a_n0和a_n1为第n层管的瑞雷阻尼系数；

(5)根据立管实际工作所处海域的海况，由以下公式计算出外层立管所受的海流作用力、顺流向涡激力的频率：

式中：C_D′为顺流向脉动拖曳力系数,取C_D′＝0.05～0.1；C_D为拖曳力系数，取C_D＝0.6～2.0；C_M为惯性力系数，取C_M＝1.0～2.0；ρ为外管密度；D为外管直径；u为海水流速；

为外管顺流向振动速度；

为外管顺流向振动加速度；

ω_s为涡泄圆频率；S_t为斯托哈尔数，取S_t＝0.18～0.2；

(6)根据步骤(1)计算得出的各层管的固有频率与步骤(5)计算得出的顺流向涡激力的频率，找出各层管与顺流向涡激力的频率相近的固有频率与阵型。

(7)外层管与内层管之间，在外层管阵型峰值点处设置扶正器，内层管与油管之间，在内层管峰值处设置扶正器；

(8)在进行立管耦合动力分析前，应确定立管两端的链接情况，由此决定剔除立管质量模型矩阵与刚度模型矩阵中的哪些元素，当为铰接时，应剔除相应节点处位移矩阵元素，当为固定链接时，应剔除相应节点处位移与转角矩阵元素；

(9)三层管之间布置有扶正器，在扶正器处，当管与管发生碰撞时，扶正器处相邻管节点的位移、速度、加速度都是相同的，扶正器只传递横向作用力，通过变换立管质量矩阵、阻尼矩阵、刚度矩阵与外力矩阵，就可以进行耦合动力分析，以下公式给出各层管所受横向等效力的大小：

式中：{F_3i}、{F_2i}、{F_1i}分别表示油管、内管与外管在单元节点i处所受等效横向作用力的大小；{f_d(2i)}表示在节点i处，油管与内管之间扶正器横向作用力的大小；{f_c(2i)}表示在节点i处，外管与内管之间扶正器横向作用力的大小；{f_x(2i)}表示在节点i处，外管所受海流横向作用力的大小；

(10)通过联立求解步骤(9)中的三个公式，可以得出最外层管的动力响应表达式：

(11)由于在扶正器处，相邻两层管的横向响应同步，由步骤(10)可知，可通过将油管与内管扶正器处单元节点i处的质量矩阵、阻尼矩阵与刚度矩阵的2i行与2i列加到外管上，再将油管与内管质量矩阵、阻尼矩阵、刚度矩阵的2i行与2i列剔除，形成耦合矩阵进行求解；

(12)根据以下公式进行三层顶张紧式立管的耦合动力分析：

式中：M₁、M₂、M₃分别为外管、内管与油管的质量矩阵；C₁、C₂、C₃分别为外管、内管与油管的阻尼矩阵；K₁、K₂、K₃分别为外管、内管与油管的刚度矩阵；{F₁}、{F₂}、{F₃}分别为外管、内管与油管的受力；

分别为外管、内管、油管的横向响应加速度；

为横向响应速度；x₁、x₂、x₃为横向响应位移；

(13)通过步骤(12)可以求得外管的节点位移，用插值法求解单元在外分布载荷作用下的最大弯曲变形，以此来判断外层立管与内层立管是否会发生碰撞；

(14)通过步骤(12)可以求得三根立管随时间的动力响应，由此来确定各层管各个单元节点的应力与弯矩，进行疲劳分析。

(15)公式3中，在计算外管的单元质量矩阵时，应考虑附加质量的影响：

式中：C_a为附加惯性力系数；ρ_w为海水密度；

(16)公式4中，在计算外管的单元刚度矩阵时，应考虑顶张紧力的影响：

式中：[K_G]^e为单元几何刚度矩阵；T为单元所受顶张紧力的大小。

本发明通过建立各层管的有限元模型，进而通过组合与变换各层管的总质量矩阵、阻尼矩阵、刚度矩阵和外载荷矩阵实现各层管之间的全耦合作用，正确的解决了现有技术存在的不足与缺陷，可以简单、准确、高效得计算出各层管的全耦合动力响应。

附图说明

图1本发明三层立管纵向结构示意图。

图2本发明三层立管力的传递示意图。

图3通过本发明方法计算得出的耦合动力响应结果图。

具体实施方式

下面通过附图和实施例对本发明作进一步详细说明，但本发明的保护范围不局限于所述内容。

实施例1：一种三层顶张式立管动力响应的有限元分析方法

具体内容如下：

海域水深300m，波高3m，波浪周期8s，海流顶部流速0.8m/s，立管底部海流流速0.2m/s。

表1立管设计参数

表2立管设计物理参数

(1)首先根据下列公式1建立立管的动力特性方程表达式：

式中：n＝1,2,3分别表示油管、内管、外管；

[M_n]、[K_n]分别表示第n层管的质量矩阵与刚度矩阵；

ω_n为第n层管的固有频率；

为第n层管的振型向量；

(2)根据下列公式2建立立管的动力响应方程表达式：

式中：n＝1,2,3分别表示油管、内管、外管；

[M_n]、[C_n]、[K_n]分别表示第n层管的质量矩阵、阻尼矩阵与刚度矩阵；

{F}_t为n层立管所受外力的大小；

(3)建立立管的有限元模型，根据下列公式求解每根立管的单元质量矩阵与刚度矩阵，进而装配形成每根立管的质量矩阵[M_n]与刚度矩阵[K_n]：

其中，N₁N₂N₃N₄的表达式如下：

N₁(x)＝1-3γ²+2γ³ N₂(x)＝(γ-2γ²+γ³)l

N₃(x)＝3γ²-2γ³ N₄(x)＝(-γ²+γ³)l

其中，γ的表达式如下：

式中：m_nk_n为单元的质量矩阵与刚度矩阵；

l_n为第n层管划分的单元长度；

ρ_n、A_n分别为单元的密度与截面积，E_n、I_n分别为单元的弹性模量与截面惯性矩，当单元内的质量密度ρ_nA_n和刚度E_nI_n为常量时，单元质量矩阵与刚度矩阵可以直接表达出来；

N_i(x)为Hermite形函数，该函数可保证相邻单元在边界上横向位移与转角连续；

x_i为单元节点的起点位置，i表示单元编号；

(4)根据下式计算每层立管的阻尼矩阵：

[C_n]＝a₀[M_n]+a₁[K_n]

式中：δ为阻尼比；ω_ni、ω_nj分别为第n层管的i阶、j阶频率；a_n0和a_n1为第n层管的瑞雷阻尼系数；

(5)根据立管实际工作所处海域的海况，由以下公式计算出外层立管所受的海流作用力、顺流向涡激力的频率：

式中：C_D′为顺流向脉动拖曳力系数,取C_D′＝0.05～0.1；C_D为拖曳力系数，取C_D＝0.6～2.0；C_M为惯性力系数，取C_M＝1.0～2.0；ρ为外管密度；D为外管直径；u为海水流速；

为外管顺流向振动速度；

为外管顺流向振动加速度；

ω_s为涡泄圆频率；S_t为斯托哈尔数，取S_t＝0.18～0.2；

(6)根据步骤(1)计算得出的各层管的固有频率与步骤(5)计算得出的顺流向涡激力的频率，找出各层管与顺流向涡激力的频率相近的固有频率与阵型。

(7)外层管与内层管之间，在外层管阵型峰值点处设置扶正器，内层管与油管之间，在内层管峰值处设置扶正器；

(8)在进行立管耦合动力分析前，应确定立管两端的链接情况，由此决定剔除立管质量模型矩阵与刚度模型矩阵中的哪些元素，当为铰接时，应剔除相应节点处位移矩阵元素，当为固定链接时，应剔除相应节点处位移与转角矩阵元素；

(9)如附图2所示，三层管之间布置有扶正器，在扶正器处，当管与管发生碰撞时，扶正器处相邻管节点的位移、速度、加速度都是相同的，扶正器只传递横向作用力，通过变换立管质量矩阵、阻尼矩阵、刚度矩阵与外力矩阵，就可以进行耦合动力分析，以下公式给出各层管所受横向等效力的大小：

式中：{F_3i}、{F_2i}、{F_1i}分别表示油管、内管与外管在单元节点i处所受等效横向作用力的大小；{f_d(2i)}表示在节点i处，油管与内管之间扶正器横向作用力的大小；{f_c(2i)}表示在节点i处，外管与内管之间扶正器横向作用力的大小；{f_x(2i)}表示在节点i处，外管所受海流横向作用力的大小；

(10)通过联立求解步骤(9)中的三个公式，可以得出最外层管的动力响应表达式：

(11)由于在扶正器处，相邻两层管的横向响应同步，由步骤(10)可知，可通过将油管与内管扶正器处单元节点i处的质量矩阵、阻尼矩阵与刚度矩阵的2i行与2i列加到外管上，再将油管与内管质量矩阵、阻尼矩阵、刚度矩阵的2i行与2i列剔除，形成耦合矩阵进行求解；

(12)根据以下公式进行三层顶张紧式立管的耦合动力分析：

式中：M₁、M₂、M₃分别为外管、内管与油管的质量矩阵；C₁、C₂、C₃分别为外管、内管与油管的阻尼矩阵；K₁、K₂、K₃分别为外管、内管与油管的刚度矩阵；{F₁}、{F₂}、{F₃}分别为外管、内管与油管的受力；

分别为外管、内管、油管的横向响应加速度；

为横向响应速度；x₁、x₂、x₃为横向响应位移；

由于上式耦合求解过程过于繁琐，所以我们应当借助数学软件求解，Matlab软件具有解复杂方程且能绘制复杂函数曲线的功能，推荐使用其对模型求解与图像绘制，求解程序如下：

dsp＝zeros(sdof,nt/dt)；

vel＝zeros(sdof,nt/dt)；

acc＝zeros(sdof,nt/dt)；

dsp(:,1)＝dsp0；

vel(:,1)＝vel0；

acc(:,1)＝acc0；

ekk＝mm*4/(dt^2)+cc*2/dt+kk；

for in＝1:nt/dt

efd＝fd(:,in)+(mm*4/(dt^2)+cc*2/dt)*dsp(:,in)+(mm*4/dt+cc)*vel(:,in)+mm*acc(:,in)；

dsp(:,in+1)＝ekk\efd；

vel(:,in+1)＝2*dsp(:,in+1)/dt-2*dsp(:,in)/dt-vel(:,in)；

acc(:,in+1)＝4*dsp(:,in+1)/(dt^2)-4*dsp(:,in)/(dt^2)-4*vel(:,in)/dt-acc(:,in)；

end

(13)通过步骤(12)可以求得外管的节点位移，用插值法求解单元在外分布载荷作用下的最大弯曲变形，以此来判断外层立管与内层立管是否会发生碰撞；

(14)如附图3所示，通过步骤(12)可以求得三根立管随时间的动力响应，由此来确定各层管各个单元节点的应力与弯矩，进行疲劳分析。

(15)公式3中，在计算外管的单元质量矩阵时，应考虑附加质量的影响：

式中：C_a为附加惯性力系数；ρ_w为海水密度；

(16)公式4中，在计算外管的单元刚度矩阵时，应考虑顶张紧力的影响：

式中：[K_G]^e为单元几何刚度矩阵；T为单元所受顶张紧力的大小。

由附图3所示，该发明方法能够准确分析出三层立管的耦合动力响应，得出三层管之间合理的扶正器布置方式与每层管的横向位移响应曲线，进而分析出立管之间是否会发生碰撞以及疲劳分析。

Claims (2)

1.三层顶张式立管动力响应的有限元分析方法，其特征在于，包括以下步骤：

(1)首先根据下列公式(1)建立立管的动力特性方程表达式：

式中：n＝1，2，3分别表示油管、内管、外管；

[M_n]、[K_n]分别表示第n层管的质量矩阵与刚度矩阵；

ω_n为第n层管的固有频率；

为第n层管的振型向量；

(2)根据下列公式(2)建立立管的动力响应方程表达式：

式中：n＝1，2，3分别表示油管、内管、外管；

[M_n]、[C_n]、[K_n]分别表示第n层管的质量矩阵、阻尼矩阵与刚度矩阵；

{F}_t为n层立管所受外力的大小；

(3)建立立管的有限元模型，根据下列公式求解每根立管的单元质量矩阵与刚度矩阵，进而装配形成每根立管的总质量矩阵[M_n]与总刚度矩阵[K_n]：

其中，N₁、N₂、N₃、N₄的表达式如下：

N₁(x)＝1-3γ²+2γ³ (5)；

N₂(x)＝(γ-2γ²+γ³)l (6)；

N₃(x)＝3γ²-2γ³ (7)；

N₄(x)＝(-γ²+γ³)l (8)；

其中，γ的表达式如下：

式中：m_n、k_n为单元的质量矩阵与刚度矩阵；

l_n为第n层管划分的单元长度；

ρ_n、A_n分别为单元的密度与截面积，E_n、I_n分别为单元的弹性模量与截面惯性矩，当单元内的质量密度ρ_nA_n和刚度E_nI_n为常量时，单元质量矩阵与刚度矩阵可以直接表达出来；

N_i(x)为Hermite形函数，该函数可保证相邻单元在边界上横向位移与转角连续；

z_i为单元节点的起点位置，i表示单元编号；

(4)根据下式计算每层立管的阻尼矩阵：

[C_n]＝a₀[M_n]+a₁[K_n] (10)；

式中：δ为阻尼比；ω_ni、ω_nj分别为第n层管的i阶、j阶频率；a_n0和a_n1为第n层管的瑞雷阻尼系数；

(5)根据立管实际工作所处海域的海况，由以下公式计算出外层立管所受的海流作用力、顺流向涡激力的频率：

式中：C_D′为顺流向脉动拖曳力系数，取C_D′＝0.05～0.1；C_D为拖曳力系数，取C_D＝0.6～2.0；C_M为惯性力系数，取C_M＝1.0～2.0；ρ为外管密度；D为外管直径；u为海水流速；

为外管顺流向振动速度；

为外管顺流向振动加速度；

ω_s为涡泄圆频率；S_t为斯托哈尔数，取S_t＝0.18～0.2；

(6)根据步骤(1)计算得出的各层管的固有频率与步骤(5)计算得出的顺流向涡激力的频率，找出各层管与顺流向涡激力的频率相近的固有频率与阵型；

(7)外层管与内层管之间，在外层管阵型峰值点处设置扶正器，内层管与油管之间，在内层管峰值处设置扶正器；

(8)在进行立管耦合动力分析前，应确定立管两端的链接情况，由此决定剔除立管质量模型矩阵与刚度模型矩阵中的哪些元素，当为铰接时，应剔除相应节点处位移矩阵元素，当为固定链接时，应剔除相应节点处位移与转角矩阵元素；

(9)三层管之间布置有扶正器，在扶正器处，当管与管发生碰撞时，扶正器处相邻管节点的位移、速度、加速度都是相同的，扶正器只传递横向作用力，通过变换立管质量矩阵、阻尼矩阵、刚度矩阵与外力矩阵，就可以进行耦合动力分析，以下公式给出各层管所受等效横向力的大小：

式中：{F_3i}、{F_2i}、{F_1i}分别表示油管、内管与外管在单元节点i处所受等效横向作用力的大小；{f_d(2i)}表示在节点i处，油管与内管之间扶正器横向作用力的大小；{f_c(2i)}表示在节点i处，外管与内管之间扶正器横向作用力的大小；{f_x(2i)}表示在节点i处，外管所受海流横向作用力的大小；

(10)通过联立求解步骤(9)中的三个公式，可以得出最外层管的动力响应表达式：

(11)由于在扶正器处，相邻两层管的横向响应同步，由步骤(10)可知，可通过将油管与内管扶正器处单元节点i处的质量矩阵、阻尼矩阵与刚度矩阵的2i行与2i列加到外管上，再将油管与内管质量矩阵、阻尼矩阵、刚度矩阵的2i行与2i列剔除，形成耦合矩阵进行求解；

(12)根据以下公式进行三层顶张紧式立管的耦合动力分析：

式中：M₁、M₂、M₃分别为外管、内管与油管的质量矩阵；C₁、C₂、C₃分别为外管、内管与油管的阻尼矩阵；K₁、K₂、K₃分别为外管、内管与油管的刚度矩阵；{F₁}、{F₂}、{F₃}分别为外管、内管与油管的受力；

分别为外管、内管、油管的横向响应加速度；

为横向响应速度；x₁、x₂、x₃为横向响应位移；

(13)通过步骤(12)可以求得外管的节点位移，用插值法求解单元在外分布载荷作用下的最大弯曲变形，以此来判断外层立管与内层立管是否会发生碰撞；

(14)通过步骤(12)可以求得三根立管随时间的动力响应，由此来确定各层管各个单元节点的应力与弯矩，进行疲劳分析。

2.根据权利要求1所述的三层顶张式立管动力响应的有限元分析方法，其特征在于：

(1)公式(3)中，在计算外管的单元质量矩阵时，应考虑附加质量的影响：

式中：C_a为附加惯性力系数；ρ_w为海水密度；

(2)公式(4)中，在计算外管的单元刚度矩阵时，应考虑顶张紧力的影响：

式中：[K_G]^e为单元几何刚度矩阵；T为单元所受顶张紧力的大小。

Images