颗粒摩擦因数对胶结砂岩力学特性的影响

许明^{1, 2}，刘先珊^{1, 2}

(1. 重庆大学 土木工程学院，重庆，400045；

2. 重庆大学 山地城镇建设与新技术教育部重点实验室，重庆，400045)

摘要：为确定胶结砂岩的力学参数，探讨岩体破坏机理，采用细观力学分析方法，建立三维颗粒流数值模型，对胶结砂岩的剪切特性进行模拟。根据应力应变曲线以及平行黏结分布和接触网络变化说明胶结砂岩不同于无黏结砂岩的破坏机理，验证数值模型的可行性，分析不同胶结半径比的砂岩模型的应力比、体应变、配位数和平行黏结破坏数的变化规律。计算结果表明：胶结半径一定时，不同摩擦因数对岩样宏观力学特性的影响不同，特别是随着胶结半径比的减小，颗粒间摩擦力在胶结砂岩受力过程中占主导作用，对砂岩抗剪能力起着决定作用，不同摩擦因数对胶结砂岩力学特性的影响较大。颗粒摩擦因数对胶结砂岩数值试验中的力学行为有着重要的影响，是准确选择三维颗粒流模型细观参数的关键。

关键词：颗粒流数值模拟；胶结砂岩；剪切试验；颗粒摩擦因数

中图分类号：TU4            文献标志码：A         文章编号：1672-7207(2014)01-0299-07

Influence of particle friction coefficient on mechanical characteristics of cemented sandstone

XU Ming^{1, 2}, LIU Xianshan^{1, 2}

(1. School of Civil Engineering, Chongqing University, Chongqing 400045, China;

2. Key Laboratory of New Technology for Construction of Cities in Mountain Area, Ministry of Education,

Chongqing University, Chongqing 400045, China)

Abstract: In order to determine the mechanical parameters and failure mechanism of the cemented sandstone, micro mechanical method based on 3-D particle flow code (PFC3D) was put forward to simulate its shearing characteristics. The stress-strain curve, variation of the parallel-bond and contact network explained the broken mechanism of the cemented sandstone which is different from that of the uncemented sandstone. And the numerical model was calibrated to analyze the stress ratio, volume strain, coordination number and number of the broken parallel bonds changing with the axial strain. The simulation results show that different friction coefficients have different influences on the micro mechanical characteristics of the cemented sandstone, especially with the decrease of the parallel bonded radius ratio, friction coefficient plays a dominated role on resistance to shear when loading on the cemented sandstone. Different friction coefficients have different influences on the mechanical characteristics of the sandstone. Therefore, friction coefficient plays a significant role on the mechanical behavior of the cemented sandstone for the shearing test and correct selection of the micro parameters of the 3-D simulation model.

Key words: numerical simulation based on particle flow code; cemented sandstone; shear testing; particle friction coefficient

颗粒流数值分析方法^[1]由于能从细观角度模拟岩土介质的宏细观力学特性，在岩土材料破坏机理研究中得到了重点应用。而在数值计算中，细观参数的选择是准确反映岩土介质力学特性的关键因素。目前，运用颗粒流程序的常规接触模型对松散砂土介质及相关课题的研究有了较多的成果，如曾远等^[2]基于PFC2D模拟不同颗粒形状、颗粒摩擦因数、孔隙率等细观参数对砂土介质双轴试验结果的影响；Cui等^[3]基于PFC2D研究了理想颗粒材料剪切试验中的颗粒细观参数对砂土宏观力学特性的影响；尹成薇等^[4-5]基于PFC2D的常规接触模型进行数值双轴试验条件下黏性土宏观强度特性与细观参数的定量分析，得出宏观强度参数与对应的细观参数之间的相关关系。上述成果针对不同细观参数对砂土介质力学特性的影响研究为松散介质本构关系提供了技术手段，然而，对于由砂粒胶结而成的砂岩体，由于胶结物限制了砂岩颗粒的转动，在砂岩胶结物质破损前后，其颗粒的受力性状将发生改变。如Potyondy等^[6-7]提出了能传递颗粒间的力和力矩的平行黏结模型来模拟颗粒间胶结力，不同于松散砂岩介质中的摩擦作用；徐金明等^[8]通过PFC2D平行黏结模型细观参数的调整来获得石灰岩的宏观力学参数；Hsieh等^[9-11]基于二维平行黏结模型初步研究了砂岩宏观力学响应，验证了该模型相较于常规的颗粒接触模型能更好地反映岩体的宏观力学特性；Wang等^[12-13]基于PFC2D模拟了不同荷载作用下具有一定胶结程度的砂岩剪切力学响应，但这两个模型胶结百分比含量设置较低，采用的黏结强度均较大，与实际胶结性状有一定的出入。为了更好地运用颗粒流数值方法描述胶结砂岩的力学特性，胶结物质采用具有一定接触面积的平行黏结接触^[9-11]来模拟，在受力过程中，必须先克服平行黏结接触的破坏。对于无黏结砂土介质，颗粒摩擦因数在岩样宏观力学特性研究中具有举足轻重的位置，而对于外力作用下砂岩胶结性状发生变化时，颗粒摩擦因数对胶结砂岩的影响如何及影响胶结砂岩力学特性的主导作用的变化，人们还不了解。为此，本文作者建立以平行黏结模拟砂岩胶结性状的三维颗粒流数值模型，考虑胶结半径比变化时不同的颗粒摩擦因数对胶结砂岩力学特性的影响，研究加载过程中黏结力与摩擦力的主导作用，其成果为选用准确的细观参数进行数值模拟来描述实际工程中胶结砂岩的力学特性(如油藏储层出砂和射孔试验等)提供了重要的理论依据和研究手段。

