基于ABAQUS的深海采矿扬矿管纵向振动性能

宋庆辉 ，肖林京 ，姜海燕 ，刘秀杰 ，烟芳萍

（1.山东科技大学机械电子工程学院，山东 青岛 266590；2.山东科技大学智能装备学院，山东 泰安 271009）

扬矿管是连接采矿船和集矿机的重要设备，在深海采矿作业时，扬矿管处于悬挂状态，由于受采矿船升沉运动的影响，扬矿管会发生纵向振动，当采矿船升沉运动频率与管道的固有频率相近时，容易引发共振，产生较大的振动位移和动载荷，从而影响集矿机的准确定位，引发扬矿管的疲劳积累甚至引起管道的突然断裂[1].因此，求解扬矿管纵向振动的固有频率，分析扬矿管的动力学响应，对深水扬矿管的安全设计和强度校核意义重大.

当前，国内外学者对深水管道已经进行了一定的研究：Wu等[2]对采矿船升沉运动时扬矿管的动力学响应进行分析，运用非线性有限元法求解了时域中含时变系数的扬矿管控制方程；Wang等[3]分析了3 000 m长钻井隔水管在安装和下井作业过程中的纵向振动，确定了隔水管的固有频率、振动位移和振动力；Wang等[4]采用蚁群优化算法对自由悬挂深水立管因上端位移而导致的底部振荡运动进行了优化，同时开展水箱试验和数值模拟，验证了优化算法的有效性；Liu等[5]建立了一种深水钻井悬挂立管非线性耦合动力学模型，运用Newmark法分析了立管的横向和纵向耦合运动，通过ANSYS验证了它的有效性；Thorsen等[6]采用一种新的扬矿管与内流耦合分析框架，对扬矿管的涡激振动进行了研究，重点分析了随时间变化的内部泥浆流和相关密度波对涡激振动的影响；邱显焱等[7]建立了附加单个动力吸振器的扬矿子系统纵向振动数学模型，运用Galerkin法求解，应用Matlab编制模拟分析程序，计算了不同工况下扬矿子系统的纵向振幅和轴向应力；肖林京等[8-9]针对波流联合作用下扬矿管的纵向振动进行研究，分析了四、六级海况时管道的动力学响应.对于上、下端铰接或固接的深海立管固有频率的分析，目前一般采用有限元法[10]、变量转换[11]、傅立叶分析[12]和动刚度法[13]等来计算，但对悬挂阶梯形管道固有频率的研究较少[14]，对扬矿管共振响应的分析也鲜有报道.

基于弹性杆振动理论，本文考虑扬矿管振动系统的刚度、质量和阻尼等因素，建立了具有复杂边界条件的阶梯扬矿管纵向振动数学模型和固有频率方程，通过ABAQUS软件对不同中间矿仓质量时管道的固有频率以及不同频率和升沉幅值时阶梯扬矿管纵向振动的动力学响应进行了仿真分析，为阶梯扬矿管纵向振动固有机制和管道稳定性的研究提供了理论基础.

1 扬矿管力学模型

深海采矿阶梯形扬矿管由4级管道组成，扬矿管顶端与采矿船铰接并与采矿船同步运动，末端通常视为自由端[15].

第i段管道的长度为Li（i=1，2，3，4），管上1 000 m处装有一个泵组，5 000 m处装有一个中间矿仓，泵组和中间矿仓被视为集中质量M1和M2.深海采矿扬矿管力学模型如图1所示.本文假设只考虑扬矿管的纵向振动，不考虑管道的偏移、弯曲与扭转.

图1 深海采矿扬矿管力学模型Fig.1 Mechanical model of lifting pipe in deep sea mining

受海浪的影响，假设扬矿管顶端随采矿船近似作简谐运动，其运动方程为

式中：u0(t) 为顶端的运动位移；r为采矿船升沉幅值；ω为波浪频率；t为时间.

采矿船的升沉幅值通过莫尔经验式[9]来估算.波浪频率由我国深海采矿矿区主要海况的P-M海浪谱确定，根据文献[16]确定P-M谱能量集中分布频带为 0.2～3.2 rad/s.

