0.   引言

月球是距离地球最近的天体，拥有丰富且独特的资源。其中，月球原位资源利用（in situ resource utilization）是中国21世纪月球探测的首要科学目标之一，同时也是保障月球科考站建设和顺利运行的重要途径^[1]，月球基地的建立将为中国月球探测提供一个良好的生存和工作环境^[2-3]。月壤是月表岩石受到小天体频繁撞击、破碎、溅射和太阳的辐射而形成的，广泛覆盖于月表，无论是月球基地建设，还是月球原位资源利用，都以月壤作为主要对象。但由于月壤有着与地球土体明显不同的颗粒形态特征和物理力学性质，因此研究其静力学特性，特别是其流动特性对月球科研站的建设及原位资源利用中可能遇到的滑坡、边坡失稳等实际工程问题提供指导性意见^[4]。

颗粒体流动性是指颗粒体在受到外部作用力时能够流动的性质及程度，比如流动性好则休止角（angle of repose）相对较小，滑动距离相对较长。转筒试验因设备结构简单，体积灵活，便于研发设计，广泛用于颗粒及粉体材料的混合、筛选及流动特性等方面研究^[5-9]。Brucks等^[8]考虑有效重力场对剪切颗粒流的影响，通过离心机进行了不同有效重力场下单一粒径的转筒试验，结果发现，随转筒转速增加或有效重力场减小，动休止角均有所增加，依据弗劳德数，可将1g重力场下的颗粒流动状态可以按比例扩展到其他重力场情况。Kleinhans等^[9]通过飞机周期抛物线飞行，进行了不同材料不同重力场下的转筒试验，结果发现对于所有材料，重力场从1g减小至0.1g时，静休止角随着重力的减小而增加大约5°，而动休止角则减小约10°。Klein等^[10]通过飞机周期抛物线飞行，进行了不同重力场下的转筒试验，发现动休止角（dynamic angle of repose）与重力的平方根呈倒数关系。由于上述试验环境难以长时间维持，且粒间接触信息难以定量化描述，难以从微观层面探究其宏观力学特性产生的机理。

自1979年Cundal等^[11]提出离散单元法（distinct element method）以来，凭借其易实现、成本低、可重复性高等优点，广泛用于岩土领域的相关模拟研究，是研究宏微观土力学的有效工具之一^[12]。Walton等^[13]进行了不同颗粒形状、粒间黏聚力及重力场的转筒DEM模拟，结果发现颗粒不规则形状和粒间黏聚力的增加会使休止角明显增加，且随着重力场的减小，在相同黏聚力情况下，材料休止角明显增大。Nakashima等^[14]通过试验和离散单元法进行了不同重力场下的休止角试验，发现重力场对休止角影响可以忽略不计。Chen等^[15]使用离散单元法对休止角和重力场关系进行了相关研究，发现重力场对休止角的影响可以忽略不计。Sunday等^[16]使用改进的开源代码CHRONO，进行了转筒试验DEM模拟，模拟结果发现当弗劳德数保持恒定时，重力场对动休止角影响不明显。Karapiperis等^[17]通过离散单元法进行推板试验，研究了不同重力环境对砂土峰值内摩擦角、残余内摩擦角和休止角的影响，上述研究发现重力场对休止角影响不明显。Antony等^[18]使用离散单元法探究了低重力环境对休止角的影响，发现随着重力场由1/6g增至2g，其休止角由30.6°增长至32.04°。Elekes等^[19]使用离散单元法研究了不同颗粒尺寸和重力环境下休止角的变化情况，发现在相同粒径条件下，休止角随着重力场的增大而减小。Jiang等^[20]使用三维完整接触模型来模拟月壤颗粒，进行了不同重力环境下的转筒DEM模拟，结果表明随着重力加速度的增加动休止角也增加，且颗粒最大速度也随重力场的增加而增加。

由于地球上真实月壤十分稀有，无法开展真实月壤的单元尺度力学特性研究，导致上述试验多使用颗粒、粉体材料或模拟月壤进行相关研究，而上述DEM模拟也没有从粒观层面考虑真实月壤的特有物理力学特性，进而难以模拟出真实月壤的宏观力学响应。因此，本文使用团队提出的可以反映颗粒形状影响的三维含抗转动和抗扭转完整接触模型^[21]，基于嫦娥五号任务所携带返回部分月壤样品的形状特征确定抗转动系数β（接触模型中描述颗粒形状的参数）与粒径的对应关系^[22]，进行了不同材料（玻璃珠、丰浦砂、月壤），不同围压（50，100，200 kPa）及不同转速（2.622，3.211，3.708 rad/s）工况下的三轴和转筒DEM模拟。随后着重分析了月壤与另两种材料在静力学及流动性方面的差异，可为提高月球原位资源利用效率、科考站的合理选址及安全建设提供理论及工程指导意见。

