陆 鑫 黄 勇 李 雯,2 赵 瑢 王 浚

(1 北京航空航天大学航空科学与工程学院,北京 100191)

(2 米兰理工大学机械工程系机器人技术实验室,意大利 米兰 20156)

1 引言

月尘是月球环境重要的组成部分,是月球探测中不可忽视的影响因素。月尘是松散的颗粒群,在微流星撞击、月面明暗交界静电效应、宇航员行走、探测车行驶、探测器着陆等一定外力作用下会产生扬尘。探测器在月面着陆过程中将与月尘颗粒发生碰撞,会导致月尘颗粒扬起悬浮。扬起和悬浮的月尘颗粒因其粘附特性会附着在与其接触的各类表面上,影响航天器和探测仪器的可靠性以及航天员的身体健康[1]。因此,研究月尘的扬起和漂浮机理及过程对于保障和提高月面探测器的可靠性具有重要意义。

国内外针对月尘做了较多的实验和仿真研究。根据公开的研究报道,内容涉及月尘颗粒的尺度分布[2]、微观形态和毒性作用[3]、粘附特性[4]、静电力作用下运动[5]、月尘防护技术和设备的研究[6-7]、月面探测车辆的行走通过特性研究[8]、月尘受外力作用漂浮的二维模拟等方面[9]。

相关的研究方法包括实验方法和数值模拟方法。由于月面环境的特殊性,且现有月尘珍贵稀少,一般使用特殊研制的模拟月尘来进行实验。数值模拟的主要方法有有限元法[10]、离散元法[11]等。离散元法是将所分析的物体看作离散颗粒集合体,根据离散物质本身的离散特性,对物质中的每个颗粒作为一个单元建立模型,进行模拟计算。在分析具有离散性质的物质方面离散元法具有较大的优势。

月球探测器在月面着陆过程中与月尘颗粒发生碰撞作用,导致月尘颗粒扬起。月尘颗粒的运动不仅受探测器的碰撞作用,同时受月尘颗粒之间的碰撞作用。因此月尘颗粒的漂浮受颗粒之间碰撞动力学的控制。本文采用离散元方法,建立月尘几何模型和接触力学模型。通过月尘的三轴压缩实验来确定月尘的细观参数。月尘颗粒的漂浮过程可分解为颗粒之间的碰撞过程和颗粒的运动过程。根据所得到的月尘参数建立了三维月尘离散元模拟的计算模型,数值模拟着陆作用下的月尘扬起和漂浮,确定漂浮特性和范围。

2 月尘三维离散元模型建立和参数确定

2.1 月尘离散元几何模型的建立

通过扫描电镜设备获得的模拟月尘颗粒形态的显微照片,分析得出模拟月尘颗粒具有复杂的颗粒形态,并以棱角形、次棱角形、长条形、次圆形颗粒形状为主[13]。

表1 模拟月尘的物理力学性质Table1 Physical and mechanical properties of lunar dust simulant

根据模拟月尘的颗粒形状特征,采用多球单元重叠法构建单个月尘颗粒的几何模型。四种颗粒单元的三维几何模型如图1所示,每个月尘颗粒是由一个基球O(x0,y0,z)和几个棱角球Bi(xi,yi,z)(i=1,2,3,…)重叠构成的。图中R0是基球半径,R1是颗粒单元外接球半径。

图1 四种颗粒单元的几何模型CAD 图Fig.1 Geometric models of the four particle unit(CAD)

月尘是覆盖在月球表面一层厚厚的由于长期环境作用而形成的微小颗粒,颗粒直径大部分小于1mm,平均直径为40～130μm[12]。月尘物理力学特性包括颗粒形态、孔隙率、内聚力、内摩擦角、容积密度、颗粒比重等。已知一种模拟月尘的基本物理力学性质[13]如表1所示。

2.2 离散单元接触力学模型的建立

