颗粒形状影响下碎石-桩界面侧摩阻力计算研究

蒋冲，丁鑫，庞里，陈陆杰

（1.中南大学资源与安全工程学院，湖南长沙 410083；2.水能资源利用关键技术湖南省重点实验室，湖南长沙 410014）

碎石强度与变形特征基本上受碎石颗粒物理属性及碎石体存在状态影响，其中前者包括矿物组成、粒径、级配及颗粒形状，后者包括碎石体结构、密实程度、有效应力等［1］.而颗粒形状作为一项重要的属性参数，在工程中越来越受到关注.

颗粒形状决定了土体材料在多尺度下的复杂力学行为［1-2］.自从Wadell［3］首先使用“球度”来描述颗粒形状以来，相关学者拓展了相关研究.随着计算机图像处理技术的应用，颗粒形状测量的精确度得以提升，大量学者使用图像处理软件获得砂土、黏土等颗粒的形态数据［4-9］，然而碎石颗粒与这些细小颗粒不同，碎石颗粒在不同角度下观察，可能表现出完全不同的颗粒形状，增加了碎石形状评价的难度［9］.在颗粒形状量化的基础上，大量学者探究了颗粒形状对材料宏观特性的影响.Yang 等［10-11］开展颗粒材料直剪试验，认为颗粒形状显著影响到临界摩擦角的大小；Han 等［12］采用不同颗粒形状的砂土与不同粗糙度的钢板进行界面剪切，结果表明界面剪切过程中，界面剪切临界摩擦角与内摩擦角的比值并不会受到颗粒形状的影响；Shinohara 等［13］通过雾化不锈钢粉末的三轴试验发现该材料的内摩擦角随棱角度及初始密实度的增大而增大；Sallam［14］通过离散元方法模拟扁平、狭长等不同形状的单元集合体的应力-应变关系，指出剪胀角受颗粒形状制约较大.直接剪切及界面剪切分别揭示了剪切过程中土体材料内部及土-结构面应力应变机制，近几年在桩基础侧摩阻力的预测中得到了广泛的应用［15-17］.

有关颗粒形状对单桩侧摩阻力影响的研究仍十分有限，Liu等［18］研究颗粒形状对桩端应力倾角的影响，认为颗粒棱角度越大，桩端土体就有越大的力链，进而能更好地支撑桩基；刘清秉等［1］借助图像处理系统将颗粒形状参数量化，并开展不同颗粒形状下的模型桩离心试验，量化分析了桩端阻力、相对密实度、应力水平和颗粒形状的相互作用关系；孙彬彬等［19］利用4 根钻孔灌注桩的现场静载试验测试火山渣地层的桩侧摩阻值，并对火山渣颗粒形状等物理力学性质进行试验，结果表明相比桩侧普通角砾介质，由于火山渣棱角度及表面粗糙度较大，桩侧土体内摩擦角及桩侧阻力都相对较大.以上研究大都通过现场或室内试验揭示了颗粒形状对桩承载性能的影响机制，但都未对此进行相应的量化，并提出相应的计算模型.对于广泛存在碎石地基，还需提出一套考虑颗粒形状影响的碎石地层中总单桩桩侧摩阻计算方法.

本文将八种碎石通过图像处理以获得形状参数，采用最小Feret直径代替传统颗粒平均粒径来反映颗粒形状的影响，并提出一种改进的整体规整度计算方法来量化颗粒形状，并通过数值模拟建立了颗粒形状与碎石-桩界面临界摩擦角及临界状态线关系计算公式.在桩侧摩阻函数中考虑颗粒形状的影响，并求解碎石地层中竖向荷载下不同深度桩基础桩侧摩阻力，与室内试验结果进行比较，验证了本文方法的适用性及合理性.

1 颗粒形状特性与数值实验结果

1.1 碎石形状量化参数OR