1.   考虑月壤粒形特征的接触模型简介

使用合理的接触模型模拟月壤颗粒是进行月壤流动特性研究的关键条件之一。月壤颗粒形状不规则，以长条状、棱角及次棱角状居多^[23]。月壤粒间接触以面-面接触为主，具有一定的抗转动能力。因此，本次离散元模拟使用团队提出的考虑抗弯转及抗扭转的三维完整接触模型^[21]，该模型可以较好的模拟由于颗粒不规则形状导致的高内摩擦角等宏观特性，在模拟不规则颗粒上具有一定优势^[12]。

1.1   三维完整接触模型

离散元经典接触模型^[11]仅考虑法向及切向力学响应，且粒间接触被视为点接触，仅适用于模拟低强度砂土。Jiang等^[21]假设粒间接触面为圆形，进一步提出了考虑抗弯转及抗扭转的三维完整接触模型，该模型包含完整4个方向力学响应，如图 1所示，其中在粒间接触处引入抗转动力矩（使用抗转动系数β来表征力矩大小）来反映颗粒形状对其宏观力学特性产生的影响，可较好的反映月壤高内摩擦角特性，具体推导及计算公式等信息详见文献[21]。

图  1  三维完整接触模型力学响应图^[21]

Figure  1.  Mechanical responses of three-dimensional complete contact model^[21]

下载: 全尺寸图片 幻灯片

1.2   嫦娥五号月壤颗粒形态特征

本次离散元模拟所使用月壤数据，基于中国探月工程嫦娥五号任务所携带返回部分月壤样品的图像信息^[24]，图像经过二值化处理并对颗粒形状参数进行捕捉，统计了粒径形状（球度）与粒径之间的关系，如图 2所示，具体信息详见文献[22]。从图 2中可以看出嫦娥五号返回样本离散程度较大，但整体上呈现出随着颗粒直径增大球度减小的趋势。

图  2  月壤颗粒球度与直径间关系^[22]

Figure  2.  Relationship between sphericity and diameter of lunar regolith ^[22]

下载: 全尺寸图片 幻灯片

1.3   月壤颗粒形态特征与抗转动系数关系

本次离散元模拟依照团队已有研究成果^[22]，依照相关公式将月壤颗粒球度转换成抗转动系数并统计了两者之间的关系，如图 3所示。从图 3中可以看出虽然总体上离散程度较高，但月壤颗粒抗转动系数仍表现出一定随着颗粒直径增加而增大的趋势。

图  3  月壤直径与抗转动系数间关系^[22]

Figure  3.  Relationship between diameter of lunar regolith and rolling resistance coefficient^[22]

下载: 全尺寸图片 幻灯片

基于上述月壤粒形统计数据，可以建立月壤粒径与粒间接触抗转动系数β之间的对应关系，基于文献[22]中统计嫦娥五号月壤粒形信息，得出抗转动系数β = -4.429P_c / P_s + 4.43312，式中P_c为与颗粒面积相同圆的周长，P_s为颗粒周长。基于上述月壤直径与抗转动系数间关系，进而从粒观层面更精细化的模拟不同粒径月壤间相互作用，使离散元模拟可以更好的反演真实月壤的独特宏观力学特性。

2.   试验试样制备、步骤及静力学分析

在离散元数值模拟中，试样的制备及模拟参数选取是保证模拟结果准确性的先决条件之一。本次研究主要进行了不同围压三轴试验及不同转速转筒试验离散元模拟，下面分别对其参数选取、试样制备及试验步骤进行介绍。

2.1   颗粒级配及接触参数选取

为研究探讨嫦娥五号月壤粒形特征对其静力学和流动特性的影响，选择丰浦砂和玻璃珠作为对比材料，以揭示月壤粒形特性与地球材料在静力学和流动特性上的差异。

