浙江大学学报(工学版)

浙江大学学报(工学版), 2023, 57(8): 1597-1606 doi: 10.3785/j.issn.1008-973X.2023.08.012

土木工程、交通工程

基于X-ray CT原位三轴剪切试验的砂土颗粒材料微观动力学

苗泽锴,, 张大任, 马刚,, 邹宇雄, 陈远, 周伟, 肖宇轩

1. 武汉大学 水资源工程与调度全国重点实验室，湖北 武汉 430072

2. 武汉大学 水工程科学研究院，湖北 武汉 430072

3. 长江设计集团有限公司，湖北 武汉 430010

4. 中铁第四勘察设计院集团有限公司，湖北 武汉 430063

Microscopic dynamics of sand particles based on X-ray computed tomography and in-situ triaxial compression

MIAO Ze-kai,, ZHANG Da-ren, MA Gang,, ZOU Yu-xiong, CHEN Yuan, ZHOU Wei, XIAO Yu-xuan

1. State Key Laboratory of Water Resources Engineering and Management, Wuhan University, Wuhan 430072, China

2. Institute of Water Engineering Sciences, Wuhan University, Wuhan 430072, China

3. CISPDR Corporation, Wuhan 430010, China

4. China Railway Siyuan Survey and Design Group Co. Ltd, Wuhan 430063, China

通讯作者: 马刚，男，教授. orcid.org/0000-0002-1865-5721. E-mail： magang630@whu.edu.cn

收稿日期: 2022-09-18  

基金资助: 国家重点研发计划资助项目(2022YFC3005503)；国家自然科学基金资助项目(51825905，U1865204)；云南省重大科技专项计划资助项目(202202AF080004)

Received: 2022-09-18  

Fund supported: 国家重点研发计划资助项目(2022YFC3005503)；国家自然科学基金资助项目(51825905，U1865204)；云南省重大科技专项计划资助项目(202202AF080004)

作者简介 About authors

苗泽锴（1998—），男，硕士生，从事高坝结构数值仿真研究.orcid.org/0000-0003-1806-4910.E-mail：2016301580064@whu.edu.cn , E-mail：2016301580064@whu.edu.cn

摘要

将X射线断层扫描技术(CT)与原位三轴剪切试验相结合，分析渥太华砂在剪切过程中的微观动力学演化规律. 在试验过程中共完成15次X射线扫描，使用图像分割算法进行颗粒分割并使用球谐函数重构颗粒的表面形貌，根据颗粒的多尺度形态指标序列实现整个加载过程中颗粒的准确匹配与追踪，并分析颗粒位移、转动、局部非仿射运动和局部孔隙率等微观动力学和微观结构指标的演化规律. 在剪切过程中颗粒体系的竖向位移分布呈现2个锥形区域，颗粒的转动分布出现明显的X型剪切带. 用于度量局部塑性变形程度的局部非仿射运动和局部体积分数呈现出较为明显的相关关系，表明颗粒微观动力学与其微观结构之间存在因果关系，局部自由体积较大的地方更易发生塑性变形.

关键词： 颗粒材料 ; X射线断层扫描(CT) ; 三轴试验 ; 微观动力学 ; 微观结构

Abstract

X-ray computed tomography (CT) and in-situ triaxial shear test were combined to analyze the microscopic dynamics evolution of Ottawa sand under triaxial compression. A total of 15 X-ray scans were taken during the experiment. The particles were separated by the image segmentation algorithm and reconstructed by spherical harmonic functions. The particles were matched exactly and tracked during the loading process based on the multi-scale morphological indicators of particles, and the evolution of microscopic dynamics and microscopic structural indicators, such as particle displacement, rotation, local non-affine motion, and local porosity were analyzed. During the shear process, the vertical displacement distribution of the particle system presents two conical regions, and the rotational distribution of the particles exhibits obvious X-shaped shear bands. The local non-affine motion, which was used to measure the local plastic deformation, was significantly correlated with the local volume fraction, suggesting a causal relationship between the microscopic dynamics and particle microstructure, i.e., plastic deformation was more likely to occur where the local free volume was large.

Keywords： granular material ; X-ray computed tomography (CT) ; triaxial compression ; microscopic dynamics ; microstructure

PDF (3189KB) 元数据 多维度评价 相关文章 导出 EndNote| Ris| Bibtex  收藏本文

本文引用格式

苗泽锴, 张大任, 马刚, 邹宇雄, 陈远, 周伟, 肖宇轩. 基于X-ray CT原位三轴剪切试验的砂土颗粒材料微观动力学. 浙江大学学报(工学版)[J], 2023, 57(8): 1597-1606 doi:10.3785/j.issn.1008-973X.2023.08.012

MIAO Ze-kai, ZHANG Da-ren, MA Gang, ZOU Yu-xiong, CHEN Yuan, ZHOU Wei, XIAO Yu-xuan. Microscopic dynamics of sand particles based on X-ray computed tomography and in-situ triaxial compression. Journal of Zhejiang University(Engineering Science)[J], 2023, 57(8): 1597-1606 doi:10.3785/j.issn.1008-973X.2023.08.012

颗粒材料在自然界、工程建设和工业生产中广泛存在，是一种由大量离散固体组成的多尺度耗散体系，其力学特性较复杂，根据颗粒的运动状态，可以表现出气体、流体和固体的力学行为[1]. 当颗粒材料以类似于固体的形式存在时，离散的颗粒在微观上相互接触、相互作用，以承受施加在它们身上的外部荷载. 在岩土工程领域，砂土地基和土工构筑物的变形和稳定问题都与颗粒物质力学密切相关[2]. 滑坡、泥石流、管涌等典型的地质灾害现象也与颗粒物质力学密切相关，其本质为碎散物质组成的地质体的固-液转变[3]. 因此，颗粒材料的宏细观多尺度力学特性一直是水利、岩土、交通、地质等领域的热门课题之一[4]. 从微观尺度研究砂土受剪切作用并逐渐失稳的动力学过程对理解颗粒材料具有重要意义[5].

当前用于颗粒体系微观特性测量的技术主要有数字图像法[6]、激光散斑法[7]、透视法[8]、声学测试技术[9]等. 大部分试验技术或受限于二维空间，或受限于发射源功率，暂无法应用于大规模三维颗粒的试验研究. 目前较为成熟、使用较多的试验技术为X射线断层扫描（X-ray computed tomography，X-ray CT）技术[10]，该技术属于透视法的一种.

