用于SPH中颗粒流模拟的应力-应变一致双应力公式

《Powder Technology》：A stress–strain consistent dual-stress formulation for granular flow simulation in SPH

　　•本文提出了一种基于应力-应变一致性的双应力公式，用于SPH中的颗粒流模拟。•该公式避免了传统应力修正方法中存在的人为压力积累问题。•在长期旋转鼓模拟中，该公式展现了更好的数值稳定性。•在不同粒子分辨率下，该公式仍能得到一致的结果。

引言
颗粒材料是一类在工业过程和自然环境中广泛存在的物质。涉及颗粒材料输送与流动的问题出现在许多领域，包括采矿工程、冶金工业、土木工程以及制药行业[1]、[2]、[3]。此外，诸如雪崩和山体滑坡等许多自然灾害也可视为典型的颗粒流现象[4]、[5]、[6]。因此，深入理解颗粒材料的流动机制对于优化工业流程、提高结构安全性以及减轻自然灾害具有重要的理论意义和工程价值。

早期关于颗粒流的研究主要依赖于实验研究[7]、[8]以及理论分析[9]、[10]。实验研究通常针对特定的颗粒流现象设计，能够提供对宏观行为及局部细节的直观观察。然而，这类研究往往成本高昂，且难以系统地扩展到复杂条件和大规模场景。而理论方法则通常基于简化假设，因此仅适用于有限类型的问题，其应用范围相比实验更为受限。随着计算技术的快速发展，数值模拟逐渐成为研究颗粒流问题的重要工具。目前，用于颗粒流研究的常见数值方法包括离散元法（DEM）[3]、[11]、[12]、[13]以及有限元法（FEM）[14]、[15]、[16]。

在DEM框架中，颗粒材料被建模为相互作用的刚性或可变形粒子的集合，各个粒子的接触力及其运动会被明确计算出来。这使得能够准确描述包括平移、旋转和碰撞在内的粒子尺度行为，因此DEM非常适合用于研究微观物理机制[13]。然而，随着粒子数量的增加，DEM的计算成本会迅速上升，这在模拟大规模系统或长期演化过程时带来了显著限制。因此，DEM通常仅限于小规模、高分辨率的局部粒子行为研究，或是作为对实验观测的补充。

相比之下，FEM将颗粒材料视为连续体，通过求解带有适当本构模型的控制方程来描述其宏观力学行为。但作为一种典型的基于网格的方法，FEM的精度在很大程度上取决于网格质量。在传统的拉格朗日FEM框架中，计算网格会随材料一同移动，这在大变形情况下可能导致严重的网格畸变，进而降低精度甚至引发数值不稳定[15]。

为了解决颗粒流中常见的较大变形和自由面演化问题，近年来基于拉格朗日框架的无网格数值方法越来越受到关注，例如平滑粒子流体动力学（SPH）方法、移动粒子半隐式（MPS）方法[17]、[18]以及物质点法（MPM）[19]、[20]。其中，SPH方法将计算域离散为携带质量、速度和密度等物理量的粒子，并利用核函数来近似场变量及其空间导数[21]、[22]。通过这种方式，无需计算网格即可构建控制方程的粒子级离散形式，从而避免了网格畸变，使得SPH特别适合处理涉及大变形和动态自由面的问题[23]、[24]。正因如此，SPH被越来越多地应用于颗粒流模拟中[25]、[26]。

为了实现可靠的数值模拟，一个能够准确描述颗粒材料力学行为的本构模型至关重要。颗粒材料通常表现出高度非线性的力学特性：在静态或低流速状态下，它们表现为具有承载能力的固体；而在强剪切或外部干扰作用下，则表现出类似流体的流动特性。这种双重特性使得颗粒材料的本构建模本质上十分复杂。在本研究中，采用了Drucker–Prager（D–P）屈服准则和μI流变模型。D–P屈服准则是对Mohr–Coulomb准则的一种平滑近似，它能够更稳定地描述土壤、岩石和颗粒介质等摩擦材料从弹性变形向塑性流动的转变过程，已被广泛应用于地质材料和颗粒流的数值研究中[25]、[27]、[28]、[29]。μI流变模型则是基于实验观察建立的，它定义了摩擦系数μI与惯性数I之间的函数关系，从而表征剪应力与压力之间的比例关系。该模型能够描述颗粒材料在准静态状态与快速流动状态之间的转变[30]、[31]。以往的研究已将μI流变学纳入D–P弹塑性框架中，成功实现了对复杂颗粒流行为的有效模拟[29]、[32]。