本文采用了八种碎石进行研究，表1 概括了这八种碎石的地质分类和简单的形状分类.试验所用的碎石尺寸均为20 mm 左右，利用圆孔分型筛将八种碎石筛分成统一级配.

表1 地质分类概述及简单的形状分类Tab.1 Summary of the geological and particle shape classification

采用三种形状参数来评价颗粒的形状：长细比AR、凸度C及球度S.AR为最小Feret直径DFmin与最大Feret 直径DFmax的比值；C为颗粒实际面积A与其凸包面积（A+B）的比值；S为颗粒实际轮廓周长Pr与等面积圆周长2πRe的比值.通过Image-Pro Plus 的形状参数测量功能，获得这八种颗粒的形状参数的频数分布，如图1 所示.图1 展示了三种形状参数的频数分布，可以看出某种碎石的形状参数集中在一个特性的区间.与Yang 等［10］测量结果不同，本文碎石凸度值范围一般都为0.9～1，长径比值及球度值范围与Yang 等［10］的结果较为相似.Yang 等［10］采用小于粒径形状50%处的形状参数值作为颗粒形状代表值，本文也采用这种方法，获得的形状数据如表2所示.

图1 三种形状参数的频数分布Fig 1 Counting number distribution for three shape parameters

表2 八种碎石的颗粒形状代表值Tab.2 Representative value of the particle shape of the eight materials

Yang 等［10］提出了整体规整度OR来评价一个颗粒在形态学上的规整程度，通过三种选取指标AR50、C50、S50的平均值来计算OR.然而，对于本文采用的碎石材料，由于试验得到的碎石凸度C与球度S分布范围相对集中（如图1），使用Yang 等的方法计算整体规整度OR，将会使凸度C与球度S占更多比重.因此本文提出了一种改进的整体规整度OR的计算方法：

式中：下标“max”“min”分别表示如图1 所示的95%置信区间的上下限；式（1）通过归一化的方法将各个形状参数的权重统一，同等地反映了三个形状参数分别对整体规整度OR的影响.

分别使用式（1）与Yang 等［10］的方法计算本文所采用的八种碎石颗粒的整体规整度OR，结果如图2所示.使用Yang 等的方法计算颗粒整体规整度OR，其值集中在0.8 左右，而使用本文的方法，OR分布在0.3～0.8，比如多孔且不规则的火山石，其整体规整度只有0.367，而较圆且光滑的珍珠石，其整体规整度为0.720.这表明对于不同种颗粒来说，本文提出的改进的方法能够明显分辨颗粒的不同形态特征.图2 也表明某种形状的碎石，其OR基本在一定范围内.

图2 Yang等［10］的方法与本文方法结果对比Fig.2 The results of the method of Yang et al［10］and the present study

工程中若采用OR，由于要测量大量的颗粒形状参数，会产生一定的不便.因此，本文提供了一种简便的方法，针对不同的碎石颗粒形状，采用相应的OR：对于颗粒接近球形的碎石可评价为“圆形的”，OR=0.8～0.9；椭球形或形状良好的可评价为“较圆的”，OR=0.6～0.8；有棱角的可评价为“较棱角的”，OR=0.4～0.6；破碎的或不规则的可评价为“棱角的”，OR=0.3～0.4.

1.2 形状参数OR对碎石-桩临界状态的影响

1.2.1 考虑颗粒形状的界面剪切数值模型

通过轮廓填充法在PFC2D软件中建立八种具有颗粒形状代表性的颗粒簇模板.颗粒簇模板建模过程：首先对碎石颗粒拍照取样，通过Photoshop 软件获得颗粒的二值图；然后在Matlab 软件中编程，筛除颗粒内部的像素点，获得颗粒的轮廓图片，读取并记录每个轮廓点的坐标；最后将剩余像素点的坐标导入PFC2D中，以像素点的坐标为圆心，一定的直径建立的圆盘单元，利用clump template 命令建立不同形状的刚性簇模板.八种碎石的实际形状和数值模型形状如图3 所示，该建模方法极大地降低了单位圆数量，且依旧能展现出颗粒的轮廓形态，既能有效降低球单元数量，又能精确模拟碎石颗粒的轮廓形态，具有较好的实用意义.

