不同重力场下月基承载特性离散元数值分析

奚邦禄，蒋明镜，莫品强，张振华，郭杨

(1. 安徽省建筑科学研究设计院 绿色建筑与装配式建造安徽省重点实验室，安徽 合肥，230031；2. 中国矿业大学 深部岩土力学与地下工程国家重点实验室，江苏 徐州，221008；3. 合肥工业大学 土木与水利工程学院，安徽 合肥，230009；4. 苏州科技大学 土木工程学院，江苏 苏州，215009；5. 同济大学 地下建筑与工程系，上海，200092)

随着嫦娥五号于2020年11月24日的成功发射及其于12月17日携带月壤样品成功着陆，我国探月工程“绕、落、回”三步走战略顺利完成。国务院公布的第五部航天白皮书《2021中国的航天》中包括月球科研站建设、近地小行星采样等一系列重大工程[1]。月球基地(月球科研站)是解决月球关键科学问题、聚焦全球优势科技资源、构建新型太空开发与治理体系的重要举措。美国于2017年联合日本、欧洲、加拿大等多个国家开展阿尔忒弥斯项目，计划于2024 年建成月球常态化驻留机制，并开展月球基地设计建设、资源原位利用的系统技术和设备研究。月球基地建设已成为我国乃至世界各国的重大战略需求。从岩土工作者角度出发，在月球基地的设计中，需要明确月面低重力环境对月壤承载特性的影响规律，进而为月球基地的科学设计与安全建设提供参考。

月壤是月球表面一层松软的散粒体状材料，是在无水、陨石撞击、宇宙射线、太阳风等环境下形成的，其颗粒大小、形态与地球上粉质砂土类似，具有黏聚力低、内摩擦角大的特性[2]。地面环境下月壤及其模拟材料的承载特性已有较多研究[2-6]。在研究低重力场对月壤承载特性的影响时，需要借助飞机飞行试验[7-8]、土工磁拟重力场模型装置[9]、落塔试验[10]或倾斜试验台[11]进行研究，成本较高。其中KOBAYASHI等[8]的飞机飞行试验结果表明Toyoura砂地基承载力随重力场的增大而增大，对于模拟月壤，当重力场小于1g时，地基承载力变化较小，当重力场大于1g时，地基承载力随重力场增大而增大。MO 等[9]基于土工磁拟重力场模型试验结果发现CUMT-1c 模拟月壤地基承载力随重力场增大而增大。中国空间技术研究院的落塔试验结果表明0g环境下的地基承载力远比1g环境下的地基承载力小[10]。也有学者利用离心机模型试验中模拟高重力场环境，分析了高重力场环境对月壤承载特性的影响。如李志刚[12]研究了1g、3g和6g下模拟月壤与人工装置的相互作用，发现不同重力场对同一尺寸的人工装置与模拟月壤的相互作用影响明显。LEE 等[13]开展了(1～60)g环境下模拟月壤承载特性的离心机模型试验，发现承载力和重力场的对数呈正比。综上所述，由于在地球实现低重力试验技术难度大、时间短，现代离心机技术已应用于高重力场环境对月壤承载特性的研究，但其成果能否推广至低重力场环境以及两者之间有何联系尚需进行进一步研究，以明确重力场对月壤承载特性的影响规律。

在理论分析方面，KOBAYASHI等[14]基于滑移线理论和极限分析法中上限定理[8]分析了低重力环境对月壤承载特性的影响，提出量纲一的参数G(重力场系数，与密度、重力场、黏聚力和承载板尺寸相关)，认为当G小于0.01 时，重力场对月壤承载特性的影响可以忽略；当G大于或等于0.01时，重力场对月壤承载特性的影响显著。LEE等[13]基于离心机试验结果提出了考虑重力场系数(G)和载荷板直径的承载力经验公式，但这些理论的适用性尚需进一步验证。也有学者采用数值模拟的方法研究重力场对月壤承载特性的影响规律。BUI等[7]采用离散单元法对模拟月壤的承载特性进行了研究，分析了重力场对模拟月壤地基承载力的影响规律，发现当重力场小于1g时，重力场对模拟月壤地基承载力无明显影响。然而，其模拟未从土力学角度进行系统分析，如未考虑边界效应影响、未检测地基初始应力等，且对月壤变形特性、破坏模式等分析较少。

离散单元法将土体作为刚性颗粒集合体，颗粒间可具有一定重叠量，可模拟土体的弹性变形。该方法可以很好地模拟土体的大变形与破坏问题，开展地基承载力分析，故本文采用离散单元法进行模拟研究，采用考虑范德华力和抗转动作用的月壤模型开展多种不同重力场下(1/100g、1/6g、1/2g、1g、2g、5g和10g)的承载力试验，分析重力场对月基承载力、月基弹性模量以及破坏模式的影响规律，以期为未来月球基地的设计提供参考。