月尘颗粒是由各种陨石和微陨石撞击、宇宙射线和太阳风持续不断轰击、月表大幅度温差变化导致月球岩石热胀冷缩等作用下形成的。因而月尘颗粒间的接触一般属于非粘性接触,月尘整体强度由摩擦强度确定,其表观粘聚力产生于颗粒间的摩擦和互锁作用。因此,将颗粒组成的单元的接触本构模型简化为弹簧-阻尼线性系统,并引入非张力连接模块和库伦摩擦模块[14],如图2所示。

两个单元间接触作用力Fc可分解为法向作用力Fn和切向作用力Fs

根据图2,法向接触作用力Fn由法向弹簧产生的弹性部分Fne和由法向阻尼器产生的非弹性部分Fnd组成。切向作用力Fs包括切向弹性力,切向阻尼力,库伦摩擦力。

图2 离散单元接触力学模型Fig.2 Discrete element contact mechanics model

当颗粒相互接触时,接触力可表示为

式中kn和ks是法向和切向刚度系数;un和us是颗粒在法向和切向上的位移;cn和cs是法向和切向阻尼系数;vn和vs是接触时速度在法向和切向上的分量,阻尼力的方向与它们的方向相反。

当接触单元发生滑动时,由于库伦摩擦力的存在,切向接触力不会超过库伦摩擦力。在模拟计算中,切向作用Fs满足下列关系

粘性阻尼力与运动方向相反,通常不直接设定阻尼系数,而由粘性阻尼的特征值临界阻尼比和临界阻尼常数共同确定,如下式

式中βn和βs分别是法向和切向上的临界阻尼比。和分别为法向和切向临界阻尼常数,这两个常数的计算公式如下

其中,ωn和ωs是无阻尼系统的固有频率,m是有效系统质量,由所接触颗粒的质量计算得到

2.3 模拟月尘的细观参数确定

离散单元模型细观参数一般无法通过物理实验直接测量,但是根据相似理论的量纲分析法中的相似第二定理(Π定理),选取月尘宏观尺度参数杨氏模量E,抗剪强度τmax,容重ρ和细观尺度参数接触刚度kn和ks,临界阻尼比cn和cs,无量纲物理量摩擦系数μ,利用FLT 量纲系统进行分析计算,得到Π关系式

可见月尘的细观参数与其宏观力学参数存在尺度定律。可以利用模拟月尘的三轴压缩实验结果,采用三维离散元三轴压缩实验数值模拟的方法,建立月尘模型细观参数与宏观力学属性的内在联系,并通过调整输入参数来确定符合月尘实际宏观力学特性的最适合离散元模型参数。

2.3.1 三轴压缩实验的离散元仿真分析

根据模拟月尘三轴压缩试样的实际实验条件,建立与实际尺寸相同的试样模型,采用柔性边界墙维持围压和刚性墙加载方式,对试样的压缩过程进行模拟。试样直径为39mm ,高度H为80mm,三轴压缩实验三维离散元仿真模型,如图3所示,左边是模型外观图,右边是截面图。

图3 三维离散元三轴压缩实验仿真模型Fig.3 DEM t riaxial test model

图中VZ为底部加载墙的加载速度,柔性边界墙由半径rb为1mm 的圆球组成的圆环,在80mm高度上分了79层。考虑到压缩时的变形,组成柔性墙的圆球数中间层中最多,并依层向两头递减。这些圆球相邻重叠并采用平行粘结的方式相连。为控制边界条件获得σ2大小的围压,在组成柔性边界墙的每个圆球上预加力Fb,可由下式近似得到

式中Rb是圆球中心距试样中心轴的距离,Ni是第i层圆球数。Fb方向从圆球心水平指向Z 轴。

采用底部刚性墙匀速加载和顶部刚性墙应力采集的方法实现三轴压缩实验的离散元仿真,试样受到的主应力差(σ1-σ2)和轴向应变ε可表示为

式中σ1为Z 轴方向主应力,Fz是顶部刚性墙受到的Z 方向的不平衡力,R是当前试样圆柱的平均半径,Swd是底部刚性墙的垂直位移。