2 理论分析

2.1 扬矿管横向偏角

在深海采矿环境中，海流是一种流向稳定且同方向的平面流，可以用指数分布函数来表征海流的速度随垂直水深的变化分布情况，根据海平面以下的深度，得到海流速度（VC，m/s）随海水深度（x，海平面为0）变化的拟合式[17]为

拟合得到的海流速度如图2中的曲线A.由于非线性复杂模型一般难以求解，通过线性最小二乘近似法将当前的VC线性化为5段[17]，如图2中的曲线B，每段的海流速度可表示为

图2 海平面下的海流速度曲线Fig.2 Curve of ocean current velocity below sea level

在海流的作用下扬矿管发生横向偏移，扬矿管属于小尺度构件，单位长度的扬矿管发生偏角 θ 后所受的法相液动力可用莫里森方程[18]求得，如式（4）.

式中：CD为阻力系数； ρw为海水密度；Di为各管段的外径；VN为扬矿管与海水质点之间的相对法向速度；V为采矿船航行速度，取0.

根据式（3）和式（4），将倾斜扬矿管所受的法向液动力做线性化处理，可得法向液动力的线性分布函数qNk（x）,k=1,2,3,4,5.

当扬矿管系统平衡时，在原点O处力矩平衡其平衡方程为

利用MATLAB软件求解式（5），可以得到扬矿管在海流作用下的偏角 θ.

2.2 扬矿管纵向振动数学模型

把扬矿管看做阶梯连续弹性杆，根据达朗贝尔原理建立四级阶梯杆在海水中的纵向振动平衡方程为

式中：ui（x,t） 为第i管段任意截面的纵向振动位移；E为弹性模量；Si为横截面积；ci为阻尼系数； ρ 为扬矿管在水中的单位体积质量.

采用分离变量法求解振动微分方程，式（6）的解为

式中： ϕi（x） 为管道纵向振动的振型函数；qi（t） 为管道的振动方式qi（t） =eλt， λ =jω，为系统的一个特征值.

式（7）代入式（6），可以得到振型函数的表达式为

式（8）的通解可表示为

式中：Ai、Bi均为待定的复常数.

扬矿管顶部与采矿船升沉运动同步，顶部纵向运动位移为u1（0,t） =rcosθsinωt；扬矿管末端与中间矿仓相连，管道末端轴向力与中间矿仓惯性力之和为0；每两阶管道在连接处位移相等；在第1级、第2级管道的交界处，上一级的轴向力与惯性力之和等于下一级的轴向力；在第2级、第3级管道以及第3级、第4级管道的交界处，两级管道的轴向力相等.将振型函数的通解式（9）和时间函数代入上述阶梯管纵向振动的边界条件和连续性条件中，可得

式中：

由式（10）可求得待定复常数Ai和Bi，将其代入式（7）中可得阶梯扬矿管的纵向振动位移.

2.3 纵向振动固有频率

固有频率与扬矿管本身结构及其振动边界条件有关，与外界载荷无关[19].固有频率方程可通过分析扬矿管的无阻尼自由振动方程得到.

求解扬矿管固有频率时，系统阻尼为0，管道顶部位移边界条件变为u1（0,t） =0 ，类似于扬矿管强迫纵向振动方程的求解过程，将分离变量的通解式（9）代入边界条件、连续性条件中可得

式（12）是一个齐次线性方程组，由于其有非零解，因此系数矩阵行列式必须为0，进而可得阶梯扬矿管的固有频率方程为

2.4 主振型的归一化处理

设ωnk和ωnj为 扬 矿 管 不 相 等 的 两 个固有频率，ϕik（x）和 ϕij（x） 分别为 第i级扬 矿 管对应ωnk和ωnj的振型函数，根据阶梯管纵向无阻尼自由振动方程可得

对于任意一级扬矿管, 式（14）乘 ϕij（x） ，并从b（i-1）到bi积分，可得

