月壤双轴试验的剪切带离散元数值分析

蒋明镜 ，郑 敏，王 闯

（1.同济大学 地下建筑与工程系，上海 200092；2.同济大学 岩土及地下工程教育部重点实验室，上海 200092； 3.云南省交通规划设计研究院，昆明 650011；4.中国空间技术研究院，北京 100094）

1 引 言

月球的表面覆盖着一层松软的散粒体状的月球土壤。月壤是在无水、陨石撞击、宇宙射线、太阳风等环境下形成，形成过程中物理风化起主导作用，且月壤所处的环境与地球环境有着很大差别（如低重力场、低气压等）。图1 为宇航员采集月壤样品的图片，月壤的颗粒大小、形态与地球上的粉质砂土类似[1－3]。

图1 月壤形态图 Fig.1 Lunar soil

随着“嫦娥二号”月球探测卫星的成功发射，我国探月计划的下一个目标就是实现“嫦娥三号”探测器在月面的软着陆。月球探测器包括着陆器和巡视器两部分，分别实现月面软着陆和巡视勘察。这两者都是人工装置与月壤相互作用的问题，都牵涉到月壤的力学特性以及月壤的破坏问题。

马克思主义认为，人是一切事物的根本。对于高校来说，思想政治工作说到底是育人工作，其出发点和落脚点都是为了育人，育人是其价值本位。毛泽东曾指出：“世间一切事物中，人是第一个可宝贵的。”[4]“课程思政”应始终围绕立德树人这一教育的根本任务，坚持大局意识，深挖各专业各学科的育人价值，形成课程整体育人的联动效应，促进学生的全面成长成才；要聚焦育人价值的本源，明确价值指向，找准目标定位，弄清楚教育教学工作是为什么而做，可以做什么，最终是为了实现什么。

对月壤力学性质的研究手段主要有：地面与绕月探测器的遥感[4－5]、登月探测器及宇航员的现场测试以及对带回地球月壤的室内试验[6－8]。根据美国Apollo 登月计划和前苏联Luna 计划中所实施的大量现场试验以及对带回地球的真实月壤样品所进行的室内试验等研究可知，真实月壤的内摩擦角φ范围在25°～50°之间，黏聚力c 的范围在0.1～ 3.8 kPa 之间。由于真实月壤极其珍贵稀少，试验的方法有其局限性，近年来，一些学者开始采用数值模拟的方法来模拟月壤的力学特性以及月球车车轮等人工装置与月壤的相互作用等问题[9－11]，主要有有限单元法与离散单元法。然而，关于月壤的破坏问题的研究鲜见报道，在登月活动中的着陆器着陆与月球车在月面上行走都牵涉到月壤的强度与破坏问题。土体的破坏问题一般都是从其剪切带入手进行研究，对剪切带问题的研究可以很好地反映材料的破坏模式与机制。此外，目前能在地球上对月壤破坏特性进行研究，这些特性在月球环境下会有什么不同，是探月工程较为关注的问题之一。本文针对真实月壤所处的环境（无水、低重力场、低气压等）和内摩擦角会很大的特点，介绍1 种考虑颗粒的抗转动作用与颗粒间的范德华力2 个因素的月壤散粒体力学接触模型，并将该模型植入离散元分析商业软件PFC2D[12]中，模拟了柔性边界条件下的双轴压缩试验，着重对月壤的剪切带问题进行了研究，并对比分析了月面环境对月壤力学特性的影响。

2 月壤散粒体力学接触模型简介

文中的数值分析采用的月壤散粒体力学接触模型是在PFC2D 商业离散元软件自带线弹性颗粒接触模型的基础上建立的，主要考虑了以下2 个因素：

（1）颗粒间抗转动作用：真实月壤的颗粒形状以棱角形为主，颗粒与颗粒间接触必然存在面接触，因此，接触处可以传递力矩，颗粒的转动必然受到相互的限制，内摩擦角较高。然而PFC2D 自带模型并不能考虑颗粒间的这种抗转动作用，用PFC 自带模型计算的试样内摩擦角普遍较低[13]，考虑颗粒间抗转动作用能够有效提高离散元数值模拟试样的内摩擦角[14]。