参考已有文献离散元模拟颗粒级配选取^[21-22]，同时考虑合适的计算时间及准确性，本次模拟中3种材料DEM模拟颗粒级配如图 4所示，其中玻璃珠由11种粒径构成，最大最小粒径分别为0.3，0.1 mm，中值粒径d₅₀为0.19 mm；丰浦砂由14种粒径构成，最大最小粒径分别为0.364，0.106 mm，中值粒径d₅₀为0.204 mm；月壤由23种粒径构成，最大、最小粒径分别为0.36，0.106 mm，中值粒径d₅₀为0.193 mm。

图  4  DEM模拟颗粒级配

Figure  4.  Grain-size distribution curves for DEM simulations

下载: 全尺寸图片 幻灯片

参考已有文献离散元模拟材料参数选取^{[21-22, 25]}，本次转筒模拟材料参数选取如表 1所示，参考已有转筒试验离散元模拟，较小量级接触模量的使用对模拟结果影响可以忽略^[26-27]，为了提高计算效率，相较于三轴DEM模拟，对不同材料的接触等效模量按相同比例缩放了10²倍。转筒底部墙体摩擦系数为1.0。在月壤抗转动系数选取上，参照团队已有研究成果^[22]，按照不同粒径对应不同抗转动系数，详见表 2，对于不同粒径颗粒间接触的抗转动系数，取接触两端颗粒抗转动系数平均值。

表  1  转筒离散元数值试样的材料参数

Table  1.  Material parameters of DEM specimens for rotating drum

参数	颗粒密度/ (kg·m-3)	接触等效模量/MPa	接触刚度比	摩擦系数	抗转动系数	破损系数	法向黏滞阻尼	切向黏滞 阻尼
玻璃珠	2500	6.5	1.5	0.15	0.01	2.1	0.13	0.13
丰浦砂	2655	7.0	5.0	0.5	0.25	4.0	0.4	0.4
月壤	2320	8.0	5.0	0.5	详见表 2	4.0	0.8	0.8

下载: 导出CSV 

| 显示表格

表  2  不同粒径月壤抗转动系数^[22]

Table  2.  Rolling resistance coefficient for lunar regolith with different diameters^[22]

组数	颗粒直径/mm	抗转动系数β		组数	颗粒直径/ mm	抗转动系数β
1	0.36	1.009		13	0.216	0.589
2	0.35	0.923		14	0.207	0.576
3	0.316	0.917		15	0.195	0.573
4	0.305	0.910		16	0.184	0.562
5	0.293	0.818		17	0.175	0.558
6	0.282	0.810		18	0.165	0.558
7	0.271	0.781		19	0.156	0.533
8	0.268	0.757		20	0.136	0.533
9	0.254	0.755		21	0.127	0.518
10	0.247	0.733		22	0.114	0.356
11	0.235	0.711		23	0.106	0.335
12	0.225	0.615

下载: 导出CSV 

| 显示表格

与转筒试验DEM模拟不同，在准静态三轴试验DEM模拟中，为提高计算效率，采用密度放大法对颗粒密度进行放大，对不同材料密度按相同比例放大了10⁵倍。在准静态DEM模拟中，采用局部阻尼代替动态模拟中的黏滞阻尼，局部阻尼取值0.7，其他材料参数与转筒模拟中保持一致，如表 1所示。

2.2   离散元试样制备

在三轴DEM试样制备中，采用Jiang等^[28]提出的分层欠压法制成长宽高分别为5.2 mm×5.2 mm×10.4 mm，初始孔隙比为0.95的松样，如图 5（a）所示。在转筒DEM试样制备中，转筒半径R为7.135 mm，宽度W为2.039 mm，且宽度方向采用周期边界来尽量消除转筒宽度方向的尺寸效应影响。转筒填充率（颗粒堆积体积与转筒体积比值）为50%，初始制样孔隙比为0.7，如图 5（b）所示。

图  5  离散元试样示意图

Figure  5.  Sketch of numerical specimen for DEM simulations

下载: 全尺寸图片 幻灯片

2.3   离散元三轴及转筒试验步骤和工况

DEM三轴压缩试验及转筒试验模拟步骤与室内试验基本保持一致，其中对于三轴试验DEM模拟可分为3个阶段：①竖向预压阶段。对初始制备完成的试样施加12.5 kPa的竖向荷载进行预压，该过程由伺服系统控制上下墙体完成；②等向固结阶段。对预压完成的试样施加等向荷载（50，100，200 kPa）进行压缩。该过程由伺服系统控制六面刚性墙体完成；③剪切阶段。对固结完成的试样的上下墙体以1.2%/min的应变速率相向移动。同时，由伺服系统控制侧墙，保持试样侧向围压恒定。本次三轴模拟针对3种不同材料分别进行了3种不同围压（50，100，200 kPa）下模拟，共计9组工况。

