基于离散元模型的扬矿管线系统动力学

分析方法研究及应用

李艳^{1, 2} , 刘少军¹，李辉¹，肖芳其¹，戴欢¹

(1. 中南大学 机电工程学院，湖南 长沙，410083；

2. 浙江大学 流体动力与机电系统国家重点实验室，浙江 杭州，310027)

摘要：扬矿管线系统的动力学分析对于整个深海采矿系统来说至关重要，为了克服有限元法、集中质量法计算工作量庞大、占用大量内存和机时的缺点，采用多刚体离散元方法对深海采矿扬矿管线系统进行三维动力学计算与分析，以提高计算速度，朝实时计算分析的方向发展。首先从静态分析和动态分析两方面给出2个典型算例，证明该方法的有效性；然后，将离散元方法用于中国1 km深海采矿系统的动力学分析，讨论不同海流方向及速度下系统的联动性能。计算结果表明：在计算精度相当的情况下，采用多刚体离散元方法比有限元方法计算速度快，所需存储空间小，对于结构动力问题的迭代求解具有实际意义。

关键词：深海采矿；扬矿管线；多刚体离散元；动力学分析

中图分类号：TD807          文献标志码：A         文章编号：1672-7207(2011)S2-0226-08

Dynamic analysis of long pipe system for deep-ocean mining by discrete element method

LI Yan^{1, 2}, LIU Shao-jun¹, LI Hui¹, XIAO Fang-qi¹, DAI Huan¹

(1. School of Mechanical and Electrical Engineering, Central South University, Changsha 410083, China;

2. The State Key Laboratory of Fluid Power Transmission and Control, Zhejiang University, Hangzhou 310027, China)

Abstract: The dynamic analysis of pipe system is one of the most crucial problems for the entire mining system. Discrete element method (DEM) is proposed for the analysis of deep-ocean mining pipe system. First two examples are given to show that the DEM model is successful. Then the DEM model is used to perform dynamic analysis of the Chinese 1 km deep-ocean mining system, which moves at different velocities under different current directions. Compared with those from finite element method, simulation results show that for long pipe system, the present DEM has faster computation speed with the equivalent precision.

Key words: deep-ocean mining; mining pipe; discrete element method; dynamic analysis

典型的深海采矿系统主要由采矿船、扬矿管线(包括扬矿硬管、中间舱、软管)和集矿机组成。深海采矿系统的扬矿管线长达数千米，其动态行为对采矿系统的整体联动性能影响巨大，扬矿管线系统的动力学分析对于整个采矿系统来说至关重要，可为系统的研制和开采作业提供理论依据和技术参考，因此，扬矿子系统的运动和力学分析在深海开采技术中具有重要作用。

1  扬矿管线动力学分析方法的比较

国内外在对深海采矿扬矿管线系统进行动力学分析时，扬矿管线的空间离散方法主要有有限元法、集中质量法和多刚体离散元法。

有限元法的分析方法基本为全拉格朗日法(TL)、更新拉格朗日法(UL) 和随体旋转法(CR) 3种。全部采用增量加载平衡迭代法求解。TL 法的应变、应力是相对于初始时刻的构形度量的，UL 法的应变、应力是相对于当前时刻的构形度量的，然后，再转换成总体坐标系。CR 法也是相对于当前时刻的，但它的坐标系是粘附于结构上，随之一起平移和旋转运动的刚性坐标系，因此，梁单元相对于该坐标系的运动便是纯粹的变形运动。这是它与 UL 法的区别所在。TL和UL 法因为结构运动中包含刚体旋转成分，在大位移、大转动中它们为有限量，因此，在采用增量法求解时，其加载步必须很小，小到可以按矢量操作，否则刚体运动的误差积累，造成发散，因而计算过程中需要花费大量的计算时间；而CR 法的坐标系是粘附于结构上，随后一起平移和旋转运动的刚性坐标系。因此，梁单元相对于该坐标系的运动为纯粹的变形运动，从而能将刚体运动从总位移中消除，使其余变形运动相对于局部刚性坐标系始终较小。故可以运用小变形梁理论进行有限元建模，不会因刚体运动的误差积累而造成发散，并大大减少增量迭代求解的程序计算量。但是，随体刚性坐标系的确定较困难，即使采用CR 法,其增量步也必须很小。尤其是在进行非线性动力响应分析时，往往需将结构分成足够多的单元，取相当小的时间步长，并用增量迭代法，在每步刚体矩阵的形成和组装过程中，数值积分的计算量很大，需花费大量的机时，导致其计算效率较低。