在三门江林场中,为了使激励发挥其本质作用,真正的实现奖惩分明、奖勤罚懒,就必须制定合理的绩效考核制度,这也是所有企业对员工工作成果评价的重要一环。对员工的工作进行绩效考评,主要体现在两个方面,一个是对工作"量"的考评,一个是对工作"质"的考评。在绩效考核中,往往是综合这两方面来进行,若人力资源管理只关注某一方面,工作就会过于片面。绩效考核,需要对员工工作的优缺点进行客观系统的评价,通过科学合理的考核制度和考核办法,将每个员工的工作考评进行量化,得到一个最终成绩,依据此成绩,对员工的工作进行奖励和惩罚决定。

（2）颗粒间范德华力：真实月壤在月球环境下表现出一定的黏聚力，其产生原因很可能是在月球近于真空环境下月壤颗粒间存在相互吸引的范德华力[15－17]。为了模拟这一现象，在颗粒接触力学关系中引入相互吸引的范德华力。

为完整起见，首先简单介绍柱面坐标系。设M（x,y,z）是空间中一点，过点M作直线和坐标面xOy垂直相交于点 M'，称点 M'为点 M在坐标面 xOy上的投影；设点M'的 极 坐 标 是 （r,θ），称 有 序 三 元 数 组 （r,θ,z） 为 点 M的柱面坐标(如图 1所示)。

（1）颗粒间抗转动作用

接触模型主要由法向接触、切向接触以及颗粒的转动接触组成。常规的接触模型只考虑颗粒的法向与切向部分，不考虑颗粒的抗转动作用[12]。下面重点对颗粒间的抗转动作用与范德华力的计算进行简要说明。

离散元软件PFC2D 中自带接触模型只考虑圆形颗粒的点接触，对颗粒间的相对转动没有约束作用，颗粒间可以自由转动。显然，这与土颗粒的真实转动情况不符。月壤散粒体力学模型中的颗粒转动接触模型由蒋明镜等[14]提出，如图2 所示，模型中假定颗粒通过一定的宽度B 接触，颗粒接触宽度B = βR，其中β 为抗转动系数，R 为两颗粒的平均半径， R =2 R1R2/( R1+ R2)， R1、 R2分别为两球的半径。模型中用β 来反应颗粒间的抗转动情况，抗转动能力的大小通过调整该参数控制。图3 为转动模型中颗粒的扭矩与颗粒间的相对转动角度的关系曲线。

首先，进一步发挥好哈萨克小说在思考人生、珍视大自然等方面独特的人文教育功能，并在小说中突出体现人性境界提升、理想人格塑造以及个人与社会价值实现的人文教育理念。如艾克拜尔·米吉提的小说《车祸》，通过对区域人物与事件的叙述描写，从细微的心理刻画上，在美心美行的浸润中探究人与人之间的美好人性。

图2 颗粒接触宽度 Fig.2 Particle contact width

图3 颗粒转动接触模型 Fig.3 Particle rotation contact model

颗粒间的扭矩通过增量法计算，与抗转动系数、颗粒接触法向力、颗粒平均半径等有关，其计算公式见式（1）～（4）。在达到临界相对转角后，扭矩维持在峰值（模型1）或为0（模型2）不再变化，如图3 所示。本文中采用模型1，具体计算公式如下：

第一，行政隶属关系导致政策不统一。传统的行政划分形成了“区位差距”和“政策壁垒”，三地的教育主管部门制定政策，高校师资培训中心按照政策执行，属地化的政策差异导致一些工作不能协同，缺乏整体统一的规划。例如，北京市教委出台政策，拥有博士学位或者具有副教授和教授职称的高校教师，可以免修岗前培训课程，直接认定教师资格证。天津仍沿袭传统的政策，认定高校教师资格需参加岗前培训。而河北省教育厅刚调整了认定高校教师资格考试的政策，规定了岗前培训的成绩等同于教师资格认定考试成绩。三地间没有形成一个统一的认定政策。