2.3.2 细观模型参数的确定

根据上述仿真方法,进行了模拟实验。三轴压缩实验仿真试样的颗粒外接球半径为1～2.5mm,密度为2 770 kg·m-3,初始孔隙率约为0.35,法向切向临界阻尼系数为0.7。需要确定的细观参数是颗粒刚度和摩擦系数。下面分别确定其中一个参数,调整另外一个参数进行试验。

图4是围压为26kPa,摩擦系数为0.3,颗粒接触刚度不同时的应力—应变关系曲线。可以看出颗粒单元的刚度越大,则主应力差也越大,峰值越高。随着刚度的增加,颗粒与压缩墙体的接触变得不稳定,接触力波动逐渐增大,更大的刚度则会致使主应力数值更不稳定。

图4 不同刚度下的应力-应变曲线Fig.4 Axial strain vs.stress curves in different stiffness

图5 不同摩擦系数下的应力-应变曲线Fig.5 Axial strain vs.stress curves in different friction coefficient

图5是围压为26kPa,法向、切向刚度均为1×105N/m,颗粒间不同摩擦系数下的应力—应变关系曲线。摩擦系数在0.4 以下时,曲线平缓上升,主应力差随着轴向应变增加而提高,呈现应变硬化。摩擦系数在0.4 以上时,出现了不是非常明显的峰值强度,即主应力差随着轴向应变增加而减小,呈现应变软化,此时颗粒群处于失稳或破坏状态。可见颗粒摩擦系数越大,则颗粒群在压缩过程中克服颗粒间的翻转、滑移所需要的能量越多,在宏观上体现为主应力差更大。

图6是用模拟月尘进行的实际三轴压缩实验结果[15]。从图中可以看出随着轴向应变ε的增大,主应力差(σ1-σ2)也逐渐增大,曲线呈上升趋势。在应变为4%～5%之间时达到峰值,约为100kPa。之后主应力差随轴向应变变化不大,曲线趋于平稳,主应力差值在100kPa 左右。

图6 实际三轴实验应力-应变曲线Fig.6 Actual axial strain vs.stress curve

根据上述实际实验结果分析,将仿真模拟得到的结果与之相对比,发现当法向和切向刚度在1×105～5×105N/m、摩擦系数在0.2～0.3 范围内时曲线趋势和峰值与实际结果最为相近。将该范围参数三轴模拟结果与实际结果详细对比,在应变为4%、主应力差曲线达到峰值时,模拟结果的主应力差值为92～100,而实际实验主应力差值为97,可见误差较小,该范围的参数可以作为仿真试样模型的参数。

3 着陆作用下月尘漂浮三维模拟

3.1 模拟方法和条件

假设模拟月面着陆器为一圆盘模拟体,并在距离月面一定高度处自由下落。一定数量的月尘颗粒在重力作用下沉降在模拟槽内。圆盘模拟体着陆时与月尘颗粒发生碰撞,使得月尘颗粒的扬起和漂浮。月尘漂浮运动过程满足颗粒动力学,圆盘模拟体与月尘颗粒以及月尘颗粒间的相互碰撞满足颗粒碰撞力学。

对圆盘模拟体和月尘颗粒运动情况进行跟踪,圆盘模拟体和月尘颗粒受微重力作用力和颗粒与圆盘模拟体以及颗粒间碰撞产生的作用力,根据牛顿第二定律,颗粒运动的控制方程为

式中mi、Ⅰi、ui和ωi分别为第i颗粒的质量、惯性距、速度和角速度,Fc为所有与第i颗粒接触的颗粒对其作用力的合力,Mc为作用在颗粒上的转矩。

圆盘模拟体和月尘颗粒的接触力学模型同月尘颗粒间的接触力学模型,见式(1)、(2)、(3)。根据式(1)和(15),可以确定不同时间颗粒的运动速度和空间位置,统计获得宏观颗粒运动特性。

