离散元法中颗粒の接触建模方法

离散元法中颗粒の接触建模方法

---

前言

本文选自《Contact Modeling for Solids and Particles》一书中的第四章。该书于2018年出版，由德国科学院院士 Peter Wriggers 和 Alexander Popp 教授两人共同编写，其中第四章的颗粒接触建模部分由 Jerzy Rojek 所写，本人对该章内容进行了整理与翻译，供研究离散元理论的同学们学习使用。离散单元法（DEM）是一类广泛用于离散不连续材料和系统建模的数值方法，通过用大量的颗粒（离散单元）来表示，并且假设离散单元之间通过接触力相互作用。本章介绍了离散元法中接触建模的基本方法。给出了离散元法的基本假设、理论公式和数值算法。本文中考虑了采用球形颗粒和软接触方法的离散元公式；综述了颗粒间相互作用的基本接触模型；讨论了基本的接触机理，包括弹性、塑性、阻尼、摩擦和黏聚力；描述了结合这些机理所选定的接触模型；分析了它们在模拟动态或准静态的单个接触事件中的性能，其分析重点是单次碰撞过程中接触力的演变。虽然本文主要讨论的是相互作用力，但也对相互作用力矩进行了介绍，还对热效应和热力耦合的离散元接触公式进行了考虑。

---

一、引言

离散单元法（DEM）是一类用于分析颗粒系统（离散单元）接触相互作用力学不连续问题的数值方法。离散元法是一种较新的数值计算方法，是在20世纪70、80年代由 Cundall（1971）、Cundall 和Strack（1979）、Walton（1982、1983）在其开创性著作中所引入的方法。离散元法在 Williams等人（1985）、Bardet 和 Proubet（1991）、Moreau（1994） 等人的著作中得到了进一步的发展。它已经成为模拟各种颗粒和非颗粒材料行为的有力工具，这些材料可以用颗粒系统来表示，如土壤（Widuli´nski 等人. 2009）、粉末（Martin 等人. 2003）、岩石（Cundall 1987；Potyondy 和 Cundall 2004；Rojek 等人. 2011；Zubelewicz 和 Mroz 1983），混凝土（Hentz 等人. 2004；Wu 等人. 2013)，陶瓷（Senapati 和 Zhang 2010），甚至金属（Fleissner 等人. 2007）。

离散单元可以是任意形状（Rothenburg 和 Bathurst 1992；Tao 等人. 2014；Cundall 1988）。然而，球形颗粒通常是一个更好的选择（Cundall 1987；Widuli´nski 等人. 2009；Plassiard 等人. 2009），因为球形颗粒接触检测算法的公式简单，计算效率高，所以在本文中主要讨论球形颗粒的离散元方程。

接触算法在离散元法中起着至关重要的作用。离散单元的运动由接触力来控制，并最终体现为颗粒集的宏观行为。在离散元中，接触问题可以分为两种不同的处理方法，即所谓的软接触法（Cundall 和 Strack 1979；Cundall 1987；Potyondy 和 Cundall 2004）和硬接触法（Hong 和 McLennan 1992；Haff 和 Werner 1987；Richardson 等人. 2011年）。软接触法采用正则化的接触约束，而硬接触法采用非光滑分析方法来解决单边接触约束问题。

在软接触法中，允许颗粒之间有少量的重叠，仅近似满足接触非穿透条件，颗粒间的接触时间比时间步长，并分析了接触力的变化。在硬接触法中，不允许颗粒穿透，假设碰撞时间很短，因此可以忽略，分析了由碰撞引起的颗粒动量的变化，而未分析接触力的变化。在本文中，只考虑了软接触法。这种方法允许我们采用一个合适的接触模型来模拟颗粒接触，以及一个合适的模型来获得所需的宏观行为。本文旨在介绍离散单元法中接触建模的基本概念。

本文的大纲如下：第2节中给出了离散元法的计算公式，简述了基本假设、运动方程和时间积分方案；第3节中专门介绍了接触建模，推导了接触条件，并引入了采用惩罚正则化后的接触约束条件，提出了接触模型中所包含的基本接触机理；第4节中介绍了更复杂的接触模型和它们在简单问题中的公式和性能，推导了相互作用力和力矩。第5节中介绍了离散元法在热力学和热力学问题上的推广；第6节中简要介绍了热力和热机械公式。

---

二、离散元法公式

2.1 基本假设

本节将考虑刚性圆形（二维）或球形（三维）颗粒相互作用系统的动力学，以及颗粒的平动和转动。假设颗粒的初始位置和速度（包括线速度和角速度）是已知的，颗粒受到外部载荷（集中力和集中力矩）、重力和接触力（颗粒与周围介质的相互作用）。假设这些颗粒相互接触或者与其他物体相互作用。所要解决的问题是求解一个常微分方程（运动方程），其初始条件由接触约束条件所定义。

2.2 运动方程

离散单元（颗粒）的运动可由刚体动力学中标准牛顿-欧拉方程表示。设第 $i$ 个圆形或球形颗粒（如图1）的平动和旋转运动的方程为：
$$\begin{array}{r}{m}_i{\ddot{\mathbf{u}}}_i={\mathbf{F}}_i\end{array}\text{\hspace{1em}(1)}$$

$$\begin{array}{r}{J}_i{\dot{\mathbf{\omega }}}_i={\mathbf{T}}_i\end{array}\text{\hspace{1em}(2)}$$

其中，$u_i$ 是在惯性坐标系 $X$ 中的颗粒质心的位移。

$$\begin{array}{r}{\mathbf{u}}_i={\mathbf{x}}_i-{\mathbf{X}}_i\end{array}\text{\hspace{1em}(3)}$$

其中，${\mathbf{\omega }}_i$ 为角速度；${\mathbf{F}}_i$ 为合力；${\mathbf{T}}_i$ 是围绕中心轴的合成力矩；$m_i$ 是颗粒质量；$J_i$ 是惯性矩，其计算公式为：

