基于嫦娥五号月壤粒形特征的离散元模拟方法

王思远; 蒋明镜

doi:10.11779/CJGE20230040

基于嫦娥五号月壤粒形特征的离散元模拟方法 English Version

王思远^1,,
蒋明镜^{2, 3, 4, ,}

1.
天津大学建筑工程学院, 天津 300350
2.
苏州科技大学土木工程学院, 江苏苏州 215009
3.
同济大学土木工程减灾国家重点实验室, 上海 200092
4.
同济大学土木工程学院岩土工程系, 上海 200092

基金项目:

海南省重点研发计划项目 ZDYF2021SHFZ264

国家重大自然灾害防控与公共安全重点专项项目 2022YFC3003403

国家自然科学基金创新研究群体项目 42221002

国家自然科学基金重大项目 51890911

土木工程防灾国家重点实验室自主研究课题基金项目 SLDRCE19-A-06

详细信息

作者简介:
王思远(1995—)，男，博士研究生，主要从事太空土数值模拟及土体宏微观力学方面研究。E-mail: wsiy@tju.edu.cn

通讯作者:
蒋明镜, E-mail: mingjing.jiang@usts.edu.cn

中图分类号: TU43
计量
- 文章访问数: 675
- HTML全文浏览量: 74
- PDF下载量: 224
出版历程
- 收稿日期: 2023-01-12
- 网络出版日期: 2024-04-09
- 刊出日期: 2024-03-31

Lunar regolith simulations with discrete element method based on Chang'E-5 mission's lunar soil particle morphology

WANG Siyuan^1,,
JIANG Mingjing^{2, 3, 4, ,}

1.
School of Civil Engineering, Tianjin University, Tianjin 300350, China
2.
School of Civil Engineering, Suzhou University of Science and Technology, Suzhou 215009, China
3.
State Key Laboratory for Disaster Reduction in Civil Engineering, Tongji University, Shanghai 200092, China
4.
Department of Geotechnical Engineering, College of Civil Engineering, Tongji University, Shanghai 200092, China

摘要

摘要: 为满足当前中国探月任务中的“堪测、采、建”战略规划，获取更加贴近真实的月壤物理力学特性，考虑月壤颗粒形状的影响，结合嫦娥五号探月任务中的月壤颗粒图像，将月壤粒形特征与颗粒级配联系起来，提出考虑中国月壤粒形特征的离散元数值模拟方法。首先，基于嫦娥五号探月任务获取的部分月壤颗粒图像，对颗粒的形状特征及尺寸信息进行提取，依据其球度特征将月壤的颗粒形态分为6类，并同月壤颗粒的粒形特征与颗粒尺寸建立对应关系；其次，结合本团队开发的三维离散元月壤接触模型，引入形状参数用以表征月壤颗粒的抗转动能力，进而在离散元中考虑月壤颗粒的粒形特征，最终建立考虑中国月壤粒形特征的离散元数值模型；通过与已知月壤试验结果的对比表明，该方法可将月壤颗粒形状的可变性直接映射到离散元中。在此基础上，通过与不考虑粒形特征的数值样进行对比，探究考虑月壤粒形特征的优势。结果表明：含嫦娥五号月壤粒形特征的离散元数值模型可以有效捕捉月壤力学行为的主要特征，为开展月球资源勘探及开发利用研究提供基础。
- 月壤 /
- 粒形特征 /
- 三轴压缩 /
- 宏微观特性
Abstract: To align with the current strategic planning of "survey, mining, and development" in China's lunar exploration mission, and to obtain more accurate physical and mechanical properties of lunar regolith, this study focuses on the influences of lunar regolith particle shape based on particle images from the lunar regolith of the Chang'E-5 mission. A discrete element numerical simulation method that considers lunar regolith particle morphology is proposed by linking particle shape characteristics with gradation. Initially, the shape characteristics and size information of the particles are extracted from the lunar regolith images. The particles are subsequently categorized into six groups based on their sphericity, establishing the corresponding relationships. Secondly, the study utilizes a three-dimensional (3D) lunar regolith contact model and calculates rolling and twisting resistances at inter-particle contact by incorporating shape parameters to account for lunar particle shape effects. Subsequently, the model considers particle size characteristics within the discrete element analysis. Ultimately, a discrete element numerical model that incorporates the particle shape characteristics of the China's lunar regolith is developed. Comparison with the results of Apollo lunar regolith laboratory tests reveals that the variability of grain shape in lunar particles can be directly incorporated into the discrete model. Additionally, the benefits of considering the grain shape characteristics of lunar regolith are discussed in comparison to numerical samples that neglect the characteristics. The results show that the proposed method can effectively capture the main characteristics of the mechanical behavior of lunar regolith, and provide a basis for the lunar resource exploration and exploitation methods.
- lunar regolith /
- particle shape characteristic /
- triaxial compressed test /
- macro- and micro-mechanical property

HTML全文

0. 引言

饱和松砂在静力作用下由孔隙水压上升而诱发的抗剪强度在临界状态前就出现软化甚至趋向于零的现象称为静态液化。松砂地层中的工程施工易引起液化失稳现象（如海底长缓坡流滑^[1]和松砂地层中盾构开挖面失稳破坏^[2]）。静态液化失稳具有明显的整体性和扩散性,有别于岩土材料的其他失稳模式（如应变局部化）^[3]。以边坡静态液化失稳为例,液化一般开始于土体表面,随即因超孔隙水压力无法立即消散,导致液化区快速向四周扩散,故也被称为扩散性失稳^[4]。区别于应变局部化,静态液化时无明显滑动面,不易被察觉,破坏具有突发性。而且,静态液化一旦触发,土体强度会迅速丧失并可导致范围巨大的土体流滑,对邻近基础设施和人身安全构造重大威胁。因此,准确判别与预测静态液化具有重要意义。