其中，

MESH网络即“无线网格网络”,它是一个无线多跳网络，是由ad hoc网络发展而来。在向下一代网络演进的过程中，无线是一个不可或缺的技术。无线MESH可以与其它网络协同通信，是一个动态的可以不断扩展的网络架构,任意的两个设备均可以保持无线互联。MESH无线网络是一项极有前途的技术，被誉为下一代无线因特网，使普通无线技术过去一直存在的可扩充能力低和传输可靠性差等问题迎刃而解。由于网络中大量终端设备能够自动通过无线连成网状结构，网络中的每个节点都具备自动路由功能，每个节点只和邻近节点进行通信，因此是一种自组织、自管理的智能网络，不需主干网即可构筑富有弹性的网络[3-4]。

式中：M 为粒间扭矩；mk 为扭转刚度；β 为由抗转系数；nk 为法向接触刚度；θ 为颗粒间相对转角；rθ 为临界相对转动角，按式（4）计算。nF 为颗粒间的法向接触力；pM 为峰值扭矩值。

（2）颗粒间范德华力

在月壤所在环境中，月壤颗粒被一定厚度的气体分子层所包围[15]，如图4 所示，图中圆圈代表被月壤颗粒吸附的气体分子，t 为所吸附气体分子层厚度。为简单起见，用2 个圆球代表月壤颗粒，如图5 所示，它们的接触面假定为一圆面，半径为a，颗粒接触面间被两颗粒吸附气体分子层隔开，距离为D。

式中：A 为Hamaker 系数，对于月壤，A 约为4.3× 10-20J 左右[15]；D 为两颗粒间所吸附分子层的厚度，该厚度由月壤环境中的气压和周围气体分子种类决定，可用Adamson[16]在1990 年提出的吸附势理论计算，在本文离散元计算中将D 作为一个输入参数。有研究表明，范德华力在距离D 小于0.5× 10-6m左右作用明显[17]，在月球几乎真空的环境中，月壤颗粒间的距离D 很小，颗粒间范德华力作用明显。而对于地球环境下的砂土颗粒，由于D 值通常非常大，其分子间作用力相对于重力来说十分微小，对砂土的宏观力学特性的影响很小，一般不需要考虑。颗粒间接触面的半径a 通过下式计算：

图4 月壤颗粒吸附气体分子示意图[15] Fig.4 Adsorbed gas thickness[15]

图5 月壤颗粒接触示意图 Fig.5 Lunar soil particle contact

两颗粒受到总的相互的范德华引力可以分为两部分：一部分来自于颗粒间两接触平面，另一部分来自于颗粒剩余部分的相互作用。范德华引力的方向为接触的两球球心的连线方向，其计算式为

《混凝土结构设计规范》(GB 50010-2010)提出，临界截面周长，取距离局部荷载或集中反力作用面积周边h0/2处板垂直截面的最不利周长，h0为板的有效厚度。临界截面周长的计算简图如图1所示。

接触面的直径2a 与两颗粒抗转动作用的接触宽度B 一致。

3.平台学习功能架构。一个完整的干部网络学习平台需包括管理模块、学习模块、资源模块、信息模块以及学习服务模块。现有各平台更加着重于学习和管理模块的建设，各平台同时还根据自己的特色设置不同的子栏目，满足学员的学习需求，但信息模块中学员之间的交流实效性较差，学习服务模块还有待完善。

3 双轴试验成样及研究方法介绍