X射线断层扫描作为无损成像技术，能够为颗粒材料性能的定量分析提供强有力的支持. 随着空间分辨率的提升，X-ray CT逐渐成为观测颗粒材料的理想工具. X射线对材料密度的敏感性意味着，如果测量的空间分辨率足够高，便可以较容易地分辨出颗粒与孔隙. 近年来，一些学者已经使用X-ray CT技术完成了颗粒形态[11-13]、孔隙结构[14]、孔隙连通性[15]等空间信息的识别，以及颗粒尺度运动学的调查研究[16-18]. 但是，这些研究使用的颗粒匹配准则较为简单，例如体积、表面积、主轴长度等[19]，对扫描精度有较高的要求.

颗粒体系内大量的离散颗粒在微小扰动下即可产生非线性相互作用，例如颗粒间的相对滑动与接触演化. 这使得颗粒材料在剪切作用下容易发生塑性变形，如发生局部化失稳[20]或分散性失稳[21]. 颗粒体系内部发生局部塑性变形意味着内部局域结构的演化，这与颗粒材料的物理力学性能密切相关. 因此，识别和量化颗粒体系的关键局部结构特征有助于预测局部塑性变形的演化，从而将其与宏观性质联系，形成对颗粒体系塑性行为跨尺度的理解. 当前，虽然诸多学者探索了颗粒材料局部塑性变形的空间和时序相关性[22]，但细观尺度塑性变形行为的产生和演化及其与宏观性质的直接联系仍有待进一步研究. 并且，目前研究偏向于单分散颗粒体系和简单球颗粒形状，对复杂颗粒形状的研究较少[23].

本研究使用X射线断层扫描技术进行渥太华砂的原位三轴剪切试验，使用图像处理算法完成渥太华砂的多尺度形态表征，并基于颗粒形貌指标序列实现大规模颗粒的准确匹配与跟踪，研究剪切颗粒体系微观动力学与宏观力学响应之间的联系，为颗粒材料的力学特性研究提供了新的思路.

1. 实验研究

1.1. X-ray CT系统及三轴试验装置

本试验所采用的X-ray CT扫描设备为天津三英精密仪器公司生产的nano Voxel-4000开管反射式高穿透CT系统. 该CT系统搭配高电压的微焦点射线源，能够获得比医学CT更高的扫描精度，适合高分辨率检测. 此外，该CT系统具备较大的内部空间及足够承载能力的样品台，满足三轴试验需求. 射线源电压为200 kV，电流为300 A，空间分辨率为16.5 µm.

本次试验使用的三轴试验装置具备体积小、重量轻的特点，满足在CT机内进行原位三轴剪切试验的需求. 如图1(a)所示为试验仪器布置图，三轴仪被固定于CT系统内部的转台上，射线源发射X射线穿透试样，探测器接收后根据射线衰减程度生成一系列二维切片图像. 如图1(b)所示为三轴试验装置的结构示意图，该装置由轴压控制系统、围压控制装置、压力室及数据采集与控制系统构成. 轴向加载装置位于机身顶部，由步进电机、力位移传感器、加载杆等部件组成. 围压由标准体积压力控制器提供，使用的液体介质为纯净水. 压力室位于机身的底部，外壳采用便于X射线穿透的聚醚酮材料，该材料具备良好的韧性和刚性，最大可承受10 MPa的围压，满足各类试验需求. 数据采集与控制系统通过数据线与三轴试验装置连接，与标准体积压力控制器一同置于CT机外.

图 1

新窗口打开| 下载原图ZIP| 生成PPT

图 1   适用于CT的三轴试验装置

Fig.1   Triaxial test device for CT


1.2. 试样制备

选用渥太华砂(见图2)作为试验材料，渥太华砂是一种由石英组成的硅质砂，形状圆润且质地坚硬. 相较于常见的福建标准砂，使用渥太华砂更有利于归纳总结颗粒材料的力学性质. 如图3所示，统计了所选渥太华砂样等效粒径的相对频率 $\varphi (d)$与累计频率 $\phi (d)$，试验砂样的等效粒径d分布为0.4~0.8 mm，累计粒度分布百分数达到50%时所对应的粒径 ${d_{50}} = 0.687$mm，可以保证高扫描精度和足够的颗粒数目. 该试验为排水三轴试验，试样无需饱和，仅须在制样前对砂粒进行干燥，减少X射线在穿透试样时的损耗，以保证颗粒轮廓清晰，方便后续图像处理.

图 2

新窗口打开| 下载原图ZIP| 生成PPT

图 2   渥太华砂及其扫描电镜图像

Fig.2   Ottawa sand and its images obtained by scanning electron microscope


图 3

新窗口打开| 下载原图ZIP| 生成PPT

图 3   渥太华砂粒度分布图

Fig.3   Ottawa sand particle size distribution chart


试样的制备流程如图4所示. 首先在三轴仪底座垫上透水石，使用橡胶膜包裹住底座和透水石(见图4(a))，随后使用制样模具包裹橡胶模. 制样模具内壁留有一孔，用于连接真空机抽真空以严格保证试样尺寸. 随后使用落雨法[24]分5层注入砂粒，每层砂粒注入完毕后轻拍制样模具以使砂粒均匀堆积. 在砂样填充完毕后，刮平试样顶部(见图4(b))，垫上透水石和试样顶帽(见图4(c)). 最后取下制样模具，在试样的顶部和底部套上密封圈(见图4(d)). 试样尺寸为12 mm×25 mm，通过CT扫描结果计算出试样初始孔隙比为0.608，属于中等密实试样.

图 4

新窗口打开| 下载原图ZIP| 生成PPT

图 4   三轴试样制备流程

Fig.4   Triaxial sample preparation process


1.3. 试验过程

在试验开始前，将三轴仪固定于转台上，将围压加载至300 kPa. 在开始加载前进行第1次扫描，之后步进电机以25 µm/min的恒定速率进行加载，每隔1%或2%轴向应变进行一次扫描，共计15次，在扫描过程中围压保持恒定. 由于扫描时间相对加载时间较短，围压对试样产生的固结影响微弱.

