李艳，陆鹏，李钰晴

(1.中南大学机电工程学院，湖南长沙，410083；2.深海矿产资源开发利用技术国家重点实验室，湖南长沙，410012)

与陆地矿相比，深海矿产资源具有极高的金属品位。目前，越来越多的国家进入工业化发展阶段，对金属的需求不断加大，开发深海矿产资源是大势所趋。海底块状硫化物(seafloor massive sulfide，SMS)由于其具有赋存水深较浅、距离陆地较近、经济价值较高等优点被国内外学者认为是最有可能成为深海采矿的首选对象[1−4]。海底块状硫化物矿床多位于1 000～3 000 m 海平面以下，承受10～30 MPa的海水围压作用。海底块状硫化物在开采过程中，由于深海环境复杂并且存在海水围压作用，其剥离破碎技术是开采过程中一大难点。目前，国内外学者认为类似于陆地采煤机螺旋滚筒的截齿切削是一种最佳的块状硫化物切削方法[5−6]，其中，端盘位于深海采矿车滚筒的端部，端盘截齿数超过滚筒截齿总数的1/2，工作环境极其恶劣，截齿负荷高，消耗的功率占工作机构总消耗的1/3 左右[7]。为了改善滚筒截割过程中的受力条件，增强采矿车工作稳定性和使端盘参数与工作过程自适应，必须掌握高海水围压条件下端盘截齿切削块状硫化物矿体过程的载荷特性。PARASIE等[8]通过离散单元法建立了高静水围压下硬岩切削过程的模型，该模型揭示了在围压作用下，岩石裂纹扩展能力变弱，切削阻力大大增加；HELMONS 等[9]使用离散单元法研究了围压和孔隙压力对于饱和岩切削进程的影响，发现围压越大，Peclet 数对切削力的影响越大，而围压和Peclet 数的增加均会使得切削力增大。VERCRUIJSSE 等[10]通过离散单元法研究海底钴结壳的破碎机理，发现在海水围压的作用下海底岩石的塑性增强，切削时裂纹的产生和扩展变得更加困难；LI等[11]基于离散单元法研究围压对于岩石破碎进程的影响，发现随着围压和岩石强度增加，裂纹沿法向的传播受到抑制，切削力增加，截割比能耗增加；李艳等[12−13]采用有限单元法研究了不同围压对块状硫化物破岩机理的影响，发现随着围压增加，截齿受到的截割阻力和进给阻力都会显著增大；HUANG等[14]使用有限单元法研究了侧压对截齿破岩破碎力的影响，发现当侧压小于临界侧压时，平均峰值切削力随侧压增加而增加，而当侧压大于临界侧压时，平均峰值切削力随侧压增加而减小。宋扬等[15]通过有限单元法研究了端盘截齿旋转截割煤岩的动态过程，发现端盘截齿更容易出现疲劳断裂，安装参数存在最优单侧截割工况；PARK等[16]采用线性切割机试验与有限元分析相结合的方法，从破碎效率、结构稳定性和耐用性研究截齿的最优截割工况，发现正倾斜角对提高切削效率和结构稳定性有利；郝志勇等[17]采用有限单元法建立了采煤机端盘截齿斜切截割模型，发现随着切向安装角的增大将会极大地影响滚筒截割力矩。由于离散单元法的二维模型无法考虑围压下端盘截齿破岩过程的侧向力，而其三维模型的计算量巨大，因此，使用较多的仿真方法是三维有限元方法[18]。影响三维有限元方法计算端盘截齿切削SMS 矿体精度的关键因素是岩石本构模型，目前采用较多的岩石本构模型是Mohr-Coulomb(M-C)模型、Drucker-Prager(D-P)模型和Holmquist-Johnson-Cook模型(HJC 模型)。M-C 模型在区分岩石拉压强度的同时，能够考虑静水压力的影响，但由于其在主应力空间中的屈服面是六棱锥形状，具有奇异点，会导致计算结果不收敛；D-P模型虽然考虑了中间主应力的作用，但不能反映π 平面上的拉伸子午线和压缩子午线的不同。通过比较，HJC 模型能较好地反映SMS 矿体高孔隙率性质特点和围压对材料的影响[12−13,19]。本文作者基于ABAQUS 平台，借助VUMAT 接口编制多金属硫化物HJC本构模型代码。自定义的HJC本构模型可以充分考虑实验数据所描述的矿体力学特性。在原始HJC 模型的基础上，考虑切削过程中岩石的损伤和强化；损伤法则定义为以等效塑性应变和塑性体积应变为变量的损伤累积，强化法则定义为各项同性强化；流动特性定义为关联塑性流动。通过对比SMS 样品的三轴压缩实验结果与实验模型结果一致的三轴压缩仿真模型，验证自定义本构模型的有效性。然后，建立三维截齿旋转切削模型，并使用经过验证的自定义本构模型作为SMS 材料本构模型，研究围压、圆周切向安装角、轴向倾斜角和二次旋转角对破岩载荷和载荷波动系数的影响规律。

