CN110046451B - 变张力细长柔性圆柱体涡激振动响应预测方法 - Google Patents

Abstract

本发明公开了一种变张力细长柔性圆柱体涡激振动响应预测方法，主要包括以下步骤：1)建立结构和流场互为耦合的振动控制方程；2)基于有限差分法对耦合的振动控制方程进行求解；3)基于分析数据，对实例进行计算分析。本发明针对变张力细长柔性圆柱体涡激振动响应预测进行了全面研究。建立一个完整的流体‑结构耦合模型，用于分析顶张式立管、钢悬链式立管等变张力柔性圆柱体涡激振动响应的预测问题。该模型实时考虑了结构与流场之间互为耦合的影响，可很好地预测变张力柔性圆柱体的振动响应特性。

Description

变张力细长柔性圆柱体涡激振动响应预测方法

技术领域

本发明涉及细长柔性圆柱体涡激振动响应领域，具体地说，特别涉及一种变张力细长柔性圆柱体涡激振动响应预测方法。

背景技术

近些年，伴随着海洋油气工程逐步向深海发展，出现了大量的超细长柔性结构，如：用于油气开发的各种形式的海洋立管(顶张式立管、钢悬链式立管等)以及用于浮式海洋平台(半潜式平台、立柱式平台等)定位的细长锚链。这种细长柔性结构在一定的洋流作用下，会在结构尾部形成交替脱落的漩涡，周期性的漩涡脱落会在结构上产生周期性的水动力载荷，从而诱发结构发生振动，称为“涡激振动”(Gao,Y,Fu,SX,Xiong,YM,Zhao,Y,andLiu,LM,Experimental study on response performance of vortex-induced vibrationof a flexible cylinder,Ships and Offshore Structures,2017,12(1):116-134)。涡激振动响应是导致结构产生疲劳的主要因素之一，因此有必要对其加以研究。

细长柔性圆柱体涡激振动响应的研究方法主要可分为实验方法以及数值方法。与数值方法相比，实验方法获取的数据更为可靠、得到的现象更为直观，但实验方法同样存在一些缺陷，比如：实验研究成本昂贵、实验结构尺寸受实验场地空间大小限制以及复杂流剖面难以模拟等。数值方法通常包括计算流体动力学(Computational fluid dynamics,CFD)方法以及经验模型方法(高云，邹丽，宗智，基于有限差分法刚性圆柱体涡激振动响应特性数值研究，计算物理，2019，36(1)：53-59)。近些年，随着计算机计算速度的快速发展，CFD方法被广泛应用于预测圆柱体涡激振动响应中。由于对整个轴线方向的全尺度细长柔性圆柱体涡激振动响应进行数值模拟非常困难，目前绝大多数CFD研究还仅局限于对刚性圆柱体的涡激振动响应展开数值模拟。在解决细长柔性圆柱体涡激振动响应问题时，由于实验方法的高实验研究成本特性以及CFD方法的高数值研究成本特性，因此有必要建立一种快速预报细长柔性圆柱体涡激振动响应的方法，而尾流振子模型方法则正是近些年广泛应用于细长柔性圆柱体涡激振动响应预报的一种半经验模型方法。

目前针对细长柔性圆柱体涡激振动响应的研究主要集中在恒张力细长柔性圆柱体，即：假定圆柱体的轴向张力沿轴线方向恒定不变。而实际海洋工程中的立管，如张紧式立管以及钢悬链式立管，其张力分布沿结构轴线方向呈变化趋势：从立管上端到立管下端其张力逐渐减小。对于这种变张力分布的细长柔性圆柱体，其涡激振动响应还有待进一步展开研究。

发明内容

为了解决现有技术的问题，本发明实施例提供了一种变张力细长柔性圆柱体涡激振动响应预测方法。所述技术方案如下：

本发明提供了一种变张力细长柔性圆柱体涡激振动响应预测方法，主要包括以下步骤：

1)建立结构和流场互为耦合的振动控制方程；

2)基于有限差分法对耦合的振动控制方程进行求解；

3)基于分析数据，对实例进行计算分析。

进一步地，所述步骤1)具体包括：

a、建立均匀流中变张力细长柔性圆柱体模型

b、在所述均匀流中变张力细长柔性圆柱体模型中，考虑一变张力细长柔性圆柱体在横流方向的涡激振动响应问题；假设结构长度为L，直径为D，结构两端采用铰接边界条件，流场为速度大小为U的均匀来流；坐标系X方向取顺流方向，Z方向为铅直方向，Y方向为横流向振动方向；柔性圆柱体上受到的张力沿着轴线Z方向发生变化，记为Θ(Z)，最大张力以及最小张力分别位于圆柱体上端以及下端处，分别记为：Θ_max以及Θ_min；柔性圆柱体的弯曲刚度为EI。

