0 引 言

随着世界范围内石油需求的不断增大，海洋油气生产从近海逐渐向深水、超深水发展，这对海洋石油装备的设计与建造带来了新的挑战。海洋立管是连接海洋平台与水下生产系统的关键环节，其内部输送高压流体等，外部受浪、潮流和冰等载荷的作用，因此，立管也是易损环节。当海流流经立管时，在一定流速范围内会产生漩涡脱落，诱使立管在与流速垂直的升力方向发生横向振动，易引起立管的疲劳破坏；立管一旦发生破坏，不仅经济损失巨大，还将造成严重的环境污染和次生灾害。目前，国内外许多学者通过试验和数值模拟的方法对立管涡激振动及其疲劳寿命进行了研究，得到了较丰富的研究成果。基于试验的涡激振动测试结果虽真实可信^{[1, 2]}，但无法对物理现象的形成过程及机制进行细致描述，由于试验尺寸与工程尺寸的巨大差异，其结论推广到工程应用中存在困难^{[3, 4, 5]}；基于数值模拟的涡激振动及其疲劳寿命的研究方法，既可以对涡激振动的形成与发展过程进行细致的刻画，也可以对工程尺寸的立管进行研究，但是，如何提高大规模数值模拟的计算效率^{[6, 7]}及计算结果的可靠性仍是该课题的难点。

笔者基于弹簧-质量模型，建立等效二维模型，通过二维数值模拟对三维问题中的涡激振动参数进行求解，对计算结果的可靠性进行验证，进而计算不同条件下立管中点的振幅和频率，以及涡激振动对立管的疲劳损伤率，对立管疲劳寿命的影响因素进行分析，从而为立管的疲劳寿命预测与抑振设计提供指导。

笔者研究的立管横截面为均匀圆环，下端铰接，上端与平台相连，可视为简支梁，竖直地处于恒定海流中。立管三维模型如图 1所示。

图 1 立管三维模型 Fig. 1 3D model of riser

图选项

静力平衡时，立管中点挠度可采用均布载荷下简支梁的挠度方程求解^[8]：

式中：q为流体作用于立管的均布载荷，N/m；L为管长，m；E为材料的弹性模量，Pa；I为管道中点处横截面对其形心的惯性矩，m⁴；Δy为立管中点挠度，m。

此时立管所受总载荷为：

联立式(1)和式(2)，可得：

由式(3)可见，对于立管中点的横截面单元来说，所受内力大小与位移(立管中点挠度)成正比，且总是指向平衡位置，满足胡克定律：

式中：负号表示力的方向与位移方向相反，k为弹性系数，表达式为：

因此，可将立管中点的振动模型简化为2自由度(或单自由度)质量-弹簧系统。立管中点涡激振动的流固耦合模型可抽象为：将刚性圆柱单元与弹簧所组成的系统置于流速为u的恒定水流中，柱体在流体力的作用下受迫振动。图 2为立管中点横截面涡激振动等效二维模型示意图。

图 2 立管横截面涡激振动等效模型示意图 Fig. 2 Equivalent model of riser cross-section VIV

图选项

作用在立管上的流体力可通过求解流体力学控制方程得到。在笛卡尔坐标下，不可压缩流体连续性方程可表示为：

对于不可压缩牛顿流体的湍流流动，采用Navier-Stokes方程进行描述。由于在立管涡激振动的二维数值模拟中，XOY面内流体受到的重力分量为0，所以，二维条件下的运动方程可表示为：

式中:p为压力，Pa；ρ为流体密度，kg/m³；ν为水的运动黏度系数，m²/s。

对于湍流模拟，采用雷诺平均(RANS)数值模型，引入重整化群(RNG)κ-ε湍流模型可较好地处理流线弯曲程度较大的流动。动量方程中的对流项采用高级精度格式进行离散，瞬态项采用二阶迎风后插值方法进行离散，可保证计算结果有较高的精度。

柱体的振动响应通过求解结构动力学控制方程得到。根据牛顿第二定律，刚性柱体X、Y方向的运动方程为：

采用时域的Newmark-β法对刚体模型求解结构动力学方程。该方法是一种加速度法，根据时间增量内假定的加速度变化规律计算结构的动力响应，通过设置γ和β2个参数来控制数值精度和稳定性^[9]。

