基于EDEM的月球极区采冰钻具参数的研究

汤厚鑫，纪玉杰

（沈阳理工大学 机械工程学院，沈阳 110159）

0 引言

水资源是各个国家在月球进行探索的最重要物质及能源原料。1961年，由美国科学家Watson等第一次提出在月球南极和北极撞击坑的永久阴影区中存在大量水冰资源[1]，他们认为在太阳赤道平面与月球赤道平面夹角小于1.6°的情况下，位于月球南极和北极的某些撞击坑底部可能永久地处于温度为40～50 K的阴影之下，在这样的低温下，位于月球撞击坑永久阴影区水资源会以冰的形式存在。自Watson等首次提出在月球两极永久阴影区存在水冰资源以来，各个国家通过遥感探测、探地雷达等多种手段对月球极区水冰资源开展了探测活动[2-3]。根据探测结果，科学家们推测在月球南极的某些区域会有水冰的存在。

本文以冰为离散元模型，选取钻头的回转转速、进给速度等重点钻具参数进行对比分析，确定各钻具参数对钻具的受力及对冰的采样量的影响趋势，为工程设计提供一定的理论依据。

1 冰离散元模型的建立

本文为了在离散元软件EDEM中模拟冰的动力作用，考虑将冰离散为具有一定质量和大小的球形颗粒单元，并且单元间具有相应的黏结作用。因此本文考虑建立弹性接触力模型及平行黏结模型来分别描述冰颗粒单元间的动力作用和黏结作用。

1.1 冰颗粒离散单元间的弹性接触力模型

冰颗粒的离散元接触力模型表达了每个颗粒单元间力与位移的主要关系；在颗粒间相互作用的过程中，颗粒单元间会产生弹性作用力和剪切力[4]；其中弹性作用力是由颗粒单元间的相对速度及弹性变形导致的，然而在计算剪切力时需要满足摩尔-库伦（Mohr-coulomb）摩擦定律；弹性作用力又可以分为线性和非线性两种类型；在理论上可以采用非线性力学模型（Hertz-Mindlin）来表示颗粒碰撞的过程[4]；但是在分析颗粒单元的动力学特性时，线性弹性接触力模型在应用时更加方便和快捷，并且可较为准确地计算出颗粒单元的动力学行为[4]。如图1所示，建立冰颗粒单元间的线性弹性接触力模型，其中：颗粒A与颗粒B的质量分别为MA和MB；法向和切向刚度系数分别为Kn和Ks；法向和切向阻尼系数分别为Cn和Cs；μ为摩擦因数。

图1 冰颗粒单元间的线性弹性接触力模型

颗粒在法线方向上的相互作用力是由摩尔-库伦（Mohr-coulomb）摩擦定律来决定，其中法向力Fn是由弹性力Fe和黏性力Fv两部分组成[4-5]，因此法向力Fn的计算公式为

式中：Fe=KnXn，Fv=-CnX˙n。因此法向力分别为颗粒的法向变形和变形率）。

颗粒在切线方向上，基于Hertz-Mindlin理论和摩尔-库伦（Mohr-coulomb）的摩擦定律[6]可知，颗粒间的切向作用力计算公式为：

式中，Xs和X˙s分别为切向变形和变形率。

两个颗粒单元A和B碰撞的法向有效刚度系数计算公式[7]为

而法向阻尼系数按下式计算：

式中:Cn为法向阻尼系数，ζn为无量纲法向阻尼系数，e为回弹系数，M为颗粒单元A或B的有效质量。法向和切向刚度阻尼系数有如下对应关系[7]：Ks=αKn，Cs=βCn。其中：α=0.5，β=0。

在线性接触模型当中，颗粒单元A和颗粒单元B从碰撞到分离的接触时间[6]可被定义为

在计算过程中，计算步长为接触时间的1/50；在线弹性的模型当中，接触时间是一个与颗粒大小和材料相关的参数。

1.2 冰颗粒离散单元间的平行黏结模型

由于冰颗粒单元间具有冻结作用，所以考虑建立冰颗粒单元间的黏结模型来表达冰颗粒间的冻结作用；黏结模型有2种类型，即接触黏结模型和平行黏结模型[6]；颗粒单元间的接触黏结模型只传递力的作用，且接触黏结的黏结作用只发生在接触点上；然而平行黏结模型不仅可以将2个冰颗粒单元胶结在一起，而且还可以传递力和力矩的作用；在本文中主要考虑建立冰颗粒间的平行黏结模型。