为探究不同材料的流动特性（如稳定倾角等），在转筒试验中合理流动状态的选取至关重要。多数研究学者^[29-30]一般用弗劳德数F_r = ω²R/g来描述转筒中颗粒的流动状态，式中ω为转筒的角速度，R为转筒的半径，g为当前环境下的重力加速度，其中在滚动状态（rolling regime）下可以观察到颗粒连续稳定的向下流动，此时其倾角基本保持稳定，是研究颗粒材料流动特性的理想状态。当弗劳德数F_r为10^-4~10^-2时^[31]，转筒中颗粒材料处于滚动状态，因此为了减少计算时间，本次离散元模拟弗劳德数F_r取值分别为5×10^-3，7.5×10^-3，1×10^-2。

对于转筒试验DEM模拟可分为两个阶段：①重力固结阶段。初始试样制备完成后，在地球重力加速度g = 9.81 m/s²下进行重力固结并平衡；②转筒转动阶段。依照已确定的弗劳德数F_r（5×10^-3，7.5×10^-3，1×10^-2）及重力加速度计算得出转筒转速ω（2.622，3.211，3.708 rad/s），使转筒沿顺时针方向匀速转动。本次转筒模拟针对3种不同材料和3种不同转筒转速（2.622，3.211，3.708 rad/s）下模拟，总计9种工况。

为研究转筒内颗粒流局部流动特性变化情况，依照材料流动所呈倾角及坐标原点进行局部坐标转换（见图 6），参照此局部坐标系并沿z′轴负方向布置测量球，其中测量球半径R_m=0.595 mm，测量球间距h_m= 0.595 mm，测量球个数为12个。

图  6  测量球位置示意图

Figure  6.  Sketch of measuring spheres for DEM simulations

下载: 全尺寸图片 幻灯片

2.4   应力-应变及体变-应变关系

图 7（a）为偏应力应变关系曲线，从图中可以看出，随着材料的改变（由玻璃珠至月壤），其材料的峰值偏应力也有显著增加，也表现出与围压正相关关系，且月壤峰值偏应力随围压增大的幅度最大，而玻璃珠最小，丰浦砂居中。随着材料的改变（由玻璃珠至月壤），其曲线逐步由应变硬化向应变软化转变且在丰浦砂及月壤工况中随着围压的增加该转变趋势也更加明显。图 7（b）为体变应变关系曲线，从图中可知，玻璃珠工况表现出明显的剪缩现象，而丰浦砂及月壤工况则表现出明显的剪胀现象，且在剪切过程中丰浦砂体变要略大于月壤工况。

图  7  偏应力及体变演化规律

Figure  7.  Evolution of deviatoric stress and volume change

下载: 全尺寸图片 幻灯片

上述现象表明随着颗粒形状不规则程度（接触抗转动系数β）的增加，会显著限制颗粒相对转动，进而导致试样的偏应力明显增加且在剪切过程中体变明显增大。

2.5   抗剪强度指标

参照图 7中不同材料的应力应变数据，由于上述曲线偏应力无明显峰值，因此本次抗剪强度指标计算参考《土工试验方法标准》^[32]选取轴向应变15%所对应的偏应力值。

图 8为不同材料的抗剪强度指标（内摩擦角及黏聚力），由图 8可知，作为典型的颗粒材料，玻璃珠的内摩擦角及表观黏聚力均明显小于丰浦砂及月壤，而相比于经典岩土材料（丰浦砂），月壤的内摩擦角及表观黏聚力均比其有一定增大。上述现象表明随着颗粒形状不规则程度（接触抗转动系数β）的增加，逐渐限制颗粒的相对转动，导致其抗剪强度指标的增大。由于本次模拟中颗粒间并未加入吸引力或黏聚力，因此抗剪强度指标中的黏聚力是一种表观黏聚力，主要由粒间“互锁效应”造成，抗转动系数的增加会提高粒间的“互锁效应”，进而导致材料表现出明显的表观黏聚力。其中较高的表观黏聚力可能是粒间接触抗转动系数和摩擦系数较大导致的。

图  8  黏聚力和内摩擦角

Figure  8.  Cohesion and internal friction angle

