基于Plesha本构的岩石节理多层结构模型研究

第 31 卷第 7 期 岩 土 力 学 Vol.31 No.7 2010 年 7 月 Rock and Soil Mechanics Jul. 2010 收稿日期：2009-06-25 第一作者简介：彭从文，男，1969 年生，博士，副教授，主要从事岩土工程教学科研工作。E-mail:ypengcw@126.com 文章编号：1000－7598 (2010) 07－2059－13 基于 Plesha 本构的岩石节理多层结构模型研究 彭从文1, 2，朱向荣1，王金昌3，许成祥2 （1.浙江大学 软弱土与环境土工教育部重点实验室，杭州 310027； 2.长江大学 城市建设学院，湖北 荆州 434023；3.浙江大学 交通工程研究所，杭州 310027） 摘 要:在 Dong 等工作的基础上[1]，研究了基于 Plesha 本构的岩石节理多层结构模型。模型将节理面分解为不同层次细观结 构面，每层结构面代一级粗糙度，力学响应发生在最底层（基本面），上层结构面受力性状由其下层结构面平均化得到， 依次进行得到节理宏观力学性状。节理破坏由下而上分层次进行，基本面破坏后，其上层结构面转化为基本面，基本面采用 Plesha 节理模型。模型能模拟弹性变形、滑动变形、磨损、剪断、压碎、分离等作用机理。定义的接触因子与绕流因子能考 虑节理接触面积及粗糙度变化对渗流的影响。通过将粗糙度定义为等效起伏角，可以模拟节理循环剪切性状。采用 ABAQUS 的用户子程序 UEL 进行了模型验证与参数。 关 键 词：岩石力学；节理；多层结构模型；有限元法；ABAQUS 中图分类号：TU 452 文献标识码：A A hierarchical model of rock joints based on Plesha’s constitutive model PENG Cong-wen1, 2, ZHU Xiang-rong1, WANG Jing-chang3, XU Cheng-xiang2 (1.Key Laboratory of Soft Soils and Environmental Engineering, Ministry of Education, Zhejiang University, Hangzhou 310027, China; 2. School of Urban Construction, Yangtze University, Jingzhou, Hubei 434023, China; 3. Research Institute of Transportation Engineering, Zhejiang University, Hangzhou 310027, China) Abstract: An improved hierarchical model of rock joints was proposed based on the work of Dong et al. (Dong J J, Pan Y W. A hierarchical model of rough rock joints based on micromechanics[J]. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts, 1996, 33(2): 111-123). According to this model, joint macro-asperity is decomposed into multi-level micro-asperity; joint behavior is decided by the performance of the lowest level asperity(also basic surface) through homogenization. Joints fail progressively from the lower level to higher level asperity, after the basic surface break down, the second lowest level asperity converts to the basic surface. The roughness of joint is represented by equivalent asperity angle. The model also include the effect of variation of contact area and roughness on the permeability of joints. The model can simulate the mechanism of elastic deformation, sliding, damage, shear-off, crush and separation of joints. Model verification and parameters sensitivity analysis were carried out by UEL of commercial soft package ABAQUS. Key words: rock mechanics; joints; hierarchical model; FEM; ABAQUS 1 引 言 自 Londe 和 Sabarly 首次进行水-力耦合试验以 来[2]，在过去的几十年中，对节理的水力性能进行 了大量的试验、理论与数值研究[3－25]。影响节理 渗透性能的因素很多，其中最关键的是节理粗糙 性[26, 27]，而剪胀与突台（asperity）损伤是造成渗流 变化的主要原因[4, 26]。 目前，计算节理渗透系数有两种。一是直 接法，即通过试验，将渗透系数表示为应力的函 数[28－30]；二是间接法，即采用立方体定理，并结合 节理粗糙性进行修正[2, 5－7, 31－39]，修正因子包括力学 开度与水力开度的关系[5, 7, 37]、粗糙度影响[2, 4, 33]、 接触面积影响[34－35]、岩粉影响[38－39]等。 节理模型可分为经验与理论两大类。经验模型 主要有双曲线模型[40－42]、指数模型[43－44]，对数模 型 [45]、 JRC-JCS 模型 [5, 46 － 48]，连续屈服模型 （CYJM）[49]等。理论模型采用弹塑性理论并考虑节 理破坏机理。最早应用的是 Hertz 接触理论， Greenwood 和 Williamson 采用 Hertz 理论研究了平 岩 土 力 学 2010 年 板与粗糙面接触特性[50]；Brown 和 Scholz 研究了两 个粗糙面间的相互作用[51]；Walsh[34]、Swan[30]、 Lespinasse 等[52]将 Hertz 理论用于节理水-力耦合分 析，Misra、Mroz 考虑了剪应力的影响和接触面的 微滑动效应[53－55]。以上模型没有考虑突台间的相互 作用，Hopkins[56]、Pyrak-Nolte 等[57]将突台拟为圆 柱体，在考虑突台间相互影响情况下，研究了节理 承受法向荷载的受力性状。Gangi 研究水-力耦合作 用时提出了钉床模型[28]。以上理论模型不能考虑节 理剪胀作用。Patton 将粗糙面拟为齿形结构，建立 了基于单个突台的双线性剪切强度模型[58]，Jeager 改进了双线性模型的不连续性 [59]，Ladanyi 和 Archambault，Saeb 和 Amadei 考虑了多个突台的情 形[60－61]。Plesha 在 Patton 模型的基础上考虑了节理 细观结构的影响，并引入了损伤[62－63]，该模型结构 简明，物理意义明确，能较好地反应节理剪胀与损 伤效应，随后被许多学者广为采用[38－39, 64－69]，其中， Selvadurai 将其应用于水-力耦合分析[38－39]。 研究表明：岩石节理具有分形特征[70－72]，节理 面可以分解为不同波长与波幅的结构面[71－73]， 也就是说，节理宏观粗糙度是不同层次粗糙度的综 合反应。Patton 首先提出了主要粗糙度（primary asperity）与次级粗糙度的概念（second asperity）。 一般认为，次级粗糙度影响峰值应力及剪切方向 性[69, 74]；小位移下剪切性状主要受次级粗糙度控 制，大位移时受主要粗糙度控制[58, 75]；低法向应力 下，剪切性状受次级粗糙度影响，法向应力增大时， 受主要粗糙度控制[73, 76]。Jing 等[77]提出了如图 1 所 示的节理多粗糙度剪切概念模型，剪切过程中，次 级粗糙度先于主要粗糙度发生破坏，并且剪切具有 方向性。Gentier 等[78]研究了节理几何形态对剪切性 能的影响，图 2 是剪切前后节理剖面变化图。按照 多级粗糙度概念，可以将该节理面简化成两级粗糙 度模型，剪切后，次级粗糙度破坏，主粗糙度也发 生一定的损伤。由此可见，要从细观角度研究节理 受力机理，模型宜考虑粗糙度的结构性。Lee 等显 式考虑了次级粗糙度的影响[69]，但假定主次粗糙度 同时损伤；Jing 等则采用了塑性力学中的软硬化参 数[67]；Dong 等将节理面分解为如图 3 所示的宏细 观多层结构模型[1]（下称 Dong 模型），将节理破坏 拟为不同层次结构面渐进破坏的过程。模型考虑了 弹性变形、滑动、突台剪断、界面磨损及分离等作 用机理。 本文在Dong模型的基础上作了进一步的工作。 基本面改用 Plesha 节理模型，增加了结构面转化准 则，定义了等效起伏角、绕流因子与接触因子。模 型不仅能够更好地模拟节理循环剪切特性，还能考 虑节理形态变化对水力开度的影响。模型通过 ABAQUS平台的用户子程序UEL进行了数值实现， 并将计算结果与试验值进行了比较，最后进行了模 型参数分析。 图 1 节理剪切概念图（摘自 Jing,1992） Fig.1 Conceptual model of joint shear behavior (after Jing, 1992) (a) 剪切前 (b) 剪切后 图 2 节理剖面（摘自 Gentier,2000） Fig.2 Joint profiles(after Gentier,2000) 2 节理多层结构面模型 2.1 Dong 模型简介 模型对象为非填充的硬性节理面，采用齿形形 态、初始状态为完全啮合。据节理粗糙度的多重性 与分形特征，将节理面分解为如图 3 所示的宏细观 多层结构模型。每层结构面代表一级粗糙度，最底 层为实际接触面，称为基本面，每层结构面可以有 不同的起伏角与方位。力学响应发生在基本面，上 一层结构面性状由其下层结构面平均化后得到，依 次进行得到节理宏观力学性状。结构面损伤与破坏 也自下而上进行，基本面破坏后，上一层结构面转 化为基本面。 模型基本假定：（1）达到基本面剪切强度之前， 节理呈弹性变形，剪切强度采用 Mohr-Coulomb 准 则；（2）达到界面剪切强度后，开始滑动；（3） 任一层平均剪切应力超过节理壁面强度后，开始损 伤。损伤包括基本面的磨损与突台剪断；（4）最低 层突台平均应力达到剪断剪应力后，该层剪断，舍 (c) 反向剪切 (d) 剪切后构形 (a) 剪切前构形 (b) 正向剪切 2060 第 7 期 彭从文等：基于 Plesha 本构的岩石节理多层结构模型研究 弃该结构面。节理破坏由下往上分层次进行，只有 下一层结构面破坏后，上一层结构面才发生损伤或 破坏。 基本面本构关系式为 c c c cd d di ij j ik u rσ = − （1） 式中： iσ c、 juc分别为基本面上接触应力与相对位移； c ijk 为基本面上接触刚度； cd ir 为考虑界面滑移与分 离时的应力放松张量，滑移时，其剪应力分量为剪 应力中超过剪应力强度部份值，分离时，其法向分 量是拉应力中超过拉伸强度部份值。 2.2 对 Dong 模型的改进与拓展 2.2.1 基本面采用 Plesha 节理模型 式（1）是带放松因子的线弹性模型[79]，从机 理上,不宜描述节理面的接触滑移，也不能考虑基本 摩擦角对剪胀的贡献，为使模型更具通用性，本文 在基本面改用基于弹塑性理论的 Plesha 节理模型。 2.2.2 增加结构面的转化准则 节理损伤或破坏与其形态及荷载有关，低压时， 表现为粗糙面磨损，高压力时为突台剪坏；起伏角 较小时表现为磨损，大角度为剪坏；接触面较小时， 还会出现局部压碎[65, 75, 80]。因此，除了剪断之外， 还应考虑接触面磨损及突台压碎导致的结构面转 化。本文中假定：（1）达到界面剪切强度后，开始 滑动，同时产生滑动磨损，采用 Plesha 损伤模型， 累计损伤达到限值后，舍弃该基本面；（2）法向应 力达到材料单轴抗压强度时，发生压碎破坏，同时 舍弃该基本面。 2.2.3 考虑节理形态变化对循环剪切性状的影响 节理粗糙度指标采用等效起伏角。图 4 为单层 结构面示意图，假定左右起伏角初始值均为 1α ，剪 切后左右起伏角分别为 2α 、 1α ，据面积等效的原则， 等效起伏角 1 2 1 2 2sin( )sin( )arctan sin( ) α αα α α ⎛ ⎞= ⎜ ⎟+⎝ ⎠ 。 2.2.4 考虑节理形态变化对水力开度的影响 考虑接触面积与绕流两个影响因素。如图 4 所 示，初始状态啮合（上突台面 BDF 与下突台面 ACE 重合），剪切发生后，上下面发生错动，接触面变 为 BC 段。采用 Walsh 的方法[34]，该层结构面接触 因子为 s 1.0 / 1.0 / k BC AEr BC AE −= + ；绕流因子为 l 1.0 kr = − ( )AC CE AE AE + − 。总体接触与绕流因子为各层因子 乘积： 1 2s s s snr r r r= " ， 1 2l l l lnr r r r= " 。粗糙度增大时， 结构层增多，绕流因子变小，接触因子变大；与此 相反，当粗糙消失时（完全光滑），变为平板结构， 绕流因子变为 1，接触因子变为 0。渗透系数采用立 方体定理计算，水力开度为 0 s le e r r e= + ∆ ， 0 , e e∆ 分 别为残余水力开度（hydraulic aperture）与力学开度 （mechanical aperture）增量。 图 3 节理粗糙度多层结构示意图（摘自 Dong,1996） Fig.3 Hierarchical structure of joint asperity (after Dong, 1996) 图 4 等效起伏角计算示意图（单层结构） Fig.4 Schematic diagram of equivalent asperity angle(single structure) 图 5 多层结构应力转换关系图 Fig.5 Stress transformation of hierarchical structure 2.3 本构关系 2.3.1 弹性矩阵 图 5 为基本面（第 1 层）与第 2 层结构面关系 示意图。基本面相对于其上层结构面的起伏角为 1φ （每个接触面可以有不同的起伏角），上下接触面初 始啮合（AC 段与 BD 段重合），滑移后，接触面变 为 BC 段。弹性变形时，第 i 个基本面应力-相对位 移关系式为 (1) (1) (1) i iK= uσ （2） 式中： (1)iσ 、 (1)iu 分别为基本面上接触应力、相对位 移列阵；刚度矩阵 (1) n s 0 0 k k ⎡ ⎤= ⎢ ⎥⎣ ⎦K ，其中刚度系数 α1 α1 α1 α2 Bα2 α2 α1 A C D F E α α α (2) nσ (2) sτ (2) nσ 1φ (2) sτ EA B C D (2) nσ (1)nσ (2) sτ (1) sτ A B C 2061 岩 土 力 学 2010 年 是材料常数，与荷载及节理形态无关。基本面合力 为 (1) (1) (1) (1) i i i A= ∑F K u （3） 式中： (1)F 为基本面合力； (1)iA 为第 i 个基本面接触 面积，即图 4 中 BC 段。采用 Huang 和 Chang 的平 均化方法[68, 79]，得到第 2 层结构第 K 个接触面平均 应力与相对位移关系为 (2) (1) (1) (1) (1) (2) (2) (2) (2) 1 ( )tk i i i k k ik A A = =∑ T K T u K uσ （4） (2) (1) (1) (1) (1) (2) 1 t i i i ik A A = ∑K T K T （5） 式中： (2)kA 为第 2 层结构第 K 个接触面的名义接触 面积，即图 4 中 AE 段。 (1)iT 与 (1)tiT 分别为第 i 个基 本面与第 2 层结构面的转换矩阵与相应的逆矩阵， 其中 1 1(1) 1 1 cos sin sin cosi φ φ φ φ −⎡ ⎤= ⎢ ⎥⎣ ⎦ T （6） 依此类推，第 l 层结构第 K 个接触面平均应力 与相对位移关系式为 ( ) ( ) ( )l l l k k= K uσ （7） ( ) ( 1) ( 1) ( 1) ( 1) ( ) 1l l t l l l i i il ik A A − − − −= ∑K T K T （8） 若节理共分解为m 层结构，其宏观刚度矩阵为 ( 1) ( 1) ( 1) ( 1)1 m t m m m i i i i A A − − − −= ∑K T K T （9） 式中： A为节理接触面积。 由式（9）可知，刚度矩阵与接触面积有关，因 此，刚度矩阵为非线性。若各层结构面消失，刚度 矩阵式（9）退化为式（2）。 2.3.2 弹塑性矩阵 基本面屈服后产生塑性流动，采用 Mohr- Coulomb 屈服准则 (1) (1) t n tancσ σ ϕ= + （10） 式中：c 、ϕ为黏聚力与基本摩擦角，是材料常数。 基本面与第 2 层结构面应力有如下关系 (2) (1) (1) (2) (1) (2) 21(1) A d A = =T Tσ σ σ （11） 式中： (1)A 为基本面接触总面积； (2)A 第 2 层结构 面名义接触总面积； 21d 定义为第 2 层结构面相对 于基本面的应力放大系数。 类似地，得到基本面与节理宏观应力关系为 (2) (3) (1) (1) (2) ( ) (1) (2) ( ) (1) (2) ( ) 21 32 0 m m m m A A A A A A d d d = ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ = ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ T T T T T T σ σ σ （12） 将上式代入式（10），得到材料屈服准则。若 取 2 层结构，应力显式为 (1) 21 02 2 1 2 1 n 2 1 2 1 (1) 21 02 2 1 2 1 2 1 2 1 n [(cos cos sin sin ) (cos sin sin cos ) ] [(cos cos sin sin ) (cos sin sin cos ) ] d d d d τ τ τ σ φ φ φ φ σ φ φ φ φ σ σ φ φ φ φ σ φ φ φ φ σ = − − + = − + + n （13） 将式（13）代入式（10），得屈服准则 n n 21 02 ( ) tan cF A B A B d dτ τ σ σ σ σ ϕ= + + − − （14） 式中： 2 1 2 1 2 1 2 1 cos cos sin sin cos sin sin cos A B φ φ φ φ φ φ φ φ = − ⎫⎬= + ⎭ （15） 式（14）大于等于 0 时，接触面发生滑移。采用非 关联流动法则，取流动势函数： nG A Bτσ σ= + （16） 损伤在滑动后产生，Plesha 认为，起伏角的损 伤与塑性剪切功有关[62]，提出了粗糙度损伤演化 式： p 0 exp( )DWφ φ= − （17） 式中： pW 为塑性剪切功。Hutson 和 dowding 验证 了 Plesha 模型[81]，并且认为，损伤参数还与法向应 力及材料抗压强度有关，将试验常数写成 n c D cσσ= ， n c,σ σ 分别为法向应力和材料单轴抗压强度， c 为 试验常数。将上述各式代入下式得弹塑性矩阵 T ep T T G F F G H ⎛ ⎞∂ ∂⎧ ⎫⎧ ⎫⎜ ⎟⎨ ⎬⎨ ⎬∂ ∂⎩ ⎭⎩ ⎭⎜ ⎟= −⎜ ⎟∂ ∂⎧ ⎫ ⎧ ⎫ −⎜ ⎟⎨ ⎬ ⎨ ⎬∂ ∂⎩ ⎭ ⎩ ⎭⎝ ⎠ K D K I K σ σ σ σ （18） p F GH W τ τ σ σ ∂ ∂= ∂∂ （19） 式中： I 为单位矩阵。 2.4 数值实现 有限元列式可参见文献 [82]。数值计算通过 ABAQUS 的用户子程序 UEL 实现。采用零厚度、 六节点、平面应变单元，应力更新采用显式算法[83]。 2062 第 7 期 彭从文等：基于 Plesha 本构的岩石节理多层结构模型研究 UEL 程序框图如图 6 所示。 图 6 UEL 计算流程图 Fig.6 Flow chart of UEL 3 模型验证及参数分析 3.1 模型验证 （1）循环剪切试验 Lee 等对 Hwangdeung 花岗岩和 Yeosan 大理岩 进行了不同压力下的循环剪切试验[69]，试件尺寸为 160 mm×120 mm×120 mm，粗糙度形态采用激光 扫描获取，采样间距 0.2 mm，法向应力为 0.2～ 5.0 MPa，剪切速率 0.05～0.08 mm/s。试验表明： 第 1 次剪切循环时有明显的峰值应力与非线性剪 胀，正反向剪切强度不相等，剪切强度与压力相关 等。试验中测量了两级起伏角，因此，计算采用双 层结构模型。参数见表 1 第 1 行。参数下标 1、2 分别代表第 1、2 层结构面， cσ 为材料单轴抗压强 度， dφ 为损伤角，当突台损伤后起伏角小于该值时， 舍弃该结构层， sφ 为剪断角，其他参数含义同前。 为了体现剪切的方向性，Lee 对正反向起伏角 初始值及损伤系数取不同的数值，过于复杂。本文 模型中每层结构只需一个初始起伏角与损伤系数， 其值与剪切方向无关，剪切过程中的各向异性通过 模型定义的等效起伏角体现。表中起伏角初始值取 Lee 文中正反向起伏角平均值，损伤系数进行了调 整。计算结果如图 7、8 所示，虚线为试验值，Ⅰ至 Ⅳ部份分别为正向剪切、正向回缩、反向剪切及反 向回缩段。由图可知，剪应力计算与试验值吻合较 好。由于次级粗糙度的作用，正向剪切呈显著的峰 值效应，随着次级粗糙度的破坏，剪切应力迅速减 小，此时，主粗糙度起控制作用，剪切应力衰减速 率取决于损伤参数；在回缩段，随着剪胀角不断增 大，剪切应力越来越小；在反向剪切段，由于正向 剪切造成等效起伏角减小，剪切峰值应力明显变小。 第 2 次循环剪切中，次级粗糙度已不存在，所以曲 线没有明显的峰值现象，仅表现为主粗糙度的损伤。 计算剪胀与试验值略有差异，这可能与计算模型有 关，本文计算采用平面应变模型，没有考虑空间效 应。由图可知，剪切位移越大，剪胀角越小；反向 剪切时剪胀角较正向剪切小。 表 1 数值模拟计算参数 Table 1 Input parameters for simulation 序 号 kn / 106 (kN/m3) kn / 106 (kN/m3) c / kPa ϕ / rad φ2 / rad φ1 / rad D2 / (kN/m) D1 / (kN/m) σn / 103 kPa σc / 105 kPa φd / rad φs / rad 1 54 5.4 0 0.60 0.22 0.28 1.8 0.100 1.0 1.50 0.001 1.0 2 20 9.0 0 0.54 0.15 0.18 1.0 0.020 2.0 2.15 0.001 1.0 3 60 12.0 0 0.54 0.20 0.13 1.0 0.012 4.0 2.15 0.001 1.0 4 60 9.0 0 0.54 0.08 0.15 1.0 0.030 2.0 2.15 0.001 1.0 5 60 9.0 0 0.54 0.09 0.18 1.0 0.030 2.0 2.15 0.001 1.0 (a) 第 1 次循环 (b) 第 2 次循环 图 7 循环剪切下剪应力-剪切位移关系曲线 Fig.7 Curves of shear stress versus shear displacement under cycle loading -16 -12 -8 -4 0 4 8 12 16 -1.5 -1.0 -0.5 0.5 1.0 1.5 剪 切 应 力 / M Pa 剪切位移/ mm 计算值 试验值 Ⅰ ⅡⅢ Ⅳ -16 -12 -8 -4 0 4 8 12 16 -1.5 -1.0 -0.5 0.5 1.0 1.5 剪 切 应 力 / M Pa 剪切位移/ mm 计算值 试验值 Ⅱ Ⅰ Ⅲ Ⅳ 2063 岩 土 力 学 2010 年 (a) 第 1 次循环 (b) 第 2 次循环 图 8 循环剪切下剪胀曲线 Fig.8 Dilation curves under cycle loading （2）剪切渗流耦合试验 Olsson 进行了 CNL 与 CNS 边界条件下的剪切 渗流耦合试验[7]，试验采用常水头，剪切试验前先 进行加卸载固结，最大固结压力等于剪切试验时的 法向应力。文中提供了基本摩擦角与节理粗糙度系 数等参数，本文采用双层结构模型与试验结果相比 较，参数见表 1 第 2～5 行。其中，主、次起伏角由 峰值附加摩擦角分解得到，采用 JRC-JCS 模型计算 峰值附加摩擦角；刚度系数参考岩块材料弹性模量， 剪切刚度随法向应力增加，原因在于物理模型与本 文模型的差异性，不同压力下试件接触面形态有所 区别，造成剪切阻力的差异。计算结果如图 9、10 所示，其中虚线为试验值。剪切应力计算与试验值 基本吻合。计算开度与试验值在曲线初始段有一定 差异，这是由于试件接触面并不完全啮合，因此， 试验过程中先剪缩，后剪胀。 3.2 参数分析 3.2.1 结构层数的影响 见图 11。计算了 1 层、2 层及 3 层结构模型的 节理剪切性状。法向应力为 1.0 MPa，法向刚度 54 GPa/m，切向刚度 5.4 GPa/m，摩擦系数为 0.5 ， 其余参数见图内数据。曲线中 A 为剪应力峰值点， B 为第 3 层结构面破坏点，C 为第 2 层结构面破坏 点，各曲线端点为 D。有以下几个特点： （1）剪切强度随结构层数的增加而显著增大， 随着剪切位移的增加，粗糙度逐渐破坏，剪应力终 值趋于一致，均受控于基本摩擦角。3 层结构模型 剪应力峰值较2层及1层结构模型分别增大约24%、 52%； （2）增加结构层数会减小整体刚度，导致剪应 力峰值位移增大。从 1 层～3 层结构模型，其峰值 位移分别是 0.5、0.8、1.9 mm； （3）开度随结构层数的增加而增大。对于力学 开度终值，3 层结构模型较 2 层及 1 层结构模型分 别增大 33%、56%；水力开度终值，分别增大 50%， 95%。从变化率上看，水力开度对结构层数的变化 更敏感； （4）力学开度与水力开度性状有差异。随着剪 切位移的增加，力学开度增大，而水力开度趋于稳 定并略有下降；另外，水力开度随结构面的破坏出 现一定的波动，结构层数越多，其波动幅度越大。 因此，与力学开度不同的是，水力开度峰值不一定 出现在剪切终值，与结构层数、粗糙度、损伤参数 等因素有关。本例中，从 3 层～1 层结构模型，水 力开度终值较峰值水力开度分别减少 9.2%、6.3%、 5.2%。 3.2.2 法向应力的影响 见图 12。计算参数见图内数据,其余同 3.2.1。 等效峰值摩擦系数定义为结构初始滑动时剪切应力 与法向应力比值，采用 3 层结构模型计算。由图可 知，达到次级粗糙度破坏压力前，峰值摩擦系数随 法向压力增加而略有增大；随着压力的增加，次级 粗糙度发生破坏，峰值摩擦系数随法向压力增加不 断减小。 图 9 剪应力-剪切位移关系曲线 Fig.9 Curves of shear stress versus shear displacement -16 -12 -8 -4 0 4 8 12 16 0.5 1.0 1.5 2.0 2.5 法 向 位 移 / m m 剪切位移/ mm 计算值 试验值 Ⅰ Ⅱ Ⅲ Ⅳ -16 -12 -8 -4 0 4 8 12 16 0.5 1.0 1.5 2.0 2.5 法 向 位 移 / m m 剪切位移/ mm 计算值 试验值 Ⅰ ⅡⅢ Ⅳ 0 2 4 6 8 10 12 14 16 1 2 3 4 5 剪 切 应 力 / M Pa 剪切位移/ mm （1）法向应力 2.0 MPa，边界刚度 0.0 kN/mm （2）法向应力 2.0 MPa，边界刚度 37 kN/mm （3）法向应力 2.0 MPa，边界刚度 75 kN/mm （4）法向应力 4.0 MPa，边界刚度 37 kN/mm ? 1? ? 2? ? 3? ? 4? 2064 第 7 期 彭从文等：基于 Plesha 本构的岩石节理多层结构模型研究 (a) CNL 边界，σn = 2.0 MPa (b) CNS 边界，σn = 4.0 MPa, K = 37 kPa/mm (c) CNS 边界，σn = 2.0 MPa, K = 37 kPa/mm (d) CNS 边界，σn = 2.0 MPa, K = 75 kPa/mm 图 10 力学开度、水力开度与剪切位移关系曲线 Fig.10 Curves of mechanical aperture and hydraulic aperture versus shear displacement (a) 剪应力-剪切位移 (b) 力学开度、水力开度与剪切位移关系曲线 图 11 多层结构对剪切性状的影响 Fig.11 The influence of hierarchical structure on the shear behavior 图 12 法向应力与等效峰值摩擦系数关系曲线 Fig.12 Curve of normal stress versus equivalent peak friction cofficient 3.2.3 起伏角、损伤系数、基本摩擦角、刚度系数的 影响 见图 13～16。计算采用 2 层结构模型，参数见 图内数据，其余同 3.2.1。有以下结论： （1）起伏角改变剪切强度，但不影响应力终值。 剪切强度及开度随起伏角增大而增加。起伏角从 0.1 rad 分别增加至 0.2 rad、0.3 rad 时，剪切强度分 别增大 28.5%、65%，力学开度分别增大 8.65%、 19.4%；水力开度分别增大 10.7%、26%。 0 1 2 3 4 5 6 0.2 0.4 0.6 0.8 1.0 1.2 1.4 力 学 开 度 和 水 力 开 度 / m m 剪切位移/ mm (1) 力学开度 (2) 水力开度 (1) (2) 0 1 2 3 4 5 6 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 力 学 开 度 与 水 力 开 度 / m m 剪切位移/ mm (1) 力学开度 (2) 水力开度 (1) (2) 0 2 4 6 0.2 0.4 0.6 0.8 1.0 1.2 1.4 力 学 开 度 与 水 力 开 度 / m m 剪切位移/ mm (1) 力学开度 (2) 水力开度 (1) (2) 1 3 5 7 0 1 2 3 4 5 6 0.2 0.4 0.6 0.8 1.0 1.2 1.4 力 学 开 度 与 水 力 开 度 / m m 剪切位移/ mm (1) 力学开度 (2) 水力开度 (1) (2) 0 2 4 6 8 10 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 φ1 = 0.28φ2 = 0.10φ3 = 0.10 剪 切 应 力 / M Pa 剪切位移/ mm (1) 3 层结构模型 (2) 2 层结构模型 (3) 1 层结构模型(1) (2) (3) A B C D C D1 = 0.1 m/kN D2 = 2.0 m/kN D3 = 2.0 m/kN A A 0 2 4 6 8 10 12 14 16 1 2 3 4 力 学 开 度 与 水 力 开 度 / m m 剪切位移/ mm (1) 3 层结构模型 (2) 2 层结构模型 (3) 1 层结构模型 (1) (2) (3) (1) (2) (3) 0 1 2 3 4 5 6 7 1.0 1.1 1.2 1.3 1.4 1.5 D1 = 0.1 m/kN D2 = 2.0 m/kN D3 = 2.0 m/kN 等 效 峰 值 摩 擦 系 数 法向应力/ MPa φ1 = 0.28 φ2 = 0.10 φ3 = 0.10 σc = 150 MPa A B C 2065 岩 土 力 学 2010 年 (a) 剪应力-剪切位移 (b) 力学开度、水力开度与剪切位移关系曲线 图 13 粗糙角对剪切性状的影响 Fig.13 The influence of asperity angles on the shear behavior (a) 剪应力-剪切位移 (b) 力学开度、水力开度与剪切位移关系曲线 图 14 损伤参数对剪切性状的影响 Fig.14 The influence of damage coefficients on the shear behavior (a) 剪应力-剪切位移 (b) 力学开度、水力开度与剪切位移关系曲线 图 15 基本摩擦角对剪切性状的影响 Fig.15 The influence of basic angle of friction on the shear behavior (a) 剪应力-剪切位移 (b) 力学开度、水力开度与剪切位移关系曲线 图 16 剪切刚度对剪切性状的影响 Fig.16 The influence of shear stiffness on the shear behavior 0 2 4 6 8 10 12 14 16 0.5 1.0 1.5 2.0 2.5 (1) φ2 = 0.3 (2) φ2 = 0.2 (3) φ2 = 0.1 剪 切 应 力 / M Pa 剪切位移/ mm (1) (2) (3) φ1 = 0.28 D1 = 0.1 m/kN D2 = 2.0 m/kN 0 2 4 6 8 10 12 14 16 1 2 3 力 学 开 度 与 水 力 开 度 / m m 剪切位移/ mm (1) (2) (3) 由上至下依次为(1), (2), (3) 0 2 4 6 8 10 12 14 16 0.5 1.0 1.5 Sn= 1.0 MPa 剪 切 应 力 / M Pa 剪切位移/ mm (1) D1 = 0 m/kN D2 = 0 m/kN (2) D1 = 0.1 m/kN D2 = 0.5 m/kN (3) D1 = 0.1 m/kN D2 = 1.0 m/kN (4) D1 = 0.1 m/kN D2 = 4.0 m/kN φ1 = 0.28φ2 = 0.10 (1) (2) (3) (4) 0 2 4 6 8 10 12 14 16 2 4 6 力 学 开 度 与 水 力 开 度 / m m 剪切位移/ mm (1) (2) (3) (4) (1) 由上至下依次为(2), (3), (4) 0 2 4 6 8 10 12 14 16 0.5 1.0 1.5 剪 切 应 力 / M Pa 剪切位移/ mm φ1 = 0.28φ2 = 0.10 D1 = 0.1 m/kN D2 = 2.0 m/kN(1) (2) (3) (1) ϕ1 = 0.4 (2) ϕ1 = 0.5 (3) ϕ1 = 0.6 0 2 4 6 8 10 12 14 16 1 2 3 力 学 开 度 与 水 力 开 度 / m m 剪切位移/ mm (1) (2) (3) 由上至下依次为(1), (2), (3) 2 4 6 8 10 12 14 160 200 400 600 800 1 000 1 200 剪 切 应 力 / k Pa 剪切位移/ mm (1) Ks = 5.4 GPa/m (2) Ks = 10.0 GPa/m (3) Ks = 20.0 GPa/m Kn = 54.0 GPa/m 从左至右依次为(3), (2), (1) 2 4 6 8 10 12 140 1 2 3 力 学 开 度 与 水 力 开 度 / m m 剪切位移/ mm 从左至右依次为(1), (2), (3) 从上至下依次为(1), (2), (3) 2066 第 7 期 彭从文等：基于 Plesha 本构的岩石节理多层结构模型研究 （2）损伤参数的变化不改变剪切强度，但改变 曲线形态，粗糙度没有损伤时，剪切性状表现为理 想塑性。随着损伤系数的增加，开度减小，水力开 度的波动性降低。相对于未损伤状态，当损伤系数 从工况 1 变化至工况 2～4 时（见图 14），力学开 度降至原来的 57.5%、50.3%、44.7%，水力开度降 至原来的 17.4%、13.6%、11.0%，水力开度终值较 水力开度峰值减小 29%、16.5%、6.2%。从数据上 看，损伤对水力开度的影响要大于力学开度。 （3）剪切强度随基本摩擦角的增加而增大，基 本摩擦角对开度的影响发生在粗糙度基本破坏之 后，开度值随基本摩擦角的增大而增加。 （4）刚度系数对剪切强度与开度影响各异。法 向刚度对力学开度与水力开度影响均较大，剪切刚 度对力学开度影响很小，对水力开度影响较大。法 向刚度由 10 kPa/m 增大至 100 kPa/m 时，剪切强度 增大 1.44%，力学开度与水力开度分别增加 61%与 56%；剪切刚度由 5.4 kPa/m 增加至 20 kPa/m 时， 剪切强度增加 1.28%，力学开度与水力开度分别减 小至原来的 6.6%和 48%。 4 结 论 本文模型是在Dong模型基础上的改进与拓展， Plesha 模型相当于单层结构模型，是本文模型的特 例，经与试验结果比较，模型能有效地描述节理剪 切性状。参数分析表明，模型力学特性与 Dong 模 型相似。剪切强度随着结构层次的增加而增大，随 基本摩擦角增加而增大；结构层次不变时，其等效 峰值摩擦系数随法向应力的增大而减小；剪应力终 值由基本摩擦角控制；剪应力-位移曲线形态取决于 粗糙度结构特性；模型中刚度系数为材料常数，刚 度对剪切强度影响很小。随着剪切位移的增加，力 学开度增大，而水力开度趋于稳定；水力开度对结 构层次与损伤系数的变化较力学开度更敏感；随着 结构面的破坏，水力开度出现波动性，其波动幅度 与结构层次及损伤系数有关。法向刚度对力学开度 与水力开度影响均较大，而剪切刚度对水力开度影 响较大，对力学开度影响很小。 随着压力的增加，次级粗糙度迅速损伤，剪切 性状主要由粗糙度及基本摩擦角控制，多层结构已 转化为单层结构，因此，本模型适用于中低压力下 的脆性材料。另外，假定模型初始状态为完全啮合， 这与实际情况不完全相符。本模型不仅适用于岩石 节理模型，还可用于其他硬性结构面受力分析。 参 考 文 献 [1] DONG J J, PAN Y W. A hierarchical model of rough rock joints based on micromechanics[J]. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts, 1996, 33(2): 111－123. [2] RUTQVIST J, TSANG C F, TSANG Y. Analysis of stress and moisture induced changes in fractured rock permeability at the Yucca mountain drift scale test[C]// Proceedings of International Conference on Coupled T-H-M-C Processes in Geosystems: GeoProc 2003. Stockholm: [s. n.], 2003: 147－152. [3] RAVEN K G.. Water flow in a natural rock fracture as a function of stress and sample size[J]. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts, 1985, 22(4): 251－261. [4] WITHSPOON P A, WANG J S Y, IWAI K, et al. Validity of cubic law for fluid flow in a deformable rock fracture[J]. Water Resources Research, 1980, 16: 1016 －1024. [5] BARTON N, BANDIS S, BAKHTAR K. Strength, deformation and conductivity coupling of rock joints[J]. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts, 1985, 22: 121－ 140. [6] 蒋宇静, 王刚, 李博, 等. 岩石节理剪切渗流耦合试验 及分析[J]. 岩石力学与工程学报, 2007