进一步地，将所述的细长柔性圆柱体看作变张力梁模型，建立如下振动方程：

式(1)中T为时间，p(Z,T)为升力，表示为：p(Z,T)＝ρU²DC_L0q(Z,T)/4，其中q(Z,T)表示尾部流场的运动，C_L0为静止圆柱体横向升力系数；式(1)中系统单位长度质量m包括结构质量m_s以及附加流体质量m_f，表示为：m＝m_s+m_f,m_f＝C_MρD²π/4，其中ρ表示流体密度，C_M表示流体附加质量系数；式(1)中阻尼R包括结构阻尼R_s和流体阻尼R_f，表示为：R＝R_s+R_f,R_f＝γΩ_fρD²，其中γ为黏滞力系数，与流体拖曳力系数C_D的关系为：γ＝C_D/(4πSt)，St数为斯脱哈尔数；对于海洋工程中的低质量阻尼比介质(如水)，通常结构阻尼与流体阻尼比起是小量，因此本发明忽略了结构阻尼的影响，只考虑了流体阻尼；Ω_f为依据斯脱哈尔关系式得到的斯脱哈尔漩涡泄放频率，写作：Ω_f＝2πStU/D；

采用改进的Van der pol方程来满足尾流振子的非线性特性，表达式如下：

式(2)中ε和A为实验测试得到的经验小参数；将方程(1)和(2)转换成无量纲形式，令:y＝Y/D,z＝Z/D,t＝T·Ω_f，并将其代入到方程(1)和(2)中，整理得到结构与尾部流场耦合振动方程的无量纲形式如下：

式(3)中，μ为质量比，表示为μ＝m/(ρD²)，无量纲参数M表示为：M＝C_L0/(16π²St²μ)；式(3)中b为无量纲弯曲刚度，表示为：

c(z)为无量纲张力，表示为

本发明假设无量纲张力c(z)沿着轴线方向呈线性变化，且假定最大无量纲张力c_max和最小无量纲张力c_min之间的关系为c_max＝2c_min，这样仅需要通过平均张力c_mean的大小便可确定整个轴线上的张力分布情况；由于平均张力c_mean＝(c_max+c_min)/2，因此c_min＝2c_mean/3，c_max＝4c_mean/3；那么任意一点z处对应的无量纲张力表示为：c(z)＝c_min+(c_max-c_min)×z/(L/D)。

进一步地，所述步骤2)具体如下：

采用具备二阶精度的有限差分法对方程(3)进行求解，将结构无量纲总长度L/D划分为M段；将无量纲总时间t_total划分为N段；因此，数值计算时空间步长Δz＝L/(D×M)；时间步长Δt＝t_total/N；被划分后的M+1个空间点记为：z＝z_i(i＝0,1,2,…,M)；被划分后的N+1时间点记为：t＝t_j(j＝0,1,2,..,N)；假设t_n时刻z_m位置处参数y和q表示为

以及

因此式(3)中各偏导数项的二阶精度差分格式表示为：

将式(4)代入式(3)，并整理得到：

y的初始条件设为在整个轴线上圆柱体振动位移以及速度均为0，即：y＝0以及

q的初始条件设为：q为一小幅扰动以及

圆柱体两端采用铰接边界条件，即结构在两个端点处位移为0且弯矩为0，表示如下：

由式(5)可以看出：当2≤m≤M-2时y直接求出；当m＝0以及m＝M时，则需结合位移边界条件对其加以求解，由两端处位移为0得到：

当m＝1以及m＝M-1时，则需结合弯矩为0边界条件，由两端弯矩为0得到：

联合式(8)以及式(5)得到m＝1以及m＝M-1时y的表达式如下：

当t_n+1时刻z_m处的y已知，便依据式(9)的第二个表达式求得t_n+1时刻z_m处的q，依此类推对式(5)进行迭代求解。

本发明实施例提供的技术方案带来的有益效果是：

本发明针对变张力细长柔性圆柱体涡激振动响应预测进行了全面研究。建立一个完整的流体-结构耦合模型，用于分析顶张式立管、钢悬链式立管等变张力柔性圆柱体涡激振动响应的预测问题。该模型实时考虑了结构与流场之间互为耦合的影响，可很好地预测变张力柔性圆柱体的振动响应特性。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明实施例的均匀流中变张力细长柔性圆柱体模型的示意图；