Hill的二阶功理论是最早的材料稳定性通用判定准则之一,认为当材料处于稳定状态时其应力增量所对应的二阶功是非负的。大量学者研究了二阶功理论对判别砂土静态液化的适用性^[4-5],发现广义剪应力达到峰值时二阶功为零。然而,二阶功判别具有一定滞后性,因其值由正转负发生在静态液化触发后一增量步。因此,严格意义上二阶功判别无法预测静态液化的触发,但可用于识别已发生液化的区域。如黄茂松等^[4]将二阶功理论应用于边值问题,分析了海底长缓坡的砂土液化范围。吕玺琳等^[6]等基于二阶功理论推导了弹塑性刚度矩阵数学特征与砂土稳定性的关系,发现弹塑性刚度矩阵对称部分行列式为零发生在静态液化触发之前,从而将此条件定义为潜在失稳条件。Andrade^[7]和Lashkari^[8]基于材料失稳时应力应变增量响应解不唯一的条件推导了液化失稳触发时所对应的塑性模量的临界值,即当材料塑性模量等于该临界值时即判定静态液化触发。但上述临界塑性模量只适用于三轴应力空间。Buscarnera等^[9]基于Nova的控制理论^[10],得到了与加载控制条件相关的材料失稳临界塑性模量的广义表达式,并应用于三轴应力路径下饱和砂土的静态液化的判定^[11]。

室内单元体试验（平面应变试验^[12]及空心圆柱扭剪试验^[13]）表明,静态液化的触发及其液化程度除了受砂土初始密实度、应力水平的影响外,与中主应力比和主应力方向旋转密切相关。然而,针对上述现有判定准则的检验往往仅是在简单三轴应力路径下进行的,关于复杂应力路径下砂土静态液化判别和预测的研究工作较为有限。其中,吕玺琳等^[14-15]基于状态相关本构模型研究了 $K_{0}$ 固结和平面应变条件下的砂土静态液化失稳特性,发现松砂中弹塑性刚度矩阵对称部分等于零对应广义剪应力达到峰值,而随着砂土密实度的增加前者会先于后者出现。Lü等^[16]进一步研究了中主应力比和主应力方向对砂土液化特性的影响,发现中主应力比b（即 $b = (σ_{2} - σ_{3}) / (σ_{1} - σ_{3})$ ）和主应力方向与沉积方向的夹角 $α$ （在下文中简称为主应力旋转角）的增大会使二阶功极小值降低,液化更易发生。

针对上述研究现状,本文的目标是探究现有准则对涉及主应力旋转和不同中主应力水平的复杂应力路径下砂土静态液化判别的适用性。为此,结合陈洲泉等^[17]提出的砂土状态相关各向异性本构模型和笔者提出的判别准则,通过与空心圆柱扭剪单元体试验对比,验证该模型模拟复杂应力路径下静态液化现象的能力。基于此模型,本文比较了二阶功、弹塑性刚度矩阵对称部分行列式和失稳模量三种判别准则预测复杂应力路径下静态液化的准确性。最后,应用失稳模量方法获得了有效平均应力-广义剪应力空间的失稳线,分析了主应力旋转、中主应力水平等因素对静态液化触发应力状态和砂土可发挥最大摩擦角的影响。

1. 砂土状态相关各向异性本构模型

本文使用的本构模型以陈洲泉等^[17]提出的砂土状态相关各向异性模型为基础并加以改进,现将上述基础模型的核心部分和本工作提出的改进做简要介绍。

基础模型通过黄茂松等^[3]提出的各向异性状态变量来考虑砂土固有各向异性：

$A = \frac{\overset{⌢}{R}}{M_{cs} g (\overset{⌢}{θ})} - {(\frac{\overset{⌢}{R}}{M_{cs} g (\overset{⌢}{θ})})}_{c}$ ,

(1)

式中, $M_{cs}$ 为三轴压缩应力路径下的临界状态应力比, $\overset{⌢}{θ}$ 为组合应力张量的Lode角, $\overset{⌢}{R} = \sqrt{3 {\overset{⌢}{r}}_{i j} {\overset{⌢}{r}}_{i j} / 2}$ 为组合应力比偏量的第二不变量, $g (\overset{⌢}{θ})$ 为关于Lode角的椭圆插值函数^[18]：

$g (\overset{⌢}{θ}) = \frac{4 c}{1 + c + (1 - c) \sin (3 \overset{⌢}{θ})} - \frac{1 + c - (1 - c) \sin (3 \overset{⌢}{θ})}{2} 。$

(2)

组合应力张量采用Li等^[19]给出的定义：

${\overset{⌢}{T}}_{i j} = \frac{1}{6} ({\overset{⌢}{σ}}_{i k} F_{k j}^{- 1} + F_{i k}^{- 1} {\overset{⌢}{σ}}_{k j})$ ,

(3)

式中, ${\overset{⌢}{σ}}_{k j}$ 为归一化应力张量, $F_{i j}$ 为组构张量, ${\overset{⌢}{r}}_{i j} = {\overset{⌢}{s}}_{i j} / \overset{⌢}{p}$ 为组合应力比偏量, $\overset{⌢}{p} = {\overset{⌢}{T}}_{i i} / 3$ 为组合应力第一不变量, ${\overset{⌢}{s}}_{i j} = {\overset{⌢}{T}}_{i j} - \overset{⌢}{p} δ_{i j}$ 为组合偏应力张量。归一化应力张量 ${\overset{⌢}{σ}}_{i j}$ 定义为

${\overset{⌢}{σ}}_{i j} = \frac{M_{cs} g (θ)}{η} r_{i j} + δ_{i j}$ ,

(4)

式中, $δ_{i j}$ 为Kronecker张量, $θ$ 为应力张量的Lode角, $η$ 为应力张量对应的应力比,组构张量 $F_{i j}$ 在其主应力空间下可采用Oda等^[20]给出的定义：