集中质量法又称集中参数法，一般采用应变-应力关系作为连续性条件。首先将管线离散为一定数量的有限个单元，将每个单元的质量集中在单元两端节点处，同时，所有的流体阻力、附加质量、重力和浮力也均集中于这些节点上。节点之间由线性无摩擦弹簧相连，至此管线被简化为弹簧-质量系统。然后，应用牛顿第二定理，形成一系列微分方程组，并通过相应的数学方法进行求解。集中质量法离散时要求节点间长度足够小，这样作用于上面的外力才可以认为是均匀的。其物理意义明确，算法简便。德国为印度设计的500 m深的采沙系统，全部采用软管输送，对软管进行分析时即采用此方法^[1-2]。但此方法在软管的抗弯、抗扭性能模拟方面均有不足，多用于锚链、拖缆等缆索在各种情况下的特性研究；且当采用显示求解算法时，为了满足稳定性条件，计算时间步长取得都比较小，需要大量的计算时间，尤其是在大幅剧烈运动时，否则将会导致数值计算发散。这是集中质量法的缺陷^[3]。

多刚体离散元法，对杆梁结构而言就是有限段方法(Finite segment method)，其基本思想是：将管线离散成有限个不发生弹性变形的刚性单元(简称刚体元)，以其表示物体的惯量特性。而其刚度向接点等效移植，相邻的刚体元之间用弹簧相连，从而建立管线结构的多刚体-弹簧系统的力学模型。然后，用多刚体动力学的方法建立该力学模型的动力学方程；最后进行数值求解。

采用多刚体离散元法所得到的模型具有简单性和力学的直观性的特点，可以用多刚体动力学的理论来描述单元的大位移、大转动，概念清晰。由于组装是对单元交界点进行的，因此，总刚度矩阵是一个高度稀疏的矩阵，只与单元编号排序有关，与结点编号无关，将显著减少计算量，节省计算机内存，这对大型结构的计算具有特别重要的意义^[4]。而且由这样的离散化模型得到的系统动力学方程，对结构的变形没有限制，而材料的非线性特性则反映在离散的非线性单元弹性系数上。因此，多刚体离散元方法可以有效地分析大位移、非线性动力学问题。

对多刚体离散元方法的深入研究是一个正在探索的课题，早期的研究多限于一、二维系统。Kawai    等^[5-6]利用多刚体离散元方法进行了平面问题和结构的地震响应分析；Connelly等^[7-10]利用该法成功地解决了链、缆索的大幅度振动问题及直梁的固有频率计算；殷学纲等^[11]解决了拱梁结构、圆环的几何非线性振动问题，随后又将该法用于杆梁结构的屈曲与后屈曲分析；王文军^[12]针对平面直梁进行了线性和非线性的静动力分析；任伟新等^[4]利用多刚体离散元法成功地进行了平面刚架的弹塑性、大位移极限承载力分析；谢先海等^[13]建立了平面柔顺机构的多刚体离散元模  型，解决了柔顺机构的特征频率和几何非线性变形  问题。

Mustoe等^[14-16]提出采用基于多刚体动力学的离散元方法计算分析扬矿硬管，针对二维深海采矿扬矿硬管进行了动态分析，并与早期的有限元计算结果进行了对比，克服了有限元方法求解速度较慢，且当边界条件变化时有可能出现不收敛的问题。随后试图开展三维深海采矿扬矿硬管动力学分析，但是，在当时的计算条件下，经过三四年的尝试之后，没有成功地推广到扬矿硬管的三维动力学分析。程保荣等^[17]采用该方法对三维深海采矿扬矿硬管进行结构模态分析，计算了立管的固有频率，同样未见进一步的三维深海采矿扬矿硬管动力学分析。