实际月面着陆器和月尘结构非常复杂,为了简化计算,将月面着陆器简化成由模拟圆球构成的圆盘模拟体,其质量接近实际着陆器的质量,初始位置定位于月尘表面,初速度按照物体在一定高度自由落体计算。月尘颗粒是按照四种颗粒几何模型等比组成。月尘分为两层,上层颗粒细小,受激漂浮。下层颗粒粒径为上层颗粒三倍,传递颗粒间的接触力。两层颗粒均在月面重力作用下沉降至平衡,月尘颗粒初始速度为零。模拟计算并统计分析月尘运动特性及月尘空间分布。

3.2 圆盘模拟体着陆和月尘颗粒漂浮过程分析

3.2.1 圆盘模拟体着陆作用下月尘漂浮分析

模拟计算中取圆盘模拟体由一组直径为0.02m模拟圆球组合构成,其直径为0.15m,密度为25 000kg/m3。模拟圆球的刚度为1×106N/m,摩擦系数为0.7。月尘颗粒刚度系数为1×105N/m,摩擦系数为0.3,密度为2 770kg/m3。圆盘模拟体和月尘颗粒在微重力条件下碰撞运动,月面重力加速度取为1.68m·s-2。

利用颗粒流程序进行模拟计算时,循环使用的时间步长取临界时间步长的一部分,取用因子为0.8。在计算循环中颗粒间接触不断发生变化,根据上式计算得到的每个循环的时间步长也随之变化,基本上每个循环均不一样,这样就使得每个循环对应的实际时间不相同。在进行数据处理时,根据变化的时间步长来确定实际时间比较麻烦,且不能较为准确地估算得到实际时间,因而这里采用循环步来代替实际时间进行处理。这样处理直接方便,对最后的结果分析也影响不大。

图8是月尘颗粒在圆盘模拟体作用下漂浮过程截面图。初始时圆盘模拟体位于月尘表面,给定初速度V为2m/s。开始计算后,随着时间的推进,圆盘模拟体在一定初速度和微重力的作用下,与月尘颗粒接触发生碰撞作用。一部分颗粒被挤压进月尘内部,随着挤压程度增加,月尘内部颗粒密集度增大,抵抗圆盘下落的能力提高,圆盘下落速度减缓。一部分颗粒被扬起,这部分颗粒在圆盘与月尘颗粒发生非弹性碰撞作用时获得能量,产生向上运动而漂浮。随着时间推进,漂浮颗粒数增多,颗粒漂浮高度增大。漂浮颗粒在微重力作用下,上升至最高处后开始逐渐回落。

统计分析漂浮月尘颗粒的平均Z 轴向速度Vpaz如图9所示。初始月尘颗粒静止,受圆盘着陆碰撞作用,被激扬漂浮,漂浮速度增大,在较短的时间内平均速度达到最大。受微重力作用的影响,月尘颗粒上浮速度不断减小。随着计算循环的增加,计算循环到100 000 步时,漂浮颗粒平均速度值为0。此时大部分漂浮颗粒停止上升,漂浮在空中。继续计算,平均速度变成负值,表明漂浮颗粒大部分开始下落,进入沉降过程,而且下落速度逐渐增大。可见,月尘颗粒的Z 轴向上升速度在受激漂浮后随着时间的推移逐渐减小。

图8 圆盘着陆和月尘漂浮过程Fig.8 Progress of simulated disc body landing and luanr dust levitation

图9 月尘颗粒平均Z 轴向速度随时间变化Fig.9 Average Z-axial velocity ofparticles vs.cycle step

图10 月尘颗粒数沿高度分布Fig.10 Number of particles along the height

图10是两个计算步时沿高度h 的月尘颗粒数变化,图中N是指漂浮的月尘颗粒数。由图可见沿高度方向漂浮的月尘颗粒数逐渐减少,表明月尘颗粒数在高度上按指数规律递减。随着计算循环的继续,漂浮的月尘颗粒数增多且漂浮高度增加。根据图中漂浮月尘沿高度上分布情况可以发现,在该计算条件下,月尘漂浮最高可达到70mm。