1  数值模型建立及模型验证

建立如图1所示的颗粒流数值模型进行三轴试验，模型高度为0.08 m，直径为0.04 m，相关参数如表1所示。数值分析过程如下：(1) 对初始的砂岩颗粒流模型逐渐加载，一旦荷载达到10 MPa，导入平行黏结模拟砂岩的胶结性状。(2) 数值模型中设置砂岩颗粒胶结为均匀分布，通过PFC3D的二次开发将平行黏结的力学参数导入以模拟具有一定胶结半径、刚度和胶结量的砂岩体。(3) 开始模型剪切试验，剪切开始前的模型孔隙率为0.6。剪切试验中，底板静止，类似于实验室的等应变加载进行模型加载；实际加载值与数值计算荷载之差控制于±2.5%范围，以保证整个数值模拟过程的准静止状态。

图1  三维颗粒流数值模型

Fig. 1  Numerical model based on 3Dimensional PFC

基于建立的数值模型模拟剪切过程中的胶结砂岩力学特性，得到对应不同胶结半径比(采用胶结半径比^[14]表示，其中，为胶结平均半径，为接触处颗粒中的较小半径)的试样模型在颗粒摩擦因数=0.25时的应力比随轴向应变的变化曲线，如图2所示。从图2可见：胶结半径较大得到的峰值应力强度和初始弹性模量均较大，与文献[12-13]所研究的不同程度的胶结砂岩在外荷载作用下的变化趋势一致，说明本文采用的数值计算模型能很好的模拟胶结砂岩的宏观力学特性。同时，还得到胶结性较好(=1.0)的岩样在不同的摩擦因数(=0.5~0.125)的平行黏结分布和接触网络图，如图3~5所示，显示的是移除破坏平行黏结之后的接触网络，即平行黏结破坏后形成的剪切带。其结果表明：当胶结半径一定时，不同对应的平行黏结破坏形成的剪切带变化不大。主要原因在于外力作用下，由于试样的胶结半径较大，颗粒间的黏结力也较大，那么剪切过程中主要是克服了颗粒之间的黏结力，最终是由于颗粒间的黏结破坏而形成剪切带^[15-16]，此时颗粒摩擦因数对试样的剪切破坏的敏感性不强，颗粒间的黏结力起主导作用，这不同于一般松散砂岩的破坏机理^[12-13]。但是，在剪切过程中，不同的颗粒摩擦因数形成的接触力有所不同，摩擦因数较大的接触力逐渐也较大，符合胶结性状岩土介质力学响应的一般变化规律^[12-13]，进一步验证了本文数值模型的可行性。

表1  模型计算参数

Table 1  Parameters of PFC3D model

图2  μ=0.250时考虑不同α的应力比-应变曲线

Fig. 2  Curve of stress ratio varying with axial strain considering μ=0.250 and different α

图3  α=1.0，μ=0.500时Y-Z平面接触网络图

Fig. 3  Contact network in plane Y-Z considering α=1.0 and μ=0.500

图4  α=1.0，μ=0.250时Y-Z平面接触网络图