$$\begin{array}{r}{J}_i=\frac{1}{2}{m}_i{R}_i^2二维圆形\end{array}\text{\hspace{1em}(4)}$$

$$\begin{array}{r}{J}_i=\frac{2}{5}{m}_i{R}_i^2三维球体\end{array}\text{\hspace{1em}(5)}$$

其中，$R_i$ 是第 $i$ 个颗粒半径。转动方程（2）是任意刚体的一般形式的简化，适用于球体和圆形（二维），其转动惯性性质由二阶张量表示。向量 ${\mathbf{F}}_i$ 和 ${\mathbf{T}}_i$ 为以下各项的和：

（Ⅰ）由外部载荷施加到第 $i$ 个颗粒上的集中力 ${\mathbf{F}}_i^{ext}$ 和集中力矩 ${\mathbf{T}}_i^{ext}$。

（Ⅱ）与邻近颗粒和其它所有元素相互作用的接触力 ${\mathbf{F}}_{ij}^{cont}$ 和接触力矩 ${\mathbf{T}}_{ij}^{cont}$， $j=1,\dots ,n_i^c$ ，$n_i^c$ 表示与第 $i$ 个颗粒接触的元素数。

（Ⅲ）由外界（背景）阻尼产生的力 ${\mathbf{F}}_i^{damp}$ 和力矩 ${\mathbf{T}}_i^{damp}$ 。

因此，向量 $F_i$ 和 $T_i$ 可以写成如下表达式：

$$\begin{array}{r}{\mathbf{F}}_i={\mathbf{F}}_i^{ext}+\sum _{j=1}^{n_i^c}{\mathbf{F}}_{ij}^{cont}+{\mathbf{F}}_i^{damp}\end{array}\text{\hspace{1em}(6)}$$

$$\begin{array}{r}{\mathbf{T}}_i={\mathbf{T}}_i^{ext}+\sum _{j=1}^{n_i^c}{\mathbf{S}}_{ij}^c\times {\mathbf{F}}_{ij}^{cont}+\sum _{j=1}^{n_i^c}{\mathbf{T}}_{ij}^{cont}+{\mathbf{T}}_i^{damp}\end{array}\text{\hspace{1em}(7)}$$

其中 ${\mathbf{S}}_{ij}^c$ 是连接第 $i$ 个颗粒质心到与第 $j$ 个元素接触点的向量（如图2）。

2.3 时间积分方案

利用显式中心差分格式对运动方程（1）和（2）进行时间积分，在第 $n$ 个时间步长上的平移运动的时间积分算子如下：

$$\begin{array}{r}{\ddot{\mathbf{u}}}_i^n=\frac{{\mathbf{F}}_i^n}{m_i}\end{array}\text{\hspace{1em}(8)}$$

$$\begin{array}{r}{\dot{\mathbf{u}}}_i^{n+\frac{1}{2}}={\dot{\mathbf{u}}}_i^{n-\frac{1}{2}}+{\ddot{\mathbf{u}}}_i^n\Delta t\end{array}\text{\hspace{1em}(9)}$$

$$\begin{array}{r}{\mathbf{u}}_i^{n+1}={\mathbf{u}}_i^n+{\dot{\mathbf{u}}}_i^{n+\frac{1}{2}}\Delta t\end{array}\text{\hspace{1em}(10)}$$

旋转运动积分方案中的前两步与方程式（8）和（9）给出的步骤相同：

$$\begin{array}{r}{\dot{\mathbf{\omega }}}_i^n=\frac{{\mathbf{T}}_i^n}{J_i}\end{array}\text{\hspace{1em}(11)}$$

$$\begin{array}{r}{\mathbf{\omega }}_i^{n+\frac{1}{2}}={\mathbf{\omega }}_i^{n-\frac{1}{2}}+{\dot{\mathbf{\omega }}}_i^n\Delta t\end{array}\text{\hspace{1em}(12)}$$

旋转增量 $\Delta \mathbf{\theta }$ 计算公式如下:

$$\begin{array}{r}\Delta {\mathbf{\theta }}_i={\mathbf{\omega }}_i^{n+\frac{1}{2}}\Delta t\end{array}\text{\hspace{1em}(13)}$$

旋转增量的计算公式可以用来更新切向接触力，如果有必要，还可以跟踪颗粒旋转位置的总变化（Argyris 1982）。进一步地，嵌入在颗粒中的移动帧和固定全局帧之间的旋转矩阵必须使用适当的乘法格式进行增量更新（Rojek等人，2001）。

在时间上的显式积分可以获得较高的计算效率。显式积分方案的缺点是其数值稳定性的条件限制了时间步长 $\Delta t$，时间步长 $\Delta t$ 不能大于临界时间步长 $\Delta t_{cr}$ ：

$$\begin{array}{r}\Delta t\le \Delta t_{cr}\end{array}\text{\hspace{1em}(14)}$$

其中，临界时间步长 $\Delta t_{cr}$ 由系统的最高固有频率 ${\nu }_{max}$ 所决定：

$$\begin{array}{r}\Delta t_{cr}=\frac{2}{{\nu }_{max}}\end{array}\text{\hspace{1em}(15)}$$

确定最高固有频率 ${\nu }_{max}$ 的精确值需要解决整个连通刚性颗粒系统所定义的特征值问题。整个系统的最高固有频率可估计为每个颗粒 $e$ 周围相接触颗粒子集的固有频率 ${\nu }_i^e$ 的最大值 ，参考 Belytschko 等人（1985）：

$$\begin{array}{r}{\nu }_{max}\le {\nu }_{max}^D其中{\nu }_{max}^D=\underset{i,e}{\max }{\nu }_i^e\end{array}\text{\hspace{1em}(16)}$$

---

三、接触模型

3.1 接触条件

当采用不同的搜索方法对离散单元的接触对进行识别时，接触颗粒应满足接触点的约束条件。接触约束条件可以用接触点的相互作用和适当的运动学参数来表示。

