廖优斌， 蒋明镜， 陈有亮， 奚邦禄, 谭亚飞鸥，3

(1.上海理工大学 环境与建筑学院，上海 200093； 2.同济大学 土木工程学院，上海 200092；3.同济大学 岩土及地下工程教育部重点实验室， 上海 200092)

1 研究背景

20世纪50年代起，美国的Survey计划、Apollo计划及前苏联的Luna计划开启了人类探月计划的第一轮高潮，揭开了月球神秘的面纱[1],发现月壤的形成机理、所处的重力场环境以及力学特性与地球土壤差别较大。21世纪，欧盟的曙光计划、印度的飞船计划、日本的月亮女神计划和中国的嫦娥探月计划掀起了第二轮探月高潮[2]，将人类对月球的认识推到了一个新的高度。通过更高分辨率的影像数据发现月球表面存在不少的高陡边坡，其坡度一般都大于30°，月海区可达48.9°，月陆区甚至可达55.7°[3]。根据Apollo15探测器对月球表面雨海东部边缘现场影像数据分析，部分地区发现了滑坡运动的证据[4]。人类活动常使自然边坡的坡度变陡，一旦进入对月球资源的大规模开发，若工程设计不合理，就极易诱发坡体滑塌，甚至整体坍塌，这些地质灾害都会严重影响未来对月球资源开采的安全。我国的嫦娥探月计划也明确指出当探测第3期工程基本完成不载人的月球探测任务后，随后将实施载人登月和月球基地建设[2]。无论是资源开采还是基地建设，这些活动不可避免地将遇到不同坡角下的坡体工程。

图1 月壤边坡

坡角的大小会直接影响地面物质的流动及能量转换的规模和强度，月壤形成坡体的坡度比常规的略大，如图1所示。另外现场的沟槽试验、深陷足迹边缘月壤不塌落及粉土结块等现象[5]也表明了月壤具有较高的自稳能力，这些现象引起了国内外学者的广泛研究[3-5]。坡角的大小和材料休止角直接相关，而休止角的大小是否会受到重力场的影响至今尚无定论。月面环境下重力场大约为1.62 m/s2,仅仅为地面环境下的1/6[1]，Horgan等[6]认为休止角随着重力减小而增大，而Nakashima等[5]认为休止角不受重力大小的影响。另一方面，月面的特殊环境(低重力、高真空及昼夜温差大)促成了颗粒之间的微距力，从而在一定程度上提高了月壤的自稳能力[7-8]。削坡试验作为研究坡体稳定性中最直接最简单的方法，研究不同切坡倾角的削坡试验对理解月面环境下坡体滑塌机理具有一定的意义。然而，月壤的削坡试验的研究仍处于起步阶段，相关的研究成果较少，主要集中在试验研究和数值模拟。在试验研究方面，Zheng等[9]采用模拟月壤物研究了不同倾角下的月壤滑坡试验，研究发现月壤滑坡过程中出现裂隙面和光滑面，然而这种模拟月壤只是级配及材料相似，对于低重力场及表面微距力无法考虑。传统室内试验研究月壤具有成本高、还原度低、且对滑坡过程微观信息难以全面观察的缺点，而数值试验作为科学研究的一种重要手段，能有效克服以上缺点，Bui等[10]采用离散单元法，将颗粒之间的自锁能力引入到颗粒接触模型中来研究月壤坡体的滑塌，从微观机理上较好地解释了月壤坡体具有较高的稳定性。虽然国内外目前对月壤滑坡研究取得了一定的成果，但能够考虑颗粒间的抗转动作用及范德华力的月壤滑坡破坏的研究报告尚少见。

离散单元法能够很好地模拟颗粒材料大变形问题，深入土体的微观领域去研究其本质规律，故本文采用离散单元法进行模拟滑坡试验。借助PFC2D软件，通过植入二维月壤接触模型，模拟不同切坡倾角下月壤坡体滑坡试验，通过观察其滑动过程的变形特性和分析其滑坡机理，建立初始坡角与工程灾害指标的关系，为设计月面环境下放坡开挖及防范和治理月面环境下滑坡灾害提供参考。

2 月壤的微观模型介绍

本文所采用的微观模型是本团队最近提出来的，主要考虑月壤颗粒间抗转动作用及范德华力，模型的有效性可以参考文献[11-12],以下对模型进行简要的介绍。

