基于计算接触力学的粗颗粒土体材料细观性质模拟

潘洪武，王 伟，张丙印

(清华大学水沙科学与水利水电重点实验室，北京 100084)

颗粒物质是离散和相互作用的大量固体颗粒所组成的复杂体系。颗粒物质广泛存在于自然界，与人类的日常生活和生产密切相关。但是，颗粒体系的很多静力学和动力学行为尚不能用一般的连续介质力学理论很好地解释。《Science》杂志于2005 年把颗粒物质与湍流并列为100 个科学难题之一。

土体由离散的颗粒组成。颗粒间以强耗散的接触摩擦为主，其宏观力学特性主要取决于细观层次的颗粒联结和排列等。近年来，人们逐渐开始重视土颗粒体系细观力学行为的研究。由于在物理试验中存在对颗粒细观行为测量技术的难题，因而各种数值模拟方法逐步成为主要的研究手段。对此，学者们已经发展出了离散元法、非连续变形方法和连续-离散耦合分析方法等多种数值计算方法。

离散元法是Cundall[1]于1971 年提出的一种基于牛顿第二运动定律的非连续变形数值方法。该法将土体看成是离散单元的集合，通过给定颗粒间的接触关系[2]，来反映颗粒材料的宏细观力学响应。离散元法具有很强的描述颗粒间不连续变形的能力，受到人们广泛的青睐和重视，目前是进行颗粒材料细观力学特性分析的最主要的方法，取得了丰硕的研究成果[3-8]。

由于计算方法的限制，离散元法对颗粒本身细观力学特性的描述(如颗粒内部应力分布等)相对较为粗糙，其细观数值模拟难以定量地再现颗粒材料的宏细观力学特性。为此，Munjiza、马刚、王永亮等[9-11]发展了连续-离散耦合分析方法(FDEM)。

有限元方法是目前应用最广泛的数值模拟计算方法。但一般认为，有限元法主要适合于连续介质的模拟计算分析，其在不连续变形的描述方面存在缺陷。近年来，基于有限元方法的计算接触力学算法[12-15]得到了迅速发展，使得有限元方法在模拟不连续变形和多体接触问题方面的能力大幅提高。

土体是由离散的土颗粒组成的，其细观力学性质的研究为一个典型的多体接触问题。从原理上讲，计算接触力学方法是适用于进行这种多体接触分析的。基于上述的思路，本文开展了基于计算接触力学的土颗粒材料细观力学特性模拟的探索性研究工作。通过对粗颗粒土体材料进行三轴试验的数值计算分析，探讨了基于有限元方法的计算接触力学算法对粗颗粒土体材料多体接触问题的适用性。

1 计算接触力学算法简介

计算接触力学方法来自固体力学领域，是非线性数值分析中最具挑战性的课题之一，提出以来得到了来自计算数学和计算力学等领域诸多学者的研究。计算接触力学研究内容繁多，涉及范围很广。主要研究的问题包括接触界面离散形式(点对点、点对面和面对面等)、接触非线性求解、接触本构方程和接触约束施加方法等。计算接触力学方法被广泛应用在机械、航空和生物力学等领域。在岩土工程领域近年来才开始得到应用[16-18]。

1.1 接触问题的一般描述

直观地讲，接触问题反映的是不同物体之间的相互作用。对于本文考虑的粗颗粒土体材料，不失一般性，考虑如图1 所示的可变形颗粒 Ω1和Ω2的接触问题。将颗粒 Ωi的边界记作 ∂Ωi(i=1,2)，每个边界可分为Dirichlet 边界、Neumann 边界和接触边界，三者组成整体边界且相互无重叠：

图1 两个颗粒的接触问题示意图Fig. 1 Schematic diagram of two particle contact problem

接触边界上，两物体被记作从面和主面。将物体表面的法向和切向单位向量分别记为 n、 τ1和τ2， 接触区域的法向间隙 gn可表示为：

式中： xA为从面点的坐标； xˆB为 点 xA在Γ上的投影点。接触区域的物体接触力 tC可记作：

式中： pn为 接触力法向分量；tτ1、为接触力的切向分量。

在计算接触力学中，一般将上述接触条件表

1.2 接触问题的弱形式

1.3 接触界面的空间离散

为了在接触界面上计算接触虚功，需要对接触界面进行离散，常用的界面离散方法有点点方法、点面方法和面面方法等，其中，面面方法具有计算精度高、能通过接触分片试验等优点，是国内外学者们目前研究的热点，也是本文数值研究采用的方法。