如图5所示为试验的宏观力学响应曲线. 图中，σ1/σ3为应力比，εV为体积应变，εa为轴向应变. 在扫描时颗粒体系会出现应力松弛现象，应力水平明显下降，这是扫描时暂停加载引起的. 在恢复加载后，应力水平会恢复至下降前的值. 本次试验的偏应力峰值出现在轴向应变6.62%下，最大应力比为4.09，且在峰值后出现了明显的应力软化现象. 在体积响应方面，试样首先体积缩小，在3%的轴向应变后试样体积开始膨胀，加载结束时体积应变为−4.83%.

图 5

新窗口打开| 下载原图ZIP| 生成PPT

图 5   渥太华砂三轴剪切试验宏观力学响应曲线

Fig.5   Macroscopic mechanical response curve of Ottawa sand triaxial shear test


2. X-ray CT图像处理

在CT扫描完成后将获得15组水平方向的试样二维切片(见图6 (a)). 初始状态包含1491张切片，在加载终止时18%的轴向应变状态下包含1241张切片.

图 6

新窗口打开| 下载原图ZIP| 生成PPT

图 6   扫描图像的分水岭算法处理过程

Fig.6   Watershed algorithm processing of scanned images


2.1. 颗粒分割

首先，15组切片须裁成相同大小，且每组切片的坐标原点须一致，以确保后续颗粒匹配的准确性. 随后，将原始的32 bit图像转化为8 bit (见图6 (b))，以缩减图片大小，提高计算效率. 之后对图像进行阈值分割，选取合适的灰度阈值对颗粒和孔隙像素进行区分(见图6 (c)). 由于CT扫描的局部体积效应，颗粒与孔隙的交界处会出现图像模糊的现象，因此在选取灰度阈值时应综合考虑多张切片图像. 如图7所示为所有切片图像的灰度值分布结果. 图中，g为灰度值，n为像素数量，在波谷范围内微调阈值，通过测算分割得到的颗粒总体积，并对比根据颗粒密度计算测得的真实颗粒体积，来确定最终的分割阈值，以保证颗粒间接触检测的准确性. 该方法是CT图像处理的常用方法，已得到多位学者的验证[19].

图 7

新窗口打开| 下载原图ZIP| 生成PPT

图 7   CT图像灰度分布图

Fig.7   Gray-scale value distribution of CT image


对阈值分割后的图片进行三维重构，提取所有被识别为颗粒的像素的三维坐标，即可得到颗粒体系的三维模型. 最后，使用分水岭算法[25]将相互连接的颗粒进行分割. 该方法计算每个像素点距其最近孔隙的距离，以局部极大值为中心向周围距离减小的方向扩散，并去除颗粒间的接触像素，由此得到相互独立的颗粒. 如图6(d)所示为渥太华砂的分割效果，不同颜色代表不同颗粒. 初始状态下共标记出10378个颗粒.

2.2. 颗粒形态特征统计

颗粒的形态学参数作为砂粒的固有特征，是颗粒匹配追踪的重要依据. 通过CT图像分割获得的颗粒表面会产生锯齿状的缺陷，不利于颗粒形态的表征，因此使用球谐函数方法恢复颗粒的光滑表面. 球谐函数重构可以精确地逼近三维颗粒的边界，准确地还原颗粒的真实形态[26-27]. 提取扫描颗粒的表面轮廓顶点，将表面顶点映射到单位球面上，由此得到一组由球坐标 $\theta $和 $\varphi $表示的曲面顶点. 通过球谐函数级数将真实砂粒表面的极半径从一个单位球面展开：

$ r\left( {\theta ,\varphi } \right) = \sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {c_n^m} } Y_n^m\left( {\theta ,\varphi } \right), $
(1)

$ Y_n^m\left( {\theta ,\varphi } \right) = \sqrt {\frac{{\left( {2n+1} \right)\left( {n - m} \right)!}}{{4{\text{π}} \left( {n+m} \right)}}} P_n^m\left( {\cos \, \theta } \right){{\rm{e}}^{{\rm{i}}m\varphi }}. $
(2)

式中：r为颗粒顶点的极半径， $r ( {\theta ,\varphi } ) = \left[{\displaystyle \sum {_{x,y,z}{{\left( {k - {k_0}} \right)}^2}} } \right]^{1/2}$， $\theta \in \left[ {0,{\text{π}} } \right]$， $\varphi \in \left[ {0,2{\text{π}} } \right]$，k为颗粒表面点坐标，k0为颗粒质心坐标； $c_n^m$为球谐系数， $n \in {\bf{N}}$为球谐级数，取值与重构精度有关； $P_n^m\left( {\cos\, \theta } \right)$为n次m阶勒让德多项式. 式(1)共有 $\left( {n+1} \right)$个未知数， $r\left( {\theta ,\varphi } \right)$作为式(1)的左侧输入，可以得到含有 ${\left( {n+1} \right)^2}$个未知数的线性方程组. 因此，可以通过最小二乘法求解这些线性方程组及其包含的所有球谐系数. $P_n^m\left( {x } \right)$可以用Rodrigues公式表示为

$ {P}_{n}^{m}\left(x\right)={\left(1-{x}^{2}\right)}^{\left|m\right|/2} \frac{{{\rm{d}}}^{\left|m\right|}}{{\rm{d}}{x}^{\left|m\right|}}\left[\frac{1}{{2}^{n}n!} \frac{{{\rm{d}}}^{n}}{{\rm{d}}{x}^{n}}{\left({x}^{2}-1\right)}^{n}\right]. $
(3)

球谐系数描述了颗粒形状在不同尺度下的特征，根据重构的精度需求确定最大球谐系数 ${n_{\max }}$，共计有 ${\left( {{n_{\max }}+1} \right)^2}$个复数描述三维颗粒的表面，即球谐级数越高，就可以获得更高的重构精度. 根据之前学者的研究[28]，当球谐级数超过15时，球谐函数足以代表砂粒的一般形态特征，故本研究选用球谐级数 ${n_{\max }} = 15$进行渥太华砂粒的球谐重构，如图8所示为某颗粒重构的结果.

图 8

新窗口打开| 下载原图ZIP| 生成PPT

图 8   颗粒球谐重构结果

Fig.8   Spherical harmonic reconstruction results of granular


