如图2-1所示，将宏观的连续型材料理想化为由一系列颗粒组成的微观离散结构系统。在这一微观离散体中，每一颗粒的质心代表一个材料点。当这种材料模型在受到外力作用时，理想化的颗粒会产生平移或旋转。在三维空间，每个颗粒包含6个自由度，即3个平移自由度和3个旋转自由度。因此，两个相邻颗粒n和p间的相对位移δi可表示为式（2-1）（Misra et al，1993）：

图1-2　连续体理想化为颗粒模型

（a）虚颗粒结构；（b）两颗粒的运动学；（c）颗粒间的相互作用示意图

式中：ui为颗粒位移；ωj为颗粒旋转；rk为连接颗粒质心与颗粒间接触点的向量；eijk为置换符号。

上标代表相互作用的颗粒；此处及本专著中以后出现的所有下标如未加特殊说明，均表示空间坐标，且遵循张量求和法则。

图2-2　颗粒接触点的局部坐标系

从严格意义上讲，反映颗粒间相互作用的虚拟键势能受其周围存在的许多颗粒的影响（Israelachvili，1985）。但为了简单起见，本专著在此借鉴原子间相互作用的势能形式定义虚拟键势能为中心间距L和切向畸变w的函数，用函数Φ（L，w）来表示。Volokh & Gao（2005）在修正的VIB模型中将势能分解为膨胀变形和剪切变形两个分量。在此，我们将势能函数分解为中心势能φ（L）与非中心势能φ（w）之和，如式（2-3）。其中，φ（L）依赖于两相邻颗粒间的中心间距L；φ（w）则依赖于两相邻颗粒相对于初始平衡位置时的切向畸变w。

此处忽略了中心运动与非中心运动的交叉影响。另外，p={pL，pw}是包含颗粒尺寸和形状的参数组合。

按照自然晶体界面的形式将虚颗粒间的界面划分为如图2-1（c）所示的结构。如果将这些颗粒看作是大量原子的聚合体，那么从统计意义上讲，颗粒间的势能代表的就是构成颗粒的大量相互作用的原子的行为。颗粒间的势能原则上可以按照类似于分子势能函数的模式进行推导（Murrell et al，1984），然而这种方法却不适用于由数量庞大（数十亿）的原子构成的大颗粒。

从理论上讲，可以构造一些函数来反映亚颗粒尺度的基本作用机理。这些函数的相关参数则可以通过对现有的实验数据进行拟合求得。鉴于这种方法的实用性，本专著采用了这种思路并以两颗粒间的相互作用为例来进行具体研究。图2-3（a）给出了虚拟键间中心势能函数φ（L）的一种典型的分布趋势图。当两颗相互作用的颗粒沿中心轴方向被逐渐拉开时，中心势能函数φ（L）会从某一平衡状态值逐渐衰减到零。从现象学角度分析，这种衰减反映的是亚颗粒尺度的黏聚力消失以及孔洞扩大造成的累积损伤过程（Misra et al，2007）。同理，当这两个颗粒受到相互挤压作用时，颗粒的中心间距会不断缩短，此时虚拟键势能函数揭示了由亚颗粒尺度的剪切局部化或孔洞崩溃造成的损伤过程。van Mier（2007）讨论了水泥和岩石这两种应变软化材料的颗粒间可能存在的可压缩性作用机制。在同等受压条件下，随着原子中心间距的递减，原子势能会迅速地趋近于无穷大，相比之下颗粒间的虚势能函数的变化趋势显得异常平缓。在受压（或排斥）时，微粒（如电子或质子）间或多原子的势能函数则比对应的原子间势能函数的变化趋势要平缓得多（Murrell et al，1984；Israelachvili，1985）。根据颗粒间的势能函数，可以将两相邻颗粒的相对运动产生的作用力表示为颗粒间势能函数的导数。为方便起见，此处直接给出了力的函数表达式（Misra et al，2010）。式（2-4）和式（2-5）分别表示两相邻颗粒在拉伸和压缩条件下产生的法向作用力：

图2-3　虚势能函数曲线

（a）法线方向；（b）切线方向

式中：A、B和α1为模型参数；下标n为法向分量。

式（2-4）表征了基于Morse势函数的一种特征性的损伤行为。此前Gao & Klein（1998）将其用来描述在受拉条件下的颗粒间因凝聚力消失出现的损伤。当两颗粒受纯粹的中心轴方向压缩作用时，亚颗粒尺度的孔洞收缩或剪切局部化会迫使法向力函数趋近于一个渐进值而非无限制的递增。式（2-5）构造的便是一个趋近于渐进值的双曲线型的力函数。除了式（2-5）外，在压缩条件下根据亚颗粒的结构和组成的不同，力函数也存在其他可能的表达式，但这些表达式尚未得到验证，需要进行深入研究，在此不作讨论。图2-4（a）绘出了在拉伸和压缩条件下式（2-4）和式（2-5）对应的法向力函数随颗粒间法向相对位移的变化曲线。(www.daowen.com)

图2-4　力函数随颗粒间相对位移的变化曲线

（a）法向分量；（b）切向分量

当颗粒在沿虚拟键的切线方向发生形变时，同样需要考虑亚颗粒尺度的变形机制。例如，压缩条件下引起沿切线方向的损伤破坏所需的势能往往比在拉伸条件下所需的势能高。这种现象可能是由拉伸和压缩时颗粒界面间的相互作用差异造成的。如果亚颗粒尺度的颗粒界面存在大量断裂和缺陷，那么对颗粒间在法向和切向共同受力的情况下的损伤过程的研究就会变得极其复杂。然而，在不同的法向加载条件下，颗粒间剪切变形的差异可以用因剪切产生的单一断裂显现出的压力相关行为来描述（Misra，2002）。因此，在受压条件下非中心键势能函数φ（w）呈现出较深的“势能井”。图2-3（b）给出了在颗粒受拉伸或压缩条件下一种可能的非中心势能函数φ（w）随切向畸变的变化趋势。非中心势能函数的变化趋势类似于中心势能函数曲线的拉伸部分[图2-3（a）]，随着切向扭曲变形从某一平衡值开始递增，颗粒间相互作用产生的能量逐渐地衰减到零。在图2-3中，势能最低点分别对应的是平衡状态下的中心间距L=L0和切向变形w=0。为简便起见，同样直接给出了构造的由相邻颗粒切线方向的相对运动产生的力的函数表达式：

式中：下标w表示切向分量。

以上两式与法向力的表达式（2-4）形式基本相同。为了反映出加载情况的影响，在切向力的函数中分别采用了不同的模型参数C、D和α2用以区分拉伸和压缩两种不同的受力情况。除这里提供的力的函数形式以外，还可以对这些力的函数加以修改，比如为了使其显示出对法向力的依赖性可以引入摩擦定律，但其他形式的力函数都有待进一步的研究方可使用。图2-4（b）绘出了以上形式的切向力函数随颗粒间的切向相对位移分量的变化趋势。从图2-4中可以看出：无论在拉伸还是压缩情况下，切向力都显现出了峰值前非线性以及峰值后的应变软化现象，这与在剪切力作用下颗粒界面会出现损伤的情况相符。

根据之前给出的法向力和切向力的函数表达式，可以求得法向刚度分量Kn和切向刚度分量Kw的表达式如下：

图2-5是这些刚度分量分别随法向和切向相对位移的变化趋势图。需要说明的是，之前所构造的在拉伸和压缩两种情况下的力的分量的函数都假设了刚度在初始零应变情况下是相同的（δn=0，δw=0）。

图2-5　刚度分量随相对位移分量的变化趋势

（a）法向；（b）切向