1 矿体本构模型

1.1 本构模型定义

图1所示为岩石的应力−应变曲线的3个阶段：弹性阶段、塑性阶段和断裂损伤阶段，其中，D为损伤因子。本构模型的定义便是要针对这3个阶段进行定义。图2所示为截齿切削块状硫化物时计算每一个矿体网格单元应力－应变状态的积分推导过程。选用强健的径向返回算法，算法基于弹性预估，然后通过塑性校正对比试算应力和屈服应力，通过更新应力等张量完成积分推导。在推导过程中，所有参与推导的变量均假定为时间增量步末尾t+∆t时对应的变量，即时间增量步末尾的应力表示为σt+Δt，而增量步起始的应力表示为σt，因此，SMS 矿体HJC 本构模型从以下3个方面描述：应力更新算法、损伤和单元失效删除。

图1 矿体应力−应变曲线Fig.1 Stress−strain curve of orebody

1.1.1 应力更新算法

在材料试算过程中，应变ε分为弹性应变εe和塑性应变εp，即

图2 VUMAT积分推导流程Fig.2 Integral derivation flow of VUMAT

在弹性阶段，试探应力张量由弹性应变增量引起，由广义Hook定律得出：

其中：λ为lame常数，；E为弹性模量；μ为泊松比；trace为张量矩阵求迹函数；I为转动惯量；Δε为应变增量。而试探应力的Mises等效应力计算式为

其中：Sij=σij-δijσkk/3，为试探应力的偏应力张量；系数为增量步t时试探应力的主应力。当各应力应变张量之间满足屈服条件时，材料进入塑性阶段，屈服面函数σy可以表示为

此时，材料已发生屈服，继续增大载荷，根据关联塑性流动法则可以得到应力变化引起的屈服面的法线方向与塑性增量的关系。图3所示为塑性流动法向法则示意图，其中，dεp为塑性应变增量。

塑性应变增量dεp和等效塑性应变张量dεˉp的计算表达式如下：

图3 塑性流动法向法则示意图Fig.3 Schematic diagram of plastic flow rule

进而得到塑性应变非负标量的计算式：

其中：αt为背应力张量；H为塑性模量“:”为2个张量的双点积运算符。

结合应力补偿更新算法，得到应力张量、背应力张量和等效塑性应变的更新表达式：

其中：G为剪切模量；G=E/[2(1+μ)]。

1.1.2 损伤和失效删除

块状硫化物矿体的破碎过程由塑性体积应变和等效塑性应变引起的损伤累积控制：

其中：D为损伤因子；Δεp和Δμp分别为1次计算中等效塑性应变增量和塑性体积应变增量；为常压下破碎的塑性应变，；D1和D2为损伤常数；P∗=p/fc为标准化静水压力，p为实际静水压力；fc为岩体矿石材料的准静态单轴抗压强度；T*=T/fc，为标准化峰值抗拉应力，T为材料可以承受的最大拉伸强度。当损伤因子累积超过材料所能承受的最大损伤因子D时，此时，矿体材料积分点传递的质量、刚度、应力和应变增量变为0，这样，其在整体结构计算中就不再发挥作用而退出工作机制。本构模型的推导过程采用Fortran 语言编制成代码，每次迭代计算时，ABAQUS通过VUMAT接口调用这段代码。

1.2 VUMAT单元验证

为验证VUMAT的合理性，采用单个单元立方体模型进行三轴压缩测试。数值模型的边界条件如图4(a)所示。载荷加载方式为轴向位移加载，等效应变速率为1.0 m/s。在ABAQUS/Explicit中实施计算，结合通用动力非线性有限元软件ANSYS/LS-DYNA 计算结果和块状硫化物样品三轴压缩实验结果[20−21]，分别提取计算得到的不同围压下单元的抗压峰值强度，计算结果如表1所示。

图4 三轴测试仿真模型和仿真结果Fig.4 Triaxial test numerical model and results

表1 抗压峰值强度仿真值与实验值Table 1 Simulation and experimental values of peak compressive strength MPa

抗压强度−围压曲线对比如图4(b)所示。从图4(b)可看出：在HJC 模型参数一致的条件下，ABAQUS/Explicit和ANSYS/LS-DYNA仿真拟合曲线均与实验拟合曲线较吻合，可认为自定义HJC材料子程序能更好地模拟块状硫化物在压缩损伤时的力学特性。

