各向异性颗粒材料的弹性模量：从松散状态到致密状态的微观结构依赖性

《Mechanics of Materials》：Elastic moduli of anisotropic granular materials: Microstructural dependence from loose to dense states

　　摘要：尽管颗粒介质的弹性在地质力学、力学和物理学领域早已被广泛研究，但由于无序堆积具有本质上的非仿射运动学特性，各向异性微观结构与宏观弹性响应之间的预测关系仍不完整。本文在有效介质框架下推导了横向各向同性颗粒材料的弹性模量解析表达式，并将弹性模量与填充分数、配位数以及结构各向异性等关键微观结构参数明确关联起来。这些表达式通过大量球形颗粒在三轴压缩下的离散元方法模拟得到了定量验证，涵盖了从松散到密实的初始状态。除在接近堵塞转变点的少数临界稳定状态下外，模拟得到的弹性模量都遵循由微观结构变量决定的统一标度关系。值得注意的是，所得模型系数与初始堆积状态无关，这表明微观结构与弹性之间的关系是由接触网络几何形状和颗粒属性决定的，而非制备历史。该框架进一步预测，在剪切临界状态下弹性模量会趋同于同一数值，为颗粒系统失去记忆效应提供了力学依据。此外，研究还表明体积模量在临界稳定状态下依然保持稳定，反映出其对非仿射重排及主导剪切响应的软模式具有较强免疫力。最后，通过对径向分布函数的分析发现，无论在松散还是密集堆积中都存在瞬时峰值，只不过前者的振幅较小，这说明这一特征并非随机密排结构的独有特征。

引言：颗粒介质是由沙子、谷物、粉末或颗粒等离散宏观粒子构成的集合体，在自然环境和工业过程中随处可见（de Gennes, 1999; Luding, 2016）。在自然界中，它们影响着从土壤力学和沉积物输运到滑坡、河床形成以及沙丘演化等众多现象，进而影响地貌演变和地质灾害（Deal et al., 2023; Andreotti et al., 2013; Kamrin et al., 2024; Vu, 2026b; Vo et al., 2020）。在工业领域，颗粒材料在农业、制药、采矿和化学工程中发挥着重要作用，需要大量进行搬运、储存、混合、研磨和运输（Vu, 2026a; Orozco et al., 2020; Gutiérrez et al., 2015）。根据外部条件不同，颗粒介质可表现出类似固体、液体或气体的行为，使其机械响应极为复杂（Vu et al., 2024b; MiDi, 2004; Jop et al., 2006; Kim and Kamrin, 2020）。这些丰富的现象使得对其研究既具有基础性挑战性，又对提升工艺效率及降低自然和工程系统中的故障风险具有重要的实际意义。

弹性响应是颗粒材料的重要性质，决定了这些系统在发生不可逆重排或失效之前如何承受和传递载荷（Farhang et al., 2017; Agnolin and Roux, 2007）。在此阶段，颗粒间的弹性接触力能够承受微小变形，这些接触力控制着刚度、体积模量和剪切模量以及波传播速度等关键宏观性能（Vu et al., 2023; La Ragione and Magnanimo, 2012b; Makse et al., 2004）。在地球物理领域，土壤和沉积物的弹性会影响地基稳定性、地震波传播以及滑坡或压实过程的早期发展（Cheng et al., 2020; Taeseo and W, 2015; Li et al., 2021）。在工业应用中，了解弹性响应有助于更准确地预测颗粒组装体对振动、压缩和操作的响应，从而改进存储系统、包装工艺和机械设备的设计（Kim et al., 2015; Laubie et al., 2017）。此外，弹性阶段还为微观接触力学与宏观材料行为之间建立了联系，为构建能够描述颗粒介质更复杂和塑性响应的本构模型奠定了基础（Kruyt and Rothenburg, 2002; Kruyt et al., 2010）。

弹性响应通常通过弹性模量来表征，对于各向同性材料而言，即为体积模量K和剪切模量G；而对于各向异性材料，则需用刚度张量来表示（Lubarda and Chen, 2008; Magnanimo et al., 2008; Li et al., 2022）。早期关于这些模量的理论预测是在有效介质理论框架下提出的，该理论将无序的接触网络替换为理想化的介质，其中每个颗粒都与平均环境相互作用。在这种方法中，颗粒位移被假设为遵循应变张量的仿射变换（即仿射假设）（Khalili et al., 2017; Kruyt, 2014; Vu et al., 2023）。然而，在无序系统中，单纯的仿射运动不足以维持机械平衡，因此会出现额外的非仿射位移。这些非仿射位移会改变经典弹性理论所预测的位移场，也是当采用仿射假设时有效介质理论中弹性模量被高估的原因（Makse et al., 1999; Kruyt and Rothenburg, 2002; Baggioli et al., 2022）。尽管存在这一局限，有效介质理论仍能成功捕捉到弹性模量与堆积微观结构之间的关键趋势（Kruyt, 2014; Vu et al., 2023; Agnolin and Kruyt, 2008）。不过，现有的有效介质理论公式大多仅适用于各向同性颗粒材料，而由于其相关变换的复杂性，将其扩展到各向异性系统的尝试仍然有限（Agnolin and Kruyt, 2008; Zaccone and Scossa-Romano, 2011）。