1 月壤微观模型

本文所采用的微观模型考虑抗转动作用及范德华力的影响，模型的有效性验证见文献[15-16]，本文仅进行简要介绍。

1) 颗粒间抗转动作用。真实月壤颗粒形状不规则，以长条状、次棱角状和棱角状较为常见(图1)[17]，即月壤颗粒与颗粒间接触以面接触或多点接触为主，当接触点所受弹性变形不同时则产生弯矩，因此，颗粒间接触具有一定的抗转动能力，而不是球形颗粒间无抗转动能力的点接触。

图1 嫦娥5号月壤颗粒图[17]Fig. 1 Images of grains in lunar soil sample of Chang'e 5[17]

2) 颗粒间范德华力。NASA试验结果表明高真空环境下月壤力学特性变大[18]，这是因为月表高真空环境会导致月壤颗粒表面吸附气体分子层厚度远比地面土壤颗粒表面吸附气体分子层厚度小，需要考虑范德华力作用[19]，而地球环境下颗粒间范德华力相对于重力来说十分小，对宏观力学特性的影响很小，一般不需要考虑。数值模拟结果表明考虑范德华力使月壤试样表现出一定的黏聚力，与真实月壤的黏聚力更接近[15-16]。

2 离散元模型的建立

2.1 月壤静载荷试验模拟参数设置

真实月壤的颗粒粒径以小于1 mm为主(大多分布于30～100 μm)，颗粒粒径跨度较大[20]。然而，在离散元数值模拟中若采用与真实月壤相同的颗粒级配，所需的颗粒数目巨大，难以用当前普通计算机模拟。为提高计算效率，采用与JIANG等[21-22]研制的TJ-1 模拟月壤大颗粒分布一致的颗粒级配。图2所示为数值模拟试样级配曲线与真实月壤级配曲线对比。模拟材料的力学特性与级配、模型参数、颗粒形状及试样的孔隙比相关，因此，在模型建立过程中通过调整微观模型参数，使数值试样的宏观力学特性与真实月壤相近，即可用于月壤工程特性分析[23]，如月球车车轮行驶性能[24]、月壤开挖性能[25]等研究。通过参数反演法可以获得模拟月壤力学特性的相关参数[16]，如表1所示。由表1可知：月基中土样表现明显的软化现象，内摩擦角φ=38.8°，黏聚力c=0 kPa，本文着重介绍低重力环境对月壤承载特性的影响规律，因此，对具体调参过程不再进行介绍。

表1 离散元法模拟中土槽中月壤的参数Table 1 Parameters for lunar soil bin in DEM simulations

图2 数值模拟中试样与真实月壤级配曲线[15-16]Fig. 2 Grading curves of sample used in DEM simulations and true lunar soil[15-16]

2.2 静载荷试验模型建立

在静载荷试验离散元建模中，有必要充分考虑边界效应的影响，以免对低重力环境影响分析造成干扰[4,26]。若基于太沙基理论考虑边界效应，则会发现在深度方向上地基底部应力明显增大[26]，边界效应显著。基于附加应力传递理论[27]，当土体深度为条形基础宽度3～6倍时，附加应力为基础应力的0.1～0.2，因此，月基深度宜为基础宽度的6倍。在宽度方向上，基础与边界距离宜大于2.82倍[28]或3 倍[29]的基础宽度。故基础宽度取0.08 m(半模型试验中为0.04 m)，月基深度固结后为0.45 m(5.6倍基础宽度)，月基宽度为0.3 m(基础与边界距离为3.25倍基础宽度)，满足边界效应要求。

静载荷试验月壤土槽离散元模型如图3 所示。月基边界采用3 道无摩擦墙模拟，采用分层欠压法[30]在模型槽中生成月壤，然后，在不同重力场下固结。图4所示为不同重力场下月基固结完成后月基中竖向应力σv随月基深度h分布情况。由图3和图4 可知：在1g环境下，数值试样中σv随深度增加呈线性增加规律，σv=20.6h，与1g环境下应力水平理论值(σv=20.29h)较符合。不同重力场下应力分布形式相同，且在不同重力场下，月基中竖向应力与重力场呈正比，表明不同重力场下月基固结良好，可用于不同重力场下静载荷试验研究。图5所示为不同重力场下固结后月基的孔隙比分布情况，由图5可知：在不同重力场下，月基孔隙比均分布在0.228，上下波动相对误差在6%以内，表明可忽略不同重力场下月基密实度对月基承载特性的影响。

