0.   引言

在天然重力沉积或人工填筑过程中，粒状土体颗粒的排布方式及复杂的应力状态会导致其在刚度^[1]、强度^[2]、渗透性^[3]等力学特性上呈现明显的各向异性特征。弹性波速是评价土体性质的重要物理力学指标之一^[4-5]。压电陶瓷弯曲元能方便快速地测定土体弹性波速，其已被广泛安装在三轴仪^[6]和固结仪^[7]等土工试验设备中用于砂土弹性波速及其各向异性特性的研究。研究结果表明，在等向应力状态下，天然砂土水平向和竖向的波速明显不同，即呈现初始各向异性。例如：砂雨法制备的Ticino砂试样的水平向波速大于竖直向^[8-10]，而同样采用砂雨法制备的Virginia海滩砂试样的竖直向波速小于水平向^[11]。同时，砂土试样的制备方法对试样波速各向异性特征也有影响，例如：湿夯法制备的Hostun砂土试样的水平向波速小于竖直向，而振动法及砂雨法制备的试样水平方向波速则大于竖直向^[12]。另一方面，在不等向应力条件下，当竖向与水平向应力比(SR=${\sigma _{\text{v}}}/{\sigma _{\text{h}}}$) 较小时，Ticino砂试样不同方向的波速与该方向上的应力呈幂函数关系，而当应力比超过一定值后，沿小主应力方向传播的波速开始衰减，波速与应力间关系不再满足经验公式^[13]。Kuwano等^[10]认为波速的显著减低及各向异性程度的改变是由于颗粒试样发生剪胀，导致粒间接触分布发生变化，但其在实验中未能获取微观接触组构信息加以验证。针对这一试验现象，Gu等^[14]采用离散元数值模拟了模量（或波速）随应力状态的演化情况，发现波速的衰减与接触配位数随应力比的增加而降低有关，而波速各向异性程度与接触法向组构各向异性相关。Gu等^[15]还进一步建立了颗粒材料试样应力修正后的波速各向异性与接触法向组构各向异性的线性关系。

目前，从本质上探求颗粒材料复杂宏观特性的微细观机理，广泛采取的手段仍为数值模拟^[14-15]，通过试验获取试样微观组构的研究很少。然而，数值模拟需以试验结果作为模型建立的依据和验证，但对真实颗粒形状和试样微观结构的精细化模拟还十分困难。随着成像技术的发展，CT扫描已经成为一种无损伤表征的新技术^[16]，可通过其重构出颗粒土试样的三维构造，从而进一步获取试样的微细观结构。国内已有研究大多利用CT扫描研究岩石及砂土内部孔隙特征^[17-18]，而砂土微观接触组构信息的研究很少。国外已有研究中，Sun等^[19]基于CT图像计算了砂土试样的颗粒长轴方向、接触法向和支矢量等组构参数，同时发现颗粒形状会影响砂土的固有组构各向异性。Imseeh等^[20]和Wiebicke等^[21]利用CT扫描研究了剪切过程中砂土微观接触组构各向异性的演化，其与宏观应力应变曲线特性的变化一致。可见，颗粒材料微观组构对其宏观响应至关重要，但利用CT扫描等技术研究颗粒材料微观结构设备要求高、难度大、费时费力。如能建立微观组构和波速各向异性的机理联系，从而利用宏观波速定量和快速表征微观组构，具有重要的理论意义和实用价值。因此，从试验手段探究粒状土宏观波速各向异性的微细观机理十分必要，而目前未见宏微观相结合的试验研究。

综上所述，本文利用配置了两对弯曲元波速传感器的圆柱扭剪仪，从宏观上研究了不同应力条件下聚氯乙烯（PVC）颗粒试样中不同方向的压缩波（P波）及剪切波（S波）的波速变化情况，探讨了波速各向异性程度随应力状态的演化规律。在微观上，结合CT扫描图像重构试样的三维结构，并提取接触配位数、颗粒长轴方向、接触法向分布等微观组构参数，进而探讨颗粒材料宏观波速各向异性特征的微细观机理。