基于离散元方法的颗粒动力学模拟为测量颗粒材料的弹性模量提供了一种有效的替代方法。以往的研究大多集中在球形颗粒的各向同性密实堆积上，此时体积模量和剪切模量与填充分数、配位数、接触刚度以及围压等微观结构参数有关，尤其是当采用赫兹-明德林接触定律时（Agnolin and Roux, 2007; Velicky and Caroli, 2002; Otsubo et al., 2025）。最近，Khalili等人（2017）研究了通过等容压缩产生的各向异性球形颗粒堆积的弹性模量，并提出了弹性各向异性与围压之间的关系。其他研究也探讨了由三维球形颗粒或二维圆盘形颗粒构成的松散堆积的弹性行为（Khalili et al., 2017; Kruyt et al., 2010）。此外，Vu等人（2023）分析了通过三轴压缩产生的多面体颗粒各向异性堆积的弹性模量。他们首次提出了将弹性模量与填充分数、约束数（对于球形颗粒而言即配位数）以及结构各向异性相关联的解析表达式。然而，这些结果是基于密实堆积的模拟得出的。因此，为了得出更普遍的结论，有必要使用松散堆积的数据来验证这些表达式。

本文首先利用有效介质理论推导了横向各向同性材料的弹性模量表达式。随后，针对球形颗粒的松散和密实堆积进行了颗粒动力学模拟，通过三轴压缩生成各向异性结构。我们研究了从初始各向同性状态到剪切临界状态这一宽范围结构各向异性下的弹性模量。通过将有效介质理论的预测结果与模拟数据相互校准，我们重构出了Vu等人（2023）针对球形颗粒堆积提出的弹性模量解析表达式。随后，我们利用松散堆积的数据对该表达式进行了分析，从而为弹性模量与堆积微观结构之间的关系提供了新的见解。本文的其余部分结构如下：第2节详细推导了横向各向同性弹性模量的有效介质理论公式；第3节介绍了数值模型和模拟流程，同时分析了初始松散和密实球形颗粒堆积的微观结构以及三轴压缩过程中机械变量的变化情况，还介绍了用于测量弹性模量的应变探针，随后分析了弹性模量随剪切应变的演变规律；第4节研究了松散和密实堆积的弹性模量与微观结构之间的关系；最后，第5节总结了主要研究结果并进行了讨论。

理论框架：有效介质理论
当受到两侧应力相等的三轴加载时，材料会在垂直于加载轴的平面内呈现各向同性行为，而在加载方向上则呈现各向异性行为，这种对称性使得独立弹性常数的数量从完全各向异性情况下的21个减少到5个。在沃伊特标度下，相应的刚度矩阵可简化为以下形式（Fast et al., 1995; Lubarda and Chen, 2008）：
δσ11δσ22δσ33δσ

颗粒动力学模拟
在本节中，我们介绍了模拟中使用的颗粒动力学方法以及样品制备的数值流程。接着，我们介绍了用于施加剪切力并研究宏观和微观量随剪切应变变化的三轴压缩试验。此外，还讨论了用于计算三轴剪切过程中产生的各向异性结构的弹性模量的应变探针。最后，分析了……

堆积微观结构与弹性模量之间的关系
在第2节中，我们基于依赖仿射变形假设的有效介质理论，给出了弹性模量的解析表达式。然而，众所周知，有效介质理论无法准确预测颗粒材料的弹性模量。这一局限性源于非仿射位移，这类位移是对仿射运动的偏离，会导致弹性模量低于有效介质理论预测的值。在第3节中，我们研究了各向异性密实堆积材料的实际弹性模量……

结论
在本研究中，我们通过将弹性模量与不断变化的堆积微观结构明确关联起来，建立了一个用于预测各向异性颗粒材料弹性的微观力学框架。通过结合有效介质理论与离散元方法模拟，我们证明了横向各向同性颗粒堆积的弹性模量可以用填充分数、配位数以及结构各向异性等关键状态变量来表示。

利益冲突声明
作者声明不存在任何可能影响本文研究结果的已知财务利益或个人关系。

致谢
作者衷心感谢Farhang Radjai教授的深入指导，感谢Lhassan Amarsid博士在Rockable软件使用方面提供的帮助，同时也感谢Xavier Frank博士通过Matplotlib库提供的示例。

Duc Chung Vu