成样方法按蒋明镜等[18]提出的分层欠压法进行，该方法能够避免成样过程中试样下密上松的现象。成样过程中分8 层成样，颗粒总数目为24 000。为了使每层试样密度均匀，新一层试样的最终孔隙比较前几层的平均孔隙比要小，使得最终生成的8层试样的平均孔隙比达到目标取值（文中取为0.20），本文成样过程中从下到上采用的各层孔隙比为 ep(1)=0.215、 ep(2)=0.213、 ep(3)=0.211、 ep(4)= 0.209、 ep(5)=0.207、 ep(6)=0.205、 ep(7)=0.203、ep(8)=0.200。双轴压缩试验离散元模拟中所用的其他参数见表1。

文中模拟了含与不含范德华力2 种试样的双轴压缩试验，分别对应于月球环境与地球环境下的月壤剪切带分析，2 种试样的峰值内摩擦角φ= 47.4°，其中含范德华力试样的黏聚力c=2.5 kPa，不含范德华力试样无黏聚力。

图6 为双轴试验成样后的示意图，将试样中颗粒颜色分为网格状是为了在试验过程中清楚地观测试样的局部变形情况。双轴试验数值模拟采用颗粒为边界施加围压，这种边界可以模拟三轴等室内试验中的橡皮膜边界，该边界可以产生变形，称为柔性边界。所成试样的颗粒级配曲线如图7 所示，考虑到要模拟与真实月壤一样的颗粒级配，离散元模拟中需要的颗粒数目会非常大，计算耗时会非常长的原因，文中选择粒径跨度较小的级配曲线。该级配与蒋明镜等所研制的TJ-1 模拟月壤的级配曲线中的大颗粒分布相一致[19－20]。

表1 双轴压缩试验模拟参数 Table 1 Parameters used for biaxial compression tests

模拟过程分为成样、固结和压缩3 个步骤。成样过程如上文所述，固结过程中水平向、竖向按相应固结压力（文中试验围压为25 kPa）双向固结，固结结束后，采用伺服系统保持试样围压不变，竖向施压，施压过程通过应变控制的方法，本文采用的应变率为5%/min，在施压过程中记录试样应力、应变、孔隙比及平均转动率（APR）等参量的变化。

4 月壤双轴试验剪切带分析

4.1 试样宏观剪切特性

图8、9 分别给出了2 种试样的应力-应变关系曲线和体变曲线。从图可以看出，双轴压缩试验过程中，强度随轴向应变的增加而增加，达到峰值后，试样出现应变软化现象，强度减小最终维持在试样的残余强度，含范德华力试样的峰值强度与残余强度较高。2 种试样的体变规律也基本一致，都是先剪缩后剪胀，到达峰值点后都有一定的剪缩，但含范德华力试样的最终体变较小。说明原位月壤强度相对于取回地球的月壤样品的室内试验强度指标要高，这是由于在月面环境下月壤具有一定的黏聚力，这对于月面上的工程建设来说是有利的因素。

图8 试样应力-应变曲线 Fig.8 Stress-strain curves of two samples

图9 试样体变曲线 Fig.9 Volume strain curves of two samples

4.2 试样局部变形

试验过程中记录1 个时间段内每个颗粒的初始与最终位置，这样可以得到这个时间段内试样中的速度分布情况，以箭头的形式画出不同轴向应变时的速度场如图12、13 所示。箭头方向代表速度的方向，箭头大小代表速度相对大小。从速度场图中可以明显看出剪切带的位置与形状，在轴向应变为1%时，试样中速度较大的区域主要集中在试样上下两端，在轴向应变为2%时，可以明显地看出2 种试样中剪切带的分布情况，随着试样剪切带的继续发展，试样中的颗粒速度逐渐以相反的方向分布在剪切带的两边，剪切带内的颗粒速度表现出一定的旋涡状，说明剪切带两侧的颗粒发生了相对滑动。与图11 类似，月球环境下与地球环境下的试样破坏模式不同。

图10 不含范德华力试样网格变形 Fig.10 Grid deformation of sample without Van der Waals forces

图11 含范德华力试样网格变形 Fig.11 Grid deformation of sample with Van der Waals forces

4.3 速度场