1.   试验材料、设备、方案

1.1   试验材料及设备

试验材料采用单一粒径聚氯乙烯（PVC）颗粒，其表面光滑，长轴均长约为4.0 mm，短轴均长约为2.1 mm，平均等效粒径约为2.85 mm。PVC材料的相对质量密度G_s=1.75，杨氏模量和泊松比分别为3.92 GPa和0.32。PVC颗粒的最大孔隙比e_max和最小孔隙比e_min分别为1.0和0.728。试样采用干燥PVC颗粒通过分层压实法制备，制样时目标孔隙比为0.730，为密实试样。试样为直径70 mm，高150 mm的实心圆柱体。

采用GDS圆柱扭剪仪控制试样的应力状态。同时，在扭剪仪中安装了竖向和水平向弯曲元，可开展多方向波速测试（见图 1），包括竖直面波速v_PZ和v_SZY，水平向波速v_PX和v_SXY，倾斜向波速v_Pθ和v_θY，其中θ为弹性波传播方向与试样竖向对称轴所成夹角，从而实现横观各向同性弹性模型参数的完整测定。

图  1  弯曲元测试中的试样及波速示意图

Figure  1.  Schematic diagram of specimen and wave velocity

下载: 全尺寸图片 幻灯片

当试样达到特定应力状态，对其充入环氧树脂进行固化，再利用CT扫描获得试样不同断面照片，并通过计算机重构试样三维结构及获取试样的微观组构信息。该CT设备及其参数分别见图 2和表 1。

图  2  CT扫描设备图

Figure  2.  X-ray computed tomography system

下载: 全尺寸图片 幻灯片

表  1  CT扫描成像系统参数

Table  1.  Settings of X-ray computed tomography system

细节辨识能力/μm	空间分辨率/(LP·mm-1)	像素尺寸/μm	像素/pixels	帧数(frames·s-1)
< 3	23.0	200×200	1024×1024	15

下载: 导出CSV 

| 显示表格

1.2   试样的不同应力状态

采用两种不同的应力固结路径，如图 3所示，包括等向压缩固结和不等向压缩固结。等向压缩固结时竖向应力与水平向应力相等，应力依次固结至50，75，100，150，200 kPa。不等向固结时，水平向应力保持为100 kPa不变，而竖向应力以间隔20 kPa逐渐从初始100 kPa增加到200 kPa。当竖向应力超过200 kPa后，试样发生破坏，无法维持平衡稳定状态。试样固结到指定应力状态后，先进行弯曲元波速测试，随后在保持应力状态条件下利用注射器通过试样的排水管缓慢充入环氧树脂进行固化，以避免对试样的颗粒状态产生扰动。

图  3  固结应力路径图

Figure  3.  Stress paths in consolidation

下载: 全尺寸图片 幻灯片

2.   宏观波速试验结果分析与讨论

2.1   信号时域分析与波速确定

图 4给出了围压100 kPa的PVC颗粒试样中弯曲元在不同激发频率下不同方向的P波及S波的输出信号。由图 4可知，输出信号频率先随激发频率的增加而增加，当激发信号频率高于10 kHz后，输出信号比较稳定，不随激发信号频率增加而明显改变，此现象与离散颗粒材料的高频过滤特性相关^[22]。由Yamashita等^[23]及Gu等^[24]的研究可知，时域初达波法不受输入、输出信号频率不一致的影响，能较准确地确定弹性波到达时间，因此本文利用时域初达波法确定P波和S波的初始到达时间，进而获取波速。

图  4  PVC颗粒试样中P波及S波信号

Figure  4.  P- and S-wave signals of PVC particle specimens

下载: 全尺寸图片 幻灯片

2.2   波速表征的横观各向同性弹性刚度矩阵