(1)颗粒间的抗转动作用：月壤的形成与地壤的形成存在较大的差异，主要是陨石和微陨石的撞击、宇宙射线和太阳风粒子的持续轰击以及大幅度的昼夜温差使岩石热胀冷缩导致的破碎[1]，颗粒多属多棱角多气孔结构，表面凹凸不平，具有很强的啮合作用且具有一定的抗转动作用。

(2)颗粒间的范德华力：月壤的土性类似砂土，传统土力学认为其不应该具有黏聚力，但由于月球表面大气层极其稀薄，白天大气中分子含量仅为104个/cm2，属于超真空状态[1]，月壤颗粒表面吸附气体分子层厚度远小于地面土壤颗粒，在地面环境下无需考虑的范德华力在月面环境下应该给予考虑。

根据文献[11-12]，模型的力学响应如图2所示，其中Fv为颗粒之间的范德华力。法向力学响应主要由弹性力和范德华力两部分组成，其中线性函数的斜率表征颗粒之间的弹性接触，而其截距表征范德华力。弹性力会随粒间重叠量增加而增加，但范德华力不受重叠量改变的影响，恒为常数，当颗粒间没有接触时弹性力和范德华力均为零。切向力学响应采用增量法计算，最大值不超过滑动摩擦力μKnUn。转动力学响应也采用增量法计算，最大值不超过Knunβr/6。其中Kn为法向接触刚度，Ks为切向接触刚度，Km为扭转刚度，un为法向重叠量，us为剪切位移，θ为相对转动角度，μ为摩擦系数。

图2 月壤微观接触模型的力学响应[11-12]

3 月壤边坡模型建立

在采用离散单元法模拟月壤边坡滑坡过程试验中，若直接采用与真实月壤颗粒相同的级配，则颗粒粒径跨度过大，所需颗粒数目非常多，考虑到计算机能力及计算效率，对月壤级配做合理简化[11-12]，如图3所示。根据Mitchell等[13]测得真实月壤的黏聚力值范围为0.1～1.0 kPa、内摩擦角范围为30°～50°，通过模拟双轴试验对月壤宏观参数进行标定，微观参数如表1所示，得峰值黏聚力为0.2 kPa，峰值内摩擦角约为42.6°[11-12]，处于真实月壤的范围之内。

月壤滑坡数值模拟模型如图4所示，地基附近及下部边界采用3道无摩擦的刚性墙模拟，顶部采用自由边界模拟，与常见的离心机模型试验一致。

为了制备相对密实均匀的地基，采用Jiang等[14]提出的分层欠压法生成地基，生成的模型地基宽度为5 m，初始高度为2.5 m，如图4(a)所示。颗粒总数为1.5×105，平面孔隙比为0.20。为了模拟月面环境真实重力场，试样需要在(1/6)g重力场下完成固结平衡，固结完成后地基中水平应力及空隙比分布如图5所示。由图5可知固结平衡后的地基初始应力数值解和理论解基本相同，孔隙比分布均匀。而后需要对地基进行切坡来模拟不同初始坡角下的边坡，以边坡斜面中点为坐标原点，不断改变倾角的同时保证坡后颗粒总体积不变，删除A-B-C-D区域内的颗粒。图4(b)所示为切坡完成后45°的坡体，在坡肩处、坡面处及坡趾处各放置1个监测点。最后，在(1/6)g重力场下完成月壤滑坡系统的第二次平衡，观察并记录滑坡信息。

表1 离散元中的月壤微观参数

4 结果分析

4.1 滑坡形态分析

图6给出本文数值模拟及相关的室内试验[9]中坡度为90°时模拟月壤滑坡过程形态示意图。不同坡角下滑坡变形规律基本一致，此处仅以90°为特例做进一步介绍。无论数值模拟还是室内试验的滑坡过程都经历3个阶段：起滑阶段、滑动阶段及稳定阶段。起滑阶段，数值模拟的坡肩处出现明显的倒圆角，数值模拟中坡趾处轮廓线比室内试验大，这主要是因为室内试验采用的月壤黏聚力为1 kPa，比数值模拟0.2 kPa要大得多，此外，室内试验中模拟月壤仅干燥了12 h，可能残留的水分增大了黏聚力。滑动阶段，数值模拟和室内试验都出现叶片状滑体，90°坡体在滑动阶段出现前凸起上细下粗的叶片状，这主要是因为在叶片中心线上下两侧含有一定黏聚力的土体颗粒运动速度不一致产生变形差异所致。

