A Hypoplastic Constitutive Model Coupled with Temperature Effects for Lunar Soil Simulant

Abstract

In this paper, the QH-E（Qinghua-E sample）lunar soil simulant was taken as the research object, and a temperature dependent cohesive tensor was introduced based on the Gudehus-Bauer hypoplastic model framework to establish a hypoplastic constitutive model coupled with temperature effects. By comparing the model results with physical experiments, the effectiveness of the model was verified, achieving a unified characterization of the strain softening and compressive hardening response of lunar soil under deep ultra-low temperature environmental conditions on the lunar surface. The deformation characteristics such as shear dilation/contraction transition and critical state approaching were effectively simulated, reflecting the significant influence of temperature on the mechanical properties of lunar soil simulant. Compared with traditional elastoplastic models, the hypoplastic model coupled with temperature effects in this paper can effectively reveal the strength and deformation characteristics of lunar soil in the lunar environment, and better reflect the inherent influence of temperature on the strength and deformation evolution mechanism of lunar soil simulant. It provides theoretical support for engineering problems such as stability analysis of lunar landers and construction of lunar base.

Graphical abstract

引　言

月壤是覆盖月球表面的松散物质层，主要由陨石撞击破碎的月岩颗粒、太阳风的注入物质及宇宙尘埃等组成^[1]。月球探测器着陆，首先需研究月球探测器与月壤之间的相互作用，以设计探测器相关参数，但真实月壤样品非常少，为研究月壤的力学性能，顺利实现月面科学探测任务，故只能用性质相近（基本物理指标、颗粒级配、分布以及颗粒形状等）的模拟月壤来代替真实月壤^[2-^6]。合理的本构模型可准确描述月壤的力学性能，为月面钻探取样、着陆器与月球车设计及月球科研站建设提供理论支撑^[7]。

早期关于模拟月壤的力学本构关系，主要基于试验测试和数值模拟数据拟合建立起来的一些经验本构模型。Hasan等^[8]采用Particle Flow Code in 3 Dimensions（PFC 3D）软件模拟了美国约翰逊航天中心（Johnson Space Center，JSC）研制的JSC-1A型模拟月壤样品三轴压缩试验，并基于试验数据，建立了峰值内摩擦角、残余摩擦角、密度和平均法向应力的经验本构模型。Hicher等^[9]基于三轴压缩试验数据，建立了一个考虑月壤颗粒静电力引起的表面能效应本构模型，分析了表面能对月壤抗剪强度的影响。Richard^[10]基于三轴试验数据，使用剑桥模型描述了月球土壤本构行为，并通过调整内摩擦角，拟合了适用月壤的修正剑桥模型。Liu等^[11-^12]基于离散元模拟分析了月壤的高内摩擦角和低黏聚力特性，并利用弹簧张力模拟颗粒胶结作用，分别建立了无黏聚力和表观黏聚力的本构关系，分析了月壤在松散和密实条件下的抗剪强度指标与弹塑性行为。张璐^[13]基于离散元计算获得的模拟月壤三轴压缩试验数据，建立了模拟月壤的接触粘结本构模型。Hou等^[14] 基于月壤形貌特征数据，建立了考虑不规则形状的月壤颗粒扭转非线性接触模型，研究了模拟月壤的力学性能。

基于离散元数值模拟或试验数据拟合得到的经验弹塑性本构模型，虽然能较好地描述月壤的强度特征，但很少关注其变形特征，而传统弹塑性模型难以精确描述复杂应力状态的力学响应。亚塑性理论作为一种新兴的理论框架摒弃了传统弹塑性理论的塑性势、流动准则、硬化准则等概念，通过非线性张量方程直接描述材料的应力−应变关系，避免了弹塑性分离假设^[15]。Wu等^[16] 对不需要屈服面和分解变形的亚塑性理论进行了深入研究，提出第一个适用于砂土力学特性的简单实用亚塑性本构模型。Gudehus^[17]提出一种用于描述由孔隙比和应力张量充分表征颗粒材料力学性能的8参数亚塑性模型。Bauer^[18]进一步讨论了该亚塑性模型在广泛压力、变形和密度变化范围内材料的非线性和应变软化特征。Xu等^[19-^20] 考虑了冻土界面胶结强度和冰胶结作用对温度的敏感性，构建了用于描述冻土力学特性的扩展亚塑性模型。

月面采样和月球科研站建设的顺利与否主要取决于对月壤本构特性的认识程度，在本构模型中应合理反映月面环境月壤的强度和变形特征，才能准确描述模拟月壤的力学性能。月球表面环境的昼夜温差大，月壤力学行为不仅受颗粒间黏聚力和密实度的影响，而且还受到温度的显著影响。因此，为深入探究钻采设备与月壤的相互作用机制及月面特殊环境下月壤的力学行为，建立准确描述月壤力学特性的本构模型至关重要。