平行黏结模型可以被想象为一对法向和切向刚度恒定的弹簧设置在两颗粒接触平面的圆盘上[6]，并且两颗粒间的圆盘在接触平面的正中央；图2所示为两冰颗粒单元间的平行黏结模型；两颗粒间的力和力矩可转化为作用在黏结材料边界上的最大法向和切向应力；其中任意一个最大的应力超过了对应的连接强度，冰颗粒间的平行黏结模型就会发生断裂[6]；其中：XA和XB分别为冰颗粒单元A和冰颗粒单元B的位置矢量，n为接触方向的矢量，在两颗粒接触点的圆盘上可以传递冰颗粒单元间的作用力及力矩，即颗粒单元间的拉力、剪力、弯矩和转矩；与平行黏结有关的总的力和力矩可以用Fb和Mb表示；每一个向量都可以分解为法向分量和切向分量：

图2 两冰颗粒单元间的平行黏结模型

式中：Fb和Mb分别为力和力矩的矢量，分别为力和力矩的法向和切向分量。

在平行黏结模型当中，根据梁的拉伸、扭转和弯曲理论，黏结圆盘上最大拉应力和最大剪应力的计算公式为：

式中，A、J、I分别为两颗粒间圆盘的面积、极惯性矩和惯性矩。通常有A=πR2，J=πR4/4，I=πR4/4，其中R为两颗粒间圆盘的半径。当两颗粒圆盘间的最大拉应力及剪应力超过其拉伸强度及剪切强度时，两颗粒之间的平行黏结模型就会发生断裂。

2 基于EDEM钻具钻取冰块的仿真试验设计

为了研究钻具的进给速度、转动速度对钻具的受力和冰采样量的影响，需要制定出不同的仿真工况进行对比分析。

2.1 钻具设计

如图3所示，本文主要设计一种钻具，该钻具的内径为40 mm，外径为60 mm，钻具的高度为275 mm，螺旋叶片的螺距为35 mm。

图3 内径为40 mm钻具的结构

2.2 钻具运动参数的设计

衡量钻具的回转钻进运动的主要运动参数有回转转速和进给速度；主要根据所查文献设计，对10、15、20 rad/s等3种转动速度进行对比分析[8]。对于进给速度，本文主要设计出一种进给速度为0.15 m/s进行对比分析。因此本文主要针对转动速度设计3组仿真工况（如表1），对钻具的受力及冰的采样量进行对比分析。

表1 仿真工况设计

2.3 冰的离散元参数以及模型的建立

在离散元仿真软件EDEM中对颗粒模型设置需要进行两个方面的处理：1）定义冰颗粒单元的几何模型；2）定义冰颗粒单元的计算参数及在前处理中生成冰颗粒单元。其生成的冰颗粒离散元模型如图4所示，主要的计算参数如表2所示[7,9]。

表2 冰颗粒离散元模型计算参数

图4 冰颗粒离散元模型

2.4 钻具工作过程的仿真

钻具与冰的相互作用过程中，钻具以一定的回转速度和进给速度运动，而冰以剪切破碎为主。本文主要截取T=0.9 s时钻具的位置，如图5所示。

图5 T=0.9 s时钻具的位置

3 仿真计算结果对比分析

取表1所示的第1、2、3组仿真数据进行对比分析；3组数据分别对应钻具的3种不同的转动速度；对于钻具的受力状况，分别记录钻具所受的力矩。对于冰颗粒的采样量，在离散元软件EDEM后处理中测量出钻具内部冰颗粒的高度。再从离散元软件的后处理中记录数据进行对比分析，研究钻具的转动速度对钻具的受力及冰采样量的影响。不同转动速度的钻具对冰的采样量如图6所示，不同转动速度的钻具力矩如图7所示。

图6 不同转动速度的钻具对冰的采样量

图7 不同转动速度的钻具力矩图

从冰颗粒的采样量分析，对于不同转动速度的钻具，钻具的转动速度分别为10、15、20 rad/s，可测量出钻杆内部冰样品的高度分别为79.05、86.01、101.11 mm。由此可以看出，随着转动速度的增大，钻杆内部冰样品的高度逐渐增加。

从钻具的受力情况分析，对于不同转动速度的钻具，钻具的转动速度分别为10、15、20 rad/s。如图7所示，随着转动速度提高，钻具所受的力矩呈现逐渐减小的趋势。

4 结语

本文首先建立了冰颗粒离散元间弹性接触力模型及冰颗粒间的平行黏结模型。通过查阅文献得到冰颗粒离散元模型的计算参数，采用此参数在离散元仿真软件EDEM中生成了冰颗粒的离散元模型，利用此模型完成了3种工况下月球极区冰钻取过程的离散元仿真。通过记录3种工况下钻具的受力情况及冰的采样量，分析钻具的转动速度对回转钻进的影响，并得到以下结论：随着钻具转动速度增大，冰的采样量逐渐增大，钻具所受的力矩逐渐减小。