Fig. 4  Contact network in plane Y-Z considering α=1.0 and μ=0.250

图5  α=1.0，μ=0.125时Y-Z平面接触网络图

Fig. 5  Contact network in plane Y-Z considering α=1.0 and μ=0.125

2  颗粒流模拟结果分析

基于上述PFC3D数值模型，考虑不同的颗粒摩擦因数在不同胶结程度的砂岩(=0.1~1.0)进行剪切试验，得到了剪切过程中胶结砂岩的宏观力学特性(应力比、体应变)和颗粒的细观力学响应规律(配位数、黏结破坏数)，不同参数组合下的岩样力学特性变化趋势如图6~17所示。

(1) 图6和图7，图10和图11以及图14和图15分别描述了在胶结半径比一定时考虑不同摩擦因数对岩样应力比和体应变随轴向应变()的影响。当=1.0时，就某一给定的摩擦因数而言，总的变化趋势为：应力比随着的增大而增大，峰值之后的应变软化很明显，其后软化率随着应变的增大而减小，符合岩体的应力-应变变化规律，随着的减小，峰值后的软化不是很明显；体应变的变化趋势说明峰值应力到达之前出现了剪缩(正的体应变)和在应变软化开始后的剪胀效应。结合图9、图13和图17的平行黏结破坏发展趋势可以看出：剪胀效应在平行黏结破坏开始之后出现，黏结破坏发生后，一部分颗粒就会发生转动，颗粒的约束力减小，模型整体的变形相应增大，那么总应变也随之增大，这一现象与实际胶结砂岩的力学行为相似，同时也说明了该数值模型的可行性。另外，由图6~7、图10~11和图14~15可以看出：对任意胶结半径比的砂岩，随着颗粒摩擦因数的增大，模型都会出现很高的初始刚度和峰值应力比，峰值应力比随着轴向应变的增大而增大，摩擦因数越大，最大峰值应力比对应的应变值越小，摩擦因数越小，峰值之后的应变软化越明显，其软化率随着应变的增大而减小；体应变的变化特性也说明摩擦因数越大，剪缩越明显。由此说明胶结砂岩胶结形状发生变化时，平行黏结和摩擦力在砂岩受力过程中所起的控制作用不同。

(2) 配位数是衡量颗粒接触的一个重要指标，从细观角度反映岩样颗粒的运动状况。图8、图12和图16描述了一定时考虑不同摩擦因数对配位数的影响。对某一特定的摩擦因数，初始配位数变化缓慢，随后配位数出现峰值，之后配位数减小的速率增加。配位数总的变化并不大，说明剪切过程中大多数颗粒能充分接触，大多数颗粒都对模型的应力状态起作用。配位数的变化趋势说明摩擦因数对初始配位数的影响不大，随着摩擦因数的增大，配位数也随着增加，由图8、图12和图16的比较可知：对于胶结半径比较小的岩样，配位数的变化趋势要小一些，主要在于胶结半径比较时时，岩样类似于松散介质，其颗粒运用和承力模式主要由常规接触分析来确定。

图6 =1.0时应力比随轴向应变变化曲线

Fig. 6 Curve of stress ratio varying with axial strain considering =1.0

图7 =1.0时体应变随轴向应变变化曲线

Fig.7  Curve of volume strain varying with axial strain considering =1.0

图8 =1.0时配位数随轴向应变变化曲线

Fig. 8  Curve of coordination number varying with axial strain considering =1.0

图9 =1.0时破坏的平行黏结数随轴向应变变化曲线

Fig. 9  Curve of broken-bond number varying with axial strain considering =1.0

图10 =0.5时应力比随轴向应变变化曲线

Fig. 10  Curve of stress ratio varying with axial strain considering =0.5

图11 =0.5时体应变随轴向应变变化曲线

Fig. 11  Curve of volume strain varying with axial strain considering =0.5

图12 =0.5时配位数随轴向应变变化曲线

Fig. 12  Curve of coordination number varying with axial strain considering =0.5

图13 =0.5时破坏的平行黏结数随轴向应变变化曲线

Fig. 13  Curve of broken-bond number varying with axial strain considering =0.5

图14 =0.1时应力比随轴向应变变化曲线

Fig. 14  Curve of stress ratio varying with axial strain considering =0.1

图15 =0.1时体应变随轴向应变变化曲线

Fig. 15  Curve of volume strain varying with axial strain considering =0.1

图16 =0.1时配位数随轴向应变变化曲线