由于大洋多金属结核主要赋存于水深3~6 km的大洋底面，未来的深海多金属结核商业开采系统扬矿管线将长达数千米，考虑到有限元法处理管线系统，计算工作量庞大，占用大量内存和机时的缺点。而集中质量法在管线的抗弯、抗扭性能模拟方面均有不足，且当采用显示求解算法时，为了满足稳定性条件，计算时间步长取得都比较小，同样需要占用大量内存和机时。本文尝试采用多刚体离散元方法对扬矿管线系统进行三维动力学计算与分析，以提高计算速度，朝实时计算分析的方向发展。

2  三维多刚体离散元模型

深海采矿系统的扬矿硬管是由若干根硬管用管接头连接而成的，每一段硬管本身不能变形，可以视为刚性构件，但是，长达数千米的管线已经不能完全作为一个刚性构件考虑，必须考虑其大位移、小变形的几何非线性特性。至于柔性软管更应采用考虑大位移和旋转的分析技术来计算。

本文基于多刚体动力学的理论，采用离散元方法(Discrete element method，简称DEM法)对上述管线系统(包括硬管、中间仓和软管)进行分析。其基本思想是：将管结构离散成有限个不发生弹性变形的刚性单元(简称刚体元)，以其表示物体的惯量特性。而其刚度向接点等效移植，相邻的刚体元之间用无质量、无长度的柔性联接元相连。

如图1所示，第i个刚体元与第j个刚体元之间用柔性联接元相连，该柔性联接元可代表相应段的伸张、弯曲和扭转特性。两刚体元质心分别为C_i和C_j。刚体元j相对于刚体元i的运动以位移变量(u_x，u_y，u_z)和转动变量(θ_x，θ_y，θ_z)表示。刚体元i作用在刚体元j上的相对力和力矩分别以f_x，f_y，f_z和m_x，m_y，m_z表示。由结构力学理论，当变形较小时力和力矩分量与位移和转动分量之间的关系可表达为：

式中：E为弹性模量；G为剪切模量；A为刚体元的截面积；l为刚体元长度；I_y和I_z分别为截面对y轴和z轴的惯性矩；I_p为极惯性矩；脚标i和j指明属于哪段刚体元。

图1  管结构的离散元模型

Fig.1  DEM model of pipe

若各刚体元皆相同，即，，，,，，，则上述方程等号右边的系数矩阵简化为1个对角矩阵，并由方程解出刚体元i作用在刚体元j上的相对力和力矩如下：

进一步具体化，该柔性联接元等效于如下6个弹簧，即

1个纵向弹簧：；

2个横向弹簧：，；

1个扭转弹簧：；

2个弯曲弹簧：，。

在动力学分析中，刚体元只有动能而无弹性势能，而柔性联接元只有弹性势能而无动能。管结构的变形能仅仅储存在柔性联接元中，刚体元中任一点的位移可由刚体元质心的刚体位移来描述。取刚体i的质心在惯性参考系中的3个直角坐标和反映刚体方位的3个欧拉角作为笛卡儿广义坐标，即。采用Lagrange乘子法建立系统动力学方程：

        (1)

完整约束方程：。其中：为系统动能；为系统广义坐标列阵，；为广义力列阵；为对应于完整约束的拉氏乘子列阵；为系统广义速度列阵。

把式(1)写成更一般的形式：

，

，

                (2)

式中：为广义速度列阵；为约束反力及作用力列阵；为系统动力学微分方程列阵；为描述广义速度的代数方程列阵；为描述约束的代数方程列阵。推导出上述常微分方程和代数方程后，进一步应用Gear刚性积分算法对式(2)进行数值求解。

3  算例考证

本文从静态分析和动态分析两方面，分别选择典型算例，检验上述多刚体离散元模型的正确性与可行性。