在接触界面上，引入拉格朗日乘子，通过节点量插值记为：

1.4 接触问题的非线性迭代

由于接触区域和接触状态的双重不确定性，接触问题具有极强的非线性。计算接触力学中常采用PDASS (primal dual active-set strategy)方法[21]对接触问题进行非线性迭代求解。其基本思路为：首先将接触界面上的从点划分为脱开、滑移和粘结等不同状态的点集，并分别引入相应的接触边界条件；然后再根据相应互补的接触条件进行校核，判别从点的接触状态是否正确；对接触状态不正确的从点重新划分其接触状态。如此迭代直至接触界面从点的状态不再发生改变为止。

1.5 接触问题有限元方程的代数表达式

本文在界面离散上采用Dual mortar 元方法，其得到的D 矩阵为对角矩阵，可以高效地进行求逆操作。在激活节点上引入相对位移：

1.6 计算程序和验证算例

基于课题组发展的多体接触有限元计算程序Contana，针对颗粒集合体的特点，发展了相应的土体颗粒多体接触有限元计算程序。

采用四球堆积算例对所发展的有限元程序进行验证。如图2 所示，四个球体的半径R 相等。球体B2、B3和B4位于平面W1上，其球心连线为边长等于l=0.105 m 的等边三角形。球体B1对称放置在三个球体的顶部，四个球体的球心连线形成一个四面体。球体半径R=0.05 m，密度ρ=7800 kg/m3，弹性模型E=100 GPa，泊松比ν=0.3。本算例中，采用六面体单元对球体进行建模，每个球体含有1734 个节点、1152 个单元。

当球体间以及球体和底面间的摩擦系数μ1和μ2足够大时，四个球体处于稳定状态。记球体之间以及球体和底面之间接触力的法向和切向分量分别为Fn1、Ft1、Fn2和Ft2。根据静力平衡条件，可以求得不考虑球体变形条件下各接触力的大小。由于在算例情况下球体变形非常小，下文称该解答为理论解。

表1 给出了有限元接触力计算结果与理论值的对比情况。其中，理论解扣除了有限元网格划分所引起的质量误差。从表1 中可以看出，接触力的计算值和相应理论解的最大相对误差为0.129%，数值解和理论解吻合较好。

图2 四球堆积算例Fig. 2 Four sphere accumulation example

表1 接触力计算结果的对比Table 1 Comparison between calculation and theoretical values

2 三轴数值试验的计算模型

利用所发展的土体颗粒多体接触有限元计算程序，进行了不同围压的常规三轴数值试验，探讨了计算接触力学方法对土体颗粒系统细观力学特性分析的适用性。

2.1 计算模型概述

三轴试样由粒径分别为21 mm、19 mm、17 mm、15 mm 和13 mm 的460 个球颗粒组成。上述不同粒径的颗粒数目分别为70、95、130、95 和70。

如图3 所示，计算采用10 cm×10 cm×20 cm的柱体试样。试样的四个侧面和两个顶面共设置6 块加载板。加载板使用六面体单元离散。在加载板之间不设任何的接触关系，它们在计算中不会发生相互作用，可独立施加所需的应力或位移。每个球颗粒均被剖分为510 个四面体单元、133个节点。整个计算模型共包含有78748 个节点、245564 个单元。

图3 计算模型示意图Fig. 3 Schematic diagram of simulation model

计算中，颗粒和加载板均采用线弹性模型模拟，弹性模量分别取E=1 GPa 和25 GPa，泊松比均取ν=0.3。不考虑重力影响。颗粒-颗粒之间以及颗粒-加载板之间均采用库仑摩擦定律模拟，摩擦系数分别取0.4 和0.2。

2.2 数值试验计算过程

为了能与土体三轴试验的结果进行定性对比，参照土工三轴物理试验的流程，将本文数值试验的计算流程也划分为制样、固结和加载三个过程。数值试验的试样参照目前离散元方法所采用的三轴试样生成方法进行生成[5]。

采用生成的同一个数值试验试样，分别进行了围压σ3=200 kPa、400 kPa 和600 kPa 的常规三轴数值压缩试验。计算时，首先模拟施加围压的固结过程。采用Δσ3=50 kPa 分级施加围压力σ3。

对每一个围压的试验，固结完成后，保持侧向围压力σ3大小不变，采用应变控制的加载方式进行剪切加载。具体方式为，每个加载步使加载顶板向下位移0.01 cm，对应试样发生0.05%的轴向应变。

