离散元法用于绞刀齿切削珊瑚礁岩过程的三维模拟

林挺 ，冯沛洪 ，马源 ，郭涛 *

（1．上海交通大学船舶海洋与建筑工程学院，上海 200030；2．中交疏浚技术装备国家工程研究中心有限公司，上海 201208）

绞吸挖泥船是一种典型的疏浚工程船，作业时绞刀通过桥架下放到水底，绞刀旋转并松散、破碎水底的密实砂土或岩体，再由管路、泥泵等输送系统抽吸、排送砂土、碎岩与水的混合浆体，实现疏浚或吹填的目的。近年来，国内疏浚行业发展迅速且不断承接国际的港口或海洋工程，沿海疏浚或岛礁吹填项目中不乏挖掘珊瑚礁石的工程。未来，水下挖岩的工程量可能进一步增加，设计专用挖掘礁岩、风化岩等岩土的绞刀、刀齿很有必要。

当前，挖岩绞刀的功率预估、结构设计、强度校核主要依靠同类施工经验，这很可能造成过量设计与匹配浪费，也对新型挖岩绞刀的优化与设计造成了困扰。但是绞刀齿与岩石接触并破碎岩石的具体受力难以通过绞刀轴上宏观表现来考察，即使通过物理实验也难以直接测量，而且物理实验具有耗时长、费用大等缺陷。

离散元方法（Discrete element method，DEM）源起上世纪70年代，可将物质处理为由不连续的颗粒组成的集合，通过求解各颗粒的运动方程实现模拟物质的运动与变形，是岩土力学领域中的典型数值研究方法[1]，已在海底岩体破碎等研究中得到良好应用[2]。本文将采用DEM三维求解绞刀挖掘珊瑚礁岩时刀齿上的受力情况，为挖岩绞刀的设计提供参考。

1DEM物理模型

处理弹性颗粒球形面的变形时，DEM广泛采用考虑切向库伦摩擦力[3]以及回复阻尼系数[4]的Hertz-mindlin模型。该模型中，相互接触的球面1、球面2因弹性变形而产生的相互作用力为：

式中：Fn、Ft分别为法向力和切向力；E*、G*分别为等效杨氏模量、等效剪切模量；R*为等效半径；δn、δt分别为两球面的法向重叠距离与切向重叠距离，即两球面的球心距离相较于初接触时刻（此时两球面点接触，球面无变形，亦无弹性接触力）的径向变化量以及两球接触面中点相对于初接触时的切向变化量；μs为颗粒间的静摩擦系数。并有：

式中：υ为泊松比；其他各字母含义同上，下标1、2分别代表不同球面。可见，两球面间剪切力的

式中：ρ为颗粒密度。实际上对于高配位数的准静态颗粒堆积（比如三维球形颗粒），通常时间步小于 0．2tRay是合适的。

Hertz-Mindlin模型没考虑颗粒间的黏聚力，并不适于表征具有黏聚力的岩石类物质。Potyondy与Cundall在此基础上提出了一种用于岩石的颗粒粘合模型[6]，该模型类似于在相近颗粒之间填充了“固体粘合剂”，可以较大程度上模仿岩石类物质。本文使用该模型来表征珊瑚礁岩中各颗粒间的相互作用。该模型假定在初始时刻，若相邻离散颗粒的球面中心距离小于接触半径Rc，则球面之间将存在一种粘合作用（bonding），本文称为键。键的两端分别为两球面中心，颗粒发生相对运动时，作用键使得球面间存在类似线性弹簧（包括力与转矩）的相互影响，直到键被破坏。键对颗粒的作用力变化关系满足式（4）：

式中：v为颗粒间相对速度；S为单位面积键刚度，其量纲为N·m-3；A为键的横截面积；t为时间间隔；下标n、t分别代表法向、切向分量。作用键在特定条件下会被破坏而导致“粘合效应”失效，这种破坏不可逆。作用键的失效准则包括：键上承受了极限法向应力（CNS）、极限切应力（CTS）或者键达到极限长度。当存在作用键时，时间步长需要考虑键力的传播。

2 珊瑚礁岩的力学性能及DEM模拟参数

珊瑚礁是由珊瑚虫钙质遗骸及其他造礁生物的石灰质长期累积而成的岩土体，主体为礁灰岩，在中国南部沿海地区分布广泛。由于珊瑚虫种类繁多，珊瑚礁石形成过程中环境不断变化，礁灰岩的岩体结构复杂，其力学性能也差异明显。中科院王新志等人[7]曾对永暑礁及渚碧礁的礁灰岩进行了采样及力学研究，得到礁灰岩力学特性如表1所示[7]。研究中还发现礁灰石内部孔隙率较高，极限不能超过最大静摩擦力，否则球面间将发生相对滑动或转动。