3.1  某矩形截面悬臂梁的静态分析

如图2所示矩形截面悬臂梁，弹性模量E=2.0×10¹¹ Pa，长度L=1.5 m，高h=0.008 m，宽b=0.006 m。

图2  悬臂梁模型

Fig.2  Model of cantilever girder

首先考察该悬臂梁竖向振动的固有频率和横向振动的固有频率。表1所示为离散元法计算的前6阶固有频率与解析值的对比。可见：计算值与解析值误差较小。

然后，分以下2种情况考察该悬臂梁自由端的静态位移：(a) 考察仅有Y向载荷F_y=5 N作用时自由端的静态位移s_y；(b) 考察仅有Z向载荷F_z=5 N作用时自由端的静态位移s_z；表2列出了离散元法计算的自由端静态位移与解析值的对比。可见，该算法计算精度较高。

表1  悬臂梁的固有频率

Table 1  Natural frequency of cantilever

表2  悬臂梁自由端的静态位移

Table 2  Static displacement of free end for cantilever

3.2  某圆形截面简支梁的动态响应分析

设圆形截面简支梁的直径为8 mm，长为1.5 m，弹性模量为2.0×10¹¹ Pa，质量密度为7.8×10³ kg/m³。表3列出了离散元法计算的简支梁前6阶固有频率与解析值的对比。在简支梁中点处分别作用2种载荷：(1) 正弦载荷F=sin(20t)；(2)正弦载荷F=sin(44t)。考察中点处的响应，由于该简支梁系统的第一阶固有频率为7.070 Hz(44.4 rad/s)，在第2种载荷作用下会出现共振。

解析值如图3和图4中实线所示。采用离散元法，计算结果如图3和图4中虚线所示。可见：计算值与解析值较吻合。

表3  简支梁的固有频率

Table 3  Natural frequency of simply supported beam

图3  在正弦载荷F=sin(20t)作用下中点的响应

Fig.3  Dynamic response of middle point of simply supported beam subjected to load F=sin(20t)

图4  在正弦载荷F=sin(44t) 作用下中点的响应

Fig.4  Dynamic response of middle point of simply supported beam subjected to load F=sin(44t)

4  基于离散元模型的1 km海试系统整体联动动力学分析

本文针对中国大洋多金属结核1 km海试系统，应用多刚体离散元法进行扬矿管线系统动力学分析。应用多刚体离散元法，900 m长的扬矿硬管离散为90个刚体元，相邻刚体元用柔性联接元相连；400 m长的软管离散为200个刚体元，相邻刚体元亦用柔性联接元相连。扬矿硬管上端通过球铰与采矿船相连，软管下端通过“十”字铰与集矿机相连。每段管线系统所受液动力采用莫里森公式计算，并施加到每个刚体元的质心。

设定海洋表面海流速度为0.772 m/s，海底海流速度为0.15 m/s。集矿机与采矿船朝+X方向联动运行，且采矿船同时做升沉运动，幅度为4 m，周期为8 s。主要对表4所示6种典型工况(包括3种海流方向、2种速度)进行计算机模拟研究。

表4  6种典型工况

Table 4  Six typical conditions

图5所示为中间舱相对采矿船的X向偏移变化曲线。仿真结果表明：在同一海流方向下，系统运行速度越大，中间舱偏移越大，如工况a和d所示。海流与系统运行方向一致(即同为顺流运行)。a工况下系统以0.5 m/s的速度运行160 m后，中间舱相对采矿船的X向偏移为-22.5 m；d工况下系统以1.0 m/s的速度运行160 m后，中间舱相对采矿船的X向偏移为-61.9 m。

仿真结果也表明：在同一运行速度下，顺流运行时中间舱偏移较小，逆流时偏移最大，横流时偏移量居中。如d，e和f 3种工况，系统运行速度均为1.0 m/s；运行160 m后，中间舱相对采矿船的X向偏移分别为-61.9，-78.5和-96.2 m。

