1. 引言
类岩堆体一般指岩土混合体,组成物质复杂,粒径跨度大,呈碎裂、松散状,节理、裂隙极其发育,自稳能力差,且容易受到地下水的影响,在我国红土高原的云南地区分布极为广泛,属于特殊软弱围岩体。类岩堆体既不属于岩体,也不属于土体,且不同的类岩堆体围岩在机械开挖或爆破开挖扰动下的力学响应不同。基于类岩堆体的特征,通常可将其划分为类土型岩堆体、类岩型岩堆体和土石混合岩堆体。其中类土型岩堆体,尤其是砂土型类岩堆体是一种典型的非连续介质,其宏观力学性能及变形特征与颗粒形状密切相关,迫切需要进一步研究。
不论是宏观还是细观层面,颗粒形状对岩土体的力学性能有着重要影响,为此,许多学者研究了如何较为简单客观地得到真实颗粒的形状参数,在将颗粒形状参数量化之后,国内外的学者围绕此点使用现场原位试验、室内模型试验以及数值实验等进行了诸多研究 [1] - [8]。
Xiao等 [1] 通过对塔城的堆石料开展了一系列大型三轴压缩试验,测定堆石料剪切过程中颗粒破碎量计算得到相对破损指数,研究其对摩擦角、模量和变形的影响。结果显示颗粒的分形维数与孔隙率和围压的对数呈线性关系。马丽琴 [2] 研究了碎石颗粒大小及形态特征对混合料强度的影响,选用不同的圆度、扁平率等形状指标进行CBR试验,结果发现集料的CBR值随碎石的扁平率增大而减小,随碎石圆度的增大而增大。赵书辉等 [3] 扫描了三种不同形状的粗砂进行室内直剪试验,探究颗粒形状对粗砂抗剪强度的影响。结果表明四种颗粒形状参数敏感性由大到小依次为:长宽比 > 扁平度 > 磨圆度 > 球度,粗砂颗粒形状愈规则,峰值内摩擦角就越小,剪切后期颗粒间滑动摩擦就愈小。Guo等 [4] 使用三维离散单元法(DEM)对各向异性颗粒材料进行了排水/不排水剪切试验,研究其细观力学特性,结果显示试件的各向异性在整体抗剪强度中占据主导地位。Luo [5] [6] 选取福建土、玻璃珠以及碎玻璃珠三种不同形状的颗粒材料,将其按不同比例混合,微观尺度上使用激光扫描技术得到颗粒的细长比、圆度和凸度,将此三个形状参数取平均作为颗粒的新形状指数,进行了宏观和微观试验,证明了颗粒形状会影响整体力学响应、应力空间和体积应变。Yang [9] 随后做了进一步的室内试验来研究颗粒级配的影响,结果发现在给定相似条件下,临界状态摩擦角与级配无关,而在给定级配下,临界状态摩擦角与颗粒形状有很大关系,其他研究还发现,相同孔隙比和围压条件下,级配良好的土比均匀级配土更容易液化。Tsomokos [10] 研究了颗粒形状和棱角性指数对砂土不排水剪切特性的影响,使用四种不同形状相同级配的均匀砂进行不排水剪切试验,结果发现圆形颗粒在达到峰值偏应力后剪应力变小,而棱角状颗粒表现出稳定的响应,在达到瞬时峰值后强度持续增加。Keramatikerman等 [11] 使用静三轴压缩试验研究了颗粒形状的影响,用圆度(R)、球度(S)和规则度(ρ)描述三种不同形状的砂颗粒,结果表明碎砂和混合砂表现出了剪胀性,而天然砂表现出应变软化收缩行为。Igwe等 [12] 制备了三种不同级配的硅砂样品,使用环形剪切仪研究了粒径分布对硅砂力学性能的影响,结果表明相同应力条件下,级配良好的试件比其他试件具有更高的峰值和稳态强度,均匀系数越高,抗剪强度越高,级配不良的砂土不仅在破坏后容易强度降低,其稳态强度也很容易降到零。
综上,国内外研究成果虽然考虑了颗粒形状,但在改变颗粒形状的室内模型试验中只得到了定性的结论,缺乏定量结果和细观力学结构的分析;且在数值模拟中采用的形状参数单一,较难反映颗粒的复杂形状特征。本论文基于室内单元体试验,考虑细长比AR和凸度C两个形状参数对砂土强度和变形特性的影响,开展数值试验,通过建立宏观应力应变参数与微观力学指标的关系,明确颗粒形状对砂性岩堆体物理力学性质的影响机制。
2. 数值模型建立方法
2.1. 形状参数选取
土体颗粒形状定量表征方法有多种,为了综合描述颗粒的不规则形状,揭示颗粒形状对砂土宏观的力学性能以及变形性能的影响规律,本文选取两个被广泛应用的形状定量表征指标:细长比AR (Aspect Ratio)和凹凸度C (Convexity)。其中细长比的定义公式为 [10]:
(1)
式中
代表与颗粒外轮廓相切平行线之间的最短距离,
代表与颗粒外轮廓相切平行线之间的最长距离。凹凸度表示为 [11]:
(2)
式中:A表示颗粒的几何面积,B表示颗粒凹陷区域的几何面积,也可以用颗粒凸包的面积与颗粒面积做差得出。考虑真实土体颗粒形状各向异性,采用Luo等 [5] 的试验作为标定参考,使用福建砂测得的AR和C两个形状参数统计如图1和图2所示。
2.2. 颗粒团簇单元构建
考虑到计算机的算力,如果用几十个圆去近似拟合一个不规则形状的砂土颗粒,那相比于圆颗粒模型,同等大小的模型盒将使用几十倍的颗粒数量。如果一个由圆颗粒组成的模型盒有5万个颗粒,那颗粒形状精细化建模得到的模型盒将有百万计的颗粒数,以现在的计算机算力还无法处理这么多的颗粒数量,以往的学者大多使用圆颗粒进行模型简化也正是因为这个原因。
为了同时考虑颗粒细长比AR和凹凸性C的影响,也兼顾计算机的运算速度,本论文构建了四颗粒模型,即一个颗粒团簇(clump)仅由4个圆颗粒粘结构成,位置分布呈矩形结构,如图3所示。
Figure 1. Cumulative distribution of shape parameter AR for Fujian sand
图1. 福建砂AR形状参数累积分布图
Figure 2. Cumulative distribution of shape parameter C for Fujian sand
图2. 福建砂C形状参数累积分布图
Figure 3. Illustration of the four-grain model
图3. 四颗粒模型示意图
建模过程为:首先确定颗粒细长比AR,构建出一个长宽比等于AR的矩形,然后在其内部四个顶点附近做四个内切圆,这样改变内切圆半径的时候,颗粒细长比AR不会改变。然后根据颗粒凹凸性C,用下式方程求解出需要的小圆半径r [11]:
(3)
式中,
,
。
设置好某形状的颗粒簇基本单元后,在PFC软件中clump建模时使用diameter命令,并设置好对应的圆颗粒直径,软件即可自动将颗粒簇的尺寸等比例放缩,最终建出的颗粒模型将与目标直径的圆颗粒有着相同的面积。
根据Luo [5] 等室内试验用到的福建砂AR以及C的累计分布图,做出5组不同形状的颗粒簇,如图4所示,分别代表20%的颗粒总数。
(a) AR = 0.59, C = 0.92, r = 0.265 (b) AR = 0.68, C = 0.94, r = 0.27 (c) AR = 0.74, C = 0.96, r = 0.285 (d) AR = 0.80, C = 0.98, r = 0.32; (e) AR = 0.87, C = 0.994, r = 0.37
Figure 4. The five grain shape parameters used in the calibration
图4. 标定试验所用五种颗粒形状图
为了得到形状参量对力学特性的影响,除了标定的砂性地层外,还改变了AR,C参数,做出了多组不同形状的颗粒单元,如表1和图5所示。
必须指出的是,四颗粒模型也有其局限性,只有AR > 0.5以及C较接近1的时候,关于内切圆半径r的方程才有解,且AR越接近0.5,让方程有解的C参数范围也越小,不过从福建土这一天然土颗粒的AR,C累积分布图可以看到,AR参数和C参数极端的情况几乎不存在,大部分的类土岩堆体颗粒都在四颗粒模型的模拟范围之内。因此,本文的研究颗粒形状参数都在
和
范围内进行变化,这也与类土岩堆体砂性地层的形状参数分布较为吻合。
Table 1. Summary of the grain shape parameters
表1. 颗粒单元形状参数汇总
(a) AR = 0.7, C = 0.99, r = 0.35 (b) AR = 0.7, C = 0.97, r = 0.302 (c) AR = 0.7, C = 0.95, r = 0.276 (d) AR = 0.85, C = 0.99, r = 0.346 (e) AR = 0.85, C = 0.97, r = 0.301 (f) AR = 0.85, C = 0.95, r = 0.275 (g) AR = 1.0, C = 0.99, r = 0.363 (h) AR = 1.0, C = 0.97, r = 0.317 (i) AR = 1.0, C = 0.95, r = 0.29
Figure 5. Grain shape adopted in the parametric sensitivity study
图5. 敏感性分析中所用不同颗粒形状图
2.3. 接触本构模型选取
本文建模使用的接触本构模型是PFC软件内置的本构模型之一:滚动阻力线性模型,这是一个以线性模型为基础拓展而来的本构模型(图6)。
Figure 6. Illustration of the linear contact model
图6. 线性接触模型示意图
线性模型假设了模型颗粒之间的接触力与相对位移之间的线性关系,相当于在接触位置放置了一个切向和一个法向两个弹簧,力–位移满足胡克定律。只要设置好颗粒的刚度属性,两颗粒接触就像弹簧串联,其接触的切向刚度和法向刚度就能通过如下公式计算获得 [13] [14] [15]:
(4)
(5)
上式中,
、
分别是颗粒1和2的法向刚度属性,
为两颗粒接触的法向刚度,
、
分别是颗粒1和2两颗粒的切向刚度属性,
为两颗粒接触的切向刚度。
滚动抵抗线性模型则是在线性模型的基础上,新设了一个滚动摩擦系数,当滚动抵抗力矩超过一定限值时,滚动抵抗力矩不再继续提升,具体公式如下 [8]:
(6)
(7)
式中
是最大滚动抵抗力矩,
是滚动摩擦系数,
是接触有效半径,通过接触两端颗粒的半径倒数和再取倒数求得,
是法向接触力,
为最终的滚动抵抗力矩。
2.4. 试样生成及加载
由于是使用二维模型模拟三轴试验,模型边界选择可平移的四个wall单元构成一个矩形区域,横向宽度71.1 mm,竖向高度142.2 mm。
试样生成方法选用重力沉积法,即让颗粒簇在整个矩形区域内随机生成,赋予颗粒刚度和重力,让颗粒簇在重力场作用下作自由落体运动,最终沉积在试样盒的底部,这样的生成方法可以让颗粒随机分布在整个模型空间内,且与实验室制样过程十分近似。
颗粒的半径按Luo等 [5] 的室内试验测得的级配曲线设置,如图7所示。
具体建模步骤如下:
Figure 7. Cumulative distribution curve of the AR parameter
图7. AR参数累积分布图
(1) 分五次填充不同形状(表1)的颗粒,每次填充(1 − 孔隙比) × 20%试样盒体积的颗粒,粒径按照级配曲线分布;
(2) 设置墙体刚度和颗粒刚度,设置接触本构模型,墙体刚度远大于颗粒刚度,无摩擦;
(3) 设置重力,添加重力场;
(4) 循环,让试样颗粒自由扩散、堆积,有些clump重复面积过大,互相嵌套,受力保持一个平衡,无法通过扩散堆积分离,需删除堆叠重复的颗粒;
(5) 设置目标围压,每个时步计算颗粒作用在墙体上的力,计算该力和目标围压之间的差,以此定义墙体运动速度,循环直至围压稳定;
(6) 计算围压稳定后的孔隙比,若当前孔隙比大于目标给定的孔隙比,下调扩散、堆积时的摩擦系数,反之则下调摩擦系数;
(7) 保持水平围压不变,撤去竖直围压,设置上下两个墙体一个恒定向中心的加载速度,加载过程中每1%遍历输出所有颗粒、接触的数据,直至整体试样达到30%轴应变。
最终建立模型如图8所示。
模型中的颗粒簇均为四颗粒模型,共36928个颗粒簇(clump),147712个圆颗粒(pebble)。
2.5. 参数标定流程
参数标定分为可调参数和不可调参数,其中颗粒和模型盒的几何尺寸、颗粒密度、固结围压、初始孔隙度、压缩速率等,主要参考Luo等 [5] 的室内试验及现有文献常识取了一个合理值,在之后标定流程中不再更改,具体见表2取值。
可调参数包括颗粒刚度、摩擦系数、转动摩擦系数、阻尼等,通过与室内试验应力应变曲线和体积应变曲线匹配。参数标定时,不可调参数不做更改,当进行一轮试样生成及加载流程后,若偏应力–轴应变曲线和体积应变–轴应变曲线无法和室内试验的结果对应时,则微调可调参数,直至两条曲线匹配上室内试验的结果。最终标定完成的参数如表3。
3. 结果分析及讨论
3.1. 应力–应变曲线与试验应力–应变曲线对比
数值模拟应力–应变曲线与试验应力–应变曲线对比如图9和图10所示。在刚进入剪切阶段时,数
Figure 8. The generated model in PFC2D
图8. PFC2D生成模型示意图
Table 2. Non-calibratable parameters
表2. 不可调标定参数表
Figure 9. The calibrated deviatoric stress-axial strain curve
图9. 标定偏应力–轴应变曲线图
Figure 10. The calibrated volumetric strain-axial strain curve
图10. 标定体积应变–轴应变曲线图
值模拟模型的弹性模量比室内试验模型的弹性模量小一些,但是当轴应变达到约8%后,数值模拟模型的偏应力应变曲线很快趋于平缓,进入了屈服阶段,数值模型的峰值强度与室内试验的峰值强度相差无几,但是曲线有着一定程度的上下波动,原因是数值模型内部在剪切过程中会不断出现局部结构的调整,每次调整都会导致偏应力–应变曲线波动。
体积应变–轴应变曲线中的体积应变是压缩应变,由于围压500 kPa较大,数值试验和室内试验都表现出了一定的剪缩性,在剪切阶段初期,数值试验的剪缩性比室内试验更显著,而当到达轴应变约8%后体积应变曲线很快趋于平缓,此状态对应了偏应力应变曲线中的屈服点。之后室内试验和数值试验的体积应变都趋近于2%,相差不大。
综上所述,标定模拟结果与室内试验结果二者十分接近,因此标定参数可作为后续细观力学性能数值试验的参数。
进一步在200 kPa至1500 kPa范围内改变围压时,通过数值试验做出类土岩堆体砂性地层的破坏包络线如图11所示。
Figure 11. The calibrated test failure envelope
图11. 标定试验的破坏包络线图
砂土颗粒的粘聚力很小,所以主要关心其内摩擦角,在破坏线包络图中体现这一性质的是拟合直线的斜率。在室内试验中,Luo [5] 给出的福建土破坏包络线斜率为1.21,数值模型的破坏包络线斜率为1.239,误差在2%以内,已十分接近。
根据破坏线的斜率可以算出福建砂的内摩擦角,用下式进行计算:
(8)
式中,M为破坏包络线的斜率,
为砂土的内摩擦角。
根据室内试验的结果算出内摩擦角为30.23˚,根据数值试验的结果算出内摩擦角为30.90˚,误差为2.2%,可见数值模拟计算出的内摩擦角和室内试验得到的内摩擦角十分接近。
经过软件导出的颗粒数据并计算后,该模型在施加上围压稳定后,剪切开始前模型盒的孔隙率为0.228。为了保证之后形状变异的砂土与数值试验中的砂土在压缩前有着除了形状因素外所有因素都相同的初始条件,可调整形状变异的砂土模型在固结阶段的摩擦系数,使所有不同形状的砂土模型在进入压缩阶段前,有着相同的孔隙比0.228,孔隙比误差控制在0.4%内。数值模型的摩擦系数将在进入剪切阶段前调回标定值。
3.2. 颗粒形状参数敏感性分析
3.2.1. 颗粒细长比影响
为了研究颗粒细长比AR对力学特性的影响,采用控制变量法,固定颗粒凹凸性C为定值,然后更改AR参数,生成试样进行剪切模拟,获得的偏应力–轴应变曲线图和体积应变–轴应变曲线如图12~图15所示。
在偏应力轴应变曲线中可以看出,改变AR后的曲线初始阶段几乎是重合的,说明AR参数不影响模型在剪切开始时的弹性模量。不同AR下曲线进入屈服阶段的时刻几乎没有改变,都在5%~10%之间,说明AR参数不影响模型的屈服点。综合图12~15,可以从中看出:C不变时,AR参数越小,试样的强度越高和压缩变形越大,且AR对体积变形的影响较对强度的影响更为明显;C在0.95~0.996范围内变化时都呈现相同的趋势。改变AR参数对体积应变曲线的影响在AR参数较接近1时,不太明显;当AR参数远离1时,影响程度越来越大。
由此可见,当颗粒凹凸性不变时,细长比越小,整体模型的变形能力越强,强度略有提高,弹性模量几乎没有影响。随着颗粒凹凸程度的增加,即C值减小时,细长比的影响也更为明显。
3.2.2. 颗粒凹凸性影响
固定颗粒细长比AR保持不变,可以获得如下偏应力–轴应变曲线图和体积应变曲线图。此外,根据试验的偏应力应变曲线和式(6)可算出不同形状颗粒组成的模型的内摩擦角,如表4所示。
(a) (b)
Figure 12. At C = 0.996 (a) deviatoric stress-axial strain curve; (b) volumetric strain-axial strain curve
图12. C = 0.996时(a) 偏应力–轴向应变曲线图;(b) 体积应变–轴向应变曲线图
(a) (b)
Figure 13. At C = 0.99 (a) deviatoric stress-axial strain curve; (b) volumetric strain-axial strain curve
图13. C = 0.99时(a) 偏应力–轴向应变曲线图;(b) 体积应变–轴向应变曲线图
(a) (b)
Figure 14. At C = 0.97 (a) deviatoric stress-axial strain curve; (b) volumetric strain-axial strain curve
图14. C = 0.97时(a) 偏应力–轴向应变曲线图;(b) 体积应变–轴向应变曲线图
(a) (b)
Figure 15. At C = 0.95 (a) deviatoric stress-axial strain curve; (b) volumetric strain-axial strain curve
图15. C = 0.95时(a) 偏应力–轴向应变曲线图;(b) 体积应变–轴向应变曲线图
在偏应力轴应变曲线图中可以看出,改变C后曲线的初始阶段不再重合,C值越小,其斜率越大,说明凹凸程度会影响模型的整体弹性模量,凹凸程度越大,C值越小,弹性模量越大。不同C下曲线进入屈服阶段的时刻几乎没有改变,都在5%~10%之间,说明C参数不影响模型的屈服点。由图12~15可以看出:AR不变时,C参数越小,抵抗变形能力越强,而强度几乎没有改变,且当C越接近1时,其抵抗变形能力的变化也越小,这一规律在AR = 0.7~1.0范围内都成立。
因此,当颗粒细长比不变时,凹凸性指数C越大,整体模型的变形能力越强,凹凸性不影响模型的整体强度(图16)。
(A-1)(A-2)(A) (B-1)(B-2)(B) (C-1) (C-2)(C)
Figure 16. (A) At AR = 0.7: (A-1) deviatoric stress-axial strain curve; (A-2) volumetric strain-axial strain curve; (B) At AR = 0.85: (B-1) deviatoric stress-axial strain curve; (B-2) volumetric strain-axial strain curve; (C) At AR = 1.0: (C-1) deviatoric stress-axial strain curve; (C-2) volumetric strain-axial strain curve
图16. (A) AR = 0.7时(A-1) 偏应力–轴向应变曲线图;(A-2) 体积应变–轴向应变曲线图;(B) AR = 0.85时(B-1) 偏应力–轴向应变曲线图;(B-2) 体积应变–轴向应变曲线图;(C) AR = 1.0时(C-1) 偏应力–轴向应变曲线图;(C-2) 体积应变–轴向应变曲线图
Table 4. Internal friction of samples with different grain shape
表4. 不同形状颗粒试样的内摩擦角
4. 结论
本文采用颗粒流软件针对类土岩堆体砂性地层颗粒形状变异对宏观力学性能及变形特征的影响进行了数值仿真试验,提出选取细长比AR和凹凸性C作为量化的颗粒形状参数,并建立了四颗粒模型,并与已有的试验结果进行了对比分析。然后通过改变颗粒的形状,对不同形状参数的颗粒在其他参数完全相同的情况下进行压缩试验,使用控制变量法,研究了两个形状参数对试验宏观力学性能的影响,获得结论如下:
(1) 在控制颗粒细长比AR介于0.6~1.0范围以及凹凸性C介于0.95~1.0范围内时,两者对压缩试验的体积应变–轴应变曲线都有较为显著的影响,其中AR还会对偏应力–轴应变曲线有不同程度的影响。
(2) 在不同的颗粒凹凸性C保持不变的情况下,试样的压缩变形随着AR的减小而增大,且C越小,此影响越显著,即颗粒的形状越扁,试样整体抵抗压缩能力越弱,在压力的作用下能够被压得更密,且整体强度也有一定程度的提高。
(3) 在颗粒细长比AR保持不变的情况下,减小凹凸性C会降低体积应变曲线,且C越小,此影响越显著,即颗粒的凹陷程度越大,试样整体抵抗压缩能力越强,在压力的作用下很难被压密,但凹凸性对整体强度几乎没有影响。
致谢
本项目得到了中电建路桥集团有限公司(编号HHZ-JGY-FW-03)的资助,在此一并致谢。