图3 八种碎石颗粒实际形状和数值模型形状Fig.3 Actual shape and numerical model shape of eight kinds of gravel particles

建立不规则碎石-桩界面剪切的离散元模型，过程如下.

1）利用PFC2D中圆盘单元建立剪切模型，将数值模拟结果与Zhang 等［20］使用颗粒-结构面循环剪切仪（TH-20t CSASSI）进行的室内剪切试验进行验证，该仪器长宽高分别为500 mm、360 mm 和300 mm.为保证模型碎石级配与Zhang 等人［20］的均质碎石级配完全相同，本文采取以下步骤设计级配：首先根据Zhang等人［20］的级配曲线，将碎石的粒径和对应累积分布分组，本模型共分10 组，对应粒径为5.4 mm、5.7 mm、6.0 mm、6.6 mm、7.0 mm、7.6 mm、8.1 mm、8.7 mm、9.4 mm、10.0 mm，粒径对应的累积分布分别为10%、20%、30%、40%、50%、60%、70%、80%、90%、100%.将此分组定义为筛分级配曲线，并将级配曲线导入PFC2D软件中；然后利用“ball distribute”命令建立碎石模型，生成八组不同粒径范围的颗粒单元.为了验证碎石模型级配的准确性，设置“测量圆”，监测数值模型的尺寸分布，并使用“measure dump”命令将颗粒模型的尺寸分布导出到表中；最后将Zhang等人［20］的粗粒砂级配曲线展示在PFC2D数值软件的表中，与建立的碎石模型的级配进行比较.最终级配设计结果如图4 所示，筛分级配、数值模型级配和Zhang 等人［20］剪切试验所用粗粒砂的级配是极其相似的，可见，颗粒流软件PFC2D可以完美建立与试验用砂级配相似的数值模型.

图4 数值模拟所用级配与试验用级配对比图Fig.4 The figure of the comparison between the simulation gradation and the test

2）在数值模拟中，为消除边界效应产生的影响，模型尺寸与试样平均颗粒半径的比值应足够大.Jamiolkowski等人［21］建议，试样尺寸与最大粒径的比值应大于5，理想比值为8.因此本研究在保证试样长高比L/H=5/3 不变的情况下，令试样长度250 mm，高度150 mm，颗粒平均粒径7 mm，孔隙率0.16，建立界面剪切数值模型如图5所示.

图5 初步建立的界面剪切数值模型Fig.5 The preliminary established numerical interface model

3）在PFC2D中利用fish 语言编写命令流，导入建立的不同形状特征的颗粒簇模板，对模型的圆盘颗粒进行替换.将八种形状的颗粒模板按照不同的百分比替换试样中的圆颗粒.替换流程如下：遍历试样中的所有圆盘颗粒，记录圆颗粒的id 号、坐标值、面积；随机替换一种形状的刚性簇模板，指定刚性簇的中心点位于替换圆颗粒的圆心处，刚性簇的面积等于圆颗粒的面积，并给刚性簇一个随机的转动角度；替换完成后，删除圆盘颗粒，对刚性簇进行id 编号，并将八种不同形状的刚性簇分成不同的组；模型建立完成，碎石试样如图6所示.

图6 碎石数值模拟模型Fig.6 The simulation model of the gravel

4）进行界面剪切和参数标定.其中剪切过程如下：首先对颗粒试样进行伺服加压，压实颗粒材料；其次固定侧边的两面墙，赋予底部混凝土结构板一个均匀速度，结构板向右移动，颗粒位移保持不变，从而实现剪切.试样模型参数展示在表3 中，图7 为剪应力-切向位移曲线碎石-混凝土界面剪切标定结果，可以看出此参数下，数值、室内试验结果拟合得很好，故可认为参数标定正确.