本文以QH-E模拟月壤为研究对象，该模拟月壤是清华大学以吉林靖宇县火山灰岩为原料，经雷蒙磨粉碎研制而成的低钛玄武质散粒材料，成分与“阿波罗”（Apollo）低钛月海样品接近^[21]。基于亚塑性模型的本构框架、考虑材料线性和非线性部分，研究模拟月壤在月面特殊环境的受力状态下表现出的非线性、非弹性行为，以及应变软化/硬化等关键力学特性。并考虑月面环境温度的影响，将温度相关的黏聚力张量引入本构方程，建立适用于考虑温度效应的模拟月壤的亚塑性本构模型，分析模拟月壤在月面深低温环境下的力学性能。

Gudehus-Bauer亚塑性模型

模型本构理论

Gudehus-Bauer（G-B）模型是亚塑性理论中引入孔隙比作为状态变量的高阶本构模型，核心假设包括^[17]：①材料基本性质假设，假设颗粒本身为刚性且不可压缩，在加载过程中不会发生破碎或磨损；②变形机制假设，变形由颗粒重分布主导，材料的宏观变形仅由颗粒的相对运动和重排列引起，不考虑颗粒本身的弹性或塑性变形；③简化孔隙流体影响，仅通过有效应力原理影响总应力，忽略其黏滞性或动态效应。

模型结合了临界状态理论，能统一描述不同密实度土的力学行为，通过因子化形式表征颗粒材料的压力与密度依赖性，实现复杂力学行为的特征描述。其本构方程可表示为^[17-^18]

(1)\dot{\sigma }=A(e,\sigma ,\dot{\epsilon })+B(e,\sigma )\Vert \dot{\epsilon }\Vert

其中：\dot{\sigma }为Jaumann应力率张量；A(e,\sigma ,\dot{\epsilon })为线性部分，与应变率\dot{\epsilon }呈线性关系；B(e,\sigma )\Vert \dot{\epsilon }\Vert为非线性部分，依赖应变率欧几里得范数\Vert \dot{\epsilon }\Vert =\sqrt{\mathrm{t}\mathrm{r}{\dot{\epsilon }}^2}，体现模型增量的非线性。

依据响应包络线将应力空间分解的定性特征可视化，即在已知某应力状态和孔隙比的基础上，由各个方向单位应变率增量引起所有应力率增量所覆盖的面。\Vert \dot{\epsilon }\Vert定义为1，则应力空间内轴对称状态下应变率的变化范围会使对应的应力率形成一个椭圆体，在{\sigma }_{11}-\sqrt{2}{\sigma }_{22}相应的截面得到一个椭圆响应包络线^[18]，如图1所示。在轴对称条件下，应力状态{\sigma }_{11},{\sigma }_{22}={\sigma }_{33}和同轴拉伸情况{\dot{\epsilon }}_{11},{\dot{\epsilon }}_{22}={\dot{\epsilon }}_{33}在主应力率空间形成的响应包络线，由于线性项的影响，应力率\dot{a}=A(e,\sigma ,D)轨迹会形成一个中心位于应力率空间原点的椭圆。对于特定方向相反但数值相同的应变率{\dot{\epsilon }}_a{,\dot{\epsilon }}_b会得到对应长度相等、方向相反的应力率{\dot{a}}_1{,\dot{a}}_2；线性部分由于方程是率无关的，故应力率\dot{b}=B(e,\sigma )\Vert \dot{\epsilon }\Vert仅拥有相同的数值，并不影响椭圆的形状和大小，仅仅根据相应值在响应包络线产生移动。{\dot{a}}_1 和\dot{b}共同作用导致本构方程产生的应力率\dot{t}=\dot{a}+\dot{b}会根据常数\dot{b}移动{\dot{a}}_1对应的响应包络线，得到位置偏移的应力轨迹，根据{\dot{\epsilon }}_a{,\dot{\epsilon }}_b影响的{\dot{t}}_1,{\dot{t}}_2拥有不同的长度和方向，恰好反映了增量非线性的特征，避免了重新定义准则从而也可确定加卸载差异。临界状态下的增量刚度会显示\dot{\sigma }=0也包含在亚塑性本构方程里，表现为应力率增量轨迹的起点即圆心会与响应包络线的某一点重合，使得线性和非线性响应大小相等但方向相反，即{\dot{a}}_1=-\dot{b}^[18]。

为合理反映散粒体材料加载时不同密度对力学响应的影响，分别考虑密度因子f_{\mathrm{e}}、f_{\mathrm{d}}作为相对孔隙比的无量纲标量参数，将线性项和非线性项分解。