下载: 全尺寸图片 幻灯片

3.   流动特性模拟结果

3.1   倾角变化规律

图 9为不同材料在转筒中倾角演化规律，由图可知，在本次研究中，将材料第一次开始向下滑落所对应的最大倾角定义为启滑角θ_l，将随后滑动过程中对应的最小角定义为休止角θ_r，最后将稳定阶段的倾角定义为稳态倾角θ_s。从图中可以看出，启滑角、休止角及稳态倾角均随着材料的改变（由玻璃珠至月壤）有所增加，上述倾角均表现出随颗粒形状不规则程度（接触抗转动系数β）增加而增加的趋势。在转动状态（rolling regime）中，随着弗劳德数的增加（对应转筒转速增加），所有材料的启滑角、休止角及稳态倾角均有略微的增加，表现出一定的速度正相关性。

图  9  转筒倾角演化规律

Figure  9.  Evolution of inclination angle in rotating drum

下载: 全尺寸图片 幻灯片

表 3为不同转筒弗劳德数及材料的不同类型倾角信息。

表  3  不同倾角变化规律

Table  3.  Variation for different types of inclination angles

转筒弗劳德数 Fr	内摩擦角/ (°)	倾角类型
		启滑角θl/(°)	休止角θr/(°)	差值 θl-θr	稳态倾角θs/(°)
1.0×10-2	16.42	42.72	26.21	16.51	28.99
		69.39	37.31	32.08	42.37
		81.23	44.06	37.17	47.69
7.5×10-3	31.7	41.12	25.80	15.32	28.48
		68.88	37.87	31.01	41.93
		81.21	41.99	39.22	47.20
5.0×10-3	38.1	39.88	25.38	14.50	28.05
		67.69	37.11	30.58	41.79
		79.58	43.17	36.41	46.89

下载: 导出CSV 

| 显示表格

从表 3中可以看出，不同材料启滑角和休止角差值均约等于材料内摩擦角。该差值随着转筒转速的增加有一定的增加趋。因此针对不同材料，如需通过转筒试验初步判断材料内摩擦角，合理转速的选取是十分必要的。针对不同倾角类型与材料内摩擦角之间的准确关系，本文只进行了初步探索，相关深入研究将在后续工作中进行。

3.2   剖面速度分布

图 10为稳定状态下转筒剖面速度图，其中速度依照图 6中坐标变换，求得每个测量球内颗粒沿x′方向的平均速度，图 10中横坐标中V表示平均速度，ω为转筒转速，R为转筒半径，其乘积为转筒外侧切向速度，纵坐标中M_p表示测量球球心在z′方向大小。从图 10中可以看出，剖面速度表现为底部颗粒沿x′负方向运动，随着位置的升高，颗粒逐渐由x′负方向运动过渡到正方向运动，且下部区域颗粒速度变化主要呈线性趋势。上述现象产生的原因可能是由于底部区域颗粒主要受到转筒壁提供摩擦力的影响，其最底部颗粒速度约等于与转筒转动的切向速度，此时颗粒间基本保持相对稳定状态，粒间基本没有相对运动。随着距转筒底部距离的增高，颗粒受底部筒壁的影响逐渐减小而受表面向下滑动影响逐渐增大，粒间相对运动逐渐增强。因此可将底部速度线性变化区域称为“被动区”，而将上部以滑动为主区域称为“主动区”。

图  10  转筒剖面速度图

Figure  10.  Normalized velocity profile of rotating drum

下载: 全尺寸图片 幻灯片

针对不同转速下转筒工况，其无量纲化速度基本一致，说明转筒转速对其速度分布影响不明显。在下部区域不同材料颗粒速度变化趋势基本相同，主要呈线性趋势，而在表面，玻璃珠颗粒速度最大，丰浦砂次之，月壤最小。上述现象表明对于表层颗粒，玻璃珠无量纲化速度衰减最快，可间接对应其所受剪切速率（单位深度内速度差值）最大，丰浦砂次之，月壤最小。对于不同材料，其主被动区域分界点位置也不尽相同，其中玻璃珠工况分界点位置靠下，丰浦砂工况居中，而月壤工况则位置最靠上。

3.3   剖面孔隙比分布