表3 PFC2D数值模型参数赋值Tab.3 Parameter assignment of PFC2D numerical model

图7 剪应力-切向位移曲线碎石-混凝土界面剪切标定结果Fig.7 Shear calibration results of gravel-concrete interface for shear stress-shear displacement

在PFC2D软件中，二维孔隙比与现场实际测量的三维孔隙比是不同的，且二维空间的孔隙比远小于三维空间下的孔隙比.针对孔隙比在二维和三维条件下的转化问题，周健等［22］、张刚等［23］在研究等粒径颗粒二维孔隙率和三维孔隙率转化中相差最大的三维六角晶格结构（HCP）和等粒径颗粒二维孔隙率和三维孔隙率转化中相差最小的正中心六面体结构（FCC）的两种排列结构后，经过密实度的修正，提出了等粒径颗粒二维孔隙率（e2D）向三维孔隙率（e3D）转化的公式：

式中：Dr为土体相对密实度.

1.2.2 数值实验结果

改变离散元模型的整体规整度OR为0.35、0.5、0.7、0.8 和1.0，在四个常法向应力100 kPa、200 kPa、400 kPa和700 kPa下，进行共20组碎石-桩界面剪切数值模拟，探究整体规整度OR对临界界面摩擦角的影响规律.

碎石-桩剪切模型在τ-σn平面应力变化如图8所示，这些数据点可以通过一条穿过原点的直线进行完美拟合.图8 总体表现为直线的斜率M随着整体规整度的增加而减小，即试样颗粒越不规整，其临界界面摩擦角越大.此外，进一步尝试探索临界界面摩擦角与整体规整度的关系，在图9 中，临界界面摩擦角和整体规整度OR之间出现了很强的相关性（R2=0.99），显示临界界面摩擦角随着整体规整度OR的增加而近似线性减小，可以用式（4）表示临界界面摩擦角对整体规整度OR的依赖性：

图8 不同整体规整度下τ-σn平面应力变化Fig.8 The variation of the stress in τ-σn plane under different overall regularity

图9 临界界面摩擦角与整体规整度的相关性Fig.9 Correlation between critical friction angle and overall regularity

界面剪切临界状态通常使用临界状态线表示，临界状态线为ecs-ln(σn/patm)平面上界面剪切临界孔隙比ecs与界面剪切正压力σn之间的关系，其中patm为大气压强，可由式（5）计算：

式中：e0、λ分别为临界状态线的截距及斜率.

由图10 可见，不同颗粒整体规整度OR下，临界状态线的斜率λ并未有较大变化，而截距e0变化略大.图11 提供了颗粒整体规整度OR与临界状态线截距e0之间的关系式：

图10 不同颗粒整体规整度OR下界面剪切临界状态线Fig.10 Critical state line of interfacial shear under different overall particle regularity OR

图11 本文提出OR-e0关系式Fig.11 OR-e0 relation provided by this article

2 碎石颗粒形状影响下桩侧摩阻函数

由数值试验可知，碎石颗粒形状显著影响到碎石-桩界面剪切临界状态，即改变了临界界面摩擦角大小与临界状态线.Lashkari［16］通过引入界面剪切临界状态参数，建立了双曲线形界面剪切应力应变模型：

将桩划分为等长的n个单元，则式（7）中τ为该单元桩侧剪应力；γ为该单元桩侧竖向位移；Ki为该单元碎石-桩界面初始剪切模量；τult为该单元桩侧极限剪应力，可通过式（8）计算：

式（8）中界面状态参数ψ定义为当前碎石-桩界面孔隙比e与临界孔隙比ecs之间的差值，即：

临界状态下，界面状态参数ψ=0，有：

式中：为临界状态下该单元碎石-桩界面法向应力.

本文1.2节中通过数值实验研究，得到了颗粒形状影响碎石-桩界面剪切临界状态的结论，即是对式（9）中ecs及式（10）中M分别通过式（6）及式（4）进行修正.有关式（7）中初始剪切模量Ki及式（9）中碎石-桩界面当前孔隙比e的计算方法均与Lashkari［16］的计算方法相同，详细计算公式见Lashkari［16］.