图6所示为扬矿硬管顶部X方向受力变化曲线。图7和图8所示为集矿机受软管作用力变化曲线。可见：6种工况下力曲线变化趋势一致，但在同一海流方向下，系统运行速度越大，系统受力越恶劣，如a和d 2种工况，海流与系统运行方向一致(即同为顺流运行)。a工况下系统以0.5 m/s的速度运行160 m后，软管对集矿机的作用力为4 031 N；d工况下系统以1.0 m/s的速度运行160 m后，软管对集矿机的作用力为-7 789 N。

此外，在同一运行速度下，顺流运行时受力较好，横流时受力状况居中，逆流运动时系统受力较恶劣。如a，b和c 3种工况，系统运行速度均为0.5 m/s，运行160 m后，集矿机X方向受软管作用力分别为4 031，500和-4 133 N，对于集矿机来说，顺流运行时软管对集矿机的作用力是拉力，而逆流时软管对集矿机的作用力是阻力。

总之，在高速逆流运动状态下，系统受力最恶劣，中间舱相对采矿船的X向偏移也最大。此时，采矿系统的形态已不能满足采矿要求，故可确定系统运行速度不能取1.0 m/s。同理，系统低速顺流运行是较理想的状态，0.5 m/s是比较适宜的速度。另外考虑到联动时逆流前进，即使在0.5 m/s的速度下中间舱偏移也有-51 m，因此，开始联动时，集矿机与采矿船的位置不能太近。

图5  中间舱相对采矿船的X向偏移变化曲线

Fig.5  Buffer X displacement relative to ship

图6  扬矿硬管顶部X方向受力变化曲线

Fig.6  Horizontal force of the pipe top

图7  集矿机X方向受软管作用力变化曲线

Fig.7  Horizontal force acted on miner by flexible hose

图8  集矿机Z方向受软管作用力变化曲线

Fig.8  Vertical force acted on miner by flexible hose

进一步，为了校验离散元方法的效率，又采用同样的参数用有限元法对表5所示的2个工况进行计算分析。计算均在相同配置的PC机(Pentium IV 1.8 GHz 内存512 M)上完成。

表5  离散元法与有限元法计算的关键结果与速度对比

Table 5  Comparison of key results and execution time

从表5可见：在关键计算结果误差仅为1%~2%的情况下，采用离散元方法计算速度较有限元方法提高近8倍。以上研究结果表明：对深海长管线这类形状较为简单的杆系结构，在不需要特别讨论应力、应变等情况时，采用离散元法进行动力学分析是一种可行的、有效的分析方法。

5  结论

(1) 通过算例证明本文采用多刚体离散元方法对管线系统进行计算分析是行之有效的。

(2) 将多刚体离散元法应用到深海采矿整体联动系统的动力学计算与分析，研究结果与有限元法分析结果，基本一致，可为系统的研制和开采作业提供理论依据和技术参考。

(3) 在计算精度相当的情况下，采用多刚体离散元方法比有限元方法计算速度快，所需存储空间小，对于结构动力问题的迭代求解具有实际意义。

参考文献：

[1] Ravindran M, Schwarz W, Jeyamani R, et al. Shallow-water sand- mining operation[J]. Proceedings of the Annual Offshore Technology Conference, 1999(1): 83-101.

[2] Deepak C R, Shajahan M A, Atmanand M A, et al. Developmental tests on the underwater mining system using flexible riser concept[C]//Proceedings of the ISOPE Ocean Mining Symposium. Szczecin, Poland, 2001: 94-98.

[3] 王飞. 海洋勘探拖曳系统运动仿真与控制技术研究[D]. 上海: 上海交通大学船舶海洋与建筑工程学院, 2006: 2-6.
WANG Fei. Simulation and control research of marine towed seismic system[D]. Shanghai: Shanghai Jiaotong University. School of Naval Architecture, Ocean and Civil Engineering, 2006: 2-6.