图2是本发明实施例的不同平均张力下圆柱体的位移响应变化云图,瞬时变化特性以及无量纲位移均方根值的示意图；

图3是本发明实施例的不同平均张力下圆柱体振动位移谱分析的示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

本发明提供了变张力细长柔性圆柱体涡激振动响应预测方法，主要包括以下步骤：

1)建立结构和流场互为耦合的振动控制方程；

2)基于有限差分法对耦合的振动控制方程进行求解；

3)基于分析数据，对实例进行计算分析。

步骤1：建立结构和流场互为耦合的振动控制方程，具体如下：

如图1所示，考虑一变张力细长柔性圆柱体在横流方向的涡激振动响应问题。假设结构长度为L，直径为D，结构两端采用铰接边界条件，流场为速度大小为U的均匀来流。坐标系X方向取顺流方向，Z方向为铅直方向，Y方向为横流向振动方向。柔性圆柱体上受到的张力沿着轴线Z方向发生变化，记为Θ(Z)，最大张力以及最小张力分别位于圆柱体上端以及下端处，分别记为：Θ_max以及Θ_min；柔性圆柱体的弯曲刚度为EI。

将图1中的细长柔性圆柱体看作变张力梁模型，建立如下振动方程：

式(1)中T为时间，p(Z,T)为升力，可表示为：p(Z,T)＝ρU²DC_L0q(Z,T)/4(Facchinetti,ML,de Langre,E,Biolley,F,Coupling of structure and wakeoscillators in vortex-induced vibrations,Journal of Fluids and Structures,2004,19:123-140)，其中q(Z,T)表示尾部流场的运动，C_L0为静止圆柱体横向升力系数；式(1)中系统单位长度质量m包括结构质量m_s以及附加流体质量m_f，可表示为：m＝m_s+m_f,m_f＝C_MρD²π/4，其中ρ表示流体密度，C_M表示流体附加质量系数；式(1)中阻尼R包括结构阻尼R_s和流体阻尼R_f，可表示为：R＝R_s+R_f,R_f＝γΩ_fρD²，其中γ为黏滞力系数，与流体拖曳力系数C_D的关系为：γ＝C_D/(4πSt)，St数为斯脱哈尔数。对于海洋工程中的低质量阻尼比介质(如水)，通常结构阻尼与流体阻尼比起是小量，因此本发明忽略了结构阻尼的影响，只考虑了流体阻尼。Ω_f为依据斯脱哈尔关系式得到的斯脱哈尔漩涡泄放频率，可写作：Ω_f＝2πStU/D。

采用改进的Van der pol方程来满足尾流振子的非线性特性，表达式如下(Nayfeh,AH,Introduction on Perturbation Techniques,New York:Wiley,1993)：

式(2)中ε和A为实验测试得到的经验小参数。将方程(1)和(2)转换成无量纲形式，令:y＝Y/D,z＝Z/D,t＝T·Ω_f，并将其代入到方程(1)和(2)中，整理得到结构与尾部流场耦合振动方程的无量纲形式如下：

式(3)中，μ为质量比，可表示为μ＝m/(ρD²)，无量纲参数M表示为：M＝C_L0/(16π²St²μ)；式(3)中b为无量纲弯曲刚度，表示为：

c(z)为无量纲张力，表示为

本发明假设无量纲张力c(z)沿着轴线方向呈线性变化，且假定最大无量纲张力c_max和最小无量纲张力c_min之间的关系为c_max＝2c_min，这样仅需要通过平均张力c_mean的大小便可确定整个轴线上的张力分布情况。由于平均张力c_mean＝(c_max+c_min)/2，因此c_min＝2c_mean/3，c_max＝4c_mean/3。那么任意一点z处对应的无量纲张力可表示为：c(z)＝c_min+(c_max-c_min)×z/(L/D)。

步骤2：基于有限差分法对耦合的振动控制方程进行求解，具体如下：

采用具备二阶精度的有限差分法对方程(3)进行求解，将结构无量纲总长度L/D划分为M段；将无量纲总时间t_total划分为N段。因此，数值计算时空间步长Δz＝L/(D×M)；时间步长Δt＝t_total/N。被划分后的M+1个空间点记为：z＝z_i(i＝0,1,2,…,M)；被划分后的N+1时间点记为：t＝t_j(j＝0,1,2,..,N)。假设t_n时刻z_m位置处参数y和q表示为