3.2.2 统计分析着陆速度变化对月尘漂浮的影响

统计在不同着陆速度的圆盘碰撞作用下月尘颗粒漂浮情况,分析圆盘着陆速度对月尘漂浮的影响。图11是不同着陆速度的圆盘作用下月尘飘浮的最大高度hmax随时间变化情况。从图中可以看出,圆盘着陆速度越大,碰撞后月尘漂浮最大高度越大,同时颗粒回落下来的计算时间更长,这表明颗粒漂浮在空中时间变长。且当速度增大时,漂浮高度从近60mm 扩大到近160mm。

图11 月尘颗粒漂浮最大高度随时间变化Fig.11 Maximum height of levitation particles changes with time

根据上述分析可以看出,圆盘着陆速度对月尘漂浮数量、漂浮时间和漂浮高度均有影响。且圆盘着陆速度越大,漂浮的月尘颗粒数越多,漂浮时间越长,漂浮高度越高。在以上三种着陆速度条件下的计算结果表明,漂浮最大高度能够达到近160mm。

着陆速度的变化除了沿竖直方向上的量值变化,还存在偏离竖直方向从而带有切向速度分量的情况。这里简单研究了速度为2m/s,水平方向上有速度分量,即速度的方向沿X 轴偏离竖直方向15°和30°时这两种条件下对月尘受激漂浮的影响。图12为速度是2m/s 偏离角度α是15°和30°时,月尘颗粒漂浮数沿高度分布情况。

图12 中明显看出偏离角度较小时,竖直方向上的速度分量较大,月尘漂浮的颗粒数也相应较大,大于飘离角度大时(即竖直方向上的速度分量较小时)的数量,这显然与之前获得的竖直方向上速度量的变化引起的月尘漂浮颗粒数变化相一致。而两种情况下的月尘漂浮最大高度相差不大,这是因为在模拟计算过程中,月尘颗粒受激漂浮时在飞扬过程中颗粒间也会相互碰撞,使得原本漂浮高度较小的获得飘向更大高度的能量,使得最大漂浮高度上升,因而最大漂浮高度相差不大。

3.2.3 统计分析月尘细观参数变化对月尘漂浮的影响

图12 月尘漂浮颗粒数沿高度分布Fig.12 Number of particles along the height

统计月尘颗粒其细观参数中颗粒刚度变化时月尘漂浮的情况。在上述确定的最合适月尘刚度范围内,即1×105～5×105N/m,取值进行计算模拟。分析月尘颗粒刚度分别为1×105N/m、2.5×105N/m和5×105N/m时,圆盘以2m/s 的速度着陆,月尘受激扬的漂浮情况。图13是月尘颗粒不同刚度时漂浮最大高度随时间变化情况。

图13 月尘颗粒漂浮最大高度随时间变化Fig.13 Maximum height of levitation particles changes with time

从图13可以发现,月尘颗粒细观参数中刚度变化时,月尘颗粒漂浮受到影响。刚度越大,颗粒间相同变形所需的接触力增大,月尘颗粒漂浮上扬需要的能量增大。因而相同条件下月尘颗粒刚度大则扬起的最大高度小,刚度小反而扬起的最大高度大。同时对月尘颗粒漂浮的时间也有影响,刚度较大时,漂浮的时间较短;刚度较小时,漂浮的时间较长。统计分析三种刚度情况下月尘颗粒漂浮达到最大高度时的漂浮颗粒数,发现漂浮颗粒数量的变化不大。这表明刚度变化对月尘漂浮的最大高度和漂浮时间有一定影响,而对漂浮颗粒数影响不大。

4 结论

着陆器在月球表面进行着陆操作时,会引起月面月尘颗粒的激扬行为,并且这种激扬特性与月尘颗粒的宏细观物理力学属性、着陆器的结构特征及其工作参数等因素密切相关。