双轴试验在竖向加载前将试样划分为13×25 个矩形网格，每个网格中的颗粒用不同颜色标记，这样便于观察试样在加载过程中的局部变形，可以明显地观察到试样剪切带的变化情况。2 种试样的变形过程如图10、11 所示。从图中可以看出，2 种试样在竖向应变为2%时，剪切带开始出现，随着竖向应变的继续增大，试样中的剪切带越来越明显。整个过程中剪切带内的网格变形很大，而带外的网格变形始终较小，说明剪切带内发生了很大的剪切变形，而带外的剪切变形很小。2 种试样剪切带形式明显不同，不含范德华力试样为1 条贯通的剪切带，而含范德华力试样则有2 条垂直的剪切带，1条主带，1 条次带，两者于试样底部交叉。这说明月球环境下试样的破坏模式与地球环境下可能不同，值得注意。

图12 不含范德华力试样内部速度场 Fig.12 Velocity field of sample without Van der Waals forces

图13 含范德华力试样内部速度场 Fig.13 Velocity field of sample with Van der Waals forces

4.4 试样内部孔隙比分布

图14、15 分别为2 种试样在不同轴向应变时试样内部孔隙比分布云图，从图可以看出，整个试样过程中试样内部孔隙比的变化。2 种试样在轴向应变为1%时，试样内部的孔隙比分布比较均匀，而在轴向应变达到2%时，试样中的孔隙比出现局部集中现象，1 个条带内的孔隙比显著增大，即试样中出现剪切带。剪切带附近的孔隙比等值线较为密集，而距离剪切带较远区域的孔隙比等值线较为疏松，说明剪切带内孔隙比梯度较大，剪切带内发生局部的剪胀。可见，2 种环境下，试样剪切带内的孔隙比都很大，都会发生应变局部化现象。

式中：ω 为APR 值；N 为接触点个数； Rk为第k个接触点处2 个颗粒半径平均值，其计算如式（8）所示；、分别为第k 接触点处颗粒1、2 的转动速度。、分别为第k 接触点处颗粒半径。

图14 不含范德华力试样内部孔隙比(单位: m) Fig.14 Void ratio of sample without Van der Waals forces (unit: m)

图15 含范德华力试样内部孔隙比(单位: m) Fig.15 Void ratio of sample with Van der Waals forces (unit: m)

4.5 试样内部转动场

文中的转动场分析采用的是蒋明镜等[21]提出的平均纯转动率APR（average pure rotation rate），它是1 个表征颗粒间相互转动的微观参量，其数值与颗粒转动及颗粒半径有关，计算表达式为

制作的样机实物图如图2所示,该样机拉压刚度不小于1×108 N/m,量程为-1 500 N/(N·m)～1 500 N/(N·m),误差要求控制在3.75%范围内。

图16、17 分别为不含与含范德华力试样在不同轴向应变时内部APR 的分布图。从图16 中可以看出，在轴向应变为1%时，不含范德华力试样内部的APR 值出现局部的集中，且分布较为分散；含范德华力试样内部APR 值在次剪切带处出现局部的集中。当轴向应变达到2%时，2 种试样试样内部APR 值的分布均出现了明显的贯通带状，带内APR绝对值较大，而带外APR 值基本为0，随着竖向应变的继续增大，APR 值的分布形式不再改变，带外的APR 值始终保持为0 左右，说明在整个试验过程

图16 不含范德华力试样内部APR(单位: m) Fig.16 APR of sample without van der Waals forces (unit: m)

图17 含范德华力试样内部APR(单位: m) Fig.17 APR of sample with Van der Waals forces (unit: m)

中，试样中发生相对转动的区域主要是剪切带区域，其他区域中的颗粒基本不发生相对转动。说明月球环境下，范德华力作用有减小颗粒转动的趋势。

4.6 剪切带性状分析