根据弹性波在连续介质中的传播定律，通过不同方向的P波和S波波速可完整确定横观各向同性弹性模型的5个独立参数^[25]：

式中：C₁₁为水平向的压缩模量M_h；C₃₃为竖向的压缩模量M_v；C₄₄为竖直面内的剪切模量G_vh；C₆₆为水平面内的剪切模量G_hh；C₁₂为C₁₁-2C₄₄；C₁₃为一独立弹性刚度系数，可由式（2）、（3）、（5）计算。

式中：ρ为试样密度；$\theta$为弹性波传播方向与试样竖直面对称轴的夹角。

当5个刚度矩阵系数确定后，即可根据式（3）~（5）计算任意方向或平面内的波速。图 5给出了等向围压100 kPa条件下竖直面内波速的分布图，其中散点为实测波速值，曲线为预测的波速。由图 5可知，倾斜向波速理论计算值与实测值吻合良好，且等向固结应力状态下P波及S波波速面受初始组构各向异性的影响较大，表现为近似椭圆状分布，其中水平向波速大于竖直向。

图  5  等向围压100 kPa下P波及S波的波速面

Figure  5.  Distributions of P- and S-wave velocities at an isotropic confining pressure of 100 kPa

下载: 全尺寸图片 幻灯片

2.3   等向围压大小对波速的影响

复杂应力状态下土体弹性波速的经验公式为^[26]：

式中：${\sigma _{{\text{pro}}}}$为弹性波传播方向上的应力；${n_{{\text{pro}}}}$为该方向上的应力指数；${\sigma _{{\text{osc}}}}$为颗粒振动方向上的应力，${n_{{\text{osc}}}}$为该方向上的应力指数；${\sigma _{{\text{out}}}}$为垂直于振动方向和传播方向所成平面的应力，${n_{{\text{out}}}}$为该方向上的应力指数；A为土体固有波速参数，反映了土体级配、颗粒形状等对波速的影响；F(e) 为土体孔隙比e的函数。图 6对比了PVC颗粒试样不同方向P波及S波波速随围压的变化情况。由图 6可知，P波及S波波速均与围压具备幂函数关系，应力指数n保持在0.21~0.22，符合已有波速与围压的经验关系公式^[8-10]。同时，PVC颗粒试样中竖向传播的波速值小于倾斜向，而水平向传播的波速值最大。由图 6拟合曲线所得竖向与水平向土体波速参数比值分别为A_PZ/A_PX = 0.85，A_SZY/A_SXY = 0.89，表现出明显的各向异性。

图  6  P波及S波波速与围压关系

Figure  6.  Relationship between P- and S-wave velocities and confining pressure

下载: 全尺寸图片 幻灯片

2.4   不等向固结应力对波速的影响

图 7给出了PVC颗粒试样在不等向应力状态下，不同方向P波及S波的波速情况。由图 7可见，${v_{{\text{P}}Z}}$，${v_{{\text{S}}ZY}}$，${v_{{\text{P}}\theta }}$，${v_{{\text{S}}\theta Y}}$起初均表现为随应力比SR(${\sigma _{\text{v}}}{\text{/}}{\sigma _{\text{h}}}$)增加而增加，与公式（7）预测一致。同时，${v_{{\text{P}}X}}$，${v_{{\text{S}}XY}}$随应力比SR的增加先保持不变，也符合理论上${v_{{\text{P}}X}}$，${v_{{\text{S}}XY}}$主要受水平向应力（本研究中保持不变）影响的预期。当应力比SR达到1.6后，${v_{{\text{P}}Z}}$，${v_{{\text{S}}ZY}}$增长速度变缓，${v_{{\text{P}}\theta }}$，${v_{{\text{S}}\theta Y}}$，${v_{{\text{P}}X}}$，${v_{{\text{S}}XY}}$开始衰减，不再满足经验式（6），此现象可能与试样微观组构发生改变有关^[15]。

图  7  P波及S波波速与应力比关系

Figure  7.  Relationship between P- and S-wave velocities and SR