稳定阶段，此时滑动区域的颗粒速度较小，主体变形基本完成，数值模拟和室内试验的滑动区域都从叶片状变为梭子状，但两者仍有区别，数值模拟滑坡稳定后坡面较光滑且坡度较小，而室内试验坡面不光滑坡度也较大，这种差异主要有两个原因：一是因为数值模拟采用的是简化级配，粒径跨度较小，形成坡面较光滑，而室内试验无需考虑计算效率问题，可以直接采用真实的月壤级配，其粒径跨度大，形成的坡面较粗糙；二是因为数值模拟和室内试验模拟月壤的黏聚力相差较大，黏聚力小的容易光滑，黏聚力大的容易形成裂纹。根据Bui等[10]对模拟月壤90°滑坡试验结果，出现黏性土坡的微裂缝并且最终坡度也较大，可归结于其模拟的月壤含有较高的黏聚力2 kPa。而此处模拟的月壤黏聚力0.2 kPa，不足以形成高坡角微裂缝。以上这些差异都体现了月壤因月面独特环境下而产生黏聚力，而这种黏聚力会直接影响破坏形态，决定着月壤偏向砂土还是偏向黏性土破坏。

图4 数值模型边界条件及切坡示意图(单位：m)

图5 1/6重力场下地基初始应力水平及孔隙比分布

图7给出数值模拟和室内试验下5种倾角保持坡后体积不变时，滑坡前后形态示意图。

从图7(a)中可知尽管5种坡体滑坡前形态不一致，但滑坡完成后的形态基本相似；无论在室内试验还是数值模拟中，月壤滑坡的过程均为土体颗粒围绕坡面中点旋转再一次达到平衡的过程，如图7(a)、7(b)所示。

图6 滑坡过程形态对比分析

4.2 滑坡过程流滑距离和最终倾角

为了减少滑坡带来的危害，对滑坡过程常见的指标进行预测是非常有必要的。根据Scheidgger[15]采用的滑坡体积、垂直滑落高度和水平滑移距离三者存在一定的联系来预测滑动距离的研究思路，本文通过控制坡后体积保持一定，简化分析滑坡距离与边坡倾角之间的关系，如图8(a)所示。流滑距离为坡趾处到滑坡舌的距离，为了便于理解及使用，将流滑距离和坡高的比值定义为归一化流滑距离Smax。随着边坡倾角的增大，归一化流滑距离也不断增大，这种增大趋势可以通过一条开口向下的抛物线进行拟合：

Smax=-0.01395α2+3.0455α-70.088

(1)

通过对坡高进行归一化处理，并且代入临界值为50%时，得到一条实用的结论：针对此类月壤，当坡角高于51.7°时，冲击到最远的土体的归一化值大于50%，即超过了边坡高度的一半；对于最危险的90°坡体其流滑距离也小于坡高。图8(b)给出了初始切坡倾角与坡体最终倾角的关系，由图8(b)可看出，初始坡角对最终稳定的坡度几乎没有影响，最终坡度在29°附近，以45°倾角时呈现微弱的峰值。根据文献[12]此类月壤在低围压的峰值内摩擦角为42.6°、残余内摩擦角为20.2°。5种倾角下滑坡失稳后的稳定边坡的角度都小于峰值内摩擦角，而大于残余内摩擦角，且基本稳定在休止角30.6°[16]。

图7 滑坡前后形态对比分析

4.3 滑坡过程最大速度分析

滑坡的速度直接决定滑坡时长，亦是人类对抗滑坡灾害反应的时间，如快速滑坡比缓慢滑坡更易造成对生命财产的破坏，因此研究滑坡过程的速度变化规律具有重要的工程意义。图9(a)给出了5种切坡倾角下最大速度Vmax随时间变化的规律。随着时间的增长，最大速度逐渐增长并且在0.2～0.8 s内相继达到峰值，之后随着时间的增长最大速度逐渐接近零。这与Laouafa等[17]从速度的角度上描述滑坡过程一致，认为滑坡的本质其实就是受到外力的触发由“固态”变为“液态”，最后受到外部阻力又从“液态”变为“固态”的过程。根据最大速度变化规律可以将不同坡度的滑坡分为类崩塌型和流滑型，类崩塌型的以60°、75°和90°为典型，出现较大时间跨度下的高速度场区，类似于突然间的崩塌；流滑型以30°和45°为典型，整个滑坡过程中不会有明显的特大速度，而是始终以稳定的低速度进行流动，更似颗粒的流动滑移。若将5种倾角在滑坡过程中出现的最大速度单独提取出来，会发现最大速度基本上随着切坡倾角呈现抛物线式的增大，如图9(b)所示，这条曲线也可以通过二次多项式进行拟合：

