在任一代表性体积元V内应变能密度W可以表示成：

式中：N表示代表性体积元V内颗粒间接触点的总数量。

依据Chang & Liao（1990）的思路，把连续位移场用泰勒级数式展开。取代表性体积元V的中心点x 0作为参考点，于是颗粒n的连续型位移可以用参考点x 0处的一系列位移梯度来预测：

式（2-16）忽略了四阶及其以上的导数项，下标中的逗号表示对于其后所跟相应的空间坐标的偏导数。在同一代表性体积元内，ui（x 0）、ui，j（x 0）、ui，jk（x 0）和ui，jkq（x 0）都是常数。二阶张量ui，j有9个分量，然而由于忽略了颗粒间的相对旋转，即假设非极性的连续体，于是9个分量减少到了6个独立分量。三阶张量ui，jk关于最后两个下标对称且含有18个独立分量。四阶张量ui，jkq关于最后3个下标对称且含有30个独立分量。将式（2-16）代入到式（2-1）中可以得到：

其中，几何张量分别为：

假设颗粒p位于局部坐标系的原点处并且虚颗粒系统内所有的颗粒尺寸大小一致，半径为r，式（2-18）～式（2-20）则可以简化为：

将式（2-8）和式（2-17）代入到应变能密度表达式（2-14）中得到：

分别定义零阶、一阶及二阶（上标表示阶数）应变张量如下：

由此可以得到对应的共轭应力张量如下：

联立式（2-17）、式（2-21）～式（2-26）可以得出如下的本构关系：(www.daowen.com)

其中，本构张量用几何量的形式表示如下：

由于引入了高阶应力张量，计算会变得异常复杂。为简单起见，本专著仅考虑零阶和一阶应力、应变张量，略去二阶张量。根据文献（Chang et al，2002），可假设研究的材料微观结构为中心对称型，于是可以得到Bijqmn=0，因此本构方程（2-27）和式（2-28）可以进一步简化为：

考虑到高阶应力和高阶应变张量的对称性，此处的四阶和六阶本构张量需要满足以下对称性条件：

由于代表性体积元通常是由数量庞大的颗粒聚合而成，要对这样一个单元内的所有颗粒集合体的任一张量求和，可以将其转换为积分形式。因此离散形式的高阶本构张量式（2-30）和式（2-32）可用积分形式重新表示如下：

式中：a20、a22和b22是表征材料对称性的结构参数。例如，a20=a22=b22=0表征各向同性材料；当a22=b22=0时，材料具有横向各向同性对称性；当a20≠a22≠0时，材料具有正交对称性。

为便于计算，本专著假设所研究的材料具有各向同性的微观结构，因此方向密度函数（2-42）进一步简化为：