岩体作为由完整岩块和结构面不连续特征组成的离散介质[1]，完整岩块通过与结构面不连续特征相互作用实现力的传递，故岩体结构面不连续特征的几何和力学属性对岩体的力学行为具有重要的影响作用[2]。大量的工程实践表明[3-4]，岩体的失效主要是沿结构面发生剪切滑移，并非由完整岩块失效引起的。一般情况下，岩体结构面处于非粘结状态，不能承受拉力，因此压剪应力作用下岩体结构面的剪切力学行为备受关注[5-6]。
由于岩体结构面形貌粗糙度是控制结构面剪切力学行为的一个关键因素，因此大多数剪切强度关系式的建立都考虑了结构面粗糙度的影响作用。通过对锯齿形凸起结构面的直剪试验研究，PATTON[7]提出了两段线性包络线的双线性剪切强度准则。法向应力小于过渡应力时，结构面以剪切滑移为主，强度符合库伦-滑移准则；法向应力大于过渡应力时，结构面发生剪切啃断，强度符合库伦-摩尔准则。尽管该准则考虑了结构面剪切滑移和啃断两个关键剪切失效特征，但天然结构面剪切时滑移和啃断是同时发生的。LADANYI et al[8]基于大量的锯齿形节理直剪试验提出了同时考虑结构面剪切滑移和啃断特征的剪切强度准则，但该准则中一些参数测量难度较大。通过对天然结构面的直剪试验研究，BARTON et al[9]提出了结构面JRC-JCS剪切强度准则，并得到广泛应用。此外，还有许多学者提出了不同的剪切强度准则[10-17]。目前，由于缺乏能够准确定量描述结构面粗糙度的方法和揭示剪切过程中凸起失效的机理，所以仍然没有一个普遍适用的岩体结构面剪切强度准则。
数值模拟方法在岩石力学领域的快速发展，尤其是颗粒流的广泛应用，为研究结构面剪切力学行为提供了一种有效手段。在颗粒流中，完整岩块由粘结在一起的刚性颗粒集合表征，其力学行为由刚性颗粒间粘结状态控制，粘结破坏形成微裂纹，微裂纹聚合成宏观裂隙[18-19]。可见，颗粒流程序具备模拟分析结构面凸起失效机理的能力，进而成为研究结构面粗糙致滑剪切过程的力学行为的有利工具，但采用颗粒流研究结构面剪切力学行为的方法及其可行性并未得到验证。基于此，本文对比分析了颗粒流模拟结构面剪切力学行为的方法，提出并解决了模拟中遇到的问题，给出了准确模拟结构面剪切力学行为的方法。
1 岩体结构面颗粒流模拟的两类方法
1.1 粘结解绑法
粘结解绑法是将结构面两侧一定范围内颗粒间的粘结设置为非粘结状态，如图1所示。图中黄色线段表示颗粒处于粘结状态，红色线段表示颗粒处于非粘结状态。早期，该方法被广泛用来模拟岩体结构面。颗粒的摩擦系数对结构面的剪切力学行为具有重要的影响，其值的确定应根据结构面的摩擦角进行标定。
图1 岩体平直结构面粘结解绑法示意图
Fig.1 Schematic diagram of bond debonding method to
simulate flat rock discontinuities
采用粘结解绑法的岩体平直结构面直剪模拟结果如图2所示。与理论结果对比分析表明，粘结解绑法不能够准确模拟岩体结构面剪切力学行为。粘结解绑法的岩体平直结构面剪应力达到峰值剪应力后发生快速下降，而实际平直结构面在达到峰值剪应力后将按照库伦准则发生滑移，剪应力保持恒定不变。此外，粘结解绑法的岩体平直结构面在剪切过程中发生了明显的剪胀，这与实际直剪无剪胀特征不符。模型中颗粒的尺寸和位置随机分布将使得平直结构面具备一定的粗糙度，且颗粒只能沿着与其相接触颗粒的边缘运动，这是造成上述模拟结果与实际岩体平直结构面剪切力学行为不符的根本原因。
图2 岩体平直结构面粘结解绑法直剪模拟与理论结果对比
Fig.2 Results of direct shear test simulation and theoretical
for bond debonding method
1.2 光滑节理接触模型
为准确模拟结构面，开发了光滑节理接触模型[20]。将光滑节理接触模型赋予结构面两侧颗粒间的接触，可模拟处于粘结和非粘结状态的结构面。施加光滑节理接触模型后，颗粒可相互嵌入和穿越，不必沿着颗粒边缘运动。
采用光滑节理接触模型，岩体平直结构面直剪结果如图3所示。其中模型中最小颗粒直径为0.5 mm，剪切位移小于最小颗粒直径0.5 mm时，模拟剪应力与理论剪应力吻合；剪切位移大于最小颗粒直径0.5 mm时，模拟剪应力背离理论剪应力而增大，不能准确模拟平直结构面直剪力学行为。
图3 岩体平直结构面光滑节理接触法直剪模拟和理论值
Fig.3 Results of direct shear test simulation and theoretical
for smooth joint contact method
光滑节理接触法剪切位移为1 mm时，岩体平直结构面力链分布并非均匀分布，存在两个应力集中部位，如图4所示。两个部位应力集中造成剪应力快速增大，并导致颗粒间平行粘结破坏产生微裂纹。岩体平直结构面剪切后两个部位发生的应力集中是由剪切后结构面上接触模型发生变化而导致的。
图4 岩体平直结构面光滑节理接触法直剪前后接触模型
Fig.4 Contact model before and after direct shear
by smooth joint contact method
如图4(a)，剪切前，橙色颗粒(按颗粒中心定位)位于结构面上部，粉色颗粒位于结构面下部，橙色颗粒与其接触且位于结构面下部的灰色颗粒间为光滑节理接触；剪切1 mm后，粉色颗粒与橙色颗粒相接触，粉色颗粒和橙色颗粒均位于结构面下部，它们之间的接触不与结构面相交，因此它们之间为平行粘结接触。由于平行粘结接触刚度大于光滑节理接触，故产生应力集中进而引起剪应力增大。
如图4(b)，剪切前，橙色、紫色、绿色颗粒均位于结构面下部，红色颗粒位于结构面上部；剪切1 mm后，绿色颗粒与红色颗粒相接触，由于绿色颗粒和红色颗粒均位于结构面下部，它们的接触不和结构面相交，故它们之间接触为平行粘结接触。此外，虽然紫色、橙色、绿色三个颗粒剪切前之间的接触为平行粘结，但剪切1 mm后，紫色颗粒位于结构面上部，橙色颗粒和绿色颗粒位于结构面下部，紫色和橙色颗粒的接触与结构面相交，紫色和绿色颗粒的接触与结构面相交，与结构面相交的接触应为光滑节理接触模型，而实际上紫色和橙色颗粒的接触、紫色和绿色颗粒的接触仍为剪切前的平行粘结接触。由于平行粘结刚度大于光滑节理接触刚度，因此上述平行粘结接触存在多部位应力集中，并将引起剪应力增大。
2 岩体结构面模拟改进方法
由1.2小节可知，造成岩体平直结构面模拟剪应力背离理论值增大的原因是结构面两侧颗粒在剪切过程中容易穿越结构面而引起颗粒间接触模型的不合理设置。因此，可采用结构面相交接触判别法和结构面两侧颗粒组接触判别法来解决光滑节理接触法存在的问题。
结构面相交接触判别法是指剪切过程中随时判断颗粒间接触与结构面的相交情况，如果颗粒间接触与结构面相交，则将该接触赋予光滑节理模型。该方法可确保颗粒穿越结构面时，结构面所在位置不存在应力集中现象，从而避免剪切过程中剪应力背离理论值而增大的问题。采用该方法模拟法向应力1 MPa，摩擦角57°的岩体平直结构面直剪结果如图5所示，模拟结果与理论值吻合。虽然该方法能准确模拟结构面剪切力学行为，但颗粒剪切过程中穿越结构面，强制将该颗粒与位于结构面另一侧相接触颗粒间的初始平行粘结设置为光滑节理接触模型的做法不符合实际岩石破坏情况。
图5 岩体平直结构面相交接触判别法直剪结果
Fig.5 Direct shear results of discrimination intersection of
contact and discontinuity method
结构面两侧颗粒组接触判别法是指将结构面两侧的颗粒分别分组(见图6中红色颗粒组和绿色颗粒组)。剪切过程中随时判定结构面两侧颗粒的接触情况，若新接触为两个颗粒组内的颗粒相接触，则设置该新接触为光滑节理接触模型。该方法可有效避免剪切过程中结构面非实际的应力集中现象，采用该方法模拟法向应力1 MPa，摩擦角57°的岩体平直结构面直剪结果见图7，模拟结果吻合理论值。
图6 结构面两侧颗粒组接触判别法示意图
Fig.6 Schematic diagram for discrimination of particle groups
on both sides of discontinuity method
图7 岩体平直结构面两侧颗粒组接触判别法直剪结果
Fig.7 Direct shear results for discrimination of particle groups
on both sides of discontinuity method
3 结构面两侧颗粒组接触判别法可行性验证
为了进一步验证结构面两侧颗粒组接触判别法的可行性，对摩擦角为57°的岩体结构面分别进行了法向应力为3 MPa、5 MPa和7 MPa的直剪模拟，模拟结果如图8所示。由图8可知，随着法向应力的增大，峰值剪切位移增大；不同法向应力峰值剪切强度符合库伦准则，模拟结构面摩擦角为57.3°，与实际57°相差甚小，以上结果表明结构面两侧颗粒组接触判别法可有效模拟再现岩体结构面的剪切力学行为。
图8 岩体平直结构面两侧颗粒组接触判别法
不同法向应力直剪结果
Fig.8 Results of direct shear under different normal stress for
discrimination of particle groups method
4 结论
本文采用颗粒流程序重点研究了岩体结构面直剪力学行为模拟方法，指出了粘结解绑法和光滑节理接触法模拟结构面剪切力学行为的弊端，提出了结构面相交接触判别法和结构面两侧颗粒组接触判别法两种准确模拟结构面剪切失效的方法，具体结论如下：
1) 采用粘结解绑法模拟岩体结构面时，模型中颗粒的尺寸和位置随机分布引起一定的粗糙度，进而造成结构面剪切过程中产生应力集中现象，不能准确模拟结构面剪切力行为。
2) 采用光滑节理接触模拟岩体结构面时，结构面两侧颗粒在剪切过程中容易穿越结构面而引起颗粒间接触模型的不合理设置，进而使结构面产生非实际应力集中现象，造成模拟失败。因此可采用结构面相交接触判别法和结构面两侧颗粒组接触判别法解决光滑节理接触法存在的弊端。
3) 虽然结构面相交接触判别法能准确模拟结构面剪切力学行为，但颗粒剪切过程中穿越结构面，强制将该颗粒与位于结构面另一侧相接触颗粒间的初始平行粘结设置为光滑节理接触模型的做法不符合实际岩石破坏情况。
4) 摩擦角57°的岩体结构面在法向应力为3 MPa、5 MPa和7 MPa的直剪模拟结果表明了结构面两侧颗粒组接触判别法可有效模拟再现岩体结构面剪切力学行为。
[1] 夏才初,孙宗颀.工程岩体节理力学[M].上海：同济大学出版社,2002.
[2] HØEG B K.Numerical modelling of block size effects and influence of joint properties in multiply jointed rock[J].Tunnelling and Underground Space Technology,1997:407-415.
[3] BARTON N,BANDIS S.Some effects of scale on the shear strength of joints[J].International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts,1980,17(1):69-73.
[4] TSE R,CRUDEN D M.Estimating joint roughness coefficients[J].International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts,1979,16(5):303-307.
[5] 吴月秀,刘泉声,刘小燕.岩体节理粗糙度系数与统计参数的相关关系研究[J].岩石力学与工程学报,2011,29(S1):2593-2598.
WU Y X,LIU Q S,LIU X Y,et al.Study on relationship between joint roughness coefficient and statistical parameters[J].Chinese Journal of Rock Mechanics and Engineering,2011,29(S1):2593-2598.
[6] LI Y,ZHANG Y.Quantitative estimation of joint roughness coefficient using statistical parameters[J].International Journal of Rock Mechanics & Mining Sciences,2015,77:27-35.
[7] PATTON F D.Multiple modes of shear failure in rock[C].Proceeding of the1st Congress of International Society of Rock Mechanics,Lisbon,Portugal,1966:509-515.
[8] LADANYI B,ARCHAMBAULT G.Direct and indirect determination of shear strength or roch mass[C].Society of Mining Engineers of AIME,AIME Annual Meeting,Las Vegas,Nevada,1980:80-95.
[9] BARTON N R,CHOUBEY V.The shear strength of rock joints in theory and practice[J].Rock Mechanics and Rock Engineering,1977,10(1):1-54.
[10] XIA C C,TANG Z C,XIAO W M,et al.New peak shear strength criterion of rock joints based on quantified surface description[J].Rock Mechanics and Rock Engineering,2014,47(2):387-400.
[11] 夏才初,唐志成,宋英龙.基于三维形貌参数的偶合节理峰值抗剪强度公式[J].岩石力学与工程学报,2013,31(S1):2833-2839.
XIA C C,TANG Z C,SONG Y L.The formula of peak shear strength of coupling joints based on three dimensional morphology parameters[J].Chinese Journal of Rock Mechanics and Engineering,2013,31(S1):2833-2839.
[12] 王小江.岩石结构面力学及水力特性实验研究[D].武汉：武汉大学,2013.
[13] 孙辅庭,佘成学,万利台,等.基于三维形貌特征的岩石节理峰值剪切强度准则研究[J].岩土工程学报,2014,36(3):529-536.
SUN F T,SHE C X,WAN L T,et al.Study on criterion of peak shear strength of rock joints based on three-dimensional topography[J].Chinese Journal of Geotechnical Engineering,2014,36(3):529-536.
[14] 唐志成,夏才初,宋英龙.粗糙节理的峰值抗剪强度准则[J].岩土工程学报,2013,35(3):177-183.
TANG Z C,XIA C C,SONG Y L.Criteria for peak shear strength of rough joints[J].Chinese Journal of Rock Mechanics and Engineering,2013,35(3):177-183.
[15] 唐志成,黄润秋,张建明,等.含坡度均方根的节理峰值剪切强度经验公式[J].岩土力学,2015,36(12):3433-3438.
TANG Z C,HUANG R Q,ZHANG J M,et al.Empirical formula for peak shear strength of joints with root mean square slope[J].Rock and Soil Mechanics,2015,36(12):3433-3438.
[16] 葛云峰.基于BAP的岩体结构面粗糙度与峰值抗剪强度研究[D].武汉：中国地质大学,2014.
[17] 蔡毅.岩体结构面粗糙度评价与峰值抗剪强度估算方法研究[D].武汉：中国地质大学,2018.
[18] CUNDALL P A,POTYONDY D O.A bonded-particle model for rock[J].International Journal of Rock Mechanics and Mining Sciences,2004,41(8):1329-1364.
[19] 耿毅德,梁卫国,刘剑,等.高温及三轴应力条件下抚顺油页岩渗透规律试验研究[J].太原理工大学学报,2019,50(3):272-278.
GENG Y D,LIANG W G,LIU J,et al.Experimental study on permeability law of Fushun oil shale under high temperature and triaxial stress[J].Journal of Taiyuan University of Technology,2019,50(3):272-278.
[20] IVARS D M,PIERCE M E,DARCEL C,et al.The syntheti crock mass approach for jointed rock mass modelling[J].International Journal of Rock Mechanics and Mining Sciences,2011,48(2):219-244.