然而，在SPH的连续体假设下，由D–P屈服准则所控制的应力演化在数值积分过程中并不总能严格保持在屈服面之内。尤其是当SPH粒子处于拉伸状态时，应力路径可能会进入与屈服面顶点相关的拉伸区域，从而违反屈服条件。为了解决这一问题，通常需要引入额外的应力修正步骤。现有的方法采用了一种相对简单的修正方式：当应力状态超出屈服面顶点时，会对主应力进行各向同性修正，将其重新投影回屈服面顶点上[25]。尽管在某些情况下，此类处理可以防止应力进入非物理的拉伸区域，但在涉及长期演化以及反复的拉伸-压缩循环的模拟中，它们可能会导致人为压力的积累，进而降低结果的精度和稳定性。

为克服这些局限性，本研究提出了一种双应力公式，其中区分了演化应力与有效应力。演化应力用于描述材料的应力-应变演化过程，确保本构响应的连续性和一致性；而有效应力则用于施加屈服约束，并被用于动量方程中。通过这种方式，既保持了应力-应变演化的连续性，又确保了进入动量方程的应力状态始终位于D–P屈服面之内。这一方法提升了SPH模拟颗粒流时的数值稳定性与一致性，尤其是在涉及周期性拉伸-压缩过程的长期模拟中。

章节节选
Drucker–Prager屈服准则
在本研究中，采用了非关联型的Drucker–Prager（D–P）屈服准则来构建颗粒材料的弹塑性本构模型，旨在描述材料在变形过程中从弹性响应向塑性流动的转变。张量分量在笛卡尔坐标系中表示，上标α、β和γ分别代表空间方向x、y和z。重复出现的指标遵循爱因斯坦求和约定。

非关联型D–P模型
颗粒流的控制方程包括体现质量守恒的连续性方程以及体现动量守恒的动量方程，具体表达式如下：
dρ/dt = ?ρ??v
dv/dt = 1/ρ??σ + g
其中ρ为密度，v为速度矢量，σ为应力张量，g则表示外力引起的加速度。需要注意的是，在密集的颗粒流中，应力主要通过粒子间的接触力传递，而单个粒子的可压缩性……

模型验证
为验证所提出的应力修正方法的有效性，本研究进行了周期性体积变形试验，以此说明为何需要应力修正，传统修正方法可能导致应力与应变之间出现不一致的原因，以及所提出的双应力公式如何在反复的压缩-伸展循环中保持应力-应变的一致性。在完成上述验证后，还研究了两个经典的基准问题：颗粒柱坍塌问题以及……

讨论
上述结果表明，所提出的双应力公式在典型的颗粒流问题中具有良好的稳定性和准确性。不过，其适用范围、局限性以及计算成本等问题仍需进一步探讨。本节将对这些方面进行简要分析。

本研究提出的双应力公式主要适用于无粘性颗粒材料的模拟。对于这类问题，从理论上来说，材料无法……

结论
本研究提出了一种基于应力-应变一致性的双应力公式，旨在解决传统应力修正方法可能引入人为压力这一问题，因为这种压力会导致非物理的体积膨胀、守恒性质的退化，以及颗粒材料长期模拟中数值稳定性的下降。该公式在SPH框架内得到了实现，并通过针对典型问题的数值模拟进行了系统验证。

RediT作者贡献声明
Fan Zhang：写作——审稿与编辑，写作——初稿，验证，软件，方法论，研究，正式分析。Chenyu Bao：软件。Hongqiang Lyu：监督，项目管理，研究。

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

致谢
本研究得到了中国国家自然科学基金的支持（项目编号：12472236）。
Fan Zhang | Chenyu Bao | Hongqiang Lyu