$F_{i j} = \frac{1}{3 + Δ} [\begin{array}{l} 1 - Δ 0 0 \\ 0 1 + Δ 0 \\ 0 0 1 + Δ \end{array}]$ ,

(5)

式中, $Δ$ 为各向异性参量,不同制样方法下所对应的 $Δ$ 取值可参见文献[21]。式（1）中 ${(\frac{\overset{⌢}{R}}{M_{cs} g (\overset{⌢}{θ})})}_{c}$ 表示三轴压缩状态下的 $\frac{\overset{⌢}{R}}{M_{cs} g (\overset{⌢}{θ})}$ 值。

为考虑砂土临界状态的各向异性^[19],将临界状态线定义为各向异性状态参数A的函数：

$e_{c} = e_{Γ} (A) - λ {(\frac{p}{p_{a}})}^{ξ}$ ,

(6)

式中, $λ$ 和 $ξ$ 为材料参数, $p$ 为有效平均应力, $p_{a}$ 为标准大气压, $e_{Γ} (A)$ 表示 $p$ 为零时的临界状态孔隙比。根据Verdugo等^[22]和Yoshimine等^[13]通过试验测得的临界状态数据（见图1）,本文建议 $e_{Γ} (A)$ 采用以下表达式：

图 1 有效平均应力为零时的临界状态孔隙比随各向异性状态参数A的变化

Figure 1. Variation of critical state void ratio with anisotropic state-parameter A

下载: 全尺寸图片幻灯片

$e_{Γ} (A) = e_{Γ} (0) - k_{Γ} {| A |}^{x}$ ,

(7)

式中, $k_{Γ}$ 和 $x$ 为材料参数, $e_{Γ} (0)$ 为三轴压缩路径下（即 $A = 0$ ） $e_{Γ}$ 的取值。基础模型的剪胀和硬化法则均是各向异性状态变量A的函数^[17]：

$D = \frac{d ε_{v}^{p}}{d ε_{s}^{p}} = \frac{d (A)}{M_{cs}} (M_{d} - \frac{q}{p g (θ)})$ ,

(8)

$\frac{\partial M}{\partial ε_{s}^{p}} = \frac{h (A) G}{p} \frac{M_{p} - M}{M}$ 。

(9)

式中 $M_{d}$ 和 $M_{p}$ 分别为相变应力比和峰值应力比,两者均为关于Been等^[23]提出的状态参数ψ的函数； $h (A)$ , $d (A)$ 的具体表达式为

$d (A) = d_{0} (1 - k_{d} A)$ ,

(10)

$h (A) = (h_{1} - h_{2} e) (1 + k_{h} A)$ ,

(11)

式中, $d_{0}$ , $k_{d}$ , $h_{1}$ , $h_{2}$ 和 $k_{h}$ 为模型参数。

2. 静态液化判定准则

根据Hill^[24]给出的稳定性条件,当材料发生失稳时二阶功需满足

$d^{2} w = d σ_{i j} d ε_{i j} \leq 0, \exists ‖ d ε ‖ \neq 0$ ,

(12)

式中, $d σ_{i j}$ 和 $d ε_{i j}$ 分别为应力和应变增量。对于饱和砂处于不排水条件时,式（12）可用以判定静态液化是否已经发生。然而,当砂土为非饱和或其排水条件为部分排水时,上述条件仅是静态液化触发的必要条件^[25]。

2.1 弹塑性刚度矩阵对称部分行列式

吕玺琳等^{[6, 14-15]}提出了二阶功准则需要两个附加条件,以全面判别不同饱和度和排水状态下砂土的静态液化：

$\partial Q / \partial σ_{i j} δ_{i j} > 0$ ,

(13a)

$\partial F / \partial σ_{i j} δ_{i j} < 0$ ,

(13b)

式中, $F$ 为屈服函数, $Q$ 为塑性势函数。上述两个附加条件分别表征了静态液化发生时土体需要处于剪缩状态以及应力状态处屈服面需朝静水压力轴正方向张开。由于二阶功由正转负发生在静态液化触发后一增量步,滞后于液化点,因此并不能严格地称其为液化判别准则。

基于二阶功失稳理论,吕玺琳等^{[6, 14-15]}进一步提出了可以根据弹塑性刚度矩阵的数学特性预判静态液化的触发。将应力应变关系代入式（12）可得

$d^{2} w = d ε_{i j} D_{i j k l}^{ep} d ε_{k l}$ ,

(14)

式中,弹塑性刚度矩阵 $D_{i j k l}^{ep}$ 可分解成对称部分 ${[D_{i j k l}^{ep}]}_{s}$ 和反对称部分 ${[D_{i j k l}^{ep}]}_{s k}$ 之和。又因为应变增量 $d ε_{i j}$ 为对称二阶张量,其与刚度矩阵反对称部分 ${[D_{i j k l}^{ep}]}_{s k}$ 所构成的二次型 $d ε_{i j} {[D_{i j k l}^{ep}]}_{s k} d ε_{k l}$ 恒等于零,因此二阶功可写为

$d^{2} w = d ε_{i j} {[D_{i j k l}^{ep}]}_{s} d ε_{k l}$ 。

(15)

式（15）表明当弹塑性刚度矩阵对称部分为半正定矩阵时,则必然能找到一条应变路径（其方向为 ${[D_{i j k l}^{ep}]}_{s}$ 的一个特征向量）使得二阶功等于零。进一步,由于半正定矩阵必为非满秩矩阵,故可将行列式

$\det ({[D_{i j k l}^{ep}]}_{s}) = 0$

(16)

作为砂土进入潜在失稳状态的判定条件^[6]。