图 11为稳定状态下转筒剖面孔隙比图，从图中可以看出，对于不同材料及转速下工况，其局部孔隙比曲线均表现为上部及底部偏大而中部偏小的“凹”字形状。上述现象产生的原因可能是由于不同尺寸颗粒在转筒中运动时会发生分选效应（segregation），导致大颗粒多集中在表层及底层，小颗粒则多集中在中部，进而造成中部孔隙比的减小和上部及底部孔隙比的增大。

图  11  转筒剖面孔隙比图

Figure  11.  Normalized porosity profile of rotating drum

下载: 全尺寸图片 幻灯片

转筒转速对下部孔隙比影响不明显，而上部孔隙比则表现为随着转速的增加略有增长。对于不同材料，玻璃珠工况孔隙比要明显小于丰浦砂及月壤工况，而对于丰浦砂及月壤工况，其上部中部孔隙比基本相同而底部孔隙比则表现为月壤工况要略大于丰浦砂工况。上述现象产生的原因可能是由于材料形状参数（抗转动系数β）的差异，即玻璃珠抗转动系数要明显小于丰浦砂和月壤，而随着颗粒抗转动系数的增大其承载相应荷载所需接触数目减小，致使颗粒局部体积改变量显著增大，进而导致其孔隙比的明显增高。

3.4   剖面配位数分布

图 12为稳定状态下转筒剖面配位数图，配位数表示单位颗粒所包含的接触数目^[33]，从图 12中可以看出，不同材料及转速工况下，转筒上层配位数要明显小于转筒中部及底部，且在上层，随着深度的增加，配位数会发生明显的增大。上述现象产生的原因与孔隙比及速度变化原因一致，在上部区域（主动区），颗粒受重力作用影响，粒间相对运动较剧烈，导致对应区域配位数明显减少；而在中下部（被动区），颗粒运动主要受转筒壁带动，粒间几乎没有相对运动，颗粒体保持相对稳定，进而其配位数稳定在较大值。

图  12  转筒剖面配位数图

Figure  12.  Normalized coordination number profile of rotating drum

下载: 全尺寸图片 幻灯片

对于不同材料，玻璃珠工况配位数变化最为显著，丰浦砂次之，而月壤工况配位数变化最小。上述现象产生的原因可能是由于材料颗粒形状（抗转动系数β）差异（抗转动系数玻璃珠最小，月壤最大），随着抗转动系数的增加，相同的力可以由较少的粒间接触承担，进而导致单位体积内接触数量的减少，进而导致孔隙比也相应增大，该规律与图 11规律较符合。

3.5   剖面平均粒径分布

图 13为稳定状态下转筒剖面平均粒径图，图中横坐标中d_mean表示测量球内颗粒考虑质量加权的平均粒径，d₅₀表示材料的中值粒径。从图中可以看出，转速对无量纲化平均粒径基本没有影响，而整体上表现为底层和表层平均粒径较大，中部平均粒径较小。上述现象产生的原因可能是在转筒转动过程中，不同粒径颗粒由于“分选作用”的影响导致不同深度颗粒会发生重分布，大尺寸颗粒会在底部和表层聚集，而小尺寸颗粒则会在中部聚集，进而导致不同深度平均粒径的变化。

图  13  转筒剖面平均粒径图

Figure  13.  Normalized mean particle diameter profile of rotating drum

下载: 全尺寸图片 幻灯片

对于不同材料工况，其不同深度平均粒径变化趋势也不相同，玻璃珠工况曲线呈现出典型的“凹”形，与已有“分选作用”研究结果类似，表现为上下部平均粒径较大，而中部较小，而由于颗粒形状不规则程度的增加（抗转动系数β的增加）丰浦砂工况则表现出不同规律，其上中部平均粒径变化不明显，中下部平均粒径有减小趋势但略小于玻璃珠工况。对于月壤工况，由于考虑了不同粒径的抗转动系数的不同，导致其不同深度平均粒径规律与上述两种材料均有所不同，其曲线在上中部呈现出凸出形状，而其“凹点”则有一定程度下移。上述现象说明颗粒形状（抗转动系数β）及不同粒径对应不同抗转动系数对颗粒在转筒中的“分选程度”会有一定程度的影响。该现象与已有孔隙比变化规律不太一致，造成该现象的可能原因是由于分选作用的影响，导致局部区域内颗粒数目和颗粒粒径分布不一致，进而使孔隙比规律与平均粒径规律不一致。以上只是对现有研究结果的初步解释，颗粒形状（抗转动系数β）和不同粒径对应不同抗转动系数对颗粒在转筒中的“分选程度”的影响的微观机理仍需后续进一步具体研究。