以及

因此式(3)中各偏导数项的二阶精度差分格式可表示为：

将式(4)代入式(3)，并整理得到：

y的初始条件设为在整个轴线上圆柱体振动位移以及速度均为0，即：y＝0以及

q的初始条件设为：q为一小幅扰动以及

圆柱体两端采用铰接边界条件，即结构在两个端点处位移为0且弯矩为0，可表示如下：

由式(5)可以看出：当2≤m≤M-2时y可直接求出；当m＝0以及m＝M时，则需结合位移边界条件对其加以求解，由两端处位移为0得到：

当m＝1以及m＝M-1时，则需结合弯矩为0边界条件，由两端弯矩为0得到：

联合式(8)以及式(5)得到m＝1以及m＝M-1时y的表达式如下：

当t_n+1时刻z_m处的y已知，便可依据式(9)的第二个表达式求得t_n+1时刻z_m处的q，依此类推对式(5)进行迭代求解。

步骤3：实施计算与分析，具体如下：

(1)无量纲基本参数，参见表1；

(2)计算结果

图2给出了不同平均张力c_mean下圆柱体振动无量纲位移随时间与空间的变化云图、某一时刻振动位移随空间变化的瞬时特性以及沿轴线方向的振动位移无量纲均方根值。由最右边的振动位移无量纲均方根值可看出：振动无量纲位移均方根值最大值随着c_mean的增大变化不明显；而对于某个特定的c_mean，随着张力的减小(即沿z轴负方向)，振动位移均方根值呈减小趋势。整个轴线方向的无量纲位移均方根值均大于0(两个端点除外)说明：圆柱体的涡激振动位移响应是由驻波和行波共同混合组成。

当c_mean较小(10和30)时，在圆柱体两端附近，位移均方根值最小值接近0说明：此时振动响应由驻波占主导；而在圆柱体中间段附近，位移均方根值最小值要明显大于0说明：此时振动响应由行波占主导。当c_mean较大(50)时，在整个圆柱体轴线长度上位移均方根值均接近0说明：在整个轴线上均由驻波占主导。

由图2中的位移随时间和空间的变化云图可以看出：对于均匀来流下的变张力细长柔性圆柱体，行波沿着从低张力区向高张力区进行传播(即沿z轴正方向)，且传播的速度随着平均张力c_mean的增加而迅速上升。为了确定圆柱体涡激振动位移响应的振动波长，图2也给出了不同平均张力c_mean下某一时刻(t＝2025)振动位移随空间变化的瞬时特性。当c_mean为10、30以及50时，圆柱体轴线长度(L/D＝2000)上分别存在20、10以及6个波长λ，说明圆柱体涡激振动位移响应的振动波长随着平均张力c_mean的增大而呈增加趋势。对于某个特定的c_mean，振动波长λ沿着轴线方向是发生变化的：振动波长λ沿着圆柱体轴线方向从低张力区向高张力区(即沿z轴正方向)呈增加趋势。

取一稳定时间段(t＝2700-3000)的振动位移响应对其进行快速傅立叶(FastFourier Transform,FFT)变换，图3给出了不同平均张力c_mean下圆柱体沿z方向的振动位移频率响应特性。由图3可以看出三种不同平均张力c_mean下圆柱体的涡激振动位移频率响应均呈现单频分布特性：即位移响应仅存在一个峰值频率，这是由发明所研究的流场为均匀流场的特性所决定的。当平均张力c_mean为10、30以及50时，圆柱体振动位移响应主导频率(图3中白色虚线)分别为0.984、0.942以及0.922，即：振动位移响应主导频率随着平均张力的增加呈缓慢下降趋势。且从低张力区向高张力区(即沿z轴正方向)，位移响应的频率带宽逐渐变大。

本发明实施例提供的技术方案带来的有益效果是：

本发明针对变张力细长柔性圆柱体涡激振动响应预测进行了全面研究。建立一个完整的流体-结构耦合模型，用于分析顶张式立管、钢悬链式立管等变张力柔性圆柱体涡激振动响应的预测问题。该模型实时考虑了结构与流场之间互为耦合的影响，可很好地预测变张力柔性圆柱体的振动响应特性。

以上仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种变张力细长柔性圆柱体涡激振动响应预测方法，其特征在于，包括以下步骤：

1)建立结构和流场互为耦合的振动控制方程；

2)基于有限差分法对耦合的振动控制方程进行求解；

3)基于分析数据，对实例进行计算分析；

所述步骤1)具体包括：