文中统计了2 种试样在轴向应变为4%、8%时剪切带的厚度及倾角。剪切带倾角主要通过试样的局部变形网格图、速度场、孔隙比分布图、转动场来统计，而剪切带厚度主要通过试样的孔隙比分布图、转动场来统计。具体统计结果见表2，其中剪切带宽度以试样的平均粒径（1.27 mm）为单位。由表中统计数据可知，不含范德华力试样的剪切带倾角在59.75°左右，剪切带宽度为试样平均粒径的23.9 倍；含范德华力试样的剪切带倾角在60.25°左右，剪切带宽度为试样平均粒径的19.65 倍。2 种试样的剪切带宽度均随着轴向应变的增大而有所增大，且含范德华力试样的剪切带宽度比不含范德华力试样小。说明月壤在地球环境与月壤环境下的剪切带性状会有所不同。

2．入库文献的选择必须遵循完整性原则。不采用传统语料库建设中随机择句或择段的方式选择语料，不论文献长短都进行全文收录，确保古籍文献的全貌；不要求语料库规模和入库语料追求“大而全”，反对不加选择地简单堆砌罗列所有文献。

图18给出了2种试样的体应变随着偏应变的变化曲线。由图可以得出，2 种试样的最大剪胀角ψ 均为26.4°左右。

表2 两种试样剪切带厚度、倾角统计表 Table 2 Thickness and inclination of the shear band

图18 两种试样体变-偏应变曲线 Fig.18 Volume strain-partial strain of two samples

对于剪切带倾角大小，国内外一般有以下几种观点：

式中：α 为试样剪切带的倾角；φ 为试样的内摩擦角；ψ 为试样的剪胀角。文中将2 种试样剪切带倾角的大小与以上几个公式的理论值作了比较，如表3 所示，其中离散元模拟的试样剪切带倾角大小取平均值。通过表3 的对比可以看出，式（9）所估计的剪切带倾角值对于2 种试样来说都过高，2 种试样剪切带倾角的离散元模拟值都在式（10）、（11）的计算值之间。说明地球环境与月壤环境下月壤试样的剪切带倾角都可以用式（10）、（11）来估计。

表3 2 种试样剪切带倾角对比 Table 3 Contrast of Shear band inclination

5 结 论

（1）月面环境下月壤试样双轴试样的剪切峰值强度比地面环境下月壤的剪切强度高，说明原位月壤的强度相对于取回地球的月壤样品的室内试验强度指标要高，这是由于在月面环境下月壤中具有一定的黏聚力，这对于未来在月面上的工程建设来说是有利的因素。

（2）月面环境下月壤颗粒间的范德华力对月壤试样的剪切带形式有较大影响，2 种试样剪切带形式明显不同，说明在月球环境下月壤的破坏模式会发生改变。

（3）月球环境下月壤试样的剪切带宽度比地面环境下的试样小，2 种环境下试样的剪切带倾角均在45°+ψ/2与45°+ (φ +ψ)/4之间。

[1] HEIKEN G H, VANIMAN D T, FRENCH B M. Lunar Sourcebook[M]. London: Cambridge University Press, 1991.

[2] 欧阳自远. 月球科学概论[M]. 北京: 中国宇航出版社, 2005.

[3] 郑永春, 欧阳自远, 王世杰, 等. 月壤的物理和机械性质[J]. 矿物岩石, 2004, 24(8): 14－19. ZHENG Yong-chun, OUYANG Zi-yuan, WANG Shi-jie, et al. Physical and mechanical properties of lunar regolith[J]. Journal of Mineralogy and Petrology, 2004, 24(8): 14－19.

[4] CHRISTENSEN E M, BATTERSON S A, BENSON H E, et al. Lunar surface mechanical properties-Surveyor 1[J]. Journal of Geophysical Research, 1967, 72(2): 801－813.

[5] JAFFE L D. Surface structure and mechanical properties of the Lunar Maria[J]. Journal of Geophysical Research, 1967, 72(6): 1727－1731.

[6] COSTES N C, MITCHELL J K. Apollo 11 soil mechanics investigation[C]//Proceedings of the Apollo 11 Lunar Science Conference. Houston: [s. n.], 1970, (3): 2025－2044.