为保证计算的准确性，需要确保时间步明显小于瑞利时间步——剪切波穿过一个颗粒的时间。瑞利时间步 tRay由式（3）[5]确定：干燥礁岩的天然密度约为（1．25～1．35）× 103kg/m3，孔隙率约为45%～55%；礁石压裂后会存在残余强度，随着压裂破坏的增加强度进一步降低，这与脆性灰岩不同。

表1 采样礁灰岩力学特性Table 1 The mechanical property of the sampled coral reef rocks

为使DEM模型能够准确模拟水下饱和珊瑚礁灰岩的力学性能，按DL/T 5368—2007《水电水利工程岩石实验规程》中的实验规程，对礁石的巴西劈裂法测抗拉强度实验以及单轴抗压强度实验进行相同的DEM数值实验，并根据实验结果调整计算参数。计算中，基本离散粒子选为球形；为保证颗粒直径总是小于试件直径的1/10，球形粒径设为4 mm；接触半径为3 mm；为克服球形颗粒易滚动的缺陷，限制颗粒的旋转自由度；基于Mohr-Coulomb强度理论，假定破碎的主要原因是岩体抗剪切强度不足。

经过多次“选择参数-计算-修改参数”的循环后，选用表2所示的参数，得到与表1中实验结果均值相近的礁岩模型。

表2 模拟礁灰岩颗粒参数Table 2 Particle parameters of the simulated coral reef rocks

巴西劈裂法测抗拉强度的实验中，圆柱形试件的周面上、下两侧各有宽约4 mm的垫条缓慢向内挤压，图1给出了实验前后的键变化。图2给出了实验中试件承受拉力随垫条位移的变化过程。图中可见，礁岩在不足1 mm的变形量下即发生突然的张性断裂，断裂面基本经过上下的垫条，其抗拉强度约为1．2 MPa。

图1 巴西劈裂实验前、后键状态Fig.1 Statues of the bonds before and after Brazilian splitting test

图2 巴西劈裂数值实验计算结果Fig.2 Numerical results of the Brazilian splitting test

图3给出了礁石单轴抗压强度实验前、后的键变化。图4给出了单轴抗压强度的数值计算结果。从图中可见，数值模拟的礁岩圆柱在轴向形变约1 mm时达到最大单轴抗压强度，约为6 MPa。此后压力随着轴向形变增大而迅速减少，这点与实验中压力阶梯式减小的真实情况不符，这可能是数值模拟中忽略了真实礁岩材料内不均匀分布的孔隙所造成的。尽管数值结果在应力衰减上存在偏差，但在极限应力、极限形变上与实验情况非常相近，仍认为当前计算参数能够表征珊瑚礁岩的主要力学性能。

图3 单轴压裂实验前、后键状态对比Fig.3 Statues of the bonds before and after uniaxial compression test

图4 单轴压裂数值实验计算结果Fig.4 Numerical results of the uniaxial compression test

3 绞刀齿切珊瑚礁岩的三维计算与分析

绞吸船工作时，桥架上的绞刀探入水下并绕轴旋转，带动刀齿、刀臂切削岩土。船体一边由横移钢缆牵动绕工作钢桩转动，一边由台车推动改变船体与工作钢桩的距离，同时桥架受岩土的反作用力变化，导致其下探角度略微波动。因考察局部刀齿挖掘过程，绞刀的变半径、变深度的弧形运动可以近似简化为绕轴旋转以及水平横移的叠加，如图5所示。

图5 简化后绞刀运动示意图Fig.5 Schematic diagram for the simplified cutter movement

DEM的计算量随粒子数增加而指数级上升。受限于计算资源，本文仅计算了某大型六臂绞刀的第4排相邻刀齿以大挖岩横移速度挖掘礁岩时的受力情况。该绞刀额定功率2 000 kW，最大旋转半径约1．74 m。由于安装半径、安装角度不同，不同刀齿的挖掘面积也不同。以本绞刀第4排刀齿为例，其相邻刀齿错开排列，齿尖半径分别为1．232 m与1．309 m，尖端半径较小的齿（后称齿A）在轴向上较邻齿（后称齿B）前突约0．1 m，相邻齿尖的周向夹角分别为67°及53°。绞刀旋转速度22 r/min、横移速度15 m/min，图6给出向上挖掘礁岩时相邻的齿A、齿B受力的计算模型。图中模型的颗粒数约为85 000个，虚线为齿尖的运动轨迹。由图可见，受刀齿错开排列的影响，齿A在横移面上的挖掘面明显小于齿B，但齿A在轴向上突前，轴向挖掘深度明显大于齿B。