Vmax=0.0018α2-0.0338α+1.2175

(2)

但这种增大趋势与归一化流滑距离增长趋势不一样，最大速度增长速率随着坡度增大不断加剧，流滑距离增大速率随着坡度增大逐渐放缓。国际地质科学联合会[18]规定对于3 m/min～5 m/s之间的滑坡可以定位为快速滑坡,而大于5 m/s的可以认为是极速滑坡，而极速滑坡易造成极大的生命财产损失；将极速滑坡的下限值5 m/s代入拟合公式(2)得到α=56°，这表明当边坡坡角高度超过56°时，有很大的可能出现月壤的极速滑坡。所以从滑坡速度方面设计月壤坡体时，尽管月壤有较强的自稳能力，仍应该减少高陡坡体。

4.4 滑坡过程S-t曲线分析

斜坡从变形的产生到最终失稳破坏，累计位移-时间曲线(以下简称S-t曲线)，一般都会经历起滑阶段、滑动阶段及稳定阶段。前人根据S-t曲线斜率的特点将斜坡变形演化分成不同的阶段，并据此建立滑坡预警判据[19-20]，所以研究不同倾角下的斜坡S-t曲线具有重要的意义。本次模拟中，选取具有代表性的3个特征点动态追踪其累计位移，3个特征点分别布置在5种边坡的坡趾处、坡面中点处及坡肩处，共计15个，主要研究滑坡过程边坡倾角及监测点布置对S-t曲线的影响，如图10所示。

由图10可知所有的特征值点的S-t曲线都分为了起滑阶段、滑动阶段及稳定阶段。边坡倾角对S-t曲线的影响在不同布置区域呈现不同的规律：

图8 月壤滑坡的流滑距离和最终坡角的预测

图9 滑坡过程最大速度随时间和边坡倾角的变化曲线

图10 特征点累计位移曲线

(1)坡肩处(图10(a))，边坡倾角越大，则同一时间监测点的累计位移越大，表明坡肩处土体运动的距离和切坡倾角是正相关的。与此同时，随着切坡倾角的增大，坡体再次达到稳定阶段的时间也在逐步延长。

(2)坡面中点处(图10(b))，坡面中点累计位移基本上仍随着切坡倾角的增大而增加，但90°例外。对于这种反常现象，通过输出全过程的变形图，发现坡面中部特征点的位移受到坡顶土体垂直下落的影响，阻碍了这部分颗粒向前运动，因此累计位移减小。

(3)破趾处，中等坡度(如45°)的位移增长得最快，坡度过陡(如90°)或者过缓(如30°)的坡体其S-t曲线都会出现累计位移增长较慢。这主要是因为30°坡趾处由于坡面较缓，滑坡呈现轻微流动，对坡趾处监测点的土体提供的推力较小；90°坡趾处在滑坡过程中，坡肩处颗粒垂直落在坡趾，使得坡趾处的颗粒压实变密，坡后的土体不但没能提供坡趾处土体足够的推力，反而成了其阻力，这部分与坡面中点特征点极为相似，只是阻力作用更为显著。而介于这两个极端的坡体即中间坡度的坡体能够由坡后土体向前推动，使得其位移最大。综上所述，采用S-t曲线做滑坡预警时，需要注意监测点布置位置差异带来对监测数据的影响。若需要根据坡角信息特征及简单时间段估算滑坡过程处于何种阶段，最好将监测点布置在坡肩顶部，此时所得信息较为可靠。

5 结 论

本文通过引入考虑范德华力及抗转动作用的月壤微观接触模型，采用离散单元法模拟了月壤切坡试验，从多个角度分析了不同切坡倾角对滑坡灾害的影响，并且建立了滑坡过程中流滑距离、最终倾角及最大速度与切坡倾角的大致关系。得到主要结论如下：