对于4级阶梯形扬矿管整体，式（17）成立.

当固有频率变为ωnj时，上式同样成立.因此可得

将x=b1和x=b4处的边界条件和连续性条件代入式（18）中，整理可得阶梯扬矿管的质量归一化条件为

求解式（19）可得对应固有频率ωnk的系数A1k，将式（19）与式（12）联立，可得无阻尼自由振动的系数值Aik和Bik，进而求得扬矿管纵向振动的正则振型函数.

3 仿真分析及验证

在国家海洋技术中心海洋动力环境试验室开展了不同波高和不同波浪频率下的深海采矿系统试验，试验模型按照1∶80的缩尺比制作，如图3所示.试验中利用拉压力传感器、6自由测量系统、姿态传感器等分别获得了采矿船的运动位移和扬矿管的顶部轴向力.

图3 深海采矿船与扬矿管试验模型Fig.3 Test model of deep sea mining vessel and lifting pipe

基于ABAQUS对扬矿管的固有特性和动力学响应进行了仿真分析.ABAQUS/Aqua模块在海洋工程结构物的有限元仿真分析中应用广泛，可有效模拟海洋环境中的风、浪、流等环境载荷[20].仿真分析时，选择ABAQUS中的线性梁单元模拟较大细长比的扬矿管，扬矿管顶端设为铰接，末端定义为自由，洋流特性和Airy线性波通过编辑关键词来定义.仿真时，扬矿管作业水深为5 000 m，海水密度为1 028 kg/m3，泵组质量为8 t，中间矿仓质量为30 t.阶梯形扬矿管的具体参数见表1.

表1 阶梯扬矿管的参数Tab.1 Parameters of stepped lifting pipe

在进行仿真分析时，定义静力学分析时间步长为1 s，动力学分析时间步长为30 s.利用不同的r和ω定义采矿船的升沉运动来表征实际作业海况表面波对采矿船的影响，进而得到不同幅值和频率时扬矿管的动态特性.其中，r=0～6.0 m，ω=0.2～3.2 rad/s，所选择的振幅和频率可以覆盖深海采矿系统在运行过程中可能遇到的大多数海洋环境条件.

3.1 固有特性分析

扬矿管纵向振动固有频率的计算结果如表2所示.

表2 扬矿管前5阶固有频率Tab.2 The first five natural frequencies of lifting pipe

对比我国深海采矿区P-M海浪谱能量集中分布频带可以发现：扬矿管的1阶固有频率处于海浪经常出现的频率范围内.

扬矿管开始作业时，中间矿仓质量会随着矿仓内矿物质的增多而逐渐增加，受中间矿仓质量变化的影响，管道的固有频率发生改变，且变化程度不同，如图4所示.

图4 不同中间矿仓质量时的固有频率Fig.4 Natural frequency of different buffer masses

由图4可知：随着矿仓质量的增加，扬矿管固有频率逐渐降低，由于中间矿仓质量相对于整个扬矿管质量而言较小，中间矿仓质量增加时扬矿管低阶固有频率值变化并不明显，但高阶固有频率值明显减小.因此，研究扬矿管高阶纵向振动时必须考虑中间矿仓质量的影响.

3.2 纵向振动位移分析

纵向振动位移仿真计算时，调整主要参数值得到扬矿管不同位置的频率-振幅曲线和升沉幅值-振幅曲线分别如图5、6所示.

图5 不同频率时扬矿管的振幅（ r =1.0 m）Fig.5 Amplitude of lifting pipe at different frequencies（ r =1.0 m）

由图5可知：频率在0.20～2.40 rad/s内，随着采矿船升沉频率的增大，扬矿管纵向振幅先增大后减小，1 000、3 000、5 000 m 处振动位移同步变化； 当升沉频率为1阶固有频率时，扬矿管纵向振幅达到最大值，此时管道末端振幅为4.72 m； 泵组处纵向振幅在2.40 rad/s时最先达到最低值0.57 m，频率在2.40～3.20 rad/s内时，随着频率的增大，泵组处的纵向振幅开始逐渐增大，但3 000、5 000 m处的纵向振幅仍在逐渐减小，因此沿管长方向纵向振幅的变化趋势与波浪频率密切相关，通过式（10）可验证该结论.

