孔隙空间分布对粗粒土强度变形特性影响研究

胡有方,袁俊平，卢 毅

(1. 河海大学 岩土力学与堤坝工程教育部重点实验室，江苏 南京 210098；2. 河海大学 江苏省岩土工程技术工程研究中心，江苏 南京 210098)

粗粒土是一种典型的颗粒材料，由大量粒径不等、形状各异、排列随机的土石颗粒组成，微观上孔隙结构的多样性决定了其宏观物理力学性质的复杂性[1]。因此可以说，粗粒土的变形过程实质上就是其孔隙结构变化的过程，如果能获得粗粒土的孔隙结构特征和变化规律，也就能掌握其变形规律[2]。近年来，许多国内外学者对粗粒土的孔隙结构特征进行了研究，主要方法有实验法和数值重构法[3-4]。实验方法最早起源于1982年，Petrovic等[5]首先在土壤的密度和含水率研究中引入CT扫描技术，之后，Warner[6]则在对CT扫描图像的分析中发现，孔隙的位置、大小和数量信息都能由图像精确地揭示。在国内，学者吕菲[7]等在扫描图像的基础上建立了孔隙网络模型，并成功预测了饱和土壤的水力学性质。相方园[8]使用切片法进行试验，将孔隙拓扑结构与形态学理论相结合，定量地获得了孔隙结构中的特征参数。与实验方法不同，数值重构法的思路是基于二维孔隙信息和不同的数学模拟算法建立三维孔隙结构模型，目前常用的方法有高斯随机场法[9]、过程法[10]和模拟退火法[11]。比较成功的案例有Pilotti[12]基于过程法，用一种主要适用于球体、椭球体、圆柱体和平行六面体的算法建立了任意形状颗粒的随机堆积体。曾建邦[13]等利用改进的模拟退火法重建了含水沉积物的三维孔隙结构，并通过孔隙结构的特征化分析揭示了分布模式对沉积物特性的影响。

综上所述，目前的研究主要针对的是粗粒土孔隙结构特征的描述方法，因此，如何将粗粒土的孔隙结构特征与强度变形特性定量地联系起来，仍是有待探讨的问题。本文在前人的研究基础上，从孔隙空间分布的不均匀性出发进行研究，定义了可以定量刻画粗粒土孔隙空间变异性的孔隙空间分布系数FSD；基于PFC3D离散元软件，分析了FSD与粗粒土强度变形参数之间的关系；同时，基于Abaqus有限元软件，分析了FSD在有限元计算中的适用性。

1 土体孔隙空间分布不均匀性的定量刻画

首先，粗粒土试样中各岩土颗粒和孔隙的位置可以通过扫描的方法得到[14]。接下来，使用最大球算法[15]对孔隙进行模拟，即对于每一部分孔隙，插入若干个以岩土颗粒为边界的最大内切球，使这些球充满孔隙空间。以各球体的质心作为孔隙的质心，以各球体的体积作为孔隙的体积，从而达到以大量规则的球体模拟不规则的孔隙体的目的，最大球算法的示意图如图1所示。

图1 最大球算法示意图Fig.1 Schematic diagram of maximum sphere algorithm

下一步，将土体划分为N个相同的区域对孔隙球体进行统计，由于粗粒土的不均匀性，每个区域内的孔隙数量和孔隙体积分布都不相同，因此引入孔隙数量变异系数CVn和孔隙体积变异系数CVv，两者的表达式分别为

(1)

(2)

CVn和CVv可以较好地将试样中孔隙数量的分布情况和孔隙体的分布情况反映出来，其值越小，则代表均匀性越好。理想情况下，孔隙分布完全均匀时，两者值为0。但要注意的是，这两个参数本身无法反映出孔隙质心点在空间中的分布情况，因为存在偏聚现象，如图2所示，a试样与b试样的CVn和CVv相同，但显然a试样的偏聚程度小于b试样。因此，还需要引入偏聚系数β来反映这种情况。

图2 孔隙空间分布状态示意图Fig.2 Diagram of pore space distribution

(3)

式中：n表示样本区域内孔隙总数量。

接着，用x0表示试样中孔隙完全规则地分布在试样空间中时孔隙之间的距离，然而对于任意的n值，x0的值都难以计算，因此本文用下式代替：

(4)

式中：v表示样本区域内孔隙总体积。如此计算时，x0表示当n个相同尺寸的球体孔隙总体积与样本区域内孔隙总体积相等时这些孔隙的直径。