图6 计算刀齿切削过程模型Fig.6 The model for simulating the cutting processes

计算时，先由齿A破环礁岩表面并留下挖掘痕迹，考察后续相邻的齿B、齿A继续切削礁岩的受力情况。珊瑚礁岩按图6中的样式构建以减少离散粒子数量，礁岩轴向厚度0．28 m，绞刀轴位于礁岩上沿0．7 m。

为了提高计算效率，根据式（1）、（3）、（4）以及相似准则，调整计算参数如表3所示。其他计算设置与前文一致，计算结果与选用表2参数的结果相似，相似比例为1∶1 000，计算耗时为后者的1/30。计算的时间步长为2×10-6s，总时长为1．2 s，使用Intel I7-4790处理器满负荷计算共耗时约30 h。图7给出了计算后礁岩的破坏情况。

表3 等效计算时选择的参数设置Table 3 Parameters used in the equivalent calculations

图7 刀齿切削后礁岩模型Fig.7 The reef model after cutting with cutter teeth

图7中，礁岩上刀齿所过处礁岩破裂，部分碎裂成粒子，部分仍保持团块状，刀齿未经过处的内部也存在着明显裂痕。在顶部，齿A在轴向上基本挖穿了礁岩；在底部，齿B的挖掘深度约为齿A的一半。

设定齿B在挖掘起点处为绞刀的0°相位角，相位角将随着绞刀旋转而增加。图8给出了计算得到的齿A与其相邻的齿B在挖掘过程中所承受的礁岩作用力随绞刀旋转的变化。

从图中可见，齿A与齿B上承受的峰值力分别约为5×104N与9×104N。齿B上承受了较高的作用力，这是由于齿B旋转半径较大，齿A切削范围内的部分礁岩颗粒已被齿B破坏了内部键。

图8 刀齿上切削合力随相位角变化Fig.8 Cutting forces on the teeth vary with rotating anlges

假设绞刀第4排的3对A、B齿在挖掘过程的受力情况相同，根据刀齿受力矢量和相对位置可以折算出该排刀齿整体所受的挖掘力矩，该值以120°为周期往复性变化。图9给出了2个周期内的估算力矩及拟合曲线。可见，绞刀第4排刀齿承受的力矩大致在（1～9）×104N·m的范围内，在齿B开始切削后又旋转了约42°时达到最大值，在齿A完成切削时达到最小值。

图9 第4排刀齿挖掘力矩随相位角变化Fig.9 The cutting torque of the fourth row teeth vary with rotating anlges

图8中，刀齿B在相位角约为40°～50°时受力最大，计算此时间段内的刀齿表面承受压力的时均值，采用插值映射的方法将该值作为刀齿载荷施加在有限元软件上，考察作用力可能引起的刀齿应力分布及变形。刀齿材料选择为某低合金耐磨钢，计算结果如图10、图11所示。图中，礁岩对刀齿的压力主要集中在齿尖上，在齿B承受最大受力时，齿尖压力峰值高达1．14×108Pa。此时，齿尖处将发生最大约0．54 mm的变形，齿背上加强筋处的von-Mises等效应力最大，约1．68×108Pa，同理此处的安全系数最小，约为7．1。

图10 刀齿最大受力时的变形Fig.10 The deformation of the cutter teeth under maximum force

图11 刀齿承受的最大等效应力Fig.11 The maximum equivalent stress on the cutter teeth

4 总结与讨论

使用DEM三维模拟了绞刀第4排刀齿挖掘珊瑚礁岩的过程，量化求解了挖岩机具的受力情况，并对计算结果进行了有限元分析，为挖岩绞刀或类似机械的设计、校核与优化提供了一种数值解决方案，得到以下结论：

1）Hertz-Mindlin接触模型与颗粒粘合模型，可以更有效地模拟珊瑚礁岩等岩体的内聚力效应，模拟结果更接近真实情况。

2）本文确定的DEM参数可以较准确地表征珊瑚礁岩的力学性能，并能较高效地处理绞刀挖掘礁岩的受力、破碎问题。

3）计算的刀齿挖掘礁岩时的最大受力约为9×104N，最大等效应力约为 1．68×108Pa，佐证了该刀齿的结构设计与安装角度满足挖掘礁岩的强度要求。

受限于资源耗费与计算效率，本文在颗粒形状、粒子尺度、粒子参数、模型范围等方面尚有不足，例如未使用不同形状、非均匀粒径的粒子以及非均匀的键值等方式来模拟礁岩内的裂隙现象，进一步研究工作可在这方面做出改善，实现更高精度的计算，并推广到整个绞刀挖岩过程的模拟。