计算在个人计算机在上进行。计算机CPU 的型号为Intel® CoreTMi7-2600K(3.40 GHz)，每个围压试验的计算耗时约为40 h。

3 三轴数值试验计算结果分析

在本文的计算结果中，应力和应变均以压为正。如前所述，采用计算接触力学的方法，可以对颗粒本身划分有限元网格。因此，所得数值试验的结果不仅能反映试样的整体宏观力学特性，并且也能反映颗粒本身的细部力学性质。

3.1 试样的宏观力学特性

图4 给出了σ3=600 kPa 数值试验过程中试样的侧视图。可以看出，在试验过程中，颗粒本身不发生大的变形，试样的变形主要由颗粒在空间相对位置的变化所致，即主要由颗粒间的滑移和颗粒转动等所导致的颗粒集合体排列方式的改变所产生。

图4 试验过程中的试样侧视图(σ3=600 kPa)Fig. 4 Side view of specimen during experiment(σ3=600 kPa)

图5 给出了不同围压的应力-应变曲线和体变计算曲线。由图5(a)所示的(σ1-σ3)-ε1曲线可以看到，各围压的应力-应变曲线均表现出一定的应变软化特征。偏差应力随轴向应变的增加逐步增大，在开始阶段增加相对较快，之后逐步变得较为平缓，并达到峰值。峰值后随轴向应变增加，偏差应力会发生一定的下降。对不同围压的试验，对应相同的轴向应变，围压大对应的偏差应力也越高，表明试样的刚度或强度随围压逐步增大，表现出了颗粒材料压硬性的特点。这些特点均符合土体常规三轴压缩试验应力-应变曲线的一般特征。

图5 不同围压数值试验的计算结果Fig. 5 Computation results of numerical experiment at different confining pressures

由图5(b)所示的εv～ε1曲线可以看出，所得体变曲线总体符合粗粒土三轴试验中的体变特性。对一个围压的数值试验，试样在初始先发生体缩。当轴向应变较大时，体积变形逐步转变为剪胀。对不同围压的数值试验，围压较小时，试验初始段的体缩较小，后段的剪胀较为显著。随着围压的提高，试样的剪缩性总体增大，剪胀性逐步减小。根据试验结果绘制莫尔圆可得试样内摩擦角约为44°。

3.2 颗粒细部应力分析

图6 给出围压σ3=600 kPa 试验试样内部的力链分布情况。由图6 可见，在试验过程中，水平方向接触力小于竖直方向接触力，荷载主要沿竖直方向传递。对比不同轴向应变下试样力链分布可知，颗粒间接触力随着轴向应变增加而增大。

图6 试样力链分布Fig. 6 Distribution of the force chain

图7 进一步给出了围压σ3=600 kPa 试验中颗粒间法向接触力大小的分布情况。由图7 可以看出，当试样固结完成尚未进行剪切时，颗粒间法向接触力总体相对较小，绝大多数分布在小于4 kN以下的范围。开始加载剪切后，颗粒间法向接触力逐步增大。随轴向应变的增大，小法向力接触对数目逐渐降低，大法向力接触对数目逐渐增加。当轴向应变达到14%时，一些颗粒的接触应力可达16 kN 以上。

如前所述，采用计算接触力学方法可以对颗粒本身划分所需的有限元的单元，从而可以计算得到其详细的细观应力和变形状态。这是计算接触力学方法相比离散元方法的一个优势。在本次计算中，每个球颗粒均包含133 个节点、510 个单元，因而可以得到任意时刻所有颗粒详细的应力状态。

本文的计算模型中，颗粒采用四面体常应变单元进行网格剖分，且单元尺寸相近，因此可根据颗粒内部单元应力的统计分布来评估颗粒的总体应力状态。图8 给出了σ3=600 kPa 围压数值试验所得某个颗粒中单元大主应力大小的分布情况。当施加围压力之后但尚未进行剪切之前，试样处于等向压缩应力状态。此时，假如试样是连续的均质材料，则试样中所有单元的三个主应力都应为600 kPa。但由于试样由土体颗粒组成，因而计算所得单元应力的大小并不相等。其中，少数单元处于受拉状态，计算所得大主应力的最大值约为1.5 MPa。

图7 颗粒间法向接触力大小统计分布Fig. 7 Distribution of normal contact forces between particles

图8 颗粒单元大主应力统计分布Fig. 8 Distribution of large principle stress in a particle