于是，偏聚系数β可以用下式计算：

(5)

在得到了孔隙数量变异系数CVn，孔隙体积变异系数CVv以及偏聚系数β之后，本文定义孔隙空间分布系数FSD作为综合性指标，从孔隙质心点空间分布的均匀性、孔隙体积空间分布的均匀性以及孔隙质心在空间上的偏聚程度三个方面综合反映孔隙空间的分布情况，其表达式定义为

FSD=β·CVn·CVv

(6)

2 孔隙空间变异系数与土体强度变形参数关系研究

在定义了孔隙空间分布系数FSD之后，为了研究其对土体宏观强度变形特性的具体影响，就需要定量分析FSD与土体各项强度参数之间的关系。因此本节中基于PFC3D软件，建立了5组三轴压缩试验试样进行数值模拟实验。试样编号为PK1—PK5，高度H=200 mm，直径D=101 mm。这些试样的总孔隙率相同，但各试样的内部分层情况不同，详细参数和示意图如图3所示。

试样模型按照以下步骤建立：(1)将试样空间按照试验方案均分为指定层数，从最底层开始生成颗粒；(2)使用半径扩大法生成试验方案中指定孔隙率的土颗粒，在每层顶部多预留5 cm的高度；(3)赋予土颗粒重力加速度，使土颗粒在重力作用下自由下落形成堆积体；(4)赋予该层顶部墙体一定速度使其下落直至获得该层指定尺寸的试样，运行一定时步使试样达到稳定状态，该层试样生成完毕；(5)重复(2)—(4)步骤直至所有层的试样均生成完毕，删除中间墙体，运行一定时步使试样达到稳定状态，试样生成完毕。模型生成过程如图4所示。

完成试样的生成后，为了防止偶然误差给数值试验结果带来影响，对各试样进行多次对照试验，结果取平均值。最终统计各试样的孔隙空间分布系数FSD如表1所示。

注：数据为各层孔隙率。图3 不同孔隙空间分布影响效应试验方案示意图Fig.3 Schematic diagram of test scheme for effect of different pore space distribution

图4 模型生成过程示意图Fig.4 Diagram of model generation process

从表1可以看出，在其他颗粒细观参数相同的条件下，各试样的孔隙空间分布系数因分层的不同而发生了明显的变化。试样PK5相对于试样PK1，空间分布系数FSD增大了46.8%，说明该方法的合理性。

为了得到上述各试样的应力和应变曲线，对表中的各组试样进行围压800 kPa的三轴固结排水试验数值模拟，结果如图5和图6所示。

从图5和图6中可以看出，在其他颗粒细观参数相同的条件下，土体的弹性模量、峰值强度和泊松比都随着孔隙空间分布系数FSD的增大而减小。试样PK5相对于试样PK1，弹性模量减小40.8%，峰值强度降低12.5%，泊松比减小39.5%。

表1 各试样孔隙空间分布系数

图5 偏应力-轴向应变曲线Fig.5 Relationship between deviatoric stress and axial strain

图6 侧向应变-轴向应变曲线Fig.6 Relationship between lateral strain and axial strain

由于加载刚进行时，试样的偏应力-轴向应变曲线近似为一条直线，因此将试样轴向应变达到1%时的割线模量作为试样的弹性模量。相关性分析结果表明，弹性模量E与孔隙空间分布系数FSD之间呈指数函数关系，如式(7)。

E=A1+B1er1FSD

(7)

用指数函数对试样弹性模量E和孔隙空间分布系数FSD进行拟合，结果如图7所示，其中A1的值为94.062 2，B1的值为-0.609 4，r1的值为-2 071.968 1。

将轴向应变1%时的切线泊松比作为试样的泊松比，相关性分析结果表明，泊松比ν与孔隙空间分布系数FSD呈指数函数关系，如式(8)。

v=A2+B2er2FSD

(8)

图7 弹性模量随空间分布系数变化曲线Fig.7 Relationship between elastic modulus and pore space distribution coefficient

图8 泊松比随空间分布系数变化曲线Fig.8 Relationship between Poisson′s ratio and pore space distribution coefficient

用指数函数对二者进行拟合，结果如图8所示，其中A2的值为0.034 3，B2的值为-21.302 8，r2的值为-3 967.103 5。