在完成颗粒的表面重构后，计算每个颗粒的几何特性，如质心、体积、表面积、主轴方向等. 本研究采用主成分分析法确定颗粒的长轴、中轴和短轴方向. 由于构成每个颗粒像素单元的权重相等，颗粒绕空间中任意一点旋转的惯性可以用惯性张量进行描述，每个颗粒的主惯性轴方向可以用3个正交的单位向量表示.

由颗粒的质心坐标 $\left( {{x_0},{y_0},{z_0}} \right)$和颗粒的边界顶点坐标 $\left( {{x_1},{y_1},{z_1}} \right)$、 $\left( {{x_2},{y_2},{z_2}} \right)$、…、 $\left( {{x_N},{y_N},{z_N}} \right)$定义一个 $3 \times N$维的矩阵：

$ {\boldsymbol{X}} = \left[ {\begin{array}{*{20}{c}} {{x_1} - {x_0}}&{{x_2} - {x_0}}& \cdots &{{x_N} - {x_0}} \\ {{y_1} - {y_0}}&{{y_2} - {y_0}}& \cdots &{{y_N} - {y_0}} \\ {{z_1} - {z_0}}&{{z_2} - {z_0}}& \cdots &{{z_N} - {z_0}} \end{array}} \right]. $
(4)

式中：N为颗粒边界顶点的个数.

其协方差矩阵的特征值对应的3个特征向量 ${{\boldsymbol{c}}_1}$、 ${{\boldsymbol{c}}_2}$、 ${{\boldsymbol{c}}_3}$分别对应颗粒的长、中、短3个主轴方向. 为了获得颗粒的主轴大小，须对颗粒空间坐标进行降维处理，将颗粒旋转至惯性主轴与笛卡尔坐标轴平行，从而得到3个主轴构建的新空间. 进而由颗粒边界在新坐标系下的3个方向的最大值与最小值坐标的差值确定颗粒的长轴a、中轴b和短轴c.

为了提高颗粒匹配的精度，对颗粒的整体形态和局部圆度进行统计[29-30]. 首先根据三主轴长度定义颗粒的伸长度和扁平度：

$ E = b/a , $
(5)

$ F = c/b. $
(6)

采用球度 $S$来描述颗粒的整体形态，从而量化粒子性状与球体的符合程度. 另一个计算颗粒紧密程度的指标是凸度 ${C}$，定义为颗粒的体积与包围颗粒的凸壳体积的比值. 圆度 $R$定义为粒子轮廓角的平均曲率半径与最大内切球半径的比值. 本研究通过计算平均曲率来计算颗粒的平均圆度.

$ S = {{\sqrt[3]{{36{\text{π}} {V^2}}}} \Big/ {{S_{\rm{A}}}}}, $
(7)

$ {C} = {V /{{V_{\rm{CH}}}}},$
(8)

$ R = {{\sum {\left( {{A_n}\frac{{{k_{\rm{s}}}}}{{{k_{\rm{M}}}}}} \right)} } \Bigg/ {\sum {\left( {{A_n}} \right)} }}. $
(9)

式中： $V$为颗粒体积， ${S_{\rm{A}}}$为颗粒的表面积， ${V_{\rm{CH}}}$为包围颗粒的凸壳的体积， ${A_n}$为颗粒表面第n个三角形网格的面积， ${k_{\rm{s}}}$为颗粒最大内接球的对应曲率， ${k_{\rm{M}}}$为颗粒的平均曲率.

如图9所示，展示了初始状态（ ${\varepsilon _{\rm{a}}} = 0\text{%} $）、峰值状态（ ${\varepsilon _{\rm{a}}} = 8\text{%} $）以及临界状态（ ${\varepsilon _{\rm{a}}} = 18\text{%} $）下颗粒的等效粒径、伸长度以及球度的计算结果. 图中，φ (·)、ϕ(·) 分别表示相应变量的相对频率、累计频率. 在各个加载状态下颗粒数量和形态学参数的频率分布有轻微的波动，可能是在颗粒分割过程中接触位置的不同对颗粒形状产生了一定影响，但没有改变整体颗粒形态学参数的分布规律，验证了颗粒基本没有发生破碎的结论.

图 9

新窗口打开| 下载原图ZIP| 生成PPT

图 9   渥太华砂形状参数的频率分布演化

Fig.9   Evolution of frequency distribution of Ottawa sand grain shape parameters


2.3. 颗粒匹配与追踪

由于同一颗粒在不同应变状态下的编号独立，须对颗粒进行匹配与追踪[17]. 如图10所示为颗粒匹配的基本流程. 对于 $\varepsilon {\text{% }}$轴向应变状态下的一个待匹配颗粒，读取该颗粒的质心坐标 $\left( {{x_0},{y_0},{z_0}} \right)$以及形态参数序列. 随后，在 $(\varepsilon {\text{+}}\Delta \varepsilon ){\text{% }}$轴向应变状态下的颗粒体系内，从潜在匹配区域内（球心为 $\left( {{x_0},{y_0},{z_0}} \right)$，半径为 $2{d_{50}}$）寻找体积差值小于10%的颗粒作为候选颗粒. 计算每个候选颗粒与待匹配颗粒的形态参数之差的平方和：

图 10

新窗口打开| 下载原图ZIP| 生成PPT

图 10   颗粒匹配过程示意图

Fig.10   Schematic diagram of particle matching process


$ \left\| {{\boldsymbol{R}}_q^k - {\boldsymbol{R}}_q^{k+1}} \right\| = \left[ {\sum\limits_{q = 1}^Q {{{\left( {\left\| {{\boldsymbol{R}}_q^k} \right\| - \left\| {{\boldsymbol{R}}_q^{k+1}} \right\|} \right)}^2}} } \right]^{1/2}. $
(10)

式中：k为轴向应变状态；q为选取的形态特征指标编号；Q为形态特征指标的数量，本研究共选取了5个指标，故Q=5. 将平方和最小的颗粒作为匹配颗粒. 若该颗粒此前已与其他颗粒匹配成功，则检查2次匹配的平方和，较小值对应的配对关系不变，另一组重新匹配.