2 仿真模型的建立

2.1 端盘截齿安装参数

端盘截齿安装角度如图5所示。端盘截齿安装参数主要包括圆周切向安装角β0、轴向倾斜角β1(一次旋转角)和二次旋转角β2。圆周切向安装角β0为齿身轴线与截齿齿尖切线正交线间的夹角，其值会影响截齿和齿尖与矿体接触形式、截齿受力。轴向倾斜角β1为齿尖截面(通过滚筒轴线并且贯穿截齿齿尖)上齿身轴线与螺旋滚筒径向线之间的夹角，若其取值合适，则可以使齿座与矿体互不干涉，避免齿座磨损。二次旋转角β2为齿尖截面正交面上齿身轴线与滚筒径向线之间的夹角，合理的二次旋转角可以减少齿座磨损。为了分析各安装参数对端盘截齿破岩破碎性能指标的影响，采用单因素变化法，每个因素按照工程施工经验取值，其中，切向安装角取3个不同值，轴向倾斜角取4个不同值，二次旋转角取2个不同值。数值模拟时，安装参数默认值为表2中括号中的数值。

图5 端盘截齿安装角度Fig.5 Installation angles of transverse flange pick

表2 截齿安装角度Table 2 Installation angles of transverse flange pick(°)

2.2 有限元模型的建立

端盘截齿采用MT/T246—2006 采掘机械用标准，滚筒直径为800 mm，截齿长为152 mm，合金头直径为18 mm，锥角为80°。块状硫化物矿体界面沿切削轨迹方向弧长为160 mm，宽度为 200 mm，厚度为100 mm，矿体本构模型调用自定义HJC 本构模型代码模拟镐形截齿旋转破碎SMS 矿体的动态过程，主要材料参数如表3所示。

表3 截齿和海底块状硫化物材料属性Table 3 Material properties of pick and SMS orebody

在仿真过程中，截割厚度固定为30 mm。为模拟围压环境，约束矿体的位移和旋转自由度，并在截齿与矿体接触的自由面上施加均布压力载荷。根据块状硫化物矿床所处海底深度，取围压范围为10～30 MPa。矿体单元设置类型为C3D8R，并且将截齿和矿体侵彻部分网格进行加密以提高计算精度。端盘截齿旋转切削矿体有限元模型如图6所示。

图6 端盘截齿旋转切削SMS有限元模型Fig.6 Finite element model of transverse flange pick rotary cutting SMS

2.3 载荷提取与处理

端盘截齿切削矿体载荷按照方向可分为：与切削轨迹相切的截割阻力；沿深度方向的进给阻力；与切削轨迹面垂直的侧向力。端盘截齿旋转切削载荷示意如图7所示，截齿三向载荷可通过换算得到，换算公式如下：

其中：FR为截割阻力；FN为进给阻力；FS为侧向阻力。

图7 端盘截齿旋转切削载荷分解示意图Fig.7 Schmatic diagram of load decomposition of transverse flange pick

在端盘截齿切削矿体过程中，其载荷特性在空间上呈现为波动无常的随机载荷，在工程作业中可能发生强烈震动，对截齿寿命产生极大影响，严重时可能发生崩齿状况。研究表明，端盘截齿工作在封闭或半封闭状态，工作条件恶劣，侧向力波动将会直接影响截齿单侧截割破碎效率，故引入载荷波动系数δ作为SMS 矿体切削效率的评价指标：

其中：为载荷均值；n为总的数据点数；Fi为第i个数据点处的载荷；σF为载荷标准差。

3 仿真结果及分析

3.1 围压对破碎性能的影响

端盘截齿在不同围压环境下破碎块状硫化物矿体时，其三向载荷和载荷波动如表4所示。在不同围压下，截齿三向阻力随时间的变化如图8所示。由表4和图8可以得到：随着围压增加，三向阻力增大，三向阻力波动系数明显减小，这主要是由于高海水围压直接导致矿体的宏观力学性能增强，抗拉及抗压强度增大，裂纹难以产生，同时，裂纹扩展也受到高海水围压抑制，切屑更难以从基岩剥离；随着海水围压增加，高压海水不但会抑制截齿的振动，而且会使块状硫化物塑性发生变化，使得截齿破岩时被切屑所包裹，周围的海水对累积在刀具前刀面的切屑有很大的压力和粘着力，最终使截齿三向阻力增大，载荷波动性减小。

表4 不同围压下破碎性能统计Table 4 Crushing performance statistics of different confining pressures

3.2 切向安装角对破碎性能的影响