下载: 全尺寸图片 幻灯片

Roesler^[27]认为波速主要受传播及振动方向上的应力影响，而${\sigma _{{\text{out}}}}$对于波速的贡献一般可忽略不计。由于等向应力固结条件下，颗粒试样不同方向的波速均与p^0.22具有正相关关系，亦可表示为(${\sigma _{{\text{pro}}}}$×${\sigma _{{\text{osc}}}}$)^0.11。因此，为区分组构及应力诱发各向异性的影响，对弹性波速利用(${\sigma _{{\text{pro}}}}$×${\sigma _{{\text{osc}}}}$)^0.11进行应力修正，修正后的波速可表示为

其中，倾斜方向的应力为

式中：${\sigma _{\text{v}}}$为竖向应力；${\sigma _{\text{h}}}$为水平向应力；${p_{\text{a}}}$为参考应力，这里取100 kPa。图 8给出了应力修正后的不同方向波速随应力比的变化规律，由图 8可知，随应力比的增加，修正后的竖向和水平向波速一开始变化较小。但当应力比超过1.6后，竖向波速的变化仍相对较小，而水平向波速则迅速衰减，且衰减幅度较大。图 9给出了应力修正后竖向波速与水平向波速的比值随应力比的变化情况。由图 9可知，修正后竖向与水平向波速的比值在应力比增加初期变化不大，表明此阶段应力对初始组构的各向异性程度影响较小。而当应力比超过1.6后，竖向与水平向波速的比值逐渐增大，并趋近于1.0，P波尤为明显，表明此阶段试样的应力诱发组构各向异性显著，具体将在后文讨论。

图  8  应力修正后P波及S波波速与应力比关系

Figure  8.  Relationship between stress normalized P- and S-wave velocities and SR

下载: 全尺寸图片 幻灯片

图  9  应力修正竖向与水平向波速比值与应力比关系

Figure  9.  Relationship between ratio of stress normalized vertical to horizontal wave velocities and SR

下载: 全尺寸图片 幻灯片

3.   微观组构试验结果分析与讨论

3.1   CT图像处理

利用法国3SR实验室开发的开源程序Spam对显微CT图像进行处理。首先，采用Otsu^[28]的方法结合视觉观察选取合适的二值化阈值将颗粒和孔隙完全分开；利用分水岭算法分割二值化后的图片，修复过分割及分割不足的颗粒，得到与原始灰度图片高度吻合的颗粒；根据分割后颗粒的体积对其进行色彩标记，具体过程见图 10。图 11显示了CT图像重构后的三维立体结构模型。计算重构试样的孔隙比为0.71，与实际制样时的孔隙比0.73相近，误差2.74%，主要源于阈值分割中孔隙像素识别的累计误差。

图  10  CT图像预处理过程

Figure  10.  Preprocessing of digital images from CT

下载: 全尺寸图片 幻灯片

图  11  PVC颗粒试样不同应力状态下三维重构模型

Figure  11.  3D reconstruction models at various stress ratios

下载: 全尺寸图片 幻灯片

3.2   单颗粒形状参数

图 12对比了实际颗粒和CT图像重构颗粒的形状及尺寸，其中等效粒径及长细比的计算见式（9）。由图可知，CT试验可实现颗粒三维形态的准确重构，并量化计算颗粒尺寸参数。表 2给出了真实颗粒与重构颗粒形状和尺寸参数的平均值。由表 2可知，重构颗粒的尺寸与实测颗粒尺寸参数相近，表明本文图像处理效果较好。

图  12  实测颗粒与图像重构颗粒的形状及尺寸

Figure  12.  Shapes and sizes of real and 3D reconstruction particles

下载: 全尺寸图片 幻灯片

表  2  实测及图像重构PVC颗粒的形状、尺寸参数

Table  2.  Shape quantities of real and 3D reconstruction particles