为了验证颗粒匹配结果的准确性，同时采取传统基于体积的颗粒匹配方法对相邻应变状态下的颗粒进行匹配，共得到1222个匹配结果不同的颗粒，计算这些颗粒与2种算法得到的匹配颗粒的形态学参数差值，即可对比2种方法匹配结果的准确性. 如图11所示为2种匹配方法的结果对比. 图中，EE、ES、EC分别为伸长度、球度、凸度偏差. 可以看出，基于形态学参数的颗粒匹配效果明显优于基于体积的匹配算法的，且匹配得到颗粒的各类形态参数误差较小，验证了匹配结果的准确性. 基于以上过程可以实现所有轴向应变状态下颗粒的精确匹配，并追踪每个颗粒在整个加载过程中的运动过程，继而进行颗粒动力学分析.

图 11

新窗口打开| 下载原图ZIP| 生成PPT

图 11   2种匹配算法的颗粒形态参数偏差

Fig.11   Deviation of particle shape parameters of two matching algorithms


3. 颗粒动力学分析

3.1. 颗粒动力学参数计算

颗粒在整个加载过程中的运动，包括平动与转动2类. 如图12所示为不同加载阶段颗粒的竖向位移增量Δz的计算结果. 图中，将处于中值范围内的颗粒做透明化处理. 可以看出，颗粒的竖向位移随加载的进行不断增大，且在峰值后在颗粒体系的顶部和底部发育了2个锥形区域，并且颗粒逐渐向侧向运动使得颗粒体系发生膨胀. 如图13所示为颗粒转动角度θ的计算结果，通过主成分分析获取颗粒在各个状态下的主惯性轴方向，计算相邻2个应变状态下颗粒的主惯性轴角度差值，得到颗粒在该应变区间内的转动角度. 由颗粒转动角度θ的计算结果可知，在达到应力峰值状态后，颗粒体系中部颗粒转动值开始增大，剪切带逐渐开始形成. 在应力峰值过后，试样中部开始出现明显的X型剪切带，剪切带的中部较厚，且越靠近剪切带中心的颗粒转动值越大，表明剪切带内部颗粒产生较多的塑性变形，是颗粒体系破坏的主要诱因.

图 12

新窗口打开| 下载原图ZIP| 生成PPT

图 12   颗粒竖直向位移增量

Fig.12   Increment of particle vertical displacement


图 13

新窗口打开| 下载原图ZIP| 生成PPT

图 13   颗粒转动角度

Fig.13   Particle rotation angle


为了量化颗粒材料的塑性行为，研究颗粒材料在外荷载作用下的破坏形式，以 $D_{\min }^2$作为局部塑性变形的量化指标，计算相邻2个轴向应变状态下的非仿射变形. $D_{\min }^2$在颗粒材料研究领域经常被用来表示颗粒的塑性变形情况[31]. $D_{\min }^2$统计了在轴向应变增量内，中心颗粒的实际位移与其邻域内颗粒的最佳拟合仿射变形的均方差[22]：

$ \begin{split} D_{\min }^2(\Delta \varepsilon ) =\;& \frac{1}{{{N_i}}}\sum\limits_j^{{N_i}} {\left\{ {{{\boldsymbol{r}}_j}\left( {\varepsilon +\Delta \varepsilon } \right) - {{\boldsymbol{r}}_i}\left( {\varepsilon +\Delta \varepsilon } \right)} \right.}- \\ \;&{\left. { {\boldsymbol{J}}\left[ {{{\boldsymbol{r}}_j}\left( \varepsilon \right) - {{\boldsymbol{r}}_i}\left( \varepsilon \right)} \right]} \right\}^2}. \end{split}$
(11)

式中： $i$为中心颗粒编号， $j$为其邻域颗粒编号， ${N_i}$为在中心颗粒半径r内的颗粒编号集合， ${r_i}\left( \varepsilon \right)$为 $\varepsilon $轴向应变状态下 $i$颗粒的位置， ${\boldsymbol{J}}$为使 $D_{\min }^2$最小的最佳拟合仿射变换张量. ${\boldsymbol{J}}$的表达式如下：

$ {\boldsymbol{X}} = \sum\limits_j^{{N_i}} {\left[ {{{\boldsymbol{r}}_j}\left( {\varepsilon +\Delta \varepsilon } \right) - {{\boldsymbol{r}}_i}\left( {\varepsilon +\Delta \varepsilon } \right)} \right] \otimes \left[ {{{\boldsymbol{r}}_j}\left( \varepsilon \right) - {{\boldsymbol{r}}_i}\left( \varepsilon \right)} \right]} , $
(12)

$ {\boldsymbol{Y}} = \sum\limits_j^{{N_i}} {\left[ {{r_j}\left( \varepsilon \right) - {r_i}\left( \varepsilon \right)} \right] \otimes \left[ {{r_j}\left( \varepsilon \right) - {r_i}\left( \varepsilon \right)} \right]} , $
(13)

$ {\boldsymbol{J}} = {\boldsymbol{X}} {{\boldsymbol{Y}}^{ - 1}}. $
(14)

如图14所示为 $D_{\min }^2$的计算结果. 可以看出， $D_{\min }^2$的分布规律与颗粒转动的分布规律类似，计算出的剪切带形状相似. 由于颗粒的 $D_{\min }^2$是根据颗粒的平动结果计算得到的，这也就证明了颗粒的转动与平动存在耦合关系[32]. 在初始加载阶段，颗粒体系整体发生重排列来适应外荷载的作用， $D_{\min }^2$大值的分布较离散，呈现随机排列状态. 在应力峰值状态下，颗粒体系中部出现明显的X型的共轭剪切带，且在随后的加载过程中不断发育，导致颗粒材料失稳破坏. 结果表明，颗粒体系内的剪切带一旦形成，后续颗粒的塑性变形基本都在X型剪切带内发展，最终导致应变局部化.

图 14

新窗口打开| 下载原图ZIP| 生成PPT

图 14   3个轴向应变状态下的颗粒非仿射运动空间分布

Fig.14   Spatial distribution of nonaffine motion of particles under three axial strain states


3.2. 颗粒结构学参数计算

为了研究颗粒的结构特性与其动力学响应的内在联系，对颗粒系统的部分结构特征量进行统计. 根据以往学者的研究[33]，颗粒的配位数CN与局部体积分数 ${\varphi _{\rm{V}}}$对于颗粒的塑性行为有着显著影响，因此首先对颗粒的配位数进行统计.