Fig. 16  Curve of coordination number varying with axial strain considering =0.1

图17 =0.1时破坏的平行黏结数随轴向应变变化曲线

Fig. 17  Curve of broken-bond number varying with axial strain considering =0.1

(3) 图9、图13和图17显示胶结半径比一定时考虑不同摩擦因数对平行黏结破坏数的影响规律。破坏的黏结数在最大峰值应力后迅速增加，说明颗粒之间的黏结破坏越多，黏结破坏的颗粒就会发生转动，运动时的约束就相应减少，与图6~8得出的结论一致。当一定时，摩擦因数减小，颗粒之间的摩擦力就越小，说明沿剪切方向的抗力减小，那么在外力作用下，剪切方向的平行黏结破坏越多。

(4) 图10和图11、图14和图15描述了胶结半径比=0.5和=0.1时的宏观力学特性变化趋势。当胶结砂岩黏结半径比较小时=0.1，颗粒摩擦因数对峰值强度和宏观弹性模量的影响很大：当颗粒内摩擦因数从0.25减小到0.125，弹性模量和峰值强度会降低约为14%和25%；随着胶结半径比的减小，配位数随着摩擦因数的减小而变大；当=0.5时，初始状态下依然是摩擦因数较大对应的配位数较大，随着剪切过程的进行，应变的增加，其配位数随着摩擦因数的减小而变大，说明在黏结较弱的情况下颗粒之间的摩擦因数越小，颗粒之间的接触越容易。由图9、图13和图17比较可知：当平行黏结破坏时，随着的减小，摩擦因数对平行黏结破坏的影响不可忽视：当黏结半径较小时，平行黏结破坏出现在最大峰值应力之前，砂岩力学特性主要由常规接触控制，摩擦因数决定了模型的抗剪能力，这也解释了该状态下摩擦因数对岩样力学响应影响较大的原因：当黏结半径较小时，砂岩胶结性状较弱，颗粒的胶结性状对岩样力学行为不再起主导作用，那么，颗粒之间的运动主要由颗粒自身特性决定，其摩擦因数的敏感性就会增强。

3  结论

(1) 以平行黏结模拟胶结砂岩的黏结性状建立PFC3D数值模型，通过考虑不同颗粒摩擦因数的平行黏结分布和接触网络图阐述了胶结砂岩的破坏机理，也验证了本文三维数值模型的可行性。同时，研究了考虑不同摩擦因数对不同程度胶结砂岩力学特性的影响规律。

(2) 随着砂岩胶结半径比的减小，颗粒摩擦因数对砂岩介质力学特性的影响起着主导作用，因而对峰值强度和宏观弹性模量的影响增加。

(3) 随着胶结半径比的减小，配位数随着摩擦因数的减小而变大，说明在弱胶结砂岩颗粒之间的摩擦因数越小，颗粒之间的接触越容易。

(4) 随着黏结半径的减小，平行黏结更容易破坏，当黏结破坏时，颗粒的抗剪力主要由摩擦力提供，此时颗粒的胶结性状对岩样力学行为不再起主导作用，颗粒之间的运动主要由颗粒自身特性决定，其摩擦因数的敏感性就会增加。

(5) 当砂岩胶结半径比减小，砂岩可看作弱胶结岩体，此时外力作用下的岩样胶结物质很容易破坏，该条件下颗粒间的摩擦力就对岩样承力起着主导作用，颗粒摩擦因数对不同胶结性状的砂岩力学特性影响不可忽略，需要在数值计算中着重考虑。因此，为了更清晰地反映摩擦力与胶结力在砂岩受力中所起的作用，采用细观数值模型模拟胶结砂岩的力学响应时，需要考虑颗粒摩擦因数和胶结程度的变异性以反映胶结砂岩的受力变化过程，以准确掌握实际工程中胶结砂岩的力学特性以及破坏机理，以满足不同需求的实际工程需求。

参考文献：

[1] Potyondy D. A bonded-disk model for rock: Relating microproperties and macroproperties[C]// Cook B K, Jensen R P. Discrete Element Methods: Numerical Modeling of Discontinua. Santa Fe, New Mexico, USA: ASCE, 2002: 340-345.

[2] 曾远, 周健. 砂土的细观参数对宏观特性的影响研究[J]. 地下空间与工程学报, 2008, 4(3): 499-504.