处于涡激振动状态的海洋立管会产生交变应力。管道若长期受到交变应力的作用，将逐渐形成疲劳裂纹并扩展，最终导致管道的疲劳断裂失效。对于特定的管道，交变应力的频率和大小决定了该管道的疲劳寿命，而交变应力的大小与管道涡激响应的振幅有关。

根据目前应用较为广泛的Miner疲劳损伤累积理论^[10]，当Miner系数D_a累积达到1时，管道将产生完全损伤即失效。Miner系数计算式为：

式中：n_i是第i种应变水平下的实际循环次数，N_i为第i种应变水平下的疲劳极限，可通过相应的S-N曲线求得，可近似表示成如下关系^[11]：

根据美国焊接协会提出的S-N曲线，通常c取6.4×10^-18，b取4。Δε为1个周期内出现的最大应变差。对于两端支承的立管，最大应变差出现在管道的中点，Δε可表示成：

式中：A为振幅，D为管径，L为管长。

响应的循环周期n_i可表示为：

式中：f_i为第i种振幅的频率，1/s；t_i是相应的时间，s。

联立式(11)~式(14)可得：

式中：A_i为管道中点的第i种振幅，m；t_i为发生振动的时间，s。

式(17)即为管道疲劳破坏的判断准则。当管道振动为单一振动形式或某一振动形式占主导时，立管疲劳寿命公式可简化为：

式中：t为立管的疲劳寿命，s；f为主导振动形式的频率。

而特定使用期限内，某种振动形式对疲劳寿命的影响可用疲劳损伤率R_i来表示：

首先对静止圆柱的二维绕流进行了数值模拟，验证流体力计算的准确性。T.K.Prasanth^[12]研究发现，当尾流区域长度L_d≥25.5D，整体区域高度H≥20D时，柱体的数值模拟将不受流体区域边界的影响。因此，选取20D×40D大小的流体计算区域，其中柱体中心距入口10D，距出口30D，距上下边界均为10D。圆柱直径取实际立管直径0.25 m。网格划分采用分块划分方法，其中圆柱体近场及尾流区域进行网格加密，其余网格比较稀疏，见图 3。

图 3 流体场的网格划分 Fig. 3 Fluid field mesh generation

图选项

边界条件：入流面采用速度入口边界；出口面采用自由出流边界；上下面采用自由滑移边界，前后面定义为对称边界；圆柱表面定义为无滑移壁面边界。

以流速u=0.4 m/s ，雷诺数Re=1×10⁵条件下的数值模拟结果为例，验证数值模拟方法的准确性。图 4和图 5分别为漩涡脱落时的速度云图和压力云图。由图可见，u=0.4 m/s时，圆柱绕流的尾流区域可以形成以一定频率交替脱落的漩涡(涡街)。

图 4 漩涡脱落时的速度云图 Fig. 4 Velocity distribution when vortex shed

图选项

图 5 漩涡脱落时的压力云图 Fig. 5 Pressure distribution when vortex shed

图选项

图 6和图 7分别为立管处于亚临界区域u=0.4 m/s、Re=10⁵时的受力特性。图 6为圆柱绕流的曳力系数C_d 和升力系数C_l 的时程曲线。从图 6可以看出，当漩涡脱落充分发展后，圆柱的升力与曳力系数均呈现周期性的震荡，且幅值保持恒定，这说明流动稳定后的漩涡脱落是周期性的；曳力系数的振荡频率约为升力系数的2倍，而曳力系数的幅值远远小于升力系数。由于曳力幅值远远小于升力幅值，其引起的振动幅值也将远远小于升力。通过式(19)可以看出，从曳力、升力引起的涡激响应对管道疲劳寿命的影响来看，可忽略曳力方向的振动对疲劳寿命的影响，而将升力方向的涡激响应作为影响立管疲劳寿命的主要振动。

图 6 Re=105时圆柱绕流升力系数和曳力系数时程曲线 Fig. 6 Time histories of drag coefficient and lift coefficient of flow around cylinder when Re=105

图选项

图 7 Re=105时圆柱升力系数频谱图 Fig. 7 The spectrum of lift coefficient of flow around cylinder when Re=105

图选项

对升力系数利用快速傅里叶变换(FFT)得到其频谱图，见图 7。图中最大幅值对应的频率即为漩涡脱落频率f_s=0.32 Hz。折算斯特劳哈尔数Sr=f_sD/U=0.2，与文献[13]中由经典试验绘制的Sr-Re包络线相吻合。通过表 1中计算结果与文献数据的对比可知，圆柱绕流的二维数值模拟结果较为可靠。