3 验证与分析

本文所提出公式适用于碎石地层中碎石-桩界面侧摩阻力的计算，然而实际工程中天然地层的情况是复杂的，大都存在多种地层.为探究单一地层中竖向荷载桩受荷响应，学者通常采用室内试验的方法，因此本文采用室内试验进行验证.

离心试验是利用离心机提供的离心力模拟重力，按相似准则研究现场工程性状的测试技术.根据相似准则，碎石-桩界面经过相似比缩小后，与砂-钢界面相似，且砂颗粒尺寸相对于钢表面粗糙轮廓尺寸与碎石颗粒尺寸相对于混凝土桩表面轮廓尺寸相似.因此本文采用Mortara［24］使用Toyoura 砂进行的砂-钢界面剪切试验及Fioravante［25］使用Toyoura 砂进行的砂-钢桩离心试验进行验证.

试验所用土体形状为“较棱角的”，根据本文所述，土体颗粒形状整体规整度OR=0.4～0.6，平均粒径D50=0.22 mm，相对密实度Dr=0.85，初始孔隙比ein=0.46，试验测定M=0.709，所需模型参数见表4.

表4 模型参数Tab.4 Model parameter

采用初始法向应力σin=200 kPa、法向刚度K=0、0.5、1、2 GPa/m 分别进行界面剪切试验.将本文计算结果与界面剪切试验结果、Lashkari［16］计算结果对比，如图12 所示.其中本文计算结果区域对应上限为OR=0.4 时计算结果、下限为OR=0.6 时计算结果及区域内部为OR=0.4～0.6时计算结果.

图12 界面剪切应力应变曲线本文计算结果对比Fig.12 Comparison of calculation results of interfacial shear stress-strain curve in this paper

可见在考虑颗粒形状影响后，在界面剪切法向刚度K=1、2 GPa/m的情况下，界面剪切计算模型能够体现“较棱角的”形状对碎石界面摩擦角的增强效应，使界面剪切有着较高的峰值剪应力，能够更好地描述界面剪切过程的应力应变软化行为.

Fioravante［25］开展Toyoura 砂中竖向荷载铝管桩离心实验.采用离心加速度Ng=30、60 进行两组试验，测定桩顶荷载-位移曲线，并测定Ng=30 情况下桩侧摩阻力-位移曲线.桩长L=24.5 cm，桩径D=1 cm.Toyoura 砂测定初始孔隙比ein=0.7，内摩擦角ϕ=33.82°，土重度γs=14 kN/m3，对于“较棱角的”碎石，eg=2.97，参数α=0.03，C=0.4，G0=125，计算所用重度γN应采用离心加速度Ng进行调整，为γN=Ng×γs.本文计算所得桩顶荷载-位移曲线与Lashkari［16］计算结果和离心实验结果如图13 所示，其中本文计算结果区域对应上限为OR=0.4 时计算结果、下限为OR=0.6 时计算结果及区域内部为OR=0.4～0.6 时计算结果.可见考虑颗粒形状后，计算模型与试验结果更为接近，能够很好地预测竖向荷载下桩受荷响应.相同竖向位移下，桩顶荷载随离心加速度增加而明显增大，说明考虑颗粒形状后，桩土体系对土重度γs的响应更为敏感.

图13 桩顶荷载-位移曲线本文计算结果对比Fig.13 Comparison of calculation results of pile top load displacement curve in this paper

β法是一种计算桩侧阻力的有效应力法，其计算公式为：

Fioravante［25］计算了离心试验中桩侧摩阻力τ对应的β值，如图14所示.其中本文计算结果区域对应上限为OR=0.4 时计算结果、下限为OR=0.6 时计算结果及区域内部为OR=0.4～0.6 时计算结果.在考虑颗粒形状后，桩侧摩阻力对应的β值与试验所测结果更为接近.在桩顶处相比于未考虑颗粒形状模型，能够较少地反应桩土体系的应力软化行为，更符合试验结果.