ZENG Yuan, ZHOU Jian. Influence of micro parameters of sandy soil on its macro properties[J]. Chinese Journal of Underground Space and Engineering, 2008, 4(3): 499-504.

[3] Cui L, O’Sullivan C. Exploring the macro- and micro-scale response of an idealized granular material in the direct shear apparatus[J]. Geotechnique, 2006, 56(7): 455-468.

[4] 尹成薇, 梁冰, 姜利国. 基于颗粒流方法的砂土宏细观参数关系分析[J]. 煤炭学报, 2011, 36(S2): 264-267.

YIN Chengwei, LIANG Bing, JIANG Liguo. Analysis of relationship between macro-micro-parameters of sandy soil based on particle flow theory[J]. Journal of China Coal Society, 2011, 36(S2): 264-267.

[5] 陈建峰, 李辉利, 周健. 黏性土宏细观参数相关性研究[J]. 力学季刊, 2010, 31(2): 304-309.

CHEN Jianfeng, LI Huili, ZHOU Jian. Study on the relevance of macro-micro parameters for clays[J]. Chinese Quarterly of Mechanics, 2010, 31(2): 304-309.

[6] Potyondy D O, Cundall P A. A bonded-particle model for rock[J]. International Journal of Rock Mechanics & Mining Sciences, 2004, 41 (8): 1329-1364.

[7] Utili S, Nova R. DEM analysis of bonded granular geomaterials[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2008, 32(17): 1997-2031.

[8] 徐金明, 谢芝蕾, 贾海涛. 石灰岩细观力学特性的颗粒流模拟[J]. 岩土力学, 2010, 31(S2): 390-395.

XU Jinming, XIE Zhilei, JIA Haitao. Simulation of mesomechanical properties of limestone using particle flow code[J]. Rock and Soil Mechanics, 2010, 31(S2): 390-395.

[9] Hsieh Y M, Li H H, Huang T H. Interpretations on how the macroscopic mechanical behavior of sandstone affected by microscopic properties revealed by bonded-particle model[J]. Engineering Geology, 2008, 99(1): 1-10.

[10] Park J W, Song J J. Numerical simulation of a direct shear test on a rock joint using a bonded-particle model[J]. International Journal of Rock Mechanics & Mining Sciences, 2009, 46(8): 1315-1328.

[11] Potyondy D O. Simulating stress corrosion with a bonded- particle model for rock[J]. International Journal of Rock Mechanics & Mining Sciences, 2007, 44(5): 677-691.

[12] Wang Y H, Leung S C. Characterization of cemented sand by experimental and numerical investigations[J]. Journal of Geotechnical Geoenvironment Engineering, 2008, 34(7): 992-1004.

[13] Jiang M J, Yan H B, Zhu H H, et al. Modeling shear behavior and strain localization in cemented sands by two-dimensional distinct element method analyses[J]. Computers and Geotechnics, 2011, 38(1): 14-29.

[14] Itasca咨询公司. PFC3D 3.10三维颗粒流程序[M]. 3版. 明尼阿波列斯: 美国明尼苏达大学, 2005: 6-34.

Itasca Consulting Group. PFC3D 3.10 particle flow code in three dimensions[M]. 3rd ed. Minneapolis: University of Minnesota, USA, 2005: 6-34.

[15] 倪骁慧, 朱珍德, 赵杰, 等. 岩石破裂全程数字化细观损伤力学实验研究[J]. 岩土力学, 2009, 30(11): 3283-3290.

NI Xiaohui, ZHU Zhende, ZHAO Jie, et al. Meso-damage mechanical digitalization test of complete process of rock failure[J]. Rock and Soil Mechanics, 2009, 30(11): 3283-3290.

[16] 许江, 陆丽丰, 杨红伟, 等. 剪切荷载作用下砂岩细观开裂扩展演化特征研究[J]. 岩石力学与工程学报, 2011, 30(5): 944-951.

XU Jiang, LU Lifeng, YANG Hongwei, et al. Study of evolution law of mocrofacturing process of sandstone under shear loading[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(5): 944-951.

(编辑  杨幼平)

收稿日期：2012-12-14；修回日期：2013-03-05

基金项目：国家自然科学基金资助项目(51109231)；国家山区公路工程技术研究中心开放基金资助项目(gsgzj-2012-10)

通信作者：许明(1975-)，男，湖北洪湖人，博士后，副教授，从事岩土工程研究；电话：15696120975；E-mail：foretech@163.com