由图6可知：随着采矿船升沉幅值的增加，扬矿管纵向振幅逐渐增大，当升沉幅值小于1.5 m时，扬矿管节点纵向振幅曲线的斜率沿管长从上到下逐渐增大，表明扬矿管节点纵向振幅的增长速度随着海水深度的增加而逐渐增大； 当升沉幅值大于1.5 m时，扬矿管节点纵向振幅的增长速度变缓，且低于采矿船升沉幅值的增长速度.

图6 不同升沉幅值时扬矿管的振幅（ ω =1.92 rad/s）Fig.6 Amplitude of the lifting pipe at different heave amplitudes （ ω =1.92 rad/s）

图7和图8为r=1.0 m，ω=1.92 rad/s时阶梯管纵向振动时历位移.从图7、8中可知：扬矿管在1阶固有频率下振动时，节点振幅较大，但由于阻尼的影响，不会无限增大, 而是经过一段时间后以增大的振幅作等幅稳态振荡，此时扬矿管的振动频率与采矿船的运动频率相同；由于重力的作用，扬矿管发生纵向拉伸变形，随海水深度的增加，管上各节点的振动平衡位置逐渐下降，在管末端中间矿仓处振动平衡位置下降6.7 m；由于扬矿管长达5 000 m，考虑外激励频率和系统阻尼，纵向振动响应时间从上到下会有延迟，产生相位差，延迟时间随着深度的增加而增大，管末端比顶端延迟0.7 s，产生的相位差约为80°；共振时沿管长方向扬矿管的纵向振幅呈现上小下大的趋势，在管末端中间矿仓处振幅达到最大值.

图7 纵向振动位移时历曲线Fig.7 Vibration displacement time-history curve

图8 沿管长方向的振动位移Fig.8 Vibration displacement along pipe length

3.3 纵向振动轴向力分析

由于扬矿管沿长度方向的纵向振幅和相位角存在差异，导致扬矿管产生动态轴向力，影响系统的安全.图9和图10给出了扬矿管的最大轴向力随频率和升沉幅值的变化规律.

图9 不同频率时扬矿管的最大轴向力（ r =1.0 m）Fig.9 Maximum axial force of lifting pipe at different frequencies （ r =1.0 m）

图10 不同升沉幅值时扬矿管的最大轴向力（ ω =1.92 rad/s）Fig.10 Maximum axial force of lifting pipe at different heave amplitudes （ ω =1.92 rad/s）

从图9可知：在0.20～3.40 rad/s范围内，随着采矿船运动频率的增大，扬矿管最大轴向力同步变化，最大轴向力值先增大后减小；1阶固有频率时，扬矿管顶端最大轴向力达到峰值11.80 MN.

从图10可知：随着采矿船升沉幅值的增加，扬矿管的最大轴向力逐渐增大；当升沉幅值大于1.5 m时，扬矿管节点最大轴向力的增长速度同样变缓，当升沉幅值由2.0 m增长到6.0 m时，扬矿管x=0处的最大轴向力由13.26 MN增长到14.35 MN，仅增长了8.2%.

在水动力实验室中开展A、B两组深海采矿系统在规则波中的试验，A组试验条件为：r= 0.1 m，ω=1.92rad/s （波向角为 90°，波高为 2.00 m）；B 组试验条件为：r= 1.0 m，ω=1.20 rad/s （波向角为 90°，波高为3.76 m）.在设定试验条件下，实时测量x= 0处扬矿管的轴向力，将试验测得的数据与ABAQUS仿真数据进行对比，结果如图11所示.由图可知：试验曲线与仿真曲线具有良好的一致性，验证了仿真模型和方法的有效性.

图11 轴向力的试验与仿真结果对比Fig.11 Comparison of test and simulation results of axial force

图12和图13为r=1.0 m，ω=1.92 rad/s时，时间和长度对管道轴向力的影响.