[7] SCOTT R F, CARRIER W D, COSTES N C, et al. Apollo 12 soil mechanics investigation[J]. Géotechnique, 1971, 21(1): 1－14.

[8] MITCHELL J K, BROMWELL L G, CARRIER W D, et al. Soil mechanical properties at the Apollo 14 site[J]. Journal of Geophysical Research, 1972, 77(29): 5641－5664.

[9] NAKASHIMA H, FUJII H, OIDA A, et al. Parametric analysis of lugged wheel performance for a lunar microrover by means of DEM[J]. Journal of Terra- mechanics, 2007, 44: 153－162.

[10] NAKASHIMA H, FUJII H, OIDA A, et al. Discrete element method analysis of single wheel performance for a small lunar rover on sloped terran[J]. Journal of Terramechanics, 2010, 47: 307－321.

[11] 杨艳静, 向树红. 月球车刚性车轮与模拟月壤相互作用有限元仿真和试验验证[J]. 强度与环境, 2010, 37(1): 47－52. YANG Yan-jing, XIANG Shu-hong. Finite element simulation and experimental validation of the interaction between lunar rover wheel and lunar soil stimulant[J]. Structure & Environment Engineering, 2010, 37(1): 47－52.

[12] Itasca Consulting Group Inc. The manuals of particle flow code in 2-dimension.Version 3.1[M]. Minneapolis: [s. n.], 2004.

[13] 郑敏, 蒋明镜, 申志福. 简化接触模型的月壤离散元数值分析[J]. 岩土力学, 2011, 31(增刊1): 766－771. ZHENG Min, JIANG Ming-jing, SHEN Zhi-fu. Discrete element numerical analysis of lunar soil with a simplified contact model[J]. Rock and Soil Mechanics, 2011, 31(Supp. 1): 766－771.

[14] JIANG M J, YU H S, HARRIS D. A novel discrete model granular material incorporating rolling resi- stance[J]. Computers and Geotechnics, 2005, 32(4): 340－357.

[15] PERKO H, NELSON J, SADEH W. Surface cleanliness effect on lunar soil shear strength[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2001, 27(4): 371－383.

[16] ADAMSON A W. Physical chemistry of surfaces[M]. New York: Wiley, 1990.

[17] NELSON J D, VEY E. Relative cleanliness as a measure of lunar soil strength[J]. Journal of Geophysical Research, 1968, 73(12): 3747－3764.

[18] JIANG M J, KONRAD J M, LEROUEIL S. An efficient technique for generating homogeneous specimens for DEM studies[J]. Computers and Geotechnics, 2003, 30(7): 579－597.

[19] 蒋明镜，李立青. TJ－1 模拟月壤的研制[J]. 岩土工程学报, 2011, 33(2): 209－214. JIANG Ming-jing, LI Li-qing. Development of TJ-1 lunar soil simulant[J]. Chinese Journal of Geotechnical Engineering, 2011, 33(2): 209－214.

[20] JIANG M J, LI L Q, SUN Y G. A new lunar soil simulant in China[C]//Earth & Space 2010, Proceedings of the 12th Biennial International Conference on Engineering, Science, Construction and Operations in Challenging Environments. Honolulu: [s. n.], 2010: 3617－3623.

[21] JIANG M J, YU H S, HARRIS D. Kinematic variables bridging discrete and continuum granular mechanics[J]. Mechanics Research Communication, 2006, 33(5): 651－666.

[22] VERMEER P A. The orientation of shear bands in biaxial tests[J]. Geotechnique, 1990, 40(2): 223－236.

[23] ROSCOE K H. The influence of strains in soil mechanics[J]. Geotechnique, 1970, 20(2): 129－170.

[24] ARTHUR J, DUNSTAN T, AL-ANI Q, et al. Plastic deformation and faliure in granular media[J]. Geotechnique, 1977, 27(1): 53－74.