图14 β-γ曲线本文计算结果对比Fig.14 Comparison of calculation results of β-γ curve

4 影响参数分析

本文采用Fioravante［25］离心试验所采用的参数进行碎石颗粒整体规整度OR参数影响分析.离心加速度Ng=30，颗粒整体规整度OR选取0.35、0.55 及0.75，分别对应“棱角的”“较棱角的”及“较圆的”.由于离心试验所采用的桩为全铝管桩，相对刚度为“刚性的”，因此在竖向荷载下，桩身可认为并不发生压缩.图15 展示了桩顶竖向荷载P0=150、300 N 两种情况下，颗粒整体规整度OR对不同深度下桩侧摩阻力的影响.由图可见在桩深度较浅处，颗粒整体规整度OR对桩侧摩阻力影响较小；而随着深度的增加，颗粒形状越趋于“棱角的”，桩侧摩阻力越大.这是由于在桩埋置深度较浅的位置，桩侧水平正压力较小，土体与桩并未充分接触，颗粒形状影响不明显；而在桩埋置深度较深的位置，桩侧水平正压力较大，碎石与桩充分接触，颗粒形状越趋于“棱角的”，则土体与桩接触形式越表现为咬合与嵌固，导致该处有着较大的界面摩擦角，界面相对运动则需要更大的侧摩阻力.

图15 整体规整度OR对桩侧摩阻力Q-深度z曲线影响分析Fig.15 Influence of overall regularity OR on pile side friction Q-depth z curve

图16 展示了桩埋置深度为z=5、10、15 cm 情况下，不同颗粒整体规整度OR对桩侧摩阻力的影响.由图可知，所有曲线均是达到侧摩阻力峰值后有近似程度的软化，说明颗粒形状变化并不能显著改变桩侧摩阻力的硬化及软化行为.而随着深度的增加，桩侧摩阻力受颗粒规整度的影响变大，表现为相邻两条曲线的差值变大，这说明在桩较为刚性的情况下，桩侧摩阻力对桩侧土体颗粒形状响应更为敏感，实际工程中应更注重深部土层的情况，以充分发挥桩侧摩阻力对桩承载力的贡献作用.

图16 整体规整度OR对桩侧摩阻力Q-竖向位移γ曲线影响分析Fig.16 Effect of overall regularity OR on pile side friction Q-vertical displacement γ curve

5 结论

本文使用能够描述应力应变硬化与软化的界面剪切函数建立了考虑碎石颗粒形状影响的单桩侧摩阻力计算方法，并基于图像处理方法，提出一种考虑三种形状参数（长细比、凸度、球度）权重的改进的颗粒整体规整度计算方法，通过离散元模拟结果建立了使用颗粒整体规整度计算碎石-桩临界界面摩擦角及临界状态线截距的方法.在桩侧摩阻函数中考虑颗粒形状的影响，并求解碎石地层中竖向荷载下不同深度桩基础桩侧摩阻力.获得了如下结论：

1）土体颗粒形状主要影响界面剪切临界摩擦角及临界状态线截距.整体规整度OR与临界界面摩擦角呈较强的线性关系，表现为土体整体规整度越大，临界界面摩擦角越小；整体规整度OR与临界状态线截距为线性分段函数关系，表现为桩土界面临界状态线截距随整体规整度的增大而先增大后减小.

2）考虑颗粒形状影响后，桩侧摩阻力计算函数能够更好地反映“较为棱角的”碎石地层中，桩-土体系荷载传递过程中的应力应变硬化与软化行为.

3）对于较为刚性的桩，桩深部侧摩阻力受该处桩侧碎石体颗粒形状影响更为明显.随颗粒整体规整度的增大，深部侧摩阻力相对于浅部侧摩阻力有着明显的增长.