3.6   剖面惯性数分布

图 14为稳定状态下转筒剖面惯性数图，惯性数（inertia number）^[34]多数用以表征局部颗粒体的运动状态激烈程度，其中惯性数为

图  14  转筒剖面惯性数图

Figure  14.  Normalized inertial number profile of rotating drum

下载: 全尺寸图片 幻灯片

式中：I为惯性数；$\dot \gamma$为剪切应变率；d为平均粒径；P为所受压力；${\rho _s}$为颗粒密度。

从图 14中可以看出，惯性数表现为从表层至底层逐步减小的趋势，且随着转筒转速的增加对应工况整体惯性数也相应增加。上述现象产生原因可能是上部颗粒相较于下部颗粒，其粒间相互运动剧烈，剪切速率大且所受压力小，进而导致上部颗粒惯性数要明显大于下部颗粒。

对于不同材料工况，上部颗粒表现为玻璃珠工况惯性数最大，丰浦砂次之，月壤最小，而对于下部颗粒，则是月壤与丰浦砂工况相近，玻璃珠工况最小。上述现象产生的原因可能是由于颗粒不规则程度增大（抗转动系数β增大），限制了颗粒间的相互运动，导致上部颗粒惯性数的减小，对于底部颗粒，可能是由于颗粒不规则程度增大导致颗粒体整体稳定倾角的增大，从而引起颗粒整体重力的法向分量减小，最终导致惯性数的减小。

3.7   剖面有效摩擦系数与惯性数关系

为探究不同材料的有效摩擦系数μ^*与惯性数I之间关系，其中有效摩擦系数μ^*为测量球内局部颗粒体切应力与正应力比值，参考Jop等^[35]研究成果，在一定惯性数范围内惯性数与有效摩擦系数展现出一定关系：

式中：${\mu _s}$为惯性数低时（准静态）对应有效摩擦系数，约为tanφ；φ为材料内摩擦角；${\mu _2}$为惯性数高时对应有效摩擦系数；I₀为常数。

图 15为不同材料有效摩擦系数与惯性数关系图，由图可知，在一定惯性数范围内均随着惯性数的增大而增大。整体上玻璃珠工况有效摩擦系数最小，丰浦砂次之，月壤工况最大，且系数${\mu _2}$及差值${\mu _2} - {\mu _s}$，均表现为类似规律。对于常数I₀则表现为玻璃珠最大，丰浦砂工况次之，月壤工况最小。上述现象产生原因可能是由于颗粒不规则程度的增大（抗转动系数β增大），致使粒间更不易发生相对滚动，使局部切向应力显著增大，进而使有效摩擦系数增大。

图  15  有效摩擦系数与惯性数关系

Figure  15.  Relationship between effective friction coefficient and inertial number

下载: 全尺寸图片 幻灯片

4.   结论与展望

基于团队提出的三维完整接触模型^[21]和嫦娥五号样品的形状特征确定的抗转动系数β（接触模型中描述颗粒形状的参数）与粒径对应关系^[22]，进行了月壤三轴和转筒三维DEM模拟，对不同工况下月壤静力学及流动特性进行了对比分析，得到以下3点主要结论。

（1）由于月壤抗转动系数随粒径的增大而增大，导致其倾角及孔隙比大于丰浦砂和玻璃珠试样，而剪切速率及配位数则小于丰浦砂和玻璃珠工况。

（2）在相同转速下，主动区内月壤惯性数要明显小于丰浦砂和玻璃珠，表明颗粒不规则程度增大会导致土体惯性数减小，运动激烈程度降低。

（3）一定惯性数范围内，月壤有效摩擦系数大于丰浦砂和玻璃珠工况，这表明，颗粒整体不规则程度（抗转动系数β）的增大，以及粒径增大导致的颗粒不规则程度增加，均会降低材料的流动性。

以上研究仅在地面环境下探究了月壤流动特性，然而由于月面环境（低重力、高真空、高温差等）的作用，月壤颗粒之间的范德华力会表现得更为明显，使月壤表现出一定黏聚力特性，进而会对月壤的相关静力学及流动特性造成一定的影响。因此，笔者后续研究拟考虑月面环境中低重力及高真空对月壤静力学及流动特性的影响，并与地面环境下结果进行对比，以期更准确的模拟出月壤流动特性，进而为月球原位资源利用、科研站选址及建设提供一定理论及工程指导。