当进行剪切后，由图8 可见，计算所得颗粒单元的大主应力总体增大，承受高应力的单元数目显著增加。少数单元仍会处于受拉状态。当轴向应变达到14%时，一些颗粒单元计算所得大主应力的最大值增大约为3.0 MPa。

图9 为σ3=600 kPa 围压的三轴数值试验，当轴向应变ε1=10%时，某典型颗粒中的大主应力分布的计算结果。由图9 可以看出，计算所得单元的大主应力在颗粒中的分布非常的不均匀。其中，在颗粒中的大部分区域，计算所得到的应力值相对均较低。但在颗粒接触点附近的局部小区域，则会计算得到相对很高的应力分布。对于颗粒材料，其所受到的作用力是通过颗粒接触点进行传递的。对于本文所采用的球形颗粒，在接触点部位会作用集中的粒间力的作用。计算所得上述的应力分布特征符合赫兹接触问题的一般规律。

图9 大主应力剖视图 /MPa Fig. 9 Section view of large principal stress

4 计算接触力学法和离散元法对比

由本文前面的分析可知，计算接触力学方法是基于有限元方法的一种进行多体接触问题的数值求解方法。与离散元方法相比，在基本原理、求解方法、颗粒本身的模拟以及接触条件的模拟等方面都存在显著的不同。表2 总结列出了两种方法在模拟颗粒材料特性方面的主要差别。

表2 两种模拟方法对比Table 2 Comparison between two simulation methods

总体而言，将计算接触力学方法引入土体颗粒集合体细观应力变形特性的模拟分析，可能会具有如下的优势：

1)在基本原理和求解方面，计算接触力学方法是一种基于变分原理的隐式数值求解方法，具有严格的理论基础。离散元方法的基本原理是牛顿第二定律，根据颗粒受到的合力计算颗粒在各时间步的位移，是一种显示求解方法。自提出至今，离散元及其相关方法的研究已有近50 年历史，众多学者发表了大量学术论文。但是，离散元方法缺乏理论的严密性。刘凯欣和高凌天[22]曾评述“研究离散元的理论和算法的文章却很少”“自诞生的那天起就带有缺乏理论严密性的先天不足”。

2)在颗粒本身的模拟方面，采用计算接触力学方法可以对颗粒本身划分所需的有限单元，从而可以得到其详细的细观应力和变形状态，可为颗粒破碎等的模拟计算提供基础。传统离散元方法的模型一般由刚性块体组成，此时无法考虑颗粒体本身变形，也无法计算得到颗粒体内部应力的分布。

3)在颗粒间的接触特性模拟方面，计算接触力学方法通过拉格朗日乘子法引入KKT 接触条件，来模拟颗粒接触界面的接触特性，是一种高精度的接触问题求解方法。拉格朗日乘子的物理意义为接触应力。该方法无须引入接触弹簧，可避免罚模量选取对人工经验的依赖等。离散元方法采用罚函数法，通过在接触点引入法向和切向弹簧模拟颗粒间的接触行为，存在弹簧刚度以及罚模量选取对人工经验的依赖等。

但是，大规模的土体颗粒集合体是一个大型强非连续颗粒系统，将海量的颗粒间接触条件引入有限元方程后，会造成系统刚度矩阵的高度病态。为此，研发针对大型复杂颗粒集合体的高效接触搜索算法、高度病态矩阵预处理技术和具有较强鲁棒性的迭代算法等，是该种方法能否成功应用于土体颗粒集合体细观力学行为模拟分析的关键。这也是该法亟待在未来研究中解决的关键技术难题。

5 结论

本文开展了基于计算接触力学的粗颗粒土体材料细观力学行为模拟分析的探索性研究工作，主要取得如下的结论：

(1)提出了一种基于计算接触力学的粗颗粒土体材料细观力学行为模拟分析新方法。该法将颗粒剖分成一定数量的有限元单元，通过计算接触力学方法模拟颗粒间的接触行为。

(2)和离散元法方法相比，本文提出的模拟计算方法是一种基于变分原理的隐式数值求解方法，具有严格的理论基础。此外，该方法在描述颗粒本身的力学特性方面具有优势，既可以计算分析粗颗粒土体材料的宏细观力学特性，也可以计算得到颗粒内部的应力分布情况，可为颗粒破碎等的细观行为的模拟计算提供依据。

(3)利用所发展的土体颗粒多体接触有限元计算程序，进行了不同围压的常规三轴数值试验，模拟计算结果符合土体三轴试验的一般规律，初步验证了所提出的计算方法的适用性。