(1)通过与试验对比，在考虑月壤之间的微距力的情况下，坡体具有一定的黏性，容易形成滑坡凸体及可能的裂缝，危害较大。初始切坡倾角对最终倾角几乎没有影响；但发现最大流滑距离随着初始坡角的增大而增大，当初始坡角超过51.7°时，最大流滑距离会超过坡高的一半。

(2)从工程危害角度分析了5种倾角下滑坡过程最大速度的规律，将其分为类崩塌型和流滑型，土体的最大速度随初始坡角的增大呈现抛物线式增长，通过抛物线拟合发现，当坡角低于56°时，可有效避免极速滑坡。

(3)通过对5种坡角的坡肩处、坡面中部及坡趾处共15个特征点的S-t曲线分析，发现S-t曲线都经历了起滑阶段、滑动阶段及稳定阶段，但只有当监测点布置在坡肩处时累计位移才会随着角度的增大而增大。若需要根据坡角信息特征及简单时间段估算滑坡过程处于何种阶段，最好将监测点布置在坡肩顶部。

参考文献：

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

[2] 欧阳自远. 我国月球探测的总体科学目标与发展战略[J]. 地球科学进展, 2004,19(3):351-358.

[3] 周增坡, 程维明, 万 丛,等.月球正面撞击坑的空间分布特征分析[J]. 地球信息科学学报, 2012,14(5):618-626.

[4] 乔 乐, 刘小倩, 赵健楠,等. 月球雨海地区三个着陆点的地质特征对比研究[J]. 中国科学(物理学 力学 天文学),2016,46(2):31-45.

[5] NAKASHIMA H, SHIOJI Y, KOBAYASHI T, et al. Determining the angle of repose of sand under low-gravity conditions using discrete element method [J]. Journal of Terramechanics, 2011,48(1):17-26.

[6] HORGAN B H N, BELL J F. Seasonally active slipface avalanches in the north polar sand sea of Mars: Evidence for a wind-related origin [J]. Geophysical Research Letters, 2015,39(9):9201-9207.

[7] CHANG C S, HICHER P Y. A constitutive model for granular materials with surface energy forces[J]. Journal of Aerospace Engineering, 2014,22(1):43-52.

[8] PERKO H A, NELSON J D, SADEH W Z. Surface cleanliness effect on lunar soil shear strength[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2001,127(4):371-383.

[9] ZHENG Hu, HUANG Yu. Model tests on flow slide of lunar regolith simulant[J]. Environmental Earth Sciences, 2014,73(8):1-7.

[10] BUI H H, KOBAYASHI T, FUKAGAWA R, et al. Numerical and experimental studies of gravity effect on the mechanism of lunar excavations[J]. Journal of Terramechanics, 2009,46(3):115-124.

[11] 郑 敏, 蒋明镜, 申志福.简化接触模型的月壤离散元数值分析[J].岩土力学, 2011,31(S1):766-771.

[12] JIANG Mingjing, SHEN Zhifu, THORNTON C. Microscopic contact model of lunar regolith for high efficiency discrete element analyses[J]. Computers and Geotechnics, 2013,54(10):104-116.

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

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

[15] SCHEIDEGGER A E. On the prediction of the reach and velocity of catastrophic landslides[J]. Rock Mechanics, 1973,5(4):231-236.

[16] 蒋明镜, 奚邦禄, 申志福,等. 不同重力下月壤水平推剪阻力离散元数值分析[J]. 岩土工程学报, 2015，37(7):1300-1306.

[17]LAOUAFA F, DARVE F. Modelling of slope failure by a material instability mechanism[J]. Computers and Geotechnics, 2002,29(4):301-325.

[18] International Union of Geological Sciences Working Group on Landslides. A suggested method for describing the rate of movement of a landslide[J]. Bulletin of the International Association of Engineering Geology, 1995,52(1):75-78.

[19] 许 强, 曾裕平, 钱江澎,等.一种改进的切线角及对应的滑坡预警判据[J]. 地质通报,2009,28(4):501-505.

[20] 许 强,汤明高,徐开祥,等.滑坡时空演化规律及预警预报研究[J].岩石力学与工程学报,2008,27(6):1104-1112.