[4] 任伟新, 郑兆昌. 平面刚架弹塑性大位移分析的多刚体离散元法[J]. 计算结构力学及其应用, 1996, 13(4): 408-417.
REN Wei-xin, ZHENG Zhao-chang. Multi-rigid body discrete element method of large displacement elasto-plastic analysis for plane frames[J]. Computational Structural Mechanics and Applications, 1996, 13(4): 408-417.

[5] Kawai T, Toi Y. A new element in discrete analysis of plane strain problems[J]. Seisan Kenkyu, 1977, 29(4): 204-207.

[6] Kawai T. New discrete models and their application to seismic response analysis of structures[J]. Nuclear Engineering and Design, 1978, 48(1): 207-229.

[7] Connelly J D, Huston R L. Dynamics of flexible multibody systems: A finite segment approach (Ⅰ). Theoretical aspects[J]. Computers and Structures, 1994, 50(2): 255-258.

[8] Connelly J D, Huston R L. Dynamics of flexible multibody systems: A finite segment approach (Ⅱ). Example problems[J]. Computers and Structures, 1994, 50(2): 259-262.

[9] Huston R L, Kamman J W. Validation of finite segment cable models[J]. Computers and Structures, 1982, 15(6): 653-660.

[10] Kamman J W, Huston R L. Modelling of submerged cable dynamics[J]. Computers and Structures, 1985, 20(1/3): 623-629.

[11] 殷学纲, 杜思义, 胡继云, 等. 应用多刚体离散化方法对梁的动力屈曲和后屈曲分析[J]. 应用数学和力学, 2005, 26(9): 1076-1082.
YIN Xue-gang, DU Si-yi, HU Ji-yun, et al. Analysis of dynamical buckling and post buckling for beams by finite segment method[J]. Applied Mathematics and Mechanics, 2005, 26(9): 1076-1082.

[12] 王文军. 直梁静动力分析的多刚体离散元法[D]. 北京: 清华大学工程力学系, 1994.
WANG Wen-jun. Multi-rigid-body discrete element method in static and dynamic analysis of straight beam[D]. Beijing: Tsinghua University. Department of Engineering Mechanics, 1994.

[13] 谢先海, 廖道训. 基于多刚体离散元模型的柔顺机构动力分析新方法[J]. 机械工程学报, 2003, 39(5): 10-15.
XIE Xian-hai, LIAO Dao-xun. New analysis method for compliant mechanisms based on multi-rigid-body discrete element model[J]. Chinese Journal of Mechanical Engineering, 2003, 39(5): 10-15.

[14] Mustoe G G W, Huttelmaier H P, Chung J S. Assessment of dynamic coupled bending-axial effects for two-dimensional deep-ocean pipes by the discrete element method[J]. International Journal of Offshore and Polar Engineering, 1992, 2(4): 289-296.

[15] Mustoe G G W. Generalized formulation of the discrete element method[J]. Engineering Computer, 1992, 9(2): 181-190.

[16] Huttelmaier H P, Chung J S, Mustoe G G W. Eigenvalues of a long vertical pipe by DEM, FEM and exact solution[C]//Proc Third Int Offshore Polar Eng Conf. Singapore, 1993: 311-314.

[17] 程保荣, 郑兆昌. 三维梁系结构刚体元—柔性连接元动力分析模型[J]. 清华大学学报: 自然科学版, 1996, 36(3): 12-17.
CHENG Bao-rong, ZHENG Zhao-chang. Discrete element with flexible connector method for dynamics analysis in the 3-D beam structure[J]. Journal of Tsinghua University: Natural Science, 1996, 36(3): 12-17.

(编辑 陈灿华)

收稿日期：2011-06-15；修回日期：2011-07-15

基金项目：浙江大学流体动力与机电系统国家重点实验室开放基金资助项目(GZKF-201001)；国际海底区域研究开发“十一五”项目(DYXM-115-04-02-01)

通信作者：李艳(1975-)，女，湖南永州人，博士，副教授，从事机电液系统控制理论与技术、深海作业装备设计与控制研究；电话：13677353895；E-mail: lylsjhome@163.com