笔者已经成功地应用式（16）对饱和砂土在三轴应力路径下^{[6, 15]}和平面应变路径下^[14]以及高饱和含气砂在三轴应力路径^[26]下的静态液化现象进行了研究。本文将进一步分析上述判定准则对复杂应力路径下静态液化失稳的预测效果。

2.2 失稳模量

当利用弹塑性刚度矩阵特性进行静态液化失稳判定时,失稳条件的具体表达式不依赖特定的应力-应变路径或加载控制条件,因此该方法具有良好的通用性。另一方面,这种通适的判定准则因没有考虑静态液化现象所蕴含的力学特征（如材料无体变）而造成对失稳触发的过于严判^[6]。另一类判定理论则显式地考虑静态液化触发时的应力-应变控制条件,进而获得了静态液化触发时弹塑性本构关系中塑性模量的临界值,即当塑性模量K_p达到这一临界值时,即判定静态液化失稳。

Andrade^[27]从应力应变增量响应解不唯一出发,推导了三轴应力路径下静态液化触发时塑性模量的临界值：

$H_{χ} = - \frac{\partial F}{\partial p} K \frac{\partial Q}{\partial p}$ ,

(17)

式中, $H_{χ}$ 为临界塑性模量, $K$ 为体变模量。当土体的塑性模量达到上述临界值时,在不排水条件下对应于任意广义剪应力增量 $d q$ 的解不唯一或不存在,即判定静态液化触发。值得注意的是,式（17）成立同时对应于Nova^[10]给出的三轴应力下不排水剪切加载路径可控性丧失的条件。

同样从加载路径是否可控的角度,Buscarnera等^[9]推导了复杂加载控制条件下岩土材料发生失稳时塑性模量所需满足的临界值的广义形式：

$H_{χ} = - {(\frac{\partial F}{\partial σ_{β}})}^{T} (D_{β β}^{e} - D_{β α}^{e} {(D_{α α}^{e})}^{- 1} D_{α β}^{e}) \frac{\partial Q}{\partial σ_{β}}$ 。

(18)

式中 F为屈服函数；Q为塑性势函数； $D_{i j}^{e}$ 为弹性刚度矩阵的子矩阵（i,j= $α$ , $β$ ）, $D_{α β}^{e}$ 联系应变增量 $d ε_{β}$ 与应力增量 $d σ_{α}$ , $D_{β α}^{e}$ 联系应变增量 $d ε_{α}$ 与应力增量 $d σ_{β}$ ,其中 $σ_{α}$ 和 $ε_{β}$ 作为加载时所控制的应力和应变分量, $σ_{β}$ 和 $ε_{α}$ 为对应的响应应力和应变分量。当塑性模量等于式（18）给出的临界值时,则在给定的加载控制路径下塑性乘子的取值不唯一或不存在。当以体应变为零和广义剪应力（或某一特定剪应力分量）单调递增作为控制条件时,式（18）具体化为静态液化触发时塑性模量的临界值。三轴条件下,以广义剪应力 $q$ 和体应变 $ε_{v}$ 为控制变量,对应的广义剪应变 $ε_{q}$ 和有效平均应力p为响应变量,相应的失稳模量就退化为式（17）的形式。对于常中主应力比和常主应力旋转角扭剪试验,其所控制的剪应力分量 $τ_{z θ}$ 与q同时达到峰值,因此可以继续使用式（17）判定静态液化的触发。对单剪路径判定静态液化时,由于 $b$ 不是常数,Lode角变化的影响不能忽略,因此需要直接以剪应力 $τ_{z θ}$ 和体应变为控制变量,此时失稳模量可表示为

$H_{χ} = - \frac{\partial F}{\partial σ_{i j}} D_{i j k l}^{e} \frac{\partial Q}{\partial σ_{k l}} + \frac{\partial F}{\partial τ_{z θ}} 2 G \frac{\partial Q}{\partial τ_{z θ}}$ ,

(19)

式中,G为弹性剪切模量, $D_{i j k l}^{e}$ 为弹性刚度矩阵。

下文中,笔者将具体分析基于塑性模量临界值的判定方法对复杂路径下静态液化失稳的预测效果,以及其与笔者提出的基于弹塑性刚度矩阵对称部分行列式进行判定的异同。

3. Toyoura砂静态液化模拟及判别

本文使用Toyoura砂单元体试验的数据来验证本构模型重现复杂应力路径下砂土静态液化的能力。为此需要标定两类材料参数。前一类描述水平沉积砂土在三轴压缩路径下的基础力学特性（表1中左侧四列）,而后一类则控制主应力方向偏离沉积方向时以及中主应力比变化时土体力学特征的改变（表1中第五列）。对于前一类参数,本文采用了Li等^[28]基于Verdugo等^[22]进行的三轴试验所标定的参数。Toyoura砂各向异性材料参数则基于Yoshimine等^[13]进行的空心圆柱扭剪试验获得（如图1所示）。具体标定过程可参考文献[17,28]。

表 1 Toyoura砂模型参数

Table 1. Model parameters of Toyoura sand

弹性参数	临界状态参数	剪胀参数	硬化参数	各向异性参数
G₀=125, $ν$ =0.05	M_cs=1.25, e_г=0.934, λ_c=0.019, ξ=0.7	d₀=0.88, m=3.5	h₁=3.15, h₂=3.05, n=1.1	Δ=0.214, K_г=0.13, k_d=1.8, k_h=0.9, c=0.75, x=0.5,

下载: 导出CSV

| 显示表格

3.1 常规不排水三轴压缩试验