计算的来流流速为0.1~1.0 m/s，对应的雷诺数为2.5×10⁴~2.5×10⁵，雷诺数较高，由文献[13]查得此范围的漩涡脱落属于亚临界区域。

利用ANSYS软件的多场耦合模块对海洋立管涡激振动进行二维数值模拟。结构动力学模块中设置流固耦合面(FSI)和弹性支撑面均为管道外表面，设定二维模型的弹性系数k以及结构的材料属性等；流体场设置同第2节；最后对结构与流场进行耦合求解。最终求得升力方向的振幅和频率，并利用式(17)计算不同海流流速、管长和壁厚条件下，20 a使用期内立管的疲劳损伤率，进而分析各因素对疲劳寿命的影响。

对长30 m，直径0.25 m，壁厚3 mm的立管，在0.1~1.0 m/s的海洋流速范围内对管道中点的涡激响应进行求解，并求出对应流速下的疲劳损伤率。

图 8为流速对立管中点振幅和频率的影响曲线。图 9为不同海流流速下的立管疲劳损伤率曲线。由图可见，随着海流流速的增加，立管振动的幅值和频率逐渐增大；疲劳损伤率增加，相应疲劳寿命缩短。这与亚临界区域内随着雷诺数的增加漩涡脱落强度增加的理论结果^[14]相一致。

图 8 海洋流速对立管中点振幅和频率的影响 Fig. 8 Effect of current velocity on vibration amplitude and frequency of riser midpoint

图选项

图 9 不同海流流速下的立管疲劳损伤率 Fig. 9 Fatigue damage rate of riser at different current velocities

图选项

对长度为15、30、45和60 m的立管的涡激响应进行计算，立管直径0.25 m，壁厚3 mm，海流流速0.5 m/s。图 10和图 11分别给出了该条件下不同长度的立管对应的振幅、频率和疲劳损伤率。由图可知，在雷诺数Re=1.25×10⁵下，随着长度的增加，立管中点的频率基本不变，而振幅上升较多；疲劳损伤率增大，相应的疲劳寿命缩短。这主要是因为单位长度的流体力学特性(升力大小、频率)基本不变，而管道长度增加引起管道中点挠度增大。

图 10 立管长度对立管中点振幅和频率的影响 Fig. 10 Effect of length of riser on vibration amplitude and frequency of riser midpoint

图选项

图 11 不同长度立管的疲劳损伤率 Fig. 11 Fatigue damage rate of risers with different lengths

图选项

对壁厚为3、4、5、6和7 mm的立管的涡激响应进行计算，立管管长30 m，直径0.25 m，海流流速0.5 m/s。图 12和图 13分别给出了该条件下不同壁厚对应的振幅、频率和疲劳损伤率。

图 12 立管壁厚对立管中点振幅和频率的影响 Fig. 12 Effect of wall thickness on vibration amplitude and frequency of riser midpoint

图选项

图 13 不同壁厚立管的疲劳损伤率 Fig. 13 Fatigue damage rate of riser with different wall thickness

图选项

由图 12和图 13可知，随着管道壁厚的增加，管道振动幅度逐渐下降，而频率变化不大；疲劳损伤率下降，相应的疲劳寿命延长。这与管道壁厚增大使管道的抗弯刚度上升，管道抵抗变形的能力增大，而流体对管道的力学特性基本保持不变的理论分析结果一致，也说明了所用方法可靠。

(1)海洋立管涡激振动的二维数值模拟结果与理论分析结果相一致，验证了等效二维模型的合理性及流体场计算的准确性。

(2)采用简化模型的二维涡激振动数值模拟可以很好地计算管道中点的振动频率和幅值，并进行疲劳寿命因素分析，较之长管道的三维数值模拟具有硬件开销小、步骤较少和求解耗时短的优点。

(3)在亚临界区域，随着海洋流速的增加，立管振幅与频率均增大，立管疲劳寿命缩短；随着立管长度增加，振幅增大，疲劳寿命延长；随着管道壁厚增加，立管振幅减小，疲劳寿命延长。因此，在海洋立管涡激振动的抑振与抗疲劳设计中，应着眼于改善近立管流场和管道结构特性。