将未分割的二值化图像减去分水岭分割后的图像数据，可以识别出相邻2个颗粒之间的接触信息[34]. 颗粒接触对的识别须检索每个接触像素的邻接像素，像素单元之间的邻接情况有3种：以面的形式连接、以边的形式连接和以角的形式连接. 一个像素共包含6个面、12条边和8个顶点，总共包含26种连接方式，只要满足一种连接即认为2个像素单位为邻接像素.

通过计算可以得到初始状态下颗粒体系的平均配位数为8.48，18%轴向应变状态下的平均配位数为6.46，颗粒体系平均配位数不断减小. 如图15所示为颗粒在剪切过程中配位数的演化过程. 可以看出，随着剪切的进行，颗粒间配位数的频率分布向配位数小的方向移动.

图 15

新窗口打开| 下载原图ZIP| 生成PPT

图 15   剪切过程中配位数频率的分布演化

Fig.15   Evolution of frequency distribution of coordination number during shearing


局部体积分数采用Set Voronoi[35-36]剖分方法进行计算. 根据颗粒的边界来划分颗粒体系的孔隙结构，计算所得到的Voronoi元胞孔隙是最靠近对应颗粒表面的像素集合. 通过计算颗粒的体积 $V$与其Voronoi元胞体积 ${V_{\rm{V}}}$的比值，可以得到每个颗粒的局部体积分数：

$ \begin{split} & \;\\& \;\\ &\qquad {\varphi _{\rm{V}}} = {V /{{V_{\rm{V}}}}}. \end{split} $
(15)

如图16所示为不同轴向应变状态下的颗粒体系局部空隙率分布演化. 可以看出，随着加载的进行，颗粒的平均局部体积分数不断减小，也验证了颗粒体系发生剪胀的现象.

图 16

新窗口打开| 下载原图ZIP| 生成PPT

图 16   剪切过程中局部体积分数频率的分布演化

Fig.16   Evolution of frequency distribution of local volume fraction during shearing


由颗粒动力与结构分析可知，在剪切初期颗粒体系的塑性区在空间中呈现弥散状态，已有的一些数值研究表明这些看似随机的塑性活动与颗粒的微观结构密切相关[37-38]. 在剪切过程中这些弥散的塑性区会逐渐集中，在应力峰值后产生明显的剪切带，此后塑性变形主要集中在剪切带内[39].

3.3. 颗粒结构学参数与动力学参数相关性分析

本研究重点关注应力峰值前颗粒的结构特性与动力学响应的相关性. 使用 $D_{\min }^2$量化局部塑性变形， $D_{\min }^2$越大的颗粒发生的塑性活动越强烈. 首先统计了3种轴向应变下的相同配位数颗粒的 $D_{\min }^2$均值（见图17），其中配位数为当前轴向应变下的统计结果， $D_{\min }^2$为下一个相邻应变窗口的非仿射变形，通过除以当前应变窗口的平均 $D_{\min }^2$，即< $D_{\min }^2 $>，对其做归一化处理. 由图17可知，3个轴向应变状态下颗粒的 $D_{\min }^2$会随着配位数的减小而不断增大，即配位数更少的颗粒更容易发生塑性破坏. 对比不同轴向应变下的计算结果，发现轴向应变越大，曲线的斜率就越大，即颗粒的动力学响应受其配位数的影响越大. 由此可以推测颗粒间的接触会制约中心颗粒的重排，从而使颗粒处于更加稳定的状态.

图 17

新窗口打开| 下载原图ZIP| 生成PPT

图 17   配位数对于颗粒非仿射运动的影响

Fig.17   Effect of coordination number on nonaffine motion of particles


颗粒局部体积分数与动力学响应的统计规律如图18所示，由于颗粒体系在加载过程中发生剪胀，颗粒的局部体积分数也在不断减小，较高轴向应变下的关系曲线会向左偏移. 由图18可知，局部体积分数较小的颗粒的塑性行为更加显著，且随着加载的进行，这种相关关系会更加明显. 由此可知，颗粒周围的孔隙越小，其位移空间就越小，就越难发生塑性破坏. 配位数和局部体积分数与 $D_{\min }^2$的相关性分析，也验证了Zhang等[33]在数值试验中的成果.

图 18

新窗口打开| 下载原图ZIP| 生成PPT

图 18   局部体积分数对于颗粒非仿射运动的影响

Fig.18   Effect of local volume fraction on nonaffine motion of particles


4. 结　　论

（1）基于图像处理技术和分水岭分割算法，获取了各个加载阶段颗粒的多尺度形态信息. 统计了颗粒整体形态和局部圆度的信息，并提出了基于颗粒多尺度形态特征的匹配算法，实现了相邻应变间隔下颗粒的精准匹配与追踪，为颗粒动力学分析奠定了基础.

（2）通过非放射变形 $D_{\min }^2$量化了试样在剪切过程中的局部塑性变形. 在剪切初期，局部塑性区在空间中呈现弥散状态，并在加载过程中逐渐集中. 当达到应力峰值后，试样中产生明显的X型剪切带，剪切带内颗粒表现出剧烈的旋转运动和非仿射重排.

（3）在当前状态下颗粒配位数和局部体积分数与后续颗粒的局部塑性活动之间有着较强的相关性，表明颗粒微观动力学与其微观结构之间存在显著因果关系. 具体而言，有较低配位数和局部体积分数的部位更易发生局部塑性变形.

（4）本研究证明了颗粒材料的局部结构对局部塑性行为存在显著影响，此结果可为探究颗粒材料局域结构对宏观力学性能的影响提供依据，但本研究未能实现对颗粒塑性变形的预测，后续可引入机器学习的方法建立颗粒局部结构与局部塑性变形的定量联系，以进一步了解颗粒材料的变形特性.

参考文献

孙其诚, 金峰

颗粒物质的多尺度结构及其研究框架

[J]. 物理, 2009, 38 (4): 225- 232

DOI:10.3321/j.issn:0379-4148.2009.04.002      [本文引用: 1]

SUN Qi-cheng, JIN Feng

The multiscale structure of granular matter and its mechanics

[J]. Physics, 2009, 38 (4): 225- 232