不排水三轴压缩试验下,实测与模拟的Toyoura砂的有效应力路径和应力-应变关系如图2（a）和图2（b）所示。从图2中可见,模拟曲线与试验曲线吻合较好。初始有效平均应力 $p_{0}$ 较大的两组试样均出现静态液化现象,但 $p_{0}$ 为100 kPa的试样的不排水强度单调上升。其原因是降低初始围压导致砂土初始状态在e-p平面上向临界状态线靠近,抑制了液化现象发生。图2（a）,（b）还表明当塑性模量等于失稳模量时（即 $K_{p} = H_{χ}$ ）,广义剪应力达到峰值,即静态液化触发。而且,静态液化触发后下一增量步对应的二阶功非正。另一方面,弹塑性刚度矩阵对称部分行列式等于零要早于（见 $p_{0}$ 为1000 kPa的试样）或同于（参见 $p_{0}$ 为2000 kPa的试样）静态液化触发。当液化现象越弱时（峰值强度后的强度软化越弱时）,两者的差别趋于增大。对 $p_{0}$ 为100 kPa的试样,弹塑性刚度矩阵方法判定存在潜在失稳的可能性。即吕玺琳等^[6]定义的潜在失稳区。这意味着在当前应力状态下存在一条潜在的应变路径使得土体丧失稳定性,但这条路径并非试验控制条件下的路径。

图 2 不排水三轴压缩路径下试验与模拟结果对比

Figure 2. Comparison between measured and simulated responses in undrained triaxial compression tests

下载: 全尺寸图片幻灯片

3.2 定中主应力比或定主应力方向的不排水空心圆柱试验

不排水定主应力方向剪切试验的实测与模拟结果如图3所示。不排水定中主应力比剪切试验的实测与模拟结果如图4所示。上述对比表明所使用的本构模型可以较好反映主应力方向和中主应力比对砂土静态液化的影响,具体表现为随着主应力方向逐渐偏离沉积方向和中主应力比的增大,静态液化更易被触发。从不同判定准则的应用效果看,失稳模量方法可以预测广义剪应力峰值（即静态液化触发点）,且失稳点后一增量步的二阶功非正。然而,弹塑性刚度矩阵对称部分行列式取值为零或同于（α=60°,b=0.5,e=0.828；α=45°,b=0.5,e=0.855和α=45°,b=0.75,e=0.849的试样）或早于（α=75°,b=0.5,e=0.823；α=45°,b=0,e=0.853和α=45°,b=1,e=0.861的试样）广义剪应力峰值点,且预判早于实际失稳的倾向随着液化现象的减弱而加强。另外,对于未出现静态液化的三组试验（α=15°,b=0.5,e=0.825；α=30°,b=0.5,e=0.824和α=45°,b=0.5,e=0.821）,行列式方法也判定出有失稳的可能。

图 3 不排水定主应力方向剪切试验的试验与模拟结果对比

Figure 3. Comparison between measured and simulated responses in undrained shear tests with fixed directions of principal stresses

下载: 全尺寸图片幻灯片

图 4 不排水定中主应力比剪切试验的试验与模拟结果对比

Figure 4. Comparison between measured and simulated stress paths in undrained shear tests with fixed intermediate principal stress ratios

下载: 全尺寸图片幻灯片

3.3 空心圆柱单剪试验

图5,6分别表示本构模型针对等向固结和 $K_{0}$ 固结（ $K_{0}$ =0.5）不排水单剪试验的模拟和试验结果。不同判定方法对静态液化触发预测的效果与之前分析的剪切路径大致相同。需要特别注意的是在单剪试验中 $σ_{1} - σ_{3}$ 的峰值与判别的失稳点不再重合,这种不重合现象在 $K_{0}$ 固结单剪中更为明显。这是因为单剪路径所控制的是剪应力分量 $τ_{z θ}$ （定义见图5）,因此该路径下静态液化触发的标志是 $τ_{z θ}$ 达到峰值,而非最大主应力差 $σ_{1} - σ_{3}$ 。后者是剪应力 $τ_{z θ}$ 和中主应力比b的函数,因而只有当 $b$ 固定时, $τ_{z θ}$ 与 $σ_{1} - σ_{3}$ 才同时取到峰值。在等向固结单剪试验中 $τ_{z θ}$ 达到峰值之前中主应力比基本稳定,所以静态液化失稳点与 $σ_{1} - σ_{3}$ 峰值近似重合。但在 $K_{0}$ 固结单剪试验中 $τ_{z θ}$ 在加载初期即达到峰值,此时中主应力比还在快速下降,因此判定失稳点与 $σ_{1} - σ_{3}$ 峰值有显著差异。

图 5 等向固结单剪试验与模拟结果

Figure 5. Comparison between measured and simulated stress paths in undrained hollow cylindrical simple shear tests following isotropic consolidation

下载: 全尺寸图片幻灯片

图 6 　K₀固结单剪试验与模拟结果（K₀=0.5）

Figure 6. Comparison between measured and simulated stress paths in undrained hollow cylindrical simple shear tests following K₀ consolidation

下载: 全尺寸图片幻灯片

上述结果表明,复杂应力路径下二阶功和失稳模量方法判定的静态液化触发对应于由加载控制条件所决定的某一特定剪应力分量达到峰值。这是因为二阶功和失稳模量的取值都是控制条件相关的,即二阶功理论内含了当前的应变路径,而失稳模量理论则是基于特定的控制条件推导得到的。相比之下,弹塑性刚度矩阵对称部分行列式判定结果则与控制条件无关,其考量的不是特定加载路径,而是是否存在任意一种加载路径使土体丧失稳定性。

4. 静态液化失稳线

失稳线作为p-q应力空间内静态液化触发所对应应力状态的集合,定义了饱和砂土（特别是松砂）在不排水状态下可发挥抗剪强度的上限。松砂剪缩特性引起的不排水剪切下孔压的急剧上升,使砂土在发挥出临界状态摩擦角之前便出现不排水强度下降,即静态液化触发（见图2（a））。因此,基于失稳线的稳定性分析较传统的基于摩擦角的稳定性分析可更准确地预判实际工程中岩土构筑物的失稳与破坏^[25]。由第3节的分析可知失稳模量理论因考虑了静态液化现象与加载控制条件的相关性,可较为准确地预测复杂路径下静态液化的触发。本文应用失稳模量理论得到了有效平均应力-广义剪应力空间内的静态液化失稳线,并归纳了砂土密实度、有效平均应力、中主应力比及主应力方向对静态液化失稳线的影响。