图3 月壤土槽离散元模型Fig. 3 DEM model of lunar soil bin

图4 不同重力场下初始月基应力水平Fig. 4 Initial stress conditions of lunar soil ground with different gravity fields

图5 不同重力场下初始月基孔隙比分布Fig. 5 Distributions of initial void ratios of lunar soil ground with different gravity fields

月基固结完成后，在月基左侧生成宽度为0.04 m 的载荷板(取半模型进行试验，以提升计算效率)，摩擦因数为0。采用位移加载[26]的方式进行试验，加载速度v恒定为0.000 5 m/s。

3 结果分析

3.1 荷载-沉降深度曲线

图6 所示为不同重力场下月基的荷载-沉降深度曲线。由图6可知：试验初始时荷载p随沉降深度增加而迅速增加，而后荷载增加速率逐渐变小，直到月基破坏；沉降深度s迅速增加，表现为典型的整体剪切破坏模式，且不同重力场下p-s曲线形式特类似；从承载力上看，极限承载力表现出明显的随重力场增大而增大的趋势，这与KOBAYASHI等[8]的飞机飞行试验结果(图6(b))类似。此外，我国计划的小行星探测需开展星壤取样[31]，尽管小行星重力场极低，接近0g(如Bennu 小行星表面加速度为20～100 μm/s2[32]))，但其地基仍具有一定承载能力，采用悬停着陆取样方案可行。

图6 不同重力场下月基荷载-沉降深度曲线Fig. 6 Load-settlement curves of lunar soil ground with different gravity fields

3.2 归一化极限承载力系数

以荷载趋近于稳定时承载力为极限承载力进行统计(沉降深度为8～12 mm)，如图7(a)所示，由图7(a)可以发现：月基承载力随重力场增大呈非线性增大，当重力场较小时，月基承载力增大较快；当重力场较大时，月基承载力增加较慢。以1g时月基承载力为基准，对不同重力场环境下月基承载力进行归一化得到承载力归一化系数nqu(不同重力场下承载力qnug与1g环境下承载力q1ug之比值)，如图7(b)所示，其中，R2为拟合度。由图7(b)可以发现：承载力归一化系数nqu随重力场增大呈非线性增长规律，但可以取1g环境为临界重力场，采用“二折线”可以较好地进行拟合。采用相同方法对MO 等[9]、KOBAYASHI 等[8]和LEE 等[13]的试验结果进行处理，试验中极限承载力归一化系数nqu与重力场系数ng关系曲线如图7(c)和图7(d)所示，其中，Dr为相对密实度，拟合参数见表2。由图7(c)、图7(d)和表2可以发现：MO等[9]低重力场试验和LEE 等[13]的高重力场试验结果均可采用线性公式较好地进行拟合，KOBAYASHI 等[8]的部分试验结果也可以较好地用“二折线”进行拟合，但部分低重力场试验结果可能由于试验误差和低重力环境不稳定等因素，拟合结果相对较差。总体来说，可采用以1g为临界重力场的“二折线”来考虑重力场对月基极限承载力的影响：

表2 模拟和试验中“二折线”拟合参数表Table 2 Fitting parameters of "binary line" in simulation and tests

图7 月基承载力-重力场关系Fig. 7 Relationships between bearing capacity and gravity field of lunar soil ground

式中：a1、a2、b1和b2均为拟合系数。

NASA 在计算月基承载力时，以承载质量(单位kg/m2)代替承载力(N/m2)来评估上覆荷载[33]，承载质量归一化系数与重力场ng的关系为

式中：a和b均为拟合系数。图8 所示为月基承载质量归一化系数与重力场ng的关系。由图8 可知：随重力场减小，重力场系数的倒数增大，承载质量归一化系数增加；当重力场小于1g时，月基可承载质量更大，而当重力场大于1g时，月基可承载质量更小。这是因为重力环境对构筑物自身重力的影响要比对月基承载力的影响大。此外，可采用以1g为界的“二折线”来拟合月基极限承载质量和重力场系数倒数间关系(式(2))。

图8 月基承载质量归一化系数-重力场关系Fig. 8 Relationship between bearing mass normalization coefficient and gravity field of lunar soil ground

3.3 月基弹性模量

由图6(a)可知：当沉降深度为1 mm 时，不同重力场环境下荷载位移曲线均处于线性增长阶段，可用于计算月基弹性模量，计算结果如图9(a)所示。由图9(a)可以发现：月基弹性模量也随重力场增大呈非线性增大，在低重力环境时增大较快，在高重力环境时增大较慢。同样以1g环境下月基弹性模量为基准，计算不同重力场环境下月基弹性模量归一化系数nE0(E0ng/E01g，E0ng为ng环境下月基弹性模量)，如图9(b)所示。由图9(b)可以发现：与月基极限承载力类似，弹性模量归一化系数与重力场系数关系也可以进行1g环境为临界重力场，采用“二折线”较好地拟合，且低重力环境下拟合直线斜率更高。