DOI:10.3321/j.issn:0379-4148.2009.04.002      [本文引用: 1]

王光谦, 孙其诚

颗粒物质及其多尺度结构统计规律

[J]. 工程力学, 2009, 26 (Suppl.2): 1- 7

[本文引用: 1]

WANG Guang-qian, SUN Qi-cheng

Granular matter and the scaling laws

[J]. Engineering Mechanics, 2009, 26 (Suppl.2): 1- 7

[本文引用: 1]

郑虎, 牛文清, 毛无卫, 等

颗粒物质力学及其在工程地质领域中的应用初探

[J]. 工程地质学报, 2021, 29 (1): 12- 24

DOI:10.13544/j.cnki.jeg.2021-0017      [本文引用: 1]

ZHENG Hu, NIU Wen-qing, MAO Wu-wei, et al

Mechanics of granular material and the application in engineering geology

[J]. Journal of Engineering Geology, 2021, 29 (1): 12- 24

DOI:10.13544/j.cnki.jeg.2021-0017      [本文引用: 1]

周伟, 马刚, 刘嘉英, 等

高堆石坝筑坝材料宏细观变形分析研究进展

[J]. 中国科学: 技术科学, 2018, 48 (10): 1068- 1080

DOI:10.1360/N092018-00279      [本文引用: 1]

ZHOU Wei, MA Gang, LIU Jia-ying, et al

Review of macro-and mesoscopic analysis on rockfill materials in high dams

[J]. SCIENTIA SINICA Technologica, 2018, 48 (10): 1068- 1080

DOI:10.1360/N092018-00279      [本文引用: 1]

CHEN Y, MA G, ZHOU W, et al

An enhanced tool for probing the microscopic behavior of granular materials based on X-ray micro-CT and FDEM

[J]. Computers and Geotechnics, 2021, 132: 103974

DOI:10.1016/j.compgeo.2020.103974      [本文引用: 1]

HAGEMEIER T, BOERNER M, BUECK A, et al

A comparative study on optical techniques for the estimation of granular flow velocities

[J]. Chemical Engineering Science, 2015, 131: 63- 75

DOI:10.1016/j.ces.2015.03.045      [本文引用: 1]

BANDYOPADHYAY R, GITTINGS A S, SUH S S, et al

Speckle-visibility spectroscopy: a tool to study time-varying dynamics

[J]. Review of Scientific Instruments, 2005, 76 (9): 093110

DOI:10.1063/1.2037987      [本文引用: 1]

KOU B Q, CAO Y X, LI J D, et al

Granular materials flow like complex fluids

[J]. Nature, 2017, 551 (7680): 360- 363

DOI:10.1038/nature24062      [本文引用: 1]

张攀, 赵雪丹, 张国华, 等

垂直载荷下颗粒物质的声波探测和非线性响应

[J]. 物理学报, 2016, (2): 210- 216

[本文引用: 1]

ZHANG Pan, ZHAO Xue-dan, ZHANG Guo-hua, et al

Acoustic detection and nonlinear response of granular materials under vertical vibrations

[J]. Acta Physica Sinica, 2016, (2): 210- 216

[本文引用: 1]

程壮, 王剑锋

用于颗粒土微观力学行为试验的微型三轴试验仪

[J]. 岩土力学, 2018, 39 (3): 1123- 1129

DOI:10.16285/j.rsm.2016.0577      [本文引用: 1]

CHENG Zhuang, WANG Jian-feng

A mini-triaxial apparatus for testing of micro-scale mechanical behavior of granular soils

[J]. Rock and Soil Mechanics, 2018, 39 (3): 1123- 1129

DOI:10.16285/j.rsm.2016.0577      [本文引用: 1]

SU D, YAN W M

3D characterization of general-shape sand particles using microfocus X-ray computed tomography and spherical harmonic functions, and particle regeneration using multivariate random vector

[J]. Powder Technology, 2018, 323: 8- 23

DOI:10.1016/j.powtec.2017.09.030      [本文引用: 1]

ZHAO B, WANG J F

3D quantitative shape analysis on form, roundness, and compactness with μCT

[J]. Powder Technology, 2016, 291: 262- 275

DOI:10.1016/j.powtec.2015.12.029     

ZHANG D R, MA G, DENG Z R, et al

A self-adaptive gradient-based particle swarm optimization algorithm with dynamic population topology

[J]. Applied Soft Computing, 2022, 130: 109660

DOI:10.1016/j.asoc.2022.109660      [本文引用: 1]

KHALILI A, MATYKA M, MOHAMMADI R M, et al

Porosity variation within a porous bed composed of multisized grains

[J]. Powder Technology, 2018, 338: 830- 835

DOI:10.1016/j.powtec.2018.07.039      [本文引用: 1]

YANG B H, WU A X, MIAO X X, et al

3D characterization and analysis of pore structure of packed ore particle beds based on computed tomography images

[J]. Transactions of Nonferrous Metals Society of China, 2014, 24 (3): 833- 838

DOI:10.1016/S1003-6326(14)63131-9      [本文引用: 1]

WIEBICKE M, ANDO E, VIGGIANI G, et al

Measuring the evolution of contact fabric in shear bands with X-ray tomography

[J]. Acta Geotechnica, 2020, 15 (1): 79- 93

DOI:10.1007/s11440-019-00869-9      [本文引用: 1]

CHENG Z, WANG J F

A particle-tracking method for experimental investigation of kinematics of sand particles under triaxial compression

[J]. Powder Technology, 2018, 328 (1): 436- 451

[本文引用: 1]

LIM K W, KAWAMOTO R, ANDò E, et al

Multiscale characterization and modeling of granular materials through a computational mechanics avatar: a case study with experiment

[J]. Acta Geotechnica, 2016, 11 (2): 243- 253

DOI:10.1007/s11440-015-0405-9      [本文引用: 1]

ANDO E, HALL S A, VIGGIANI G, et al

Grain-scale experimental investigation of localised deformation in sand: a discrete particle tracking approach

[J]. Acta Geotechnica, 2012, 7 (1): 1- 13

DOI:10.1007/s11440-011-0151-6      [本文引用: 2]

杨忠平, 刘浩宇, 李进, 等

土石混合料-基岩接触面剪切力学特性及剪切带变形特征研究