图8 不同围压时端盘截齿三向阻力曲线Fig.8 Curves of three-direction resistance with time under different confining pressures

在20 MPa 围压下，端盘截齿以不同圆周切向安装角进行硫化物矿体破碎时，其三向载荷和载荷波动情况如表5所示。在不同切向安装角下，切削过程中三向阻力随时间的变化见图9。结合表5和图9可以看出：随着切向安装角增大，截齿截割阻力FR均值呈现先减小后增大的趋势，而截齿进给阻力FN和侧向阻力FS变化不大，说明圆周切向安装角主要影响截齿的截割阻力；当圆周切向安装角β0=45°时，截齿截割力均值有极小值且存在较小的波动，易达到最高效率。这是因为切向安装角越大，越不利于截齿截入矿体，切向安装角越小，截齿后刀面与矿体的摩擦越严重；当切向安装角取中值时，截齿垂直楔入岩石和横向楔入岩石的力分配更加合理，这样，不但使矿体楔入部分更容易破碎，而且使矿体更容易剥落。

图9 不同切向安装角时端盘截齿三向阻力曲线Fig.9 Curves of three-direction resistance with time under different tangential installation angles

3.3 轴向倾斜角对破碎性能的影响

当围压为20 MPa 时，在不同轴向倾斜角下，端盘截齿三向载荷和载荷波动见表6。在不同轴向倾斜角下，切削过程中三向阻力随时间的变化曲线见图10。对比表5和表6、图9和图10 可以得到：当轴向倾斜角为0°时，侧向力FS均值很小，基本沿z轴两侧波动；而当轴向倾斜角取5°～15°时，侧向力FS波动倾向于沿z轴负向变化，截齿出现单侧受力现象，说明截割时两侧煤岩不同时崩落，产生了侧向力差值，并随着轴向倾斜角的增大而更加明显；当端盘截齿β1=10°时，截割阻力均值出现极小值，此时，矿壁剥落使截齿单侧受力明显，截割效率最高。

3.4 二次旋转角对破碎性能的影响

二次旋转角通常取0°或5°，当围压为20 MPa时，在不同端盘截齿二次旋转角下，截齿三向载荷和载荷波动情况见图10。在不同二次旋转角下，切削过程中三向阻力随时间的变化曲线见图11。从表7和图11 可以看出：随着二次旋转角增加，端盘截齿沿二次旋转方向逐渐倾斜，侧向力FS均值和峰值沿z轴负向增大明显，说明在端盘截齿截割过程中，齿尖侧矿石剥落更加容易，可以更好地形成截割壁面，有利于单侧截割工况；而当二次旋转角β1=5°时，截割力FR均值和峰值变化幅度较小且载荷波动较小，此时，端盘截齿截割效率较高。

表5 不同切向安装角破碎性能统计量Table 5 Crushing performance statistics of different tangential installation angles

表6 不同轴向倾斜安装角破碎性能统计量Table 6 Crushing performance statistics of different axial tilt angles

图10 不同切向安装角时端盘截齿三向阻力曲线Fig.10 Curves of three-direction resistance with time under different tangential installation angles

图11 不同二次旋转角时端盘截齿三向阻力曲线Fig.11 Curves of three-direction resistance with time under different secondary rotating angles

表7 不同二次旋转角破碎性能统计量Table 7 Crushing performance statics of different secondary rotating angles

4 结论

1)借助VUMAT用户材料子程序接口，自定义HJC多金属硫化物本构模型代码，并将本构模型嵌入ABAQUS/Explicit 分析模块，通过将SMS 矿体样品三轴压缩实验和单个单元下的数值仿真进行对比，验证了VUMAT自定义HJC本构模型方法的可靠性。同时，与ANSYS/LS-DYNA 中内置的HJC本构模型数值仿真结果进行对比，说明了自定义HJC本构模型的优越性。

2)围压对于三向阻力有显著影响，随着围压增加，三向载荷增加，载荷波动减小。

3)侧向力波动是造成截齿疲劳断裂的主要原因，同时，轴向倾斜角和二次旋转角对侧向力的影响较大。当圆周切向安装角β0=45°，轴向倾斜角β1=10°，二次旋转角β2=5°时，截齿的截割效率较高，载荷波动较小。

4)端盘截齿安装参数对SMS 矿体的切削破碎有较大的影响。在进行海底多金属硫化物采掘头设计时，不能忽视圆周切向安装角、轴向倾斜角，二次旋转角这3个重要的因素。本文的仿真模型充分利用了数值分析工具的优势，在一定程度上定量分析了圆周切向安装角、轴向倾斜角和二次旋转角对截齿切削破碎阻力的影响，可为端盘截齿结构改良和采掘头参数优化提供参考。