假设两颗粒的相互作用集中在一个称为接触点的点上，两个颗粒 $i$ 和 $j$ 之间的相互作用由集中力 ${\mathbf{F}}_{ij}^{cont}$ 和等效在接触点的集中力矩 ${\mathbf{T}}_{ij}^{cont}$ 组成（如图2）。

其中，接触力 ${\mathbf{F}}_{ij}^{cont}$ 可以分解为法向分量 $({\mathbf{F}}_n^{cont}{)}_{ij}$ 和切向分量 $({\mathbf{F}}_t^{cont}{)}_{ij}$ ：

$$\begin{array}{r}{\mathbf{F}}_{ij}^{cont}=({\mathbf{F}}_n^{cont}{)}_{ij}+({\mathbf{F}}_t^{cont}{)}_{ij}=（F_n^{cont}{)}_{ij}{\mathbf{n}}_{ij}+({\mathbf{F}}_t^{cont}{)}_{ij}\end{array}\text{\hspace{1em}(17)}$$

其中 ${\mathbf{n}}_{ij}$ 是接触点法向上的单位向量，定义如下：

$$\begin{array}{r}{\mathbf{n}}_{ij}=\frac{{\mathbf{x}}_j-{\mathbf{x}}_i}{||{\mathbf{x}}_j-{\mathbf{x}}_i||}\end{array}\text{\hspace{1em}(18)}$$

假设法向接触和切向接触可以解耦，并且可以单独考虑。法向上单边接触（无粘连接触）的基本条件如下：

$$\begin{array}{r}{\mathbf{F}}_n^{cont}\le 0,g\ge 0,{\mathbf{F}}_n^{cont}g=0\end{array}\text{\hspace{1em}(19)}$$

其中，$g$ 为颗粒之间的间隙（见图3）：

$$\begin{array}{r}g=d_{ij}-R_i-R_j\end{array}\text{\hspace{1em}(20)}$$

其中，$d_{ij}$ 是颗粒质心之间的距离

$$\begin{array}{r}{d}_{ij}=||{\mathbf{x}}_j-{\mathbf{x}}_i||\end{array}\text{\hspace{1em}(21)}$$

公式（19）中的第一个不等式表示应力条件（不允许有拉应力），第二个不等式定义了不穿透性条件，第三个等式称为补充条件。令 ${\mathbf{F}}_n^{cont}\le 0$ 和 $g=0$ 或 ${\mathbf{F}}_n^{cont}=0$ 和 $g>0$，则单侧法向接触定律如图4a所示。在双边（有粘连）接触中，允许有拉应力，并将（19）中的几何不等式约束 $g>0$ 替换为等式约束 $g=0$（Curnier 1999）。对双边接触的严格数学推导要比单边接触复杂得多。

切向上的相互作用通常是由颗粒间的摩擦引起，切向上的摩擦滑动接触的补充条件（Klarbring 1999）如下：

$$\begin{array}{r}{\varphi }_t\le 0,{\lambda }_t\ge 0,{\varphi }_t{\lambda }_t=0\end{array}\text{\hspace{1em}(22)}$$

其中，${\varphi }_t$ 为滑动准则，非负参数 ${\lambda }_t$ 由滑移定律所定义：

$$\begin{array}{r}{\mathbf{v}}_{rt}={\lambda }_t\frac{{\mathbf{F}}_t^{cont}}{||{\mathbf{F}}_t^{cont}||}\end{array}\text{\hspace{1em}(23)}$$

其中，${\mathbf{v}}_{rt}$ 为接触点处的切向相对速度：

$$\begin{array}{r}{\mathbf{v}}_{rt}={\mathbf{v}}_r-v_{rn}{\mathbf{n}}_{ij}\end{array}\text{\hspace{1em}(24)}$$

其中，${\mathbf{v}}_r$ 和 $v_{rn}$ 为接触点处的总相对速度和法向相对速度：