为了得到各试样的内摩擦角，在围压200、800和2 000 kPa的条件下，对各试样进行三轴试验数值模拟，得到了各试样在不同围压下的峰值强度。根据莫尔-库伦强度理论，各试样的内摩擦角如表2所示。

表2 各试样内摩擦角

为了研究试样内摩擦角与孔隙空间分布系数之间的定量关系，对内摩擦角与孔隙空间分布系数进行相关性分析，发现两者之间呈指数函数关系，如式(9)。

φ=A3+B3er3FSD

(9)

用指数函数对试样内摩擦角和空间分布系数进行拟合，结果如图9所示，其中A3的值为38.469 2，B3的值为-0.156 7，r3的值为1 705.456 1。

3 孔隙空间分布系数的准确性验证与分析

本节中使用Abaqus有限元软件对孔隙空间分布系数FSD计算方法的准确性进行验证。相比于离散元方法，有限元方法更适合大规模的工程计算，这是由于有限元将粗粒土复杂的几何结构简化为了具有简单形状的单元，单元内的材料性质和控制方程通过单元节点的未知量来进行表达，从而使得计算的效率大大提高，然而这也使得有限元计算时忽略或低估了孔隙空间分布系数FSD的影响。那么孔隙空间分布系数FSD是否可以被忽略，是否可以使用Abaqus有限元软件进行孔隙空间变异性的模拟，这是本节中要讨论的问题。

为了解决此问题，建立编号分别为A1—A5的5种三轴固结排水实验试样，试样的高度、直径均与PK1—PK5相同。同时，通过上节中的结果推算出试样的弹性模量、泊松比、内摩擦角，使得两组试样的等效强度参数全部相同，具体数值如表3所示。最后，设置试样的孔隙率均为35%，使其为均匀试样。

表3 A1—A5试样强度参数

对试样施加800 kPa的围压时，各试样的偏应力-轴向应变曲线如图10所示。

图10 试样A1—A5偏应力-轴向应变曲线Fig.10 Relationship between deviatoric stress and axial strain of sample A1—A5

从图10中可以看出，由于试样A1—A5忽略了孔隙空间分布系数FSD，消除了空间分布差异对强度的不利影响，使得虽然试样的各项强度参数均与离散元计算时相同，但试样破坏时的偏应力均偏大，并且随着离散元试样中FSD的增加，这种差异更加明显。试样A1相比试样PK1，破坏时的偏应力增大了1.6%，试样A5相比试样PK5，破坏时的偏应力增大了12.0%。

为了降低有限元单元的均匀性，增加孔隙空间分布系数的影响，采用上节图3中的分层方法建立5种三轴固结排水实验试样B1—B5，在前者的基础上额外考虑不同土层孔隙率的差异。同样对试样施加800 kPa的围压时，各试样的偏应力-轴向应变曲线如图11所示。以试样5为例，将三次试验的结果进行比较，如图12所示。

图11 试样B1—B5偏应力-轴向应变曲线Fig.11 Relationship between deviatoric stress and axial strain of sample B1—B5

图12 试样A5,B5,PK5偏应力-轴向应变曲线Fig.12 Relationship between deviatoric stress and axial strain of sample A5,B5,PK5

从图11中可以看出，对于考虑分层的有限元试样，其峰值强度相比均匀的有限元试样有所降低，但仍然比离散元试样大。这是因为，虽然分层增加了层与层之间的孔隙空间分布差异，但单独每个层内的孔隙仍然是均匀的，其FSD虽然大于0，但仍小于离散元试样。

从图12中可以看出，对试样5来说，虽然三次试验的等效强度完全相同，但三次试验的峰值偏应力分别为2.74，3.07，2.99 MPa。再次验证了试样的峰值强度随着孔隙空间分布系数FSD的增大而减小的规律，并进一步论证了孔隙空间分布系数FSD计算方法的合理性。同时也说明，在试样离散型较大时，不应当忽略孔隙空间分布变异性的影响。

4 结论

1)孔隙空间分布系数FSD可以较好地模拟粗粒土中孔隙空间分布的不均匀性。这种不均匀性包括孔隙质心点空间分布的不均匀性、孔隙体积空间分布的不均匀性以及孔隙质心在空间上不同的偏聚程度，因此FSD是一项比较系统的综合性指标。

2)孔隙空间分布系数FSD在离散元和有限元分析中均能得到较好的运用。在其他颗粒细观参数相同的条件下，土体的弹性模量、峰值强度和泊松比都随着孔隙空间分布系数FSD的增大呈指数函数形式减小。