项目	长轴长Fmax/mm	短轴长Fmin/mm	等效粒径Deq/mm	长细比Ar
实测	4.00	2.10	2.85	1.90
重构	4.04	2.17	2.94	1.86

下载: 导出CSV 

| 显示表格

式中：${D_{{\text{eq}}}}$为颗粒等效粒径；${V_{{\text{3D}}}}$为三维颗粒体积；${A_{\text{r}}}$为颗粒长细比；F_max为长轴长；F_min为短轴长。

3.3   试样微观组构参数

为反映微观组构参数的演化规律，在不同应力状态下的试样中间位置选取相同尺寸的表征单元体（400像素×400像素×400像素）进行分析。图 13为3个应力比下试样中选取的单元体。采用Andò等^[29]及Wiebicke等^[30]提出的接触检测算法，提取单元体内的接触组构信息。图 14给出了不同应力状态下，接触配位数CN（即单个颗粒周围的接触数）的概率密度分布情况。由图 4可知，在应力比较大时，配位数出现衰减，与宏观波速的变化趋势相一致。表明宏观波速的衰减与微观接触配位数的减少有关，与Gu等^[31]的研究结论相符。

图  13  PVC颗粒试样不同应力状态下三维表征单元体

Figure  13.  3D representative element volumes at various SR

下载: 全尺寸图片 幻灯片

图  14  不同应力状态下配位数的概率密度分布

Figure  14.  PDF of CN at various SR

下载: 全尺寸图片 幻灯片

宏观波速各向异性与微观组构各向异性的演化密切相关^[15]。图 15，16分别给出了不同应力状态下，颗粒长轴方向、接触法向的3D概率密度分布图及其在竖直面（XZ面和YZ面）和水平面（XY面）的分布玫瑰图。由图可知，在等向固结（应力比SR=1.0）时，XZ面内的颗粒长轴主方向与水平面约呈30°夹角，而接触法向的主方向与竖直面约呈15°夹角，总体表现为竖向接触多于水平向。通常认为弹性波沿颗粒间的接触传播，波速会随接触数的增加而增加。然而，Otsubo等^[32]的研究发现，对于椭球形颗粒试样，尽管竖向接触数多于水平向，但沿着颗粒的长轴方向（即水平方向）传播的波速最大。这表明颗粒长轴分布对弹性波速具有重要影响。压缩波的理论计算公式为^[32]：

图  15  不同应力状态下颗粒长轴方向向量的分布

Figure  15.  Distribution of long axes of particles at various stress ratios

下载: 全尺寸图片 幻灯片

图  16  不同应力状态下试样颗粒接触法向的分布

Figure  16.  Normal distribution of contact at various stress ratios

下载: 全尺寸图片 幻灯片

式中：l为颗粒间距离；k_n为颗粒间接触法向刚度；m为颗粒质量。由式（10）可知，颗粒间距离及接触刚度均会影响弹性波速值。本研究等向固结应力下，虽竖直向接触数也多于水平向，由于受沿水平向分布的颗粒长轴影响，使得水平向波速更大，说明颗粒长轴对波速的影响大于接触数量。

由图 15，16还可知，竖向加载过程中，颗粒长轴在水平方向以及接触法向在竖直方向的概率密度均增加。然而，由于PVC颗粒表面光滑，颗粒之间的咬合作用较弱，因此颗粒容易发生偏转。在XZ面内，颗粒长轴分布的主方向逐渐从倾斜趋于水平方向，而接触法向分布的主方向逐渐从倾斜趋于竖直加载方向。当应力比增加至2.0时，在YZ面内的部分颗粒长轴趋向于竖直方向，导致水平方向颗粒长轴分布相对减少，而接触数分布相对增多。通过对XY面的组构分布观察，可以发现随着应力比的增加，X方向颗粒长轴的分布逐渐减少，接触数量增多，而Y方向颗粒长轴的分布增多，接触数量减少。由图 1可知，本研究水平向弯曲元沿X方向传播，当应力比较大时，X方向的颗粒长轴的减少加剧了水平向波速${v_{{\text{P}}X}}$的衰减，而Y方向颗粒长轴分布的相对增多使得${v_{{\text{S}}XY}}$的衰减幅度相对${v_{{\text{P}}X}}$较小。