[J]. 岩石力学与工程学报, 2023, 42 (2): 292- 306

[本文引用: 1]

YANG Zhong-ping, LIU Hao-yu, LI Jin, et al

Study on shear mechanical properties and deformation characteristics of shear zone of soil-rock mixture-bedrock interface

[J]. Chinese Journal of Rock Mechanics and Engineering, 2023, 42 (2): 292- 306

[本文引用: 1]

杨晓娟, 马刚, 周恒, 等

基于复杂网络的岩土颗粒材料分散性失稳先兆研究

[J]. 岩土力学, 2022, 43 (7): 1978- 1988

[本文引用: 1]

YANG Xiao-juan, MA Gang, ZHOU Heng, et al

Study on precursors of diffuse instability of granular materials based on complex network theory

[J]. Rock and Soil Mechanics, 2022, 43 (7): 1978- 1988

[本文引用: 1]

MA G, ZOU Y X, CHEN Y, et al

Spatial correlation and temporal evolution of plastic heterogeneity in sheared granular materials

[J]. Powder Technology, 2021, 378: 263- 273

DOI:10.1016/j.powtec.2020.09.053      [本文引用: 2]

MA G, ZHOU W, REGUEIRO R A, et al

Modeling the fragmentation of rock grains using computed tomography and combined FDEM

[J]. Powder Technology, 2017, 308: 388- 397

DOI:10.1016/j.powtec.2016.11.046      [本文引用: 1]

CHENG Z, WANG J F, MATTHEW R C, et al

A miniature triaxial apparatus for investigating the micromechanics of granular soils with in situ X-ray micro-tomography scanning

[J]. Frontiers of Structural and Civil Engineering, 2020, 14 (2): 357- 373

DOI:10.1007/s11709-019-0599-2      [本文引用: 1]

ZHAO B D, WANG J F, GIOACCHINO V, et al

An investigation of single sand particle fracture using X-ray micro-tomography

[J]. Géotechnique, 2015, 65 (8): 625- 641

[本文引用: 1]

SHEN L, FARID H, MCPEEK M A

Modeling three-dimensional morphological structures using spherical harmonics

[J]. Evolution: International Journal of Organic Evolution, 2009, 63 (4): 1003- 1016

DOI:10.1111/j.1558-5646.2008.00557.x      [本文引用: 1]

ZHOU B D, WANG J F, WANG H

A novel particle tracking method for granular sands based on spherical harmonic rotational invariants

[J]. Géotechnique, 2018, 68 (12): 1116- 1123

[本文引用: 1]

ZHOU B, WANG J F, ZHAO B D

Micromorphology characterization and reconstruction of sand particles using micro X-ray tomography and spherical harmonics

[J]. Engineering Geology, 2015, 184: 126- 137

DOI:10.1016/j.enggeo.2014.11.009      [本文引用: 1]

ZHOU B, WANG J F, WANG H

Three-dimensional sphericity, roundness and fractal dimension of sand particles

[J]. Géotechnique, 2018, 68 (1): 18- 30

[本文引用: 1]

付茹, 胡新丽, 周博, 等

砂土颗粒三维形态的定量表征方法

[J]. 岩土力学, 2018, 39 (2): 483- 490

DOI:10.16285/j.rsm.2017.1825      [本文引用: 1]

FU Ru, HU Xin-li, ZHOU Bo, et al

A quantitative characterization method of 3D morphology of sand particles

[J]. Rock and Soil Mechanics, 2018, 39 (2): 483- 490

DOI:10.16285/j.rsm.2017.1825      [本文引用: 1]

MEI J Z, MA G, WANG Q, et al. Micro- and macroscopic aspects of the intermittent behaviors of granular materials related by graph neural network [J]. International Journal of Solids and Structures, 2022: 111763.

[本文引用: 1]

ZOU Y X, MA G, MEI J Z, et al

Microscopic origin of shape-dependent shear strength of granular materials: a granular dynamics perspective

[J]. Acta Geotechnica, 2022, 17 (7): 2697- 2710

DOI:10.1007/s11440-021-01403-6      [本文引用: 1]

ZHANG Y B, ZHOU W, MA G, et al

The structure-property relationship of granular materials with different friction coefficients: insight from machine learning

[J]. Extreme Mechanics Letters, 2022, 54: 101759

DOI:10.1016/j.eml.2022.101759      [本文引用: 2]

CHENG Z, WANG J F

Experimental investigation of inter-particle contact evolution of sheared granular materials using X-ray micro-tomography

[J]. Soils and Foundations, 2018, 58 (6): 1492- 1510

DOI:10.1016/j.sandf.2018.08.008      [本文引用: 1]

SCHALLER F M, KAPFER S C, EVANS M E, et al

Set Voronoi diagrams of 3D assemblies of aspherical particles

[J]. Philosophical Magazine, 2013, 93 (31-33): 3993- 4017

DOI:10.1080/14786435.2013.834389      [本文引用: 1]

邹宇雄, 马刚, 李易奥, 等

椭球颗粒体系剪切过程中自由体积的分布与演化

[J]. 力学学报, 2021, 53 (9): 2374- 2383

DOI:10.6052/0459-1879-21-255      [本文引用: 1]

ZOU Yu-xiong, MA Gang, LI Yi-ao, et al

Distribution and evolution of free volume of ellipsoidal particle systems during shearing

[J]. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53 (9): 2374- 2383

DOI:10.6052/0459-1879-21-255      [本文引用: 1]

SCHALLER F M, NEUDECKER M, SAADATFAR M, et al

Local origin of global contact numbers in frictional ellipsoid packings

[J]. Physical Review Letters, 2015, 114 (15): 158001

DOI:10.1103/PhysRevLett.114.158001      [本文引用: 1]

XIAO S, LIU H, BAO E, et al

Finding defects in disorder: strain-dependent structural fingerprint of plasticity in granular materials

[J]. Applied Physics Letters, 2021, 119 (24): 241904

DOI:10.1063/5.0068508      [本文引用: 1]

MA G, REGUEIRO R A, ZHOU W, et al

Spatiotemporal analysis of strain localization in dense granular materials

[J]. Acta Geotechnica, 2019, 14 (4): 973- 990

DOI:10.1007/s11440-018-0685-y      [本文引用: 1]

/

〈 〉