静态液化失稳线是通过令塑性模量（可通过式（9）联合屈服面一致性条件得到）等于失稳模量（即式（17））得到。如图7所示,失稳线在有效平均应力小于一定值时是不存在的（如 $e$ 等于0.88的例子）,即因砂土初始状态接近临界状态抑制了静态液化的发生。另外,随着密实度的下降、 $b$ 值或 $α$ 角度的增大,失稳线会发生下移,且失稳线下限边界所对应的有效平均应力降低,导致液化在低围压条件下也可发生。除较难液化的砂的失稳线在起始段有较明显的弯曲之外,失稳线可被近似考虑为直线。因而,失稳线的斜率,即砂土发生液化失稳前可发挥的最大摩擦角,从工程应用角度可考虑为不受有效围压影响,而仅随密实度降低、主应力旋转及中主应力比增大而降低。

图 7 有效平均应力-广义剪应力空间内静态液化失稳线

Figure 7. Instability lines in space of effective mean normal stress and deviatoric stress

下载: 全尺寸图片幻灯片

5. 结论

基于状态相关各向异性本构模型,本文研究了涉及主应力方向旋转和中主应力比变化的复杂应力路径下的砂土静态液化,重点讨论了现有判定准则对复杂三维应力路径下静态液化的预测效果,并分析了密实度、围压、主应力方向和中主应力等因素与静态液化触发时应力状态的相关性,得到以下4点结论。

（1）要真实模拟复杂应力路径下的静态液化现象,需要在本构模型中考虑砂土状态相关特性和固有各向异性。

（2）基于二阶功和失稳模量的静态液化判别都显式地考虑了加载过程的控制条件,因此在复杂应力路径下失稳模量理论判别的液化失稳点与二阶功给出的失稳位置相同,并与土体真实失稳点一致。

（3）弹塑性刚度矩阵对称部分行列式判定方法与控制条件无关,其考量的不是当前的加载路径,而是是否存在一种加载路径触发失稳。故其判别结果同于或早于真实静态液化失稳点。

（4）失稳时已发挥的最大土体摩擦角主要受材料密实度、中主应力比和主应力方向影响,而对有效围压不敏感,材料密实度越小、中主应力比越大或主应力方向旋转角越大,液化时土体已发挥的摩擦角越小。

致谢: 感谢中国科学院地质与地球物理研究所杨蔚研究员对本文提供的巨大帮助！感谢团队硕士生管成良、李天赐、石宵宵、杨越群在本文中提供的帮助。

图 1 部分月球表面颗粒^[22]

Figure 1. Multiple lunar regolith particles^[22]

下载: 全尺寸图片幻灯片

图 2 基于图像识别月壤粒形特征捕获流程图

Figure 2. Implementation of lunar soil particle shape feature capture based on image processing

下载: 全尺寸图片幻灯片

图 3 嫦娥五号月球样品的颗粒形状和尺寸分布

Figure 3. Particle shape and size distributions of Chang'E lunar sample

下载: 全尺寸图片幻灯片

图 4 月壤颗粒球度和抗转动的散点分布

Figure 4. Scatter distribution of sphericity and rolling resistances of lunar soil particles

下载: 全尺寸图片幻灯片

图 5 月壤粒间接触示意图

Figure 5. Schematic diagram of lunar soil particle contact

下载: 全尺寸图片幻灯片

图 6 考虑抗转动影响的粒间演化规律

Figure 6. Contact evolution between particles considering resistance rotation coefficient

下载: 全尺寸图片幻灯片

图 7 DEM模拟颗粒级配曲线及抗转动参数

Figure 7. DEM simulation of grain grading curve and corresponding resistance rotation coefficient

下载: 全尺寸图片幻灯片

图 8 Apollo 12/16月壤颗粒形态特性

Figure 8. Morphological characteristics of soil particles in Apollo 12/16 mission

下载: 全尺寸图片幻灯片

图 9 数值模拟与室内试验对比

Figure 9. Comparison between numerical simulations and laboratory tests

下载: 全尺寸图片幻灯片

图 10 不同围压下试样三轴试验结果

Figure 10. Results of conventional triaxial tests under different confining pressures

下载: 全尺寸图片幻灯片

图 11 剪切过程中力学配位数及APR演化规律

Figure 11. Evolution of mechanical coordination number and porosity during shearing

下载: 全尺寸图片幻灯片

表 1 月壤的颗粒形态分布表

Table 1 Representative particle morphologies of lunar particles

球度范围	颗粒占比/%	代表性颗粒形态
${S_{\text{p}}}$ < 0.6	1.3
0.6 < ${S_{\text{p}}}$ < 0.7	5.8
0.7 < ${S_{\text{p}}}$ < 0.8	18.1
0.8 < ${S_{\text{p}}}$ < 0.9	39.27
0.9 < ${S_{\text{p}}}$ < 0.94	18.3
0.94 < ${S_{\text{p}}}$ < 1.0	17.18

下载: 导出CSV

表 2 离散元模拟参数^[21]

Table 2 Parameters used in DEM simulations

接触类型	细观参数		取值
球-球		颗粒密度/(kg·m^-3)	2320
		有效模量E_p/MPa	800
		法向与切向刚度比 $\zeta$	5
		抗转动系数 $\beta$	见图 7
		局部破碎系数	4
		颗粒摩擦系数	0.5
		局部阻尼	0.7

下载: 导出CSV

参考文献(41)