通常利用组构张量对颗粒长轴方向、接触法向分布进行定量描述，组构张量R_ij的表达式为^[33]：

式中：$E(\mathit{\Omega } )$为组构概率密度函数；${\boldsymbol{n}_i}$，${\boldsymbol{n}_j}$为颗粒长轴方向或接触法向单位向量；N为颗粒或接触总数量。组构的空间分布概率密度函数可表示为

式中：$\gamma$为与Z轴夹角；$\beta$为与Y轴夹角；${\alpha _{20}}$，${\alpha _{22}}$，b₂₂为各向异性的系数，与材料自身性质有关。当${\alpha _{22}}$，b₂₂均为0时，该材料为横观各向同性。

利用Gu等^[15]提出的考虑不同波速方向组构概率密度的比值来表征颗粒长轴及接触法向组构各向异性程度对波速各向异性的影响。该概率密度比值为

式中：$E(\alpha {}_{{\text{1pro}}})$，$E(\alpha {}_{{\text{2pro}}})$为宏观波速1和2传播方向上的组构概率密度值，$E(\alpha {}_{{\text{1pro}}})$，$E(\alpha {}_{{\text{2pro}}})$为宏观波速1和2振动方向的组构概率密度值。图 17分别给出了不同应力状态下，${v_{{\text{P}}Z}}$，${v_{{\text{P}}X}}$，${v_{{\text{S}}ZY}}$，${v_{{\text{S}}XY}}$对应传播及振动方向上的颗粒长轴及接触法向分布的概率密度比值变化情况，由图可知，较大应力比时，竖向与水平向接触法向的概率密度比值减小，颗粒长轴的概率密度比值增大，但二者均逐渐趋近于1.0，与宏观应力修正后竖向波速与水平向波速随应力比的变化趋势相一致，具体数量关系有待后续进一步研究。

图  17  组构概率密度比值随应力比的变化

Figure  17.  Evolution of fabric probability density ratio with SR

下载: 全尺寸图片 幻灯片

4.   结论

本文基于PVC颗粒试样的多方向弯曲元波速试验，研究了不同应力状态下宏观弹性波速各向异性的变化规律。同时，采用X射线CT扫描重构试样的三维结构，利用微观组构参数随不同应力状态的变化揭示了宏观波速各向异性的演化机理。得到以下3点结论。

（1）在等向应力固结下，分层压实法生成的PVC颗粒试样中不同方向的波速均与围压呈幂函数关系，应力指数n保持在0.21~0.22。试样存在明显的初始各向异性，表现为水平向波速大于倾斜向，竖向波速最小，与颗粒长轴方向主要分布在水平方向有关。

（2）在不等向应力固结下，PVC颗粒试样中竖向及倾斜向波速随应力比的增加先增加，水平向波速保持不变，当应力比增加至1.6后，竖向波速的增速变缓，而倾斜向和水平向波速随应力比的增加开始衰减。经过应力修正后的竖向与水平向波速比值在应力比增加初期变化不大，这一阶段应力比对初始组构的各向异性程度影响较小。当应力比超过1.6后，竖向与水平向波速的比值逐渐增大，并趋近于1.0，其中P波尤为明显，试样的应力诱发组构各向异性显著增强。

（3）通过对CT重构的样本进行微观组构参数提取，发现应力比较大时，配位数有所减少，导致宏观波速衰减。同时较大应力比下，竖向与水平向分布的接触法向概率密度比值减小，而竖向与水平向分布颗粒长轴的概率密度比值增大，两者逐渐趋近于1.0。这与宏观应力修正后竖向与水平向波速比值随应力比的变化趋势相一致，具体定量关系有待后续进一步研究。