采用一种三维离散元计算机仿真方法来模拟着陆过程中的月尘激扬问题。通过三维三轴压缩实验仿真的方法,对比模拟结果和实际结果得到最适合的月尘模拟参数。利用圆盘模拟体-月尘的离散元三维分析模型,进行了圆盘着陆碰撞作用下月尘漂浮过程模拟。

对圆盘模拟体着陆时引起的月尘浮扬情况进行了统计分析,得出漂浮的月尘先随时间增加而增多,然后在微重力作用下逐渐回落。漂浮的月尘颗粒数沿高度按指数规律递减。在给定其他计算条件下,圆盘初速度为2m/s时,月尘漂浮高度范围为70mm。统计分析了着陆速度和月尘颗粒刚度的变化对月尘漂浮的影响。圆盘着陆速度变化对月尘漂浮的数量、时间和高度的影响较大,且随着陆速度的增大,月尘漂浮的数量、时间和高度的值也增大。当着陆速度在3m/s时,月尘漂浮高度范围可达160mm。当存在水平速度分量时,水平速度分量增大,扬起的月尘颗粒数明显减少。月尘颗粒刚度变化对月尘颗粒漂浮的高度和时间有影响。刚度越大时,月尘颗粒漂浮的最大高度和漂浮时间越小,但是刚度变化对月尘漂浮颗粒数影响不大,三种刚度下漂浮的月尘颗粒数相近。

References)

[1]Gaier J R.The effects of lunar dust on EVAsystem s during the Apollo missions [R].NASA/TM-2005-213610,2005

[2]David W C III.Lunar soil grain size dist ribution[J].Planetary and Earth Sciences Division,1972,10

[3]Park J S,Liu Y,Kihm K D,et al.Micro-morphology and toxicological effects of lunar dust[J].Lunar and Planetary Science,2006,XXXVII

[4]Walton O R.Adhesion of lunar dust[R],NASA/CR-2007-214685,2007

[5]Mazumder M K,Horenstein M N,Trigw ell S,et al.Electrostatic and gravitational transport of lunar dust in the airless atmosphere of the moon[J].IEEE 978-1-4244-2279-1/08,2008

[6]Calle C I,Buhler C R,McFall J L,et al.Particle removal by electrostatic and dielectrophoretic forces for dust control during lunar exploration missions [J].Journal of electrostatics,2009,67:89-92

[7]Clark P E,Curtis S A,Minetto F A,et al.Characterizing physical and electrostatic properties of lunar dust as a basis for developing dust removal tools[C]// NLSI Lunar Science Conference,2008

[8]邹猛,李建桥,张金换,等.月球车驱动轮牵引性能研究[J].宇航学报,2009,1:98-103

[9]王淑彦,何玉荣,刘国栋,等.二维拱形模拟体作用下月尘颗粒悬浮过程研究[C]//大庆:2007 多相流学术会议,中国工程热物理学会,2007

[10]David W C III.Lunar soil mechanics:distribution of contact stress beneath a rigid plate resting on sand[D].Massachusetss Institute of Technology,1968

[11]Koji Wada,H iroki Senshu,Takafumi Matsui.Numerical simulation of impact cratering on granular material[J].Icarus 2006,180:528-545

[12]郑永春,欧阳自远.月壤的物理和机械性质[J].矿物岩石,2004,24(4):14-19

[13]Li Wen,Gao Feng,Jia Yang,Tractive performance analysis on radially deployable wheel configuration of lunar rover vehicle by discrete element method[J].Chinese Journal of Mechanical Engineering 2008,21(4):66-72

[14]Li Wen,Huang Yong,Cui Yi,et al.Trafficability analysis of lunar mare terrain by means of the discrete element method for w heeled rover locomotion[J].Journal of Terramechanics,2010,47(3):161-172

[15]高峰,李雯,孙刚,等.模拟月壤可行驶性的离散元数值分析[J].北京航空航天大学学报,2009,35(4):501-504