图12 轴向力时历曲线Fig.12 Axial force time-history curve

由图12和图13可知：扬矿管在1阶固有频率下振动时，管上节点轴向力幅值先增大然后作与采矿船同频的等幅稳态振荡；扬矿管顶端轴向力在起振时发生瞬态阶跃响应（最大值达39.30 MN），然后进入稳态过程，末端轴向力不发生振动，力值恒为0.29 MN（即中间矿仓的重量）；不同时刻沿管长方向扬矿管的轴向载荷随深度的增加逐渐降低，在管末端到达最小值0.29 MN，同时由于扬矿管横截面积从上到下阶梯式减小，随着水深的增加，各级扬矿管轴向力的减小速度逐渐变缓；在1 000 m处由于泵组质量的影响，轴向力会发生突变.

图13 沿管长方向的轴向力Fig.13 Axial force along pipe length

3.4 纵向振动轴向应力分析

轴向应力是影响管道安全的重要因素，最大应力点和应力值是管道设计和强度计算的重要依据.图14～16分别是频率、升沉幅值和管道长度对最大轴向应力的影响.

图14 不同频率时扬矿管的最大轴向应力（ r =1.0 m）Fig.14 Maximum axial stress of lifting pipe at different frequencies （ r =1.0 m）

分析图14和图15可知：在0.20～0.40 rad/s 范围内，随着采矿船升沉频率的增大，扬矿管最大轴向应力值先增大后减小，其最大值发生在1阶固频1 000 m管道处，应力值为740.00 MPa；升沉频率和升沉幅值对扬矿管轴向应力最大值发生的位置影响较大，当频率较低时，应力最大值一般发生在扬矿管顶部，当升沉幅值小于1.5 m，1阶共振时应力最大值发生在扬矿管1 000 m处，随着采矿船升沉频率和幅值的增大，最大轴向应力发生的位置不断变化.

图15 不同升沉幅值时管道最大轴向应力（ ω =1.92 rad/s）Fig.15 Maximum axial stress of lifting pipe at different heave amplitudes （ ω =1.92 rad/s）

由图16可知：扬矿管在单位升沉幅值和1阶固有频率下振动时，由于各级阶梯管之间横截面积的变化，轴向应力在1 000、2 000、3 500 m处发生突变变大，同时受中间矿仓质量的影响，管道轴向应力在5 000 m处到达最小值43.00 MPa.同一级管道随海水深度的增加，沿管长方向的轴向应力均逐渐降低.

图16 沿管长方向的轴向应力Fig.16 Axial stress along pipe length

4 结 论

针对深海采矿场景，建立了阶梯扬矿管纵向振动模型，获得了固有频率方程，并进行ABAQUS仿真分析，得到了扬矿管的纵向振动动力学响应，结果如下：

1） 扬矿管纵向振动的固有频率随中间矿仓质量的增加而减小，与低阶固有频率相比，扬矿管的高阶固有频率对中间矿仓质量的变化更加敏感.

2） 波浪频率对扬矿管的动态响应影响显著.随着频率的增大，扬矿管的动态响应先增大后减小，在1阶共振频率时，达到最大值.最大振幅出现在管道末端，其值为4.72 m，最大轴向力出现在管道顶端，其值为11.80 MPa，最大轴向应力发生的位置随波浪频率的变化而改变.

3） 采矿船的升沉幅值与扬矿管的动态响应呈正相关.当升沉幅值大于1.50 m时，扬矿管各节点动态响应的增长速度变缓，随着升沉幅值的增加，最大轴向应力发生的位置具有不确定性.

4） 1阶共振时，扬矿管各节点的动态响应在时域内先增大后作等幅稳态振荡.沿管长方向从上到下扬矿管的振动响应发生延迟，管末端比管顶端延迟0.7 s，产生约80° 的相位差.同时节点振幅呈现上小下大的趋势，节点轴向载荷和轴向应力随深度的增加而逐渐降低，其中轴向应力在每两阶管道间急剧变大.