[1]	COSTES N C, CARRIER W D, MITCHELL J K, et al. Apollo 11 soil mechanics investigation[J]. Science, 1970, 167(3918): 739-741. doi: 10.1126/science.167.3918.739
[2]	SCOTT R F, CARRIER W D, COSTES N C, et al. Apollo 12 soil mechanics investingation[J]. Géotechnique, 1971, 21(1): 1-14. doi: 10.1680/geot.1971.21.1.1
[3]	JAFFE L D. Shear strength of lunar soil from Oceanus Procellarum[J]. The Moon, 1973, 8(1): 58-72.
[4]	SCOTT R F. Failure[J]. Géotechnique, 1987, 37(4): 423-466. doi: 10.1680/geot.1987.37.4.423
[5]	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. doi: 10.1061/(ASCE)1090-0241(2001)127:4(371)
[6]	CARRIER III W D. Particle size distribution of lunar soil[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2003, 129(10): 956-959. doi: 10.1061/(ASCE)1090-0241(2003)129:10(956)
[7]	WILLMAN B M, BOLES W W, MCKAY D S, et al. Properties of lunar soil simulant JSC-1[J]. Journal of Aerospace Engineering, 1995, 8(2): 77-87. doi: 10.1061/(ASCE)0893-1321(1995)8:2(77)
[8]	KLOSKY J L, STURE S, KO H, et al. Geotechnical behavior of jsc-1 lunar soil simulant[J]. Journal of Aerospace Engineering, 2000, 13(4): 133-138. doi: 10.1061/(ASCE)0893-1321(2000)13:4(133)
[9]	ALSHIBLI K A, HASAN A. Strength properties of JSC-1A lunar regolith simulant[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2009, 135(5): 673-679. doi: 10.1061/(ASCE)GT.1943-5606.0000068
[10]	KANAMORI H, UDAGAWA S, YOSHIDA T, et al. Properties of Lunar Soil Simulant Manufactured in Japan[M]. SPACE 98, 1998: 462-468.
[11]	JIANG M J, LI L Q, SUN Y G. Properties of TJ-1 lunar soil simulant[J]. Journal of Aerospace Engineering, 2012, 25(3): 463-469. doi: 10.1061/(ASCE)AS.1943-5525.0000129
[12]	ZHENG Y C, WANG S J, OUYANG Z Y, et al. CAS-1 lunar soil simulant[J]. Advances In Space Research, 2009, 43(3): 448-454. doi: 10.1016/j.asr.2008.07.006
[13]	李建桥, 邹猛, 贾阳, 等. 用于月面车辆力学试验的模拟月壤研究[J]. 岩土力学, 2008, 29(6): 1557-1561. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200806025.htm LI Jianqiao, ZOU Meng, JIA Yang, et al. Lunar soil simulant for vehicle-terramechanics research in labtory[J]. Rock and Soil Mechanics, 2008, 29(6): 1557-1561. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX200806025.htm
[14]	ZHOU S, YANG Z, ZHANG R, et al. Preparation and evaluation of geopolymer based on BH-2 lunar regolith simulant under lunar surface temperature and vacuum condition[J]. Acta Astronautica, 2021, 189: 90-98. doi: 10.1016/j.actaastro.2021.08.039
[15]	LI R, ZHOU G, YAN K, et al. Preparation and characterization of a specialized lunar regolith simulant for use in lunar low gravity simulation[J]. International Journal of Mining Science and Technology, 2022, 32(1): 1-15. doi: 10.1016/j.ijmst.2021.09.003
[16]	HUANG Y, LU X, ZHAO R, et al. Three dimensional simulation of lunar dust levitation under the effect of simulated sphere body[J]. Journal of Terramechanics, 2011, 48(4): 297-306. doi: 10.1016/j.jterra.2011.06.005
[17]	SOLTANBEIGI B, PODLOZHNYUK A, PAPANICOLOPULOS S. et al. DEM study of mechanical characteristics of multi-spherical and superquadric particles at micro and macro scales[J]. Powder Technology, 2018(329): 288-303.
[18]	KATAGIRI J, MATSUSHIMA T, YAMADA Y, et al. Investigation of 3d grain shape characteristics of lunar soil retrieved in Apollo 16 using image-based discrete-element modeling[J]. Journal of Aerospace Engineering, 2014, 28(4): 4014092.
[19]	JIANG M J, SHEN Z F, WANG J. A novel three-dimensional contact model for granulates incorporating rolling and twisting resistances[J]. Computers and Geotechnics, 2015, 65: 147-163. doi: 10.1016/j.compgeo.2014.12.011
[20]	RORATO R, ARROYO M, GENS A, et al. Image-based calibration of rolling resistance in discrete element models of sand[J]. Computers and Geotechnics, 2021, 131: 103929. doi: 10.1016/j.compgeo.2020.103929
[21]	JIANG M J, SHEN Z F, THORNTON C. Microscopic contact model of lunar regolith for high efficiency discrete element analyses[J]. Computers and Geotechnics, 2013, 54: 104-116. doi: 10.1016/j.compgeo.2013.07.006
[22]	月球与深空探测科学数据与样品发布系统[EB/OL]. http://202.106.152.98:8081/moondata/web/datainfo/main.action#,2007-10-24/2023-01-02. Lunar and Deep Space Exploration Scientific Data and Sample Release Syslem[EB/OL]. http://202.106.152.98:8081/moondata/web/datainfo/main.action#,2007-10-24/2023-01-02. (in Chinese)
[23]	YANG W, WANG Y, GAO L, et al. Sci-tech arts on Chang'e-5 lunar soil[J]. Innovation (Camb), 2022, 3(5): 100300.
[24]	RORATO R, ARROYO M, GENS A, et al. Particle shape distribution effects on the triaxial response of sands: a DEM study[C]//Micro to MACRO Mathematical Modelling in Soil Mechanics, Cham, 2018.
[25]	RORATO R, ARROYO M, ANDÒ E, et al. Sphericity measures of sand grains[J]. Engineering Geology, 2019, 254: 43-53. doi: 10.1016/j.enggeo.2019.04.006
[26]	LI Q, ZHOU Q, LIU Y, et al. Two-billion-year-old volcanism on the Moon from Chang'E-5 basalts[J]. Nature, 2021, 600(7887): 54-58. doi: 10.1038/s41586-021-04100-2
[27]	HE Q, LI Y, BAZIOTIS I, et al. Detailed petrogenesis of the unsampled Oceanus Procellarum: the case of the Chang'E-5 mare basalts[J]. Icarus, 2022, 383: 115082. doi: 10.1016/j.icarus.2022.115082
[28]	YANG Y, JIANG T, LIU Y, et al. A micro mid-infrared spectroscopic study of Chang'E-5 sample[J]. Journal of Geophysical Research: Planets, 2022, 127(8): e2022J-e7453J.
[29]	LI J, LI Q, ZHAO L, et al. Rapid screening of Zr-containing particles from Chang'E-5 lunar soil samples for isotope geochronology: Technical roadmap for future study[J]. Geoscience Frontiers, 2022, 13(3): 101367. doi: 10.1016/j.gsf.2022.101367
[30]	JIA L, CHEN Y, MAO Q, et al. Simultaneous in-situ determination of major, trace elements and Fe 3+ /∑Fe in spinel using EPMA[J]. Atomic Spectroscopy, 2022, 43(1): 42-52.
[31]	ZHANG D, SU B, CHEN Y, et al. Titanium in olivine reveals low-Ti origin of the Chang'E-5 lunar basalts[J]. Lithos, 2022(414/415): 106639.
[32]	LI C, LI Y, WEI K X. Study on surface characteristics of Chang'E-5 fine grained lunar soil[J]. Scientia SInica Physica, Mechanica & Astronomica, 2022: 50: 1-10.
[33]	JI J, HE H, HU S, et al. Magmatic chlorine isotope fractionation recorded in apatite from Chang'e-5 basalts[J]. Earth and Planetary Science Letters, 2022(591): 117636.
[34]	LI C, HU H, YANG M F, et al. Characteristics of the lunar samples returned by the Chang'E-5 mission[J]. Natl Sci Rev, 2022, 9(2): 1-13.
[35]	ZHANG H, ZHANG X, ZHANG G, et al. Size, morphology, and composition of lunar samples returned by Chang'E-5 mission[J]. Science China (Physics, Mechanics & Astronomy), 2022, 65(2): 106-113.
[36]	BARDET J P. Observations on the effects of particle rotations on the failure of idealized granular materials[J]. Mechanics of materials, 1994, 18(2): 159-182. doi: 10.1016/0167-6636(94)00006-9
[37]	蒋明镜. 现代土力学研究的新视野: 宏微观土力学[J]. 岩土工程学报, 2019, 41(2): 195-254. doi: 10.11779/CJGE201902001 JIANG Mingjing. New paradigm for modern soil mechanics: Geomechanics from micro to macro[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(2): 195-254. (in Chinese) doi: 10.11779/CJGE201902001
[38]	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. doi: 10.1016/S0266-352X(03)00064-8
[39]	HEYWOOD H. Particle size and shape distribution for lunar fines sample 12057, 72[C]//Proceedings of the Lunar Science Conference, Houston, 1971.
[40]	JIANG M J, YU H, HARRIS D. Kinematic variables bridging discrete and continuum granular mechanics[J]. Mechanics Research Communications, 2006, 33(5): 651-666. doi: 10.1016/j.mechrescom.2005.06.013
[41]	THORNTON C. Numerical simulations of deviatoric shear deformation of granular media[J]. Géotechnique, 2000, 50(1): 43-53. doi: 10.1680/geot.2000.50.1.43

施引文献(17)

期刊类型引用(10)

1.	钟小春，易斌斌，竺维彬，祝思然，罗淑仪，王奇. 可液化土层地震液化引起盾构隧道上浮位移模式. 长安大学学报(自然科学版). 2024(01): 80-90 . 百度学术
2.	姚方正，刘伯成，慈伟，杨尚川，李睿. 考虑土体空间变异性的可液化土层圆形隧道地震响应分析. 兰州理工大学学报. 2023(04): 117-121 . 百度学术
3.	张益瑄. 液化区盾构隧道位置对上浮效应影响研究. 山西建筑. 2022(04): 144-146 . 百度学术
4.	申玉生，雷龙，闵鹏，朱双燕，甘雨航，王彬光. 基于降饱和度原理强震作用盾构隧道抗液化机理研究. 地震工程与工程振动. 2022(04): 17-24 . 百度学术
5.	王彦臻，赵丁凤，陈国兴，梁珂. 一维场地地震反应非线性有效应力分析法及其验证. 岩土工程学报. 2021(03): 502-510 . 本站查看
6.	张稳军，曹文振. 地震下大断面盾构隧道接缝力学及防水性能研究. 岩土工程学报. 2021(04): 653-660 . 本站查看
7.	邱明兵，高文生，秋仁东. 考虑惯性和运动相互作用桩基础抗震计算实例. 建筑科学. 2021(03): 39-45 . 百度学术
8.	程选生，周欣海，王平，柴少峰，冯欢，杨莉萍. 黄土隧道结构的振动台模型试验研究. 中国公路学报. 2021(06): 136-146 . 百度学术
9.	罗驰恒. 超大直径过江盾构隧道口字件不同连接型式地震响应特征研究. 公路交通科技. 2021(06): 39-45 . 百度学术
10.	吴悦，霍永鹏，巫志农，孔令威，李河山，张益瑄. 砂土液化对越江盾构隧道上浮行为的影响研究. 交通科技. 2021(06): 87-91 . 百度学术