a、建立均匀流中变张力细长柔性圆柱体模型

b、在所述均匀流中变张力细长柔性圆柱体模型中，考虑一变张力细长柔性圆柱体在横流方向的涡激振动响应问题；假设结构长度为L，直径为D，结构两端采用铰接边界条件，流场为速度大小为U的均匀来流；坐标系X方向取顺流方向，Z方向为铅直方向，Y方向为横流向振动方向；柔性圆柱体上受到的张力沿着轴线Z方向发生变化，记为Θ(Z)，最大张力以及最小张力分别位于圆柱体上端以及下端处，分别记为：Θ_max以及Θ_min；柔性圆柱体的弯曲刚度为EI；

将所述的细长柔性圆柱体看作变张力梁模型，建立如下振动方程：

式(1)中T为时间，p(Z,T)为升力，表示为：p(Z,T)＝ρU²DC_L0q(Z,T)/4，其中q(Z,T)表示尾部流场的运动，C_L0为静止圆柱体横向升力系数；式(1)中系统单位长度质量m包括结构质量m_s以及附加流体质量m_f，表示为：m＝m_s+m_f,m_f＝C_MρD²π/4，其中ρ表示流体密度，C_M表示流体附加质量系数；式(1)中阻尼R包括结构阻尼R_s和流体阻尼R_f，表示为：R＝R_s+R_f,R_f＝γΩ_fρD²，其中γ为黏滞力系数，与流体拖曳力系数C_D的关系为：γ＝C_D/(4πSt)，St数为斯脱哈尔数；对于海洋工程中的低质量阻尼比介质，结构阻尼与流体阻尼比起是小量，因此忽略了结构阻尼的影响，只考虑了流体阻尼；Ω_f为依据斯脱哈尔关系式得到的斯脱哈尔漩涡泄放频率，写作：Ω_f＝2πStU/D；

采用改进的Van der pol方程来满足尾流振子的非线性特性，表达式如下：

式(2)中ε和A为实验测试得到的经验小参数；将方程(1)和(2)转换成无量纲形式，令:y＝Y/D,z＝Z/D,t＝T·Ω_f，并将其代入到方程(1)和(2)中，整理得到结构与尾部流场耦合振动方程的无量纲形式如下：

式(3)中，μ为质量比，表示为μ＝m/(ρD²)，无量纲参数M表示为：M＝C_L0/(16π²St²μ)；式(3)中b为无量纲弯曲刚度，表示为：

c(z)为无量纲张力，表示为

假设无量纲张力c(z)沿着轴线方向呈线性变化，且假定最大无量纲张力c_max和最小无量纲张力c_min之间的关系为c_max＝2c_min，这样仅需要通过平均张力c_mean的大小便可确定整个轴线上的张力分布情况；由于平均张力c_mean＝(c_max+c_min)/2，因此c_min＝2c_mean/3，c_max＝4c_mean/3；那么任意一点z处对应的无量纲张力表示为：c(z)＝c_min+(c_max-c_min)×z/(L/D)。

2.如权利要求1所述的方法，其特征在于，所述步骤2)具体如下：

采用具备二阶精度的有限差分法对方程(3)进行求解，将结构无量纲总长度L/D划分为M段；将无量纲总时间t_total划分为N段；因此，数值计算时空间步长Δz＝L/(D×M)；时间步长Δt＝t_total/N；被划分后的M+1个空间点记为：z＝z_i(i＝0,1,2,…,M)；被划分后的N+1时间点记为：t＝t_j(j＝0,1,2,..,N)；假设t_n时刻z_m位置处参数y和q表示为

以及

因此式(3)中各偏导数项的二阶精度差分格式表示为：

将式(4)代入式(3)，并整理得到：

y的初始条件设为在整个轴线上圆柱体振动位移以及速度均为0，即：y＝0以及

q的初始条件设为：q为一小幅扰动以及

圆柱体两端采用铰接边界条件，即结构在两个端点处位移为0且弯矩为0，表示如下：

由式(5)可以看出：当2≤m≤M-2时y直接求出；当m＝0以及m＝M时，则需结合位移边界条件对其加以求解，由两端处位移为0得到：

当m＝1以及m＝M-1时，则需结合弯矩为0边界条件，由两端弯矩为0得到：

联合式(8)以及式(5)得到m＝1以及m＝M-1时y的表达式如下：

当t_n+1时刻z_m处的y已知，便依据式(9)的第二个表达式求得t_n+1时刻z_m处的q，依此类推对式(5)进行迭代求解。

Images