图9 月基弹性模量-重力场关系Fig. 9 Relationships between elastic modulus and gravity field of lunar soil ground

3.4 位移场分析

重力场对月基破坏模式的影响规律是建立考虑低重力环境承载力理论模型的关键因素，但在试验中较难观测到月基的破坏模式。图10 所示为不同重力场下月基中颗粒位移场图，黑线围成的区域表示根据位移变化梯度和方向确定的滑移区。由图10 可以发现：尽管不同重力场下月基内扰动范围存在显著不同，但颗粒位移方向规律较为类似，均呈现整体剪切破坏模式，即基础下方存在三角形刚性核区域(位移较大且位移方向垂直向下)，下移刚性核挤压两侧土体，使得地基破坏，形成贯通至月基表面的滑移线。此外，随重力场增大，三角形刚性核区域逐渐减小，侧向挤压区域宽度逐渐扩大。

图10 不同重力场下月基中颗粒位移场图(s=12.5 mm)Fig. 10 Particle displacement fields of particles in lunar soil ground with different gravity fields(s=12.5 mm)

月基中位移场在3/7s和4/7s梯度内变化较快，以此为代表绘制不同重力场下月基内受扰动区域，如图11 所示。由图11 可以发现：随重力场增大，月基中受扰动范围三角形刚性核区域在深度方向逐渐减小，侧向扰动区域在宽度方向逐渐增大。对于月基中位移大于4/7s和3/7s的颗粒数目，尽管满足该位移要求的颗粒并不一定被扰动，但选用同一标准进行对比，可以判断出重力场对受扰动范围大小的影响规律。由此可知，当重力场小于1g时，宽度方向扰动范围增加较快，受扰动颗粒数目逐渐增加；当重力场大于1g时，深度方向受扰动范围降低较快，受扰动颗粒数目逐渐减小。因此，当重力场小于1g时，地基初始应力和受扰动范围均随重力场增大而增大，月基承载力增大更快。

图11 月基内受扰动范围Fig. 11 Affected zone in lunar soil ground

3.5 力链分析

图12 所示为不同重力场环境下月基中力链传递模式。由图12 可知：不同重力环境下基础底部下存在明显的压力链，并向侧下方逐渐延伸(黑线为压力链主要延伸方向)，且随距基础距离增大，颗粒间力链逐渐变小，在低重力环境下和浅层尤为明显。这表明不同重力环境下地基底部应力向月基中传递模式较为一致，均表现为基础下方梯形传递，与FU等[26]模拟中整体剪切破坏模式下力链传递模式较为接近。

图12 月基内力链Fig. 12 Force chains in lunar soil ground

综上所述，在不同重力场下，p-s曲线、月基位移分布规律和力链传递模式表明不同重力场下月基破坏模式呈整体剪切破坏模式特征。随重力场增大，基础下方刚性核区域深度逐渐减小，侧向挤压区域逐渐扩大；受扰动颗粒数目在小于1g时随重力场增大而增加，大于1g时随重力场增大而减小。因此，在归一化承载力系数和月基弹性模量系数随重力场关系曲线中，可采用以1g为临界重力场的“二折线”进行拟合，低于1g时月基承载力和月基弹性模量增加更快。即采用离心机模拟高重力来预测低重力环境下月基承载力特性时不能正确反映低重力场时月基内受扰动范围的演化规律，将使月基的承载特性偏大。

4 结论

1) 当重力场小于1g时，随重力场增大，月基中受扰动范围和应力水平增大；当重力场大于1g时，随重力场增加，月基中初始应力增大，但受扰动范围逐渐减小，承载力和地基弹性模量随重力场增加呈现非线性增长趋势。

2) 归一化承载力系数(qnug/q1ug)和地基弹性模量系数(E0ng/E01g)与重力场系数ng关系可采用以1g为临界重力场的“二折线”进行拟合，低重力时，归一化承载力系数和地基弹性模量系数增大更快。

3) 随重力场减小，承载质量归一化系数增大，即低重力环境下，月基可承载质量更大，且可采用以1g为临界重力场的“二折线”来拟合月基极限归一化承载质量系数和重力场系数倒数之间的关系。

4) 采用离心机模拟高重力来预测低重力环境下月基承载力特性时，不能正确反映低重力场时月基内受扰动范围的演化规律，将使月基的承载特性偏大。