(2)\{\begin{array}{rl}& A(e,\sigma ,\dot{\epsilon })=f_{\mathrm{e}}L(\sigma ,\dot{\epsilon })\\ & B(e,\sigma )=f_{\mathrm{e}}{f}_{\mathrm{d}}N(\sigma )\end{array}

其中：L(\sigma ,\dot{\epsilon })为受f_{\mathrm{e}}调控的线性项；N(\sigma )为受f_{\mathrm{e}}、f_{\mathrm{d}}调控的非线性项。

依据应力路径的相关性，当应力路径成比例时，所得到的应变路径也应成比例关系，即同等拉伸速率得到的应力率比应该等于应力比。

(3)\frac{\dot{\sigma }}{\begin{array}{c}\Vert \dot{\sigma }\Vert \end{array}}=\frac{\sigma }{\Vert \sigma \Vert }

当且仅当拉伸速率和f_{\mathrm{d}}恒定时可满足条件，且f_{\mathrm{e}}不影响应力路径的方向，假设本构方程在应力是齐次的。

(4)L(\lambda \sigma ,\dot{\epsilon })+f_{\mathrm{d}}N(\lambda \sigma )\Vert \dot{\epsilon }\Vert ={\lambda }^m[L(\sigma ,\dot{\epsilon }+f_{\mathrm{d}}N(\sigma )\Vert \dot{\epsilon }\Vert ]

其中：\lambda为一个任意正标量；m为齐次度，可进一步将两部分方程分解。

(5)\{\begin{array}{rl}& A\left(e,\sigma ,\dot{\epsilon }\right)={\left(\mathrm{t}\mathrm{r}\sigma \right)}^m{f}_{\mathrm{e}}L(\hat{\sigma },\dot{\epsilon })\\ & B\left(e,\sigma \right)={\left(\mathrm{t}\mathrm{r}\sigma \right)}^m{f}_{\mathrm{e}}{f}_{\mathrm{d}}N(\hat{\sigma })\end{array}

其中：\hat{\sigma }=\sigma /\mathrm{t}\mathrm{r}\sigma为无量纲颗粒应力比张量，根据以上非线性项和线性项齐次，响应包络线的展开与压力水平\mathrm{t}\mathrm{r}\sigma成正比。

将压力因子f_{\mathrm{b}}=f_{\mathrm{b}}(\mathrm{t}\mathrm{r}\sigma )及刚度因子f_{\mathrm{s}}=f_{\mathrm{b}}{f}_{\mathrm{e}}代入式（5）和（1），模型的本构方程可表示为

(6)\dot{\sigma }=f_{\mathrm{s}}\left[L(\hat{\sigma },\dot{\epsilon })+f_{\mathrm{d}}N(\hat{\sigma })\Vert \dot{\epsilon }\Vert \right]

刚度因子和密度因子对响应包络线的影响如图2所示，可以看出当刚度因子f_{\mathrm{s}}>1时，线性项会影响响应包络线，增加材料刚度，且该因子关于f_{\mathrm{b}}的关系能反映压力作用导致密度变化产生对刚度的影响。因此，f_{\mathrm{s}}能综合反映刚度对于压力与密度效应的影响。当密度因子f_{\mathrm{d}}>0时，则非线性项N(\hat{\sigma })会产生相应变化进而反映到应力空间响应包络线在该方向的偏移。因此，f_{\mathrm{d}}表征了密实度对压力相关剪切阻力的影响^[17]。

模型因子和参数的计算

密度因子f_{\mathrm{d}}是瞬时密度和平均压力的函数，且考虑临界状态的相关概念，即当e=e_{\mathrm{c}}时，f_{\mathrm{d}}=1。可根据当的情况下，峰值摩擦角的变化来描述密实度对抗剪强度的调控。e_{\mathrm{c}}为临界孔隙比，e_{\mathrm{d}}为下限孔隙比。根据f_{\mathrm{d}}和相对孔隙比r_e=(e-e_{\mathrm{d}})/(e_{\mathrm{c}}-e_{\mathrm{d}})(r_e⩾0)之间的线性关系，进行优化整合后可得到f_{\mathrm{d}}的非线性表达函数^[18]为

(7)f_{\mathrm{d}}={\left(\frac{e-e_{\mathrm{d}}}{e_{\mathrm{c}}-e_{\mathrm{d}}}\right)}^{\alpha }

其中：\alpha为无量纲的正常数，^[18]。

基于临界状态理论的密度因子效应分析可发现，当材料孔隙比小于当前围压临界值即当时，材料相对密实，在剪切过程中会产生剪胀效应并伴随着应变软化的特性；反之当孔隙比超出临界值f_{\mathrm{d}}>1时，材料呈现松散的结构状态，并产生剪缩变形伴随着应变硬化特性，而且无论初始结构状态如何，随着剪切变形的发展，材料均表现出向临界状态趋近的演化规律e=(1+e_0)\exp (-{\epsilon }_{\text{v}})-1。在临界状态阶段，呈现f_{\mathrm{d}}=1的稳定特征，此时材料的力学行为与初始孔隙比及应力历史无关，充分验证了f_{\mathrm{d}}作为剪胀/剪缩行为表征参数的有效性。当f_{\mathrm{d}}=0时，非线性项消失，仅体现线性部分特征，基本满足亚弹性本构关系。为了简化，同时将当时也设定为同种情况，避免超出孔隙比范围的扩展。理论上，当材料处于特定临界状态时，该模型可退化为经典的亚弹性模型。这种函数关系既保持了传统模型的物理基础，又通过引入密度因子驱动的非线性张量项实现了本构关系的拓展，增强了模型对复杂应力路径的适应性。 另一密度因子f_{\mathrm{e}}是当前孔隙比与临界孔隙比相关联的函数，可表示为^[18]

(8)f_{\mathrm{e}}={\left(\frac{e_{\mathrm{c}}}{e}\right)}^{\beta }

其中：\beta略大于1，反映孔隙比降低（密实度增加）导致的刚度增强。

由式（6），可将线性项和非线性项部分分别表示为

(9)f_{\mathrm{s}}\left[L(\hat{\sigma },\dot{\epsilon })\right]=f_{\mathrm{s}}[a_1^2\dot{\epsilon }+\hat{\sigma }\mathrm{t}\mathrm{r}(\hat{\sigma }\dot{\epsilon })]

(10)f_{\mathrm{s}}{f}_{\mathrm{d}}N(\hat{\sigma })\Vert \dot{\epsilon }\Vert =f_{\mathrm{s}}{f}_{\mathrm{d}}{a}_1(\hat{\sigma }+{\hat{\sigma }}^{\ast })\Vert \dot{\epsilon }\Vert

式（9）和（10）同规范化应力比张量\hat{\sigma }相关，{\hat{\sigma }}^{\ast }为\hat{\sigma }的偏应力张量{\hat{\sigma }}^{\ast }=\hat{\sigma }-1/3I，描述的是剪胀/剪缩特性，线性项部分为4阶张量，通过欧拉定理可得。a_1为基于临界摩擦角和应力状态定义的无量纲标量因子，可理解为临界状态下应力空间内偏平面临界应力点的半径，以此修正因Drucker-Prager圆形偏平面摩擦角差异引起的不足，使极限应力面更符合实际颗粒材料的破坏特征，通过罗德（Lode）角\theta引入应力方向依赖性；偏应力方向的函数可具体表达为^[17-^18]

(11)a_1^{-1}=c_1+c_2\Vert {\hat{\sigma }}^{\ast }\Vert [1+\cos (3\theta )]

其中：c_1、c_2分别为与临界摩擦角相关的常数，通过假设三轴压缩和拉伸的临界摩擦角相同，可确定它们的值。

(12)c_1=\sqrt{\frac{3}{8}\frac{3-\sin {\phi }_{\mathrm{c}}}{\sin {\phi }_{\mathrm{c}}}},c_2=\frac{3}{8}\frac{3+\sin {\phi }_{\mathrm{c}}}{\sin {\phi }_{\mathrm{c}}}

其中：{\phi }_{\mathrm{c}}为材料的临界摩擦角。

\Vert {\hat{\sigma }}^{\ast }\Vert为偏应力张量的范数，表征偏应力的强度，约束a_1的变化规律，避免应力路径突变符合渐进破坏，均匀变形下数值稳健；Lode角θ反映应力状态的各向异性，定义为

(13)\cos (3\theta )=-\sqrt{6}\frac{\mathrm{t}\mathrm{r}\left({\hat{\sigma }}^{\ast 3}\right)}{{\left[\mathrm{t}\mathrm{r}\left({\hat{\sigma }}^{\ast 2}\right)\right]}^{3/2}}

对于散粒体材料，其上限孔隙比e{}_{\mathrm{i}}为松散状态下孔隙比的最大值，下限孔隙比e_{\mathrm{d}}为密实状态下孔隙比的最小值以及临界孔隙比都会随平均压力p_{\mathrm{s}}=-\mathrm{t}\mathrm{r}\sigma /3的增大而减小，临界孔隙比与变形历史无关，且根据相应满足条件可推导路径的致密化行为表达式，即

(14)\frac{e_i}{e_{\mathrm{i}0}}=\frac{e_{\mathrm{d}}}{e_{\mathrm{d}0}}=\frac{e_{\mathrm{c}}}{e_{\mathrm{c}0}}=\exp \left[-{\left(\frac{3p_{\mathrm{s}}}{h_{\mathrm{s}}}\right)}^n\right]

其中：h_{\mathrm{s}}为颗粒硬度，表征材料抵抗压缩的能力，具有应力量纲，与颗粒强度相关；n为压缩指数，用来描述孔隙比随压力衰减的非线性程度。两个参数都可通过等向压缩试验的拟合上限孔隙比曲线分析确定。

压力因子f_{\mathrm{b}}基于等向压缩试验数据与理论公式的结合，在各向同性比例压缩下可以反映刚度与压力的非线性关系，表达式为

(15)f_{\mathrm{b}}=\frac{h_{\mathrm{s}}}{n{h}_i}{\left(\frac{3p_{\mathrm{s}}}{h_{\mathrm{s}}}\right)}^{1-n}\frac{1+e_i}{e_i}

其中：归一化参数h_i为确保压力因子在零下存在有限值，保证临界状态的连续性。

(16)h_i=\frac{1}{c_1^2}+\frac{1}{3}{\left(\frac{e_{\mathrm{i}0}-e_{\mathrm{d}0}}{e_{\mathrm{c}0}-e_{\mathrm{d}0}}\right)}^{\alpha }\frac{1}{c_1\sqrt{3}}

因此，刚度因子f_{\mathrm{s}}可表示为

(17)f_{\mathrm{s}}=f_{\mathrm{b}}{f}_{\mathrm{e}}=\frac{h_{\mathrm{s}}}{n{h}_i}{\left(\frac{e_i}{e}\right)}^{\beta }\frac{1+e_i}{e_i}{\left(\frac{3p_{\mathrm{s}}}{h_{\mathrm{s}}}\right)}^{1-n}

G-B亚塑性本构模型的8个参数，即：颗粒硬度h_{\mathrm{s}}、指数n、临界摩擦角{\phi }_{\text{c}}、指数\alpha\beta、下限孔隙比e_{\mathrm{d}0}、临界孔隙比e_{\text{c}0}、上限孔隙比e_{\mathrm{i}0}。

临界摩擦角{\phi }_{\text{c}}可通过直剪试验及常规三轴试验得到；压缩指数n和颗粒硬度h_{\mathrm{s}}反映了平均压力与散粒体材料的孔隙比之间关系的参数，通过等向压缩试验的上限孔隙比曲线确定，还可通过特定压力下体积模量、压缩指数和孔隙比计算推导^[17]为

(18)n=\frac{\ln (e_1{c}_{\mathrm{c}2}/e_2{c}_{\mathrm{c}1})}{\ln (p_2/p_1)}

(19)h_{\mathrm{s}}=3p{\left(\frac{ne}{c_{\mathrm{c}}}\right)}^{\frac{1}{n}}

其中：(p_1,e_1)、(p_2,e_2)为压缩曲线上不同压力点的值，对应压缩指数为c_{\mathrm{c}1},c_{\mathrm{c}2}；(p,e)为两点之间一压力点值，对应压缩指数c_{\mathrm{c}}。

3个与压力相关的孔隙比：初始上限孔隙比e_{\mathrm{i}0}为0压力条件下土体的最大孔隙比，可以想象为简单颗粒骨架的上限孔隙比，由压缩曲线上限孔隙比曲线外推得到。对于颗粒是均匀的球状体，有e_{\mathrm{i}0}/e_{\text{max}}\approx 1.2；对于颗粒为均匀的立方体，有e_{\mathrm{i}0}/e_{\text{max}}\approx 1.3。初始下限孔隙比e_{\mathrm{d}0}为0压力条件下土体的最小孔隙比，可用最小孔隙比e_{\text{min}}近似代替e_{\mathrm{d}0}\approx e_{\text{min}}。临界孔隙比e_{\text{c}0}也可通过普通的三轴压缩试验，在土体达到临界状态时得到，由于与最大孔隙比接近，可以替代为e_{\mathrm{c}0}\approx e_{max}^[22]。

指数\alpha、\beta两个材料参数作为刚度因子和密度因子的幂指数并没有具体的物理意义，由对应的公式计算得到^[18,^22]。

考虑温度效应的改进模型

基于G-B亚塑性本构模型引入考虑温度效应导致的温度相关的黏聚力张量，G-B模型临界应力面同样符合Drucker-Prager极限条件，主应力空间中是一个以原点为顶点的圆锥。将温度相关的黏聚力张量\mathit{S}作为背压应力引入G-B模型。

(20)\dot{\sigma }=A(e,\mathit{S},\sigma ,\dot{\epsilon })+B(e,\mathit{S},\sigma )\Vert \dot{\epsilon }\Vert

其中：\mathit{S}为考虑胶结作用的黏聚力二阶张量，\mathit{S}=s(T)\mathit{I}，类似于应力张量球面部分的各向同性张量，这里\mathit{I}为单位张量，s标量为温度相关的胶结作用或颗粒黏聚力大小，且忽略含水率对胶结作用的影响，仅考虑与温度相关的函数。

将温度相关的黏聚力张量\mathit{S}引入G-B模型后，可采用应力修正项将原模型的\sigma替换为\sigma -\mathit{S}，体现黏聚力对有效应力的贡献；应力比张量\hat{\sigma }=(\sigma -\mathit{S})/\mathrm{t}\mathrm{r}(\sigma -\mathit{S})，其中\mathrm{t}\mathrm{r}(\sigma -\mathit{S})作为分母项，使材料刚度特性直接与温度相关联。

黏聚力标量s(T)表达式依据温度和黏聚力的关系曲线拟合确定。比如在负温条件下，冻土的抗剪强度与温度大致可视为线性关系，提出黏聚力标量s(T)的表达形式为

(21)s(T)=aT+b

其中：T为温度；a和b为材料相关系数，根据负温下材料的黏聚力与温度线性拟合确认。

为确定式（21）的参数a和b，按照周兆曦等^[23]计算所得负温条件下模拟月壤内摩擦角和黏聚力相关的数值（表1），计算拟合出线的斜率和黏聚力轴线截距得到a=-83，b=-280，如图3 所示。

以上主要讨论了负温条件下s(T)的表达式，考虑到月表的温差范围–183～127 ℃，对于正温条件下s(T)表达式同样可依据正温条件下温度和黏聚力的关系曲线拟合确定。如果在整个温度（–183～127 ℃）范围内，s(T)可拟合成一个连续的直线，则用一个线性函数表示；如果在正负温度区间内不能拟合成一个连续的直线，则分别在正温和负温区间采用分段函数的形式。总之，在整个月面环境温度（–183～127 ℃），引入温度效应改进模型的表达形式不变，只是黏聚力标量s(T)表达式中参数a和b的数值发生变化。依据文献提供的不同负温条件模拟月壤黏聚力的相关试验数据^[23]，本文重点讨论负温条件下模拟月壤的力学行为。

模型计算结果与分析

模型验证

为验证模型的准确性，以QH-E模拟月壤为研究对象，依据先前的研究^[1,^6]假定模拟月壤的相对密度为83%，初始孔隙比为0.65。 模型采用在恒定围压下轴向加载的方式，此应力路径能有效模拟颗粒材料在不同围压条件下的力学行为。采用位移边界条件，通过固定试样位移，施加一定应变，观察材料的响应，这种边界条件可更好地控制加载过程，确保模拟结果的准确性。模型材料参数见表2。

常规围压试验与模型曲线的对比如图4所示，模型可较好地体现加载情况下的非弹性特性及应力增长非线性。体应变−轴应变曲线表明先出现少量的收缩，然后有明显的剪胀。其强度和变形特征均同试验结果基本一致。在常规围压下，随着围压增加，应力曲线的峰值强度明显提高，且峰后软化趋势不太明显，可能是压力因子受平均压力的影响，围压增加导致刚度增强，到达峰值后软化不明显。由体应变−轴应变曲线可知随着围压增加，剪缩特性也更加明显，曲线下降趋势更为平缓，围压越低，剪胀特性明显增强。通过对比模型和试验结果，验证了模型的准确性。

考虑温度效应的G-B模型

为将温度效应耦合到本构方程中，引入各向同性黏聚力张量\mathit{S}，实现模型在负温条件下对模拟月壤力学性能的描述。基于拟合的曲线s(T)=-83T-280，\mathit{S}=s(T)\mathit{I}，分别计算150 kPa围压下，在–5～–15 ℃以及–80～–120 ℃时，模型的偏应力−轴应变曲线和体应变−轴应变曲线如图5和图6所示。从不同负温水平下的偏应力−轴应变曲线可以看出，随着温度的降低，峰值应力明显增加，并且峰值应力对应的轴向应变也相应增加；偏应力曲线随温度降低逐渐变得更为陡峭，斜率增大。因为胶结作用随着温度下降而逐渐增大，模型的破坏线半径在π平面逐渐增大，符合低温胶结强化效应，胶结强度高，初始刚度增大。体应变−轴应变曲线随着温度的增加逐渐平缓，斜率增大，体积变化收缩阶段延长，符合胶结对于颗粒位移的抑制作用，体积变化也体现了更多压缩的特性，剪胀幅度减小。当温度在深低温水平（图6所示），体应变曲线更凸显剪缩特性。随着温度的降低，剪缩量增加，体积逐渐缩小。

模拟月壤在–10 ℃条件下，不同围压的偏应力−轴应变曲线如图7 所示。在相同温度条件下，随着围压的增大，偏应力曲线斜率增大，峰值偏应力增加，符合颗粒之间摩擦与胶结效应协同作用，材料更快密实，咬合力增强，导致峰值偏应力增大。体应变−轴应变曲线显示，体应变先剪缩后剪胀，随着围压的增加，剪缩量略有增加，但剪胀量降低，表明围压增加导致材料刚度增强，抑制体积膨胀变形。

模型主要参数的影响

G-B亚塑性本构模型8个参数中，指数\alpha和\beta材料参数作为刚度因子和密度因子的幂指数并没有具体的物理意义。n和h_{\mathrm{s}}是反映平均压力与散粒体材料孔隙比之间关系的参数，通过等向压缩试验的上限孔隙比曲线确定，即e_{\mathrm{d}0}、e_{\mathrm{i}0}、e_{\text{c}0}的影响可通过n和h_{\mathrm{s}}体现。因此，本文重点分析颗粒硬度h_{\mathrm{s}}、压缩指数n和临界摩擦角{\phi }_{\text{c}}的影响。

颗粒硬度h_{\mathrm{s}}对模拟月壤偏应力−轴应变和体应变−轴应变曲线的影响如图8所示。作为唯一具有应力量纲的参数，伴随着h_{\mathrm{s}}的逐渐增大，峰值偏应力增加，峰后软化量也随之增大。由于亚塑性模型在受力作用下弹性和塑性变形同时发生，并不需要区分弹性和塑性阶段，可以看出h_{\mathrm{s}}的增加直接导致曲线斜率增大；体应变曲线的斜率随着h_{\mathrm{s}}增大逐渐减小，体积变化的剪缩阶段也逐渐减小，更快进入剪胀阶段，并且剪胀量随h_{\mathrm{s}}的增加而增加。

压缩指数n的影响曲线如图9所示，可以看出压缩指数n对偏应力和体应变的影响也较大。随着n的增加，峰值偏应力逐渐增加，峰后软化程度也逐渐增大。因为n增大，平均压力增长更快，材料会更快变得密实，在较低应变就能达到更高的应力水平，导致峰值偏应力更高。对于体应变的变化，n增加使得剪缩阶段逐渐减小，加载初期体应变曲线更陡峭，后期体积变化逐渐平缓甚至体应变上升。n越大，体积变化的剪缩阶段逐渐减小，更快进入剪胀阶段，并且剪胀量随n的增加而增加。

偏应力和体应变随临界摩擦角{\phi }_{\text{c}}的变化关系如图10所示。峰值偏应力随{\phi }_{\text{c}}的增加而增加，相应到达峰值强度的轴向应变也逐渐增大。体应变−轴应变变化曲线显示模拟月壤先出现少量剪缩后剪胀。随着{\phi }_{\text{c}}增加，会延长体积收缩，延缓剪缩阶段，进入剪胀阶段后，随着{\phi }_{\text{c}}增加，剪胀量增加。

基于以上参数影响分析，可以看出模型的颗粒硬度h_{\mathrm{s}}、压缩指数n、临界摩擦角{\phi }_{\text{c}}对模拟月壤的强度和变形有较大影响。通过调整模型的相关参数，可改变模拟月壤试样的抗剪强度和变形特征。

结　论

本文基于Gudehus-Bauer模型框架引入温度相关的黏聚力张量，建立耦合温度效应的亚塑性本构模型，较好地描述了模拟月壤的强度和变形特征，并通过对比试验结果，验证了模型的准确性。模型实现了在不同温度条件下模拟月壤应变软化与压硬性响应的统一表征，有效模拟了剪胀/剪缩转换及临界状态趋近等复杂现象，比传统弹塑性模型更能揭示月面环境下月壤的强度和变形特征，以及温度对强度和变形演化机制的内在影响。通过对负温条件下模拟月壤力学性能的分析，发现低温导致的黏聚力增加，不仅显著增加了模拟月壤的峰值强度和初始刚度，还有效抑制了体积剪胀行为。

本研究忽略了含水率对模拟月壤力学行为的影响。由于孔隙冰的胶结作用会显著增强颗粒间的胶结强度，含水率的增加可能会提升月壤的抗剪强度。未来将进一步开展不同含水率和更大温度范围内（–183～127 ℃）耦合条件下，三轴剪切、直剪等试验，并在此基础上构建能够定量描述含水率−温度−力学响应之间耦合关系的本构模型，反映含水率–温度相关的月壤力学行为，为月球探测任务提供更准确的设计依据。

References

[1] Heiken G H, Vaniman D T, French, B M. Lunar sourcebook: a user’s guide to the Moon[M]. Cambridge: Cambridge University Press, 1991.

[2] Jia Y Shen Z & Dang Z L. Lunar soil simulant and its engineering application in lunar exploration program. Space Environment Engineering 31(3), 241–247 (2014).

[3] Long-fox J M Landsman Z & Easter P B. Geomechanical properties of lunar regolith simulants LHS-1 and LMS-1. Advances in Space Research 71(12), 5400–5412 (2023).

[4] Anbazhagan S Venugopal & Arivazhagan S. A lunar soil simulant（LSS-ISAC-1）for the lunar exploration programme of the Indian Space Research Organisation. Icarus 366(10), 1–15 (2021).

[5] Ruan R H Yang & Zhang D. A moderate-Ti lunar mare soil simulant: IGG-01. Acta Astronautica 224, 148–160 (2024).

[6] Zou W L Li Y & Chen L. Mechanical properties of QH-E lunar soil simulant at low confining stresses. Journal of Aerospace Engineering 29(2), 04015036 (2016).

[7] Li Y L, Wu W P, Wu H T, et al. Research progress on mechanical properties and constitutive models of lunar soil simulants[J]. China Power Science and Technology, 2025, 31（4）: 75-92.

[8] Hasan & Alshibli K A. Discrete element modeling of strength properties of Johnson Space Center（JSC-1A）lunar regolith simulant. Journal of Aerospace Engineering 23(3), 157–165 (2010).

[9] Hicher P Y, Chang C S. A constitutive model for lunar soil[C]//Proceedings of the 11th Biennial ASCE Aerospace Division International Conference on Engineering, Science, Construction, and Operations in Challenging Environments. Long Beach, California: ASCE, 2008.

[10] Richard A S. Cam cap models for lunar soil: a first look[C]//Proceedings of the 10th Biennial ASCE Aerospace Division International Conference on Engineering, Construction, and Operations in Challenging Environments. Houston, Texas: ASCE, 2006.

[11] Liu T X Wei & Liang L. Simulation and analysis of the lunar regolith sampling process based on the discrete element method. Transactions of the Japan Society for Aeronautical and Space Sciences 57(6), 309–316 (2014).

[12] Liu T X Liang & Zhao Y. An alterable constitutive law of high-accuracy DEM model of lunar soil. Advances in Space Research 66(6), 1286–1302 (2020).

[13] Zhang L. Lunar soil simulants’ constitutive relations based on PFC3D[D]. Beijing: China University of Geosciences（Beijing）, 2014.

[14] Hou X Y, Ding T X, Chen T, et al. Constitutive properties of irregularly shaped lunar soil simulant particles[J]. Powder Technology, 2019, 346: 137-149.

[15] Dafalias Y F Bounding surface plasticity. I: mathematical foundation and hypoplasticity. Journal of Engineering Mechanics 112(9), 966–987 (1986).

[16] Wu & Bauer E. A simple hypoplastic constitutive model for sand. International Journal for Numerical and Analytical Methods in Geomechanics 18(12), 833–862 (1994).

[17] Gudehus G A comprehensive constitutive equation for granular materials. Soils and Foundations 36(1), 1–12 (1996).

[18] Bauer E Calibration of a comprehensive hypoplastic model for granular materials. Soils and Foundations 36(1), 13–26 (1996).

[19] Xu G F Wu & Qi J L. An extended hypoplastic constitutive model for frozen sand. Soils and Foundations 56(4), 704–711 (2016).

[20] Xu G F Wu & Kong L W. Hypoplastic modeling for the mechanical behavior of frozen soil in stress path testing. International Journal of Geomechanics 18(6), 1–8 (2018).

[21] Yin C Chen & Dong Z F. Luanr soil simulant for drilling tool research. Geology and Exploration 48(1), 165–169 (2012).

[22] Xiong B L. The research of hypoplastic constitutive model of cohesionless soil[D]. Dalian: Dalian University of Technology, 2008.

[23] Zhou Z & Ma Q Y. Triaxial compression test and analysis of simulated lunar soil at negative temperature. Hefei University（Comprehensive Ed） 39(2), 26–31 (2022).

[24] Li Y L. Experimental test and discrete element modeling of triaxial compressive mechanical properties of QH-E lunar soil simulant under low stresses[D]. Wuhan: Wuhan University, 2018.

[25] Zou W L Chen & Zhang J F. Techniques for triaxial compression tests on simulant lunar soil QH-E and its mechanical behaviors under low confining stress. Chinese Journal of Geotechnical Engineering 37(8), 1418–1425 (2015).

Rights & permissions

The Author(s) 2026. Published by Higher Education Press. This is an Open Access article distributed under the terms of the CC BY license (https://creativecommons.org/licenses/by/4.0/).