KaTeX parse error: Unexpected end of input in a macro argument, expected '}' at position 31: …} \bm v_r=(\dot\̲b̲m̲ ̲u_j+\bm\omega_j…

$$\begin{array}{r}{v}_{rn}={\mathbf{v}}_r\cdot {\mathbf{n}}_{ij}\end{array}\text{\hspace{1em}(26)}$$

其中，KaTeX parse error: Unexpected end of input in a macro argument, expected '}' at position 5: \dot\̲b̲m̲ ̲u_i 和 KaTeX parse error: Unexpected end of input in a macro argument, expected '}' at position 5: \dot\̲b̲m̲ ̲u_j 是颗粒质心的平移速度，${\mathbf{\omega }}_i$ 和 ${\mathbf{\omega }}_j$ 颗粒的角速度， ${\mathbf{s}}_{ij}^c$ 和 ${\mathbf{s}}_{jj}^c$ 是连接颗粒质心到接触点的向量。

对于滑动摩擦的临界值有多种模型（Raous 1999）。最常用的是库仑摩擦力模型，其滑动准则为：

$$\begin{array}{r}{\varphi }_t=||{\mathbf{F}}_t^{cont}||-\mu ||{\mathbf{F}}_n^{cont}||\le 0\end{array}\text{\hspace{1em}(27)}$$

其中，$\mu$ 为库仑摩擦系数。库仑摩擦接触对应的图如图4b所示。库仑摩擦系数通常是假定的常数，但它也可以作为一个变量，例如，其大小与滑动速度速度有关。

3.2 接触条件的惩罚正则化

基于软接触离散元法对单边接触的法向和切向的约束条件进行了惩罚正则化。惩罚正则化的法向接触条件为：

$$\begin{array}{r}{\mathbf{F}}_t^{cont}=k_{n}g当g<0\end{array}\text{\hspace{1em}(28)}$$

其中，$k_n$ 是一个特定的惩罚参数。接触条件 ${\mathbf{F}}_n^{cont}\le 0$ 和 ${\mathbf{F}}_n^{cont}g=0$ 仍然有效。仅近似满足不穿透性条件 $g\ge 0$。因为此时的接触颗粒之间允许有一定的重叠量 $h$：

$$\begin{array}{r}h=-g>0\end{array}\text{\hspace{1em}(29)}$$

当 $k_n\to \infty$，达到了惩罚的目的，然而应该注意的是，较大的惩罚值会导致由等式给出的临界时间步长值 $\Delta t_{cr}$ 较小，对接触点法向上的惩罚情况如图5a所示。

通过在式（23）中引入切向惩罚项 $k_t$ 来实现切向上摩擦约束的正则化：

KaTeX parse error: Unexpected end of input in a macro argument, expected '}' at position 90: …t||}=\frac{\dot\̲b̲m̲ ̲F^{cont}_t}{k_t…

式（30）表明，库仑摩擦约束的惩罚正则化可将总滑动速度 ${\mathbf{v}}_{rt}$ 分解为可逆和不可逆部分，分别为 ${\mathbf{v}}_{rt}^r$ 和 ${\mathbf{v}}_{rt}^{ir}$：

$$\begin{array}{r}{\mathbf{v}}_{rt}={\mathbf{v}}_{rt}^r+{\mathbf{v}}_{rt}^{ir}\end{array}\text{\hspace{1em}(31)}$$

其中

KaTeX parse error: Unexpected end of input in a macro argument, expected '}' at position 41: …{rt}=\frac{\dot\̲b̲m̲ ̲F^{cont}_t}{k_t…

$$\begin{array}{r}{\mathbf{v}}_{rt}^{ir}={\lambda }_t\frac{{\mathbf{F}}_t^{cont}}{||{\mathbf{F}}_t^{cont}||}\end{array}\text{\hspace{1em}(33)}$$

库仑摩擦正则化模型对应的图如图5b所示，其中切向接触力 $F_t^{cont}$ 被绘制为与通过对切向相对速度积分得到的切向相对位移有关的函数：

$$\begin{array}{r}{u}_{rt}=\int {\mathbf{v}}_{rt}\mathrm{d}t\end{array}\text{\hspace{1em}(34)}$$

3.3 惩罚正则化的物理解释

法向和切向接触约束的惩罚相当于在界面上指定附加的本构关系。因此，上文介绍了法向接触中弹性压缩行为，并将摩擦接触重新表述为类似于弹塑性的问题。

假设由式（29）定义的颗粒接触的重叠量为由于接触相互作用导致的接触点处颗粒局部变形的影响（如图6）。该假设提供了定义接触力与重叠量之间关系的可能性，以便更好地表示接触区域内的各种变形机制。类似地，对于切向接触，可以定义接触力与滑移量之间的关系。

必须指出的是，如果真实颗粒的变形量相对较小，便可以合理假设由接触引起的颗粒变形是局部的，它不影响其他颗粒的接触。

3.4 简单接触的变形机制

通过建立接触模型可以解释接触中所涉及的各种变形机制和物理现象。下面总结了典型的简单接触变形机制和相关物理效应。

（Ⅰ）弹性

在线性弹性模型中，法向或切向上的接触力是一个与位移量 $u$ 有关的线性或非线性函数：

$$\begin{array}{r}F=ku\end{array}\text{\hspace{1em}(35)}$$

其中 $k$ 是一个刚度常数，位移量 $u$ 可以表示颗粒的重叠量或切向相对位移。而非线性弹性模型的特点是具有可变的刚度。弹性接触模型中的位移是完全可逆的。

（Ⅱ）塑性

在由接触压力引起的颗粒塑性变形的响应中，接触力是一个与位移量 $u$ 有关的线性或非线性函数。在理想的塑性接触模型中，位移是完全不可逆的，塑性理论为摩擦建模提供了可行的依据。

（Ⅲ）粘性

在接触相互作用中，由于粘性响应而产生的接触力与速度之间的关系定义为：

$$\begin{array}{r}F=\eta \dot{u}\end{array}\text{\hspace{1em}(36)}$$

其中 $\eta$ 表示接触界面的粘性特性，在线性粘性模型中 $\eta$ 为常数，而在非线性粘性模型中 $\eta$ 为变量。粘性模型可以用来表示物理现象，如阻尼或蠕变。

（Ⅳ）摩擦

摩擦是一种与接触颗粒的切向相对运动相反的耗散机制，有时被称为干阻尼，（Zonetti 等人. 1999），而不是上面提到的与速度相关的粘性阻尼。库仑模型是最常用的摩擦模型。库仑模型中的接触力图如图4b所示。

（Ⅴ）粘结力

用离散单元对岩石或混凝土等粘性材料进行建模，需要在接触模型中考虑粘结力。接触颗粒之间引入了粘结键。这些键传递接触力，阻止颗粒在法线方向上的分离，以及在切向方向上的相对运动。

（Ⅵ）损伤

损伤是指由于内部裂纹的扩张而导致的机械材料性能（如刚度和强度）的恶化。考虑到接触中的损伤效应，损伤可以表示为内聚粘结键所代表的机械性能的逐渐恶化。

（Ⅶ）破裂

损伤累积到一定程度会导致粘结键完全退化。当超过粘结强度时，可以假定粘结键以脆性方式断裂。这样，材料中裂纹的萌生和发展就可以用离散元法进行模拟。

（Ⅷ）热效应

颗粒间的摩擦接触会伴随着热量的产生，摩擦产生的热量通过颗粒吸收并传导，从而导致颗粒温度升高，可能会影响机械接触性能，如接触刚度、粘度或摩擦系数（Shillor 等人. 2004）。之后将进一步介绍考虑热效应的接触模型公式以及离散元法的热力和热机械公式。

离散元法中的接触模型通常包含上述不同的机制和效应，这使我们能够模拟真实材料的复杂行为。接触模型与本构材料模型相似，通常用流变结构图来表示。流变结构图是由代表基本机理的流变元件组成的，典型流变元件如图7所示。

线性和非线性弹簧（图7a，b）表示流变模型中的弹性特性，线性阻尼器（图7c）对应于等式（36）中描述的粘性特性，滑块（图7d）表示摩擦特性，当滑块单独表示时，为非正则库仑摩擦模型；当滑块与弹簧串联时，则表示正则化的库仑摩擦模型。

---

四、接触模型的选择

4.1 考虑库仑摩擦的线性粘弹性接触模型

4.1.1 模型的建立

这里提出的模型与 Cundall 和 Strack（1979）在开创性研究中所提出的模型相似。模型的流变结构图如图8所示。法向接触力由粘弹性 Kelvin–Voigt 单元表示，该单元由线性弹簧和线性阻尼器通过并联构成。切向接触力对应的单元由一个弹簧与滑块串联构成。

Kelvin–Voigt 单元传递的法向接触力 $F_n^{cont}$ 由弹簧传递的弹力 $F_n^e$ 和阻尼器传递的粘性阻尼力 $F_n^d$ 组成：

$$\begin{array}{r}{F}_n^{cont}=F_n^e+F_n^d\end{array}\text{\hspace{1em}(37)}$$

弹力 $F_n^e$ 根据类似于式（28）的线性关系来计算：

$$\begin{array}{r}{F}_n^e=k_{n}g\end{array}\text{\hspace{1em}(38)}$$

其中，$k_n$ 为法向接触刚度，$g$ 由式（20）定义。当 $g<0$ 时，表示颗粒重叠，式（38）成立。当 $g\ge 0$ 时，则设弹力为零 $F_n^e=0$ 。

阻尼力 $F_n^d$ 根据类似于式（36）的线性关系来计算：

$$\begin{array}{r}{F}_n^d=c_n{v}_{rn}\end{array}\text{\hspace{1em}(39)}$$

式中，$c_n$ 为法向粘性阻尼系数，$v_{rn}$ 为由式（26）所定义的接触处法向相对速度。

图8所示接触模型的切向部分对应于第3.2节所描述的正则化库仑摩擦模型。线性弹簧的刚度 $k_t$ 对应于式（30）中引入的惩罚参数。滑块阻碍切向运动直至到达到滑动准则（27）。切向接触力与切向相对位移之间的关系如图5所示。

考虑正则化摩擦接触模型与弹塑性模型具有相似性，我们可以使用类似于弹塑性模型中的径向返回算法来计算摩擦力。首先计算切向接触力的试算值：

$$\begin{array}{r}{\mathbf{F}}_t^{trial}={\mathbf{F}}_t^{old}-k_t{\mathbf{v}}_{rt}\Delta t\end{array}\text{\hspace{1em}(40)}$$

然后判断是否满足滑动准则：

$$\begin{array}{r}{\varphi }^{trial}=||{\mathbf{F}}_t^{trial}||-\mu |F_n|\end{array}\text{\hspace{1em}(41)}$$

如果 ${\varphi }^{trial}\le 0$ 则为静摩擦接触，摩擦力为试算值：

$$\begin{array}{r}{\mathbf{F}}_t^{new}={\mathbf{F}}_t^{trial}\end{array}\text{\hspace{1em}(42)}$$

否则，为滑动接触，将执行返回映射：

$$\begin{array}{r}{\mathbf{F}}_t^{new}=\mu |F_n|\frac{{\mathbf{F}}_t^{trial}}{||{\mathbf{F}}_t^{trial}||}\end{array}\text{\hspace{1em}(43)}$$

4.1.2 模型参数计算

在离散元法中评估接触刚度 $k_n$ 有不同的方法，它可以在整个离散元的集合体中被看作是均匀的（Rojek 等人. 2005年），或者在局部中看作是均匀的，通常假设它取决于接触颗粒的尺寸（Potyondy 和 Cundall 2004），可以通过颗粒半径 $R_i$ 和 $R_j$ 来计算：

$$\begin{array}{r}{k}_n=f_k(R_i,R_j)\end{array}\text{\hspace{1em}(44)}$$

Rojek 等人. (2012) 讨论了对函数 $f_k(R_i,R_j)$ 形式的不同假设。在这里，我们将介绍其中的一个。

弹簧元件可以等效为一个非均匀截面积的弹性杆（见图9），由两个段组成：

其长度分别为：

$$\begin{array}{r}{L}_i=R_i,L_j=R_j\end{array}\text{\hspace{1em}(45)}$$

横截面积分别为：

$$\begin{array}{r}{A}_i={\alpha }_i\pi (R_i{)}^2,A_j={\alpha }_j\pi (R_j{)}^2\end{array}\text{\hspace{1em}(46)}$$

其中，$0\le {\alpha }_i,{\alpha }_j\le 1$ 是颗粒横截面积的系数，将杆段的面积定义为颗粒横截面积的分数。

两个杆段组成的系统可以看作两个弹簧串联。整个系统传递的轴向力 $F_e$ 等于 $i$ 、$j$ 段轴向力 $F_i^e$ 和 $F_j^e$ 之和：

$$\begin{array}{r}{F}_e=F_i^e+F_j^e\end{array}\text{\hspace{1em}(47)}$$

假设接触系统的整体轴向变形等于重叠量 $g$ （g < 0），则可以将其分解为两个分段的变形，即 $g_i$ 和 $g_j$：

$$\begin{array}{r}g=g_i+g_j\end{array}\text{\hspace{1em}(48)}$$

各段杆的力与位移关系可表示为：

$$\begin{array}{r}{F}_i^e=k_n^i{g}_i\end{array}\text{\hspace{1em}(49)}$$

$$\begin{array}{r}{F}_j^e=k_n^j{g}_j\end{array}\text{\hspace{1em}(50)}$$

其中，和 $k_i^n$ 和 $k_j^n$ 是 $i$、$j$ 段的刚度。联立（38）（47）（48）（49）（50）得到：

$$\begin{array}{r}\frac{1}{k_n}=\frac{1}{k_n^i}+\frac{1}{k_n^j}\end{array}\text{\hspace{1em}(51)}$$

对等式两边取倒数得：

$$\begin{array}{r}{k}_n=\frac{k_n^i{k}_n^j}{k_n^i+k_n^j}\end{array}\text{\hspace{1em}(52)}$$

表达式（52）与 Potyondy 和 Cundall（2004）所使用的相同。

利用（45）和（46），$i$ 、$j$ 段的刚度可表示为：

$$\begin{array}{r}{k}_n^i=\frac{E_i{A}_i}{L_i}={\alpha }_i\pi E_i{R}_i\end{array}\text{\hspace{1em}(53)}$$

$$\begin{array}{r}{k}_n^j=\frac{E_j{A}_j}{L_j}={\alpha }_j\pi E_j{R}_j\end{array}\text{\hspace{1em}(54)}$$

其中，$E_i$ 和 $E_j$ 是杆段（或颗粒）$i$ 和 $j$ 材料的杨氏模量。将（53）和（54）代入（52）中，并假设 $E_i=E_j=E$ 和 ${\alpha }_i={\alpha }_j=\alpha$，得到等效刚度 $k_n$ 的表达式如下：

$$\begin{array}{r}{k}_n=\alpha \pi E{R}^{\ast }\end{array}\text{\hspace{1em}(55)}$$

其中，$R^{\ast }$ 是根据颗粒半径 $R_i$ 和 $R_j$ 定义的有效半径：

$$\begin{array}{r}\frac{1}{R^{\ast }}=\frac{1}{R_i}+\frac{1}{R_j}\end{array}\text{\hspace{1em}(56)}$$

对于等尺寸颗粒接触 $R_i=R_j=R$，式（55）的形式为：

$$\begin{array}{r}{k}_n=\alpha \pi ER\end{array}\text{\hspace{1em}(57)}$$

切向刚度参数 $k_t$ 的值原则上与法向刚度参数无关，但通常假设法向刚度与切向刚度成一定比例 $\beta$：

$$\begin{array}{r}\beta =\frac{k_t}{k_n}\end{array}\text{\hspace{1em}(58)}$$

$\beta$ 值非常重要，因为它对离散元模型的宏观行为有很大的影响。如杨氏模量或泊松比可以表示为 $k_t/k_n$ 比值的函数（Marczewska 等人. 2016）。

阻尼系数 $c_n$ 与系统的临界阻尼 $C_{cr}$ 有关：

$$\begin{array}{r}\zeta =\frac{c_n}{C_{cr}}\end{array}\text{\hspace{1em}(59)}$$

其中，$\zeta$ 被称为阻尼比。它是一个非负的无量纲参数（$\zeta \ge 0$）。零阻尼比 $\zeta =0$，表示无阻尼，$0<\zeta <1$ 表示欠阻尼，$\zeta =1$ 表示临界阻尼， $\zeta >1$ 表示过阻尼。对于质量为 $m_i$ 和 $m_j$ 的两个刚体组成的系统，其临界阻尼 $C_{cr}$，与相连弹簧的刚度 $k_n$ 有关（Taylor 和 Preece 1992）。

$$\begin{array}{r}{C}_{cr}=2\sqrt{m^{\ast }{k}_n}\end{array}\text{\hspace{1em}(60)}$$

其中，有效质量 $m^{\ast }$ 定义为：

$$\begin{array}{r}\frac{1}{m^{\ast }}=\frac{1}{m_i}+\frac{1}{m_j}\end{array}\text{\hspace{1em}(61)}$$

阻尼比 $\zeta$ 可以用恢复系数（COR）$e$ 来表示 （Nagurka 和 Huang 2006）：

$$\begin{array}{r}\zeta =-\frac{\ln e}{\sqrt{{\pi }^2+(\ln e{)}^2}}\end{array}\text{\hspace{1em}(62)}$$

其中，法向恢复系数 $e$ 是反映碰撞时物体变形恢复能力的参数，它只与碰撞物体的材料有关。其定义为碰撞前后两物体接触点的法向相对分离速度与法向相对接近速度之比：

$$\begin{array}{r}e=\frac{|v_{rn}^{end}|}{|v_{rn}^0|}\end{array}\text{\hspace{1em}(63)}$$

阻尼比 $\zeta$ 与恢复系数 $e$ 的关系如图（10）所示。

在离散单元法中，粘性阻尼作为一种机制，允许在颗粒碰撞中耗散能量，实现系统在动态载荷下作出不同响应，包括选取当适当阻尼与缓慢施加的载荷相结合时的准静态响应。然而，必须注意的是，在接触模型中引入粘性阻尼会导致在某些方面上与实际情况相比存在不一致性，这将在下面的第一个数值示例中讨论。

4.1.3 数值分析示例

4.1.3.1 给定初始速度下两个球的碰撞分析

采用粘弹性 Kelvin-Voigt 接触模型，模拟半径为 $10mm$ 的两个相同球的碰撞过程，其初始速度大小相等，为 $v=10m/s$，密度 $\rho =8000kg/m^3$，取杨氏模量 $E=200GPa$，系数 $\alpha =0.04$，接触刚度 $k_n$ 根据式（57）计算得到 $k_n=1.26\times 10^{8}N/m$。研究在不同的恢复系数 $e=0.1、0.5、0.8、1.0$ 和对应不同阻尼下的碰撞响应，其中 $e=1$ 对应理想弹性碰撞。

对于不同的阻尼，在碰撞过程中，球之间的间隙和其中一个球的速度被绘制成时间的函数，如图11所示。碰撞的持续时间对应于间隙为负值的时间间隔（见图11a）。从图11b中可以看出，碰撞后球体相互反弹，其速度与阻尼有关。在弹性碰撞中，球在碰撞后的速度与碰撞前的速度相同。在非弹性碰撞中，球在碰撞后的速度比碰撞前的速度要低。阻尼越高（即 COR 越低），回弹速度越低。

在不同阻尼下，总接触力及其分量的演化过程如图12所示。弹力、阻尼力和总接触力分别在图12a、b和c中绘制为时间的函数。通过对比图11和图12可以看出，弹力与间隙成正比，符合（38）的条件，阻尼力与速度成正比，符合（39）的条件。从图12b中可以看出，阻尼力在碰撞开始时具有一定的非零值，在碰撞结束时具有非零值。在碰撞前后，阻尼力为零的情况下，会导致阻尼元件在分析期间的不连续。由于粘滞阻尼力的不连续，图12c中显示的总接触力 $F_n^{cont}$ 在碰撞开始和结束时也是不连续的，而实际中接触力是连续的。

此外，由于粘性阻尼，在碰撞的最后阶段 $g<0,v_{rn}>0$，总接触力 $F_n^{cont}$ 为内聚力，而无粘结颗粒系统中的接触力始终是排斥的。

如 Hunt 和 Crossley 在1975年所提出的，用适当的非线性弹簧和阻尼器替换线性弹簧和阻尼器，在一定程度上可以缓解线性粘弹性接触模型的不一致性。

不同阻尼下总接触力与间隙的关系曲线如图12d所示。在理想碰撞情况下（COR=1），装卸力与间隙关系一致。在阻尼碰撞中，加载和卸载曲线不重合，形成一个迟滞回线，回线围成的面积作为碰撞过程中能量耗散的测量标准（Lin 和 Hui 2002）。

4.1.3.2 阶跃载荷下两个球的接触分析

使用线性粘弹性 Kelvin-Voigt 模型模拟了在阶跃压缩载荷为 $F=20kN$ 下两个半径为 $R=10mm$ 的质量相等的两个相同球的接触过程。假设材料特性与前一示例中相同，初始间隙和初始速度均为零，研究在不同的恢复系数 $e=0.05、0.5、0.8、1$ 和对应不同阻尼下的接触响应。

图13分别为在不同阻尼下球之间的间隙、其中一个球的速度和总接触力的时间响应。可以看出，在零阻尼（COR = 1）系统中，球的振荡周期与理论值完全一致。

$$\begin{array}{r}T=2\pi \sqrt{\frac{m^{\ast }}{k_n}}=72.6\mu s\end{array}\text{\hspace{1em}(64)}$$

间隙振荡的均值（实际上是重叠振荡）与静载荷作用下的静平衡态一致

$$\begin{array}{r}g=\frac{F_n}{k_n}=-159.1\mu m\end{array}\text{\hspace{1em}(65)}$$

通过计算得到的间隙，使得阻尼系统的振动衰减，达到了准静态平衡状态。

这表明可以使用动态公式来解决静态问题，这是离散元法（Rojek 等人. 2013）和有限元法（Joldes 等人. 2011）中采用的动态松弛法的基本原理。在用动态松弛法求解的弹性线性问题中，瞬态周期的解并不重要，不同的阻尼值使我们可以得到相同的静态解。动态松弛方法也可以谨慎地应用于路径相关问题。

4.1.3.3 线性递增载荷下两个球的接触分析

与前面的例子相同，初始条件相同，在从 $0$ 到 $0.5s$ 的时间间隔内，承受从 $0$ 线性增加到 $20kN$ 的压缩载荷。然后，对于 $t>0.5s$ 的时间段，载荷保持不变。采用线性粘弹性 Kelvin-Voigt 模型，研究在不同的恢复系数 $e=0.1、0.5、1$ 和对应不同阻尼下的接触响应。

图14显示了在不同阻尼下两球之间的间隙、其中一个球的速度和总接触力的时间响应。当阻尼为零（COR = 1）时，解具有振动的特征。振动在有阻尼的解中衰减。对于足够高的阻尼水平（COR 取较小值），间隙和接触力的时间响应实际上是线性的，与施加力的线性增加一致。在图14中可以观察到接触力与所施加的力完全一致，这证实了加载和响应可以被认为是准静态的。这证明了以增量形式再现准静态条件的可能性（参见图15），这对于分析非线性和路径相关问题具有十分重要的意义。

4.2 粘弹性 Hertz-Mindlin-Deresiewicz 接触模型

4.2.1 模型的建立

该模型将法向相互作用的 Hertz 粘弹性接触模型与切向摩擦的 Mindlin-Deresiewicz 接触模型相结合。模型的流变方案与图8相似，不同之处在于将线性弹簧替换为非线性弹簧。

Hertz 模型基于弹性球体之间接触问题的解析解，采用非线性关系评估弹性接触力（Hertz 1882；Johnson 1985）：

$$\begin{array}{r}{F}_n^e=-K_{nHz}{h}^{\frac{3}{2}}\end{array}\text{\hspace{1em}(64)}$$

其中，$h(h=-g)$ 为颗粒间的重叠量，接触刚度参数 $K_{nHz}$ 计算公式如下：

$$\begin{array}{r}{K}_{nHz}=\frac{4}{3}{E}^{\ast }\sqrt{R^{\ast }}\end{array}\text{\hspace{1em}(65)}$$

其中，$E^{\ast }$ 是由两个接触颗粒的杨氏模量 $E_i$ 和 $E_j$ 以及泊松比 ${\nu }_i$ 和 ${\nu }_j$ 定义的有效弹性模量：

$$\begin{array}{r}\frac{1}{E^{\ast }}=\frac{1-{\nu }_i^2}{E_i}+\frac{1-{\nu }_j^2}{E_j}\end{array}\text{\hspace{1em}(66)}$$

和 $R^{\ast }$ 是由等式（56）定义的有效半径。请注意，式（64）中根据重叠量 $h$ 而不是间隙 $g$ 定义了接触力，并且引入了负号，以便与之前研究者所使用的符号惯例保持一致，将压缩接触力视为负值。

在DEM的框架中，通常将粘性阻尼添加到弹性 Hertz 力中，以便在颗粒碰撞时耗散能量。在计算中有时采用公式（39）给出的线性阻尼项，参见 Lee（1994）。然而，最新的模型在赫兹弹性接触中使用非线性阻尼项。Hunt 和 Crossley（1975）推导了以下一般形式的非线性阻尼：

$$\begin{array}{r}{F}_n^d={\eta }_n{h}^p{v}_{rn}^q\end{array}\text{\hspace{1em}(67)}$$

Tsuji 等人（1992）提出了上述关于 $p=1/4$ 和 $q=1$ 的阻尼项：

$$\begin{array}{r}{F}_n^d={\eta }_n{h}^{\frac{1}{4}}{v}_{rn}\end{array}\text{\hspace{1em}(68)}$$

式（68）中使用的阻尼耗散系数 ${\eta }_n$ 与非线性弹簧刚度 $K_{nHz}$ 和恢复系数 $e$ 有关，如下所示（Navarro 和 de Souza Braun 2013）：

$$\begin{array}{r}{\eta }_n=\sqrt{5}\sqrt{m^{\ast }{K}_{nHz}}\frac{\ln e}{\sqrt{{\pi }^2+(\ln e{)}^2}}\end{array}\text{\hspace{1em}(69)}$$

Jankowski（2006）推导了非线性粘弹性接触模型中阻尼比和恢复系数之间的不同解析关系。改善 Hertz 弹性接触模型中粘性阻尼的可能性仍在研究中（Zdanceviˇcius 等人. 2017）。

Hertz 法向接触模型通常与 Mindlin 和 Deresiewicz （1953）提出的切向接触模型相结合。完全采用 Mindlin-Deresiewicz 理论会产生复杂的算法，参考（Renzo 和 Maio 2004；kruggel-Emdenl 等人. 2008），因此后人提出了不同的简化方法。Tsuji 等人（1992）取 Mindlin 和 Deresiewicz 的法向力的恒定值，推导出切向力的公式：

$$\begin{array}{r}{F}_t=k_t{u}_{rt}\end{array}\text{\hspace{1em}(70)}$$

其中，切向刚度 $k_t$ 的计算方法如下：

$$\begin{array}{r}{k}_t=8G^{\ast }\sqrt{R^{\ast }h}\end{array}\text{\hspace{1em}(71)}$$

其中，$G^{\ast }$ 是由两个接触颗粒的剪切模量 $G_i$ 和 $G_j$ 以及泊松比 ${\nu }_i$ 和 ${\nu }_j$ 定义的有效剪切模量：

$$\begin{array}{r}\frac{1}{G^{\ast }}=\frac{1-{\nu }_i^2}{G_i}+\frac{1-{\nu }_j^2}{G_j}\end{array}\text{\hspace{1em}(72)}$$

有效半径 $R^{\ast }$ 由（56）得到。$h$ 为粒子重叠量，接触点处的相对切向位移 $u_{rt}$ 通过相对切向速度的积分得到：

$$\begin{array}{r}{u}_{rt}={\int }_0^t{\mathbf{v}}_{rt}\mathrm{d}t\end{array}\text{\hspace{1em}(73)}$$

切向力 $F_t$ 受到库仑条件的限制：

$$\begin{array}{r}{F}_t\le \mu |F_n|\end{array}\text{\hspace{1em}(74)}$$

为了使得模型与 Mindlin-Deresiewicz 理论的具有一致性，Renzo 和 Maio（2004年）；Maio 和 Renzo（2005年）对 Tsuji 等人（1992年）的模型进行改进，包括将公式（71）给出的刚度按系数2/3进行缩放。

$$\begin{array}{r}{F}_t\le \mu |F_n|\end{array}\text{\hspace{1em}(75)}$$

4.2.2 数值分析示例

4.2.2.1 给定初始速度下两个球的碰撞分析

基于非线性 Hertz 弹性分量和阻尼分量的粘弹性 Kelvin-Voigt 接触模型，（68）和（69）被用来模拟先前用线性 Kelvin-Voigt 接触模型分析过的两个相同球的碰撞。假设球尺寸相同，半径 $R=10mm$，密度$\rho =8000kg/m^3$，杨氏模量 $E=200GPa$ 和初始速度 $v=10m/s$，泊松比 $\nu =0.3$，研究在不同的恢复系数 $e=0.1、0.5、0.8、1.0$ 和对应不同阻尼下的碰撞响应。

图16（a）和（b）分别绘制了不同阻尼下球之间的间隙和一个球的速度随时间的变化。可以看出，COR越低，冲击时间越长，反弹速度也越低。

弹性、阻尼和总接触力分别在图17（a）（b）（c）中绘制为时间的函数。在图17（b）中可以观察到，与图12（b）中显示的线性 Kelvin–Voigt 接触模型中的阻尼力相反，本模型中的阻尼接触力不再是不连续的。因此，图17（c）中的总接触力也不是不连续的。线性Kelvin-Voigt模型的另一个缺陷，表现在碰撞最后阶段相互作用中，在图17（c）（d）注意到非线性 Kelvin-Voigt 模型中没有消除内聚力。图12（d）显示了不同阻尼力与重叠量（负间隙）的关系曲线。可以看出，与线性模型类似，理想弹性情况下的加载和卸载过程中力与间隙的关系是一致的，并且在阻尼碰撞中，加载和卸载曲线形成滞后回路。但本模型与线性模型不同，加载和卸载函数是零间隙的平滑曲线。

4.3 Walton-Braun 弹性模型

4.3.1 模型的建立

粘性阻尼可以看作是颗粒碰撞过程中非弹性变形的一种机制。塑性变形取决于位移类型的变量，因此，与速率无关的滞回接触模型，如 Walton 和 Braun（1986）提出的模型，似乎更适合模拟与塑性变形相关的碰撞。Walton–Braun模型假设了一种线性力-重叠关系，但卸载坡度（刚度）高于加载坡度（刚度），当力降到零时，这导致一定的残余不可逆重叠。这使我们能够将该模型视为具有弹性卸载的弹塑性模型。力作为粒子重叠的函数绘制在图18中。请注意，尽管该图位于图的第一象限，但接触力符号（压缩接触力-负数）的约定是通过取垂直轴的力的负数来保持的。力的计算公式为：