恒定法向刚度边界条件下岩石节理剪胀特性数值模拟

吴家红，刘 杰，宋彦彬

( 1. 绍兴文理学院 土木工程学院，浙江 绍兴 312000；2. 浙江省岩石力学与地质灾害重点实验室，浙江 绍兴 312000；3. 浙江省山体地质灾害防治协同创新中心，浙江 绍兴 312000 )

深部地下工程岩体中通常包含大量的节理面，节理在强度上较弱，在应力作用下，岩体会沿着节理面错动而发生剪胀行为，因此研究岩石节理剪胀特性对准确评估岩体稳定性具有重要的工程和实践意义[1]。在深部地下工程中，岩石节理剪切破坏受围岩的约束，其法向应力随着剪胀位移的增大而增大，而法向应力的增大又约束节理剪胀的发展，因此恒定法向荷载( Constant Normal Load，CNL )边界条件不再适用，应采用恒定法向刚度( Constant Normal Stiffness，CNS )边界条件[2-3]。

针对CNL边界条件下的粗糙节理剪胀特性，国内外学者开展了大量的研究工作，如周辉[4]等研究了CNL边界条件下不同起伏高度的大理石规则齿形节理的剪胀效应，提出了评价剪胀效应的经验公式；金磊磊[5]等通过对上下盘相同抗压强度和上下盘不同抗压强度的人工节理试样在CNL边界条件下的剪切试验，分析了法向应力、三维形貌特征及节理强度比对软-硬节理的抗剪强度和剪胀角的影响；BARTON N[6]通过峰值剪切位移处的剪胀角建立了Barton剪胀模型；赵延林[7]等利用UDEC软件建立了随机形貌岩石节理剪切模型，通过对剪胀数值曲线的分析，提出了同时考虑粗糙度和法向应力的非线性剪胀本构模型；肖维民[8]等开展了不同法向应力下充填节理的直剪试验，建立了薄层充填岩石节理峰值剪胀角随充填度呈指数变化的经验关系式。对于CNS边界条件，GU X F[9]等发现结构面磨损与结构面粗糙度和法向边界条件有关；崔国建[2]等开展了CNL和CNS边界条件下的剪切试验，分析了边界条件和粗糙度对结构面峰值剪胀位移的影响；韩观胜[10]开展了不同法向刚度下的剪切试验，研究了法向刚度对节理剪胀位移的影响，研究发现法向刚度越大，其剪胀位移越小；尹乾[1]等研究了CNS边界条件下初始法向应力和粗糙度系数对岩石节理剪胀位移的影响，研究表明剪胀位移随着粗糙度的增大而增大，随着初始法向应力的增大而减小。上述研究主要以室内试验为主，表明CNS边界条件下，岩石节理剪胀位移具有强烈的非线性特征。在数值模拟方面，SHRIVASTAVA A K[11]等利用UDEC软件开展了CNL和CNS边界条件下的剪切试验，发现数值模拟结果与试验结果非常吻合；WANG X[12]等通过PFC软件研究了加载方向及法向应力对完整及节理煤岩剪切各向异性和剪胀变形的影响；PACKULAK T[13]利用UDEC软件研究了CNL及CNS边界条件下剪胀对岩石节理屈服后剪切特性的影响。上述的数值模拟手段主要采用离散元数值模拟方法，为笔者顺利开展CNS边界条件下岩石节理剪胀特性数值模拟提供了重要参考。

迄今为止，关于岩石节理CNS边界条件下剪胀特性的研究成果较少，且对剪胀本构模型的研究鲜有报道。此外，数值模拟可考虑室内试验无法考虑的因素，在时间、费用、复杂程度和重复度上具有显著优势[7]。为此，本文利用UDEC软件，开展CNS边界条件下随机形貌岩石节理剪胀效应的数值模拟研究，旨在建立适用于CNS边界条件下的岩石节理非线性剪胀本构模型。

1 随机形貌岩石节理剖面构建

采用Hurst指数独立分割法实现随机形貌岩石节理模拟，具体实现过程[14]为

( 1 ) 设初始剖面为一直线；

( 2 ) 用一随机点P分割该直线；

( 3 ) 在点P两侧，用0.5 mm等分距等分直线；

( 4 ) 在Y方向上，分割点两侧各等分点相对于前一等分点的偏移为P(x)，则

式中，M为振幅相关系数；B为正态分布随机变量，其均值为0，方差为1；x为等分点至P点的距离；H为Hurst指数。

采用表1中的节理剖面参数，利用MATLAB语言编写计算程序构建5种不同粗糙度的节理剖面( 图1 )，粗糙度系数JRC根据TSE R[15]等提出的JRC与均方根的经验关系式进行计算。

表1 随机形貌岩石节理剖面参数Table 1 Rock joints profile parameters with random morphology

图1 不同JRC节理剖面Fig. 1 Joints profiles with different JRCs

2 随机形貌岩石节理剪胀特性的数值试验

2.1 数值模型

为研究CNS边界条件下粗糙岩石节理的剪胀特性，建立如图2所示的CNS边界条件下岩石粗糙节理剪切计算模型，模型由弹簧块、上部岩块和下部岩块组成，节理选用图1中的5种不同粗糙度节理，采用UDEC自带的fish语言将粗糙节理导入数值模型中。节理本构模型选用连续屈服( Continuously-Yielding，CY )微段节理模型[16]，该模型与库仑摩擦线性模型和Barton-Bandis模型相比更为真实，考虑了在物理试验中观察到的非线性行为、剪切过程中节理的损伤、剪胀角随塑性剪切位移的增加而减小等情况。岩石材料参数和CY微段节理模型的细观参数见表2( 计算参数参考MASIVARS D[17]的研究成果 )。为研究5种不同的法向刚度对岩石节理剪胀的影响，弹簧块所对应的法向刚度nk取值分别为2，4，6，8，10 GPa/m，其刚度大小主要由弹簧块的体积模量和剪切模量决定，具体参数见表3。为研究5种初始法向应力对岩石节理剪胀的影响，模型垂直方向施加初始内应力分别为2，4，6，8，10 MPa。弹簧块左右边界水平位移固定，上部边界垂直位移固定，上部岩块左右边界水平位移固定，下部岩块下边界垂直位移固定，左边界水平方向施加0.001 m/s的剪切速率。

图2 CNS边界条件下岩石粗糙节理剪切计算模型Fig. 2 Shear calculation model of rock rough joint under CNS boundary condition

表2 岩石及CY 微段节理模型计算参数Table 2 Calculation parameters of rock and CY microsegment joints model

表3 弹簧块材料参数Table 3 Material parameters of spring blocks

2.2 试验方案

为减少试验次数，采用正交试验设计试验方案，试验选取JRC、初始法向应力σn0、法向刚度nk共3种不同因素，记为影响因素i( 取值1～3 )，每个因素设立5个水平，设计得到3因素5水平正交试验方案，见表4。

表4 正交试验方案Table 4 Orthogonal test scheme

2.3 试验能力验证

CNS边界条件下，法向应力增量与剪胀位移增量呈线性关系，其比值为法向刚度。为验证法向刚度在剪切过程中是否保持恒定，以初始法向应力2 MPa的试样为例，绘制了剪胀位移与法向应力的关系曲线，如图3所示。由图3可知，法向应力随剪胀位移变化，两者通过线性函数较好地拟合，其斜率与截距分别近似等于设定的法向刚度与初始法向应力，决定系数R2=0.999 5～0.999 7。以上结果表明，该模型能较好地实现研究所需的CNS边界条件。

图3 CNS边界条件下法向应力与剪胀位移关系曲线Fig. 3 Relation curves between normal stress and dilation displacement under CNS boundary conditions

2.4 试验结果

图4为CNS边界条件下岩石节理剪切位移与剪胀位移关系曲线。由图4可知，岩石节理剪切初始阶段，节理在法向应力作用下闭合压密及节理面两侧试样的弹性压缩变形而经历了剪缩阶段；随后，由于各微段节理的爬坡效应，剪胀位移迅速增加；随着剪切的发展，由于初始应力和法向刚度的抑制作用，节理面的凸起体发生剪切破坏，部分节理剪胀位移增长速率逐渐减小，剪胀位移开始趋于平缓；另有部分节理的剪胀位移缓慢下降而呈现软化特征，这是由于该部分节理承受的法向应力更大，凸起体剪切破坏更加严重。

图4 不同粗糙度节理剪切位移与剪胀位移关系曲线Fig. 4 Relation curves between shear displacement and shear dilation displacement of different roughness joints

3 岩石节理剪胀本构模型

3.1 剪胀的数学表达

图4中节理剪胀位移增长阶段的变化趋势符合自然负指数增长方程：

式中，ξ为剪胀位移；ξult为最大剪胀位移；cR为剪胀收敛速率；sδ为剪切位移；＊sδ为剪胀位移为0时的剪切位移。

该方程只考虑了剪胀位移从0到最大值的剪胀阶段，而未考虑岩石节理的初始剪缩阶段。当δs趋于时，ξ趋于0；当δs趋于无穷大时，ξ趋于ξult。由于该方程仅描述剪胀响应的快速启动阶段、平稳阶 段 和 收 敛 阶 段，其 有 效 阶 段 是 从δs=δs＊到δs=δs1(δs1为在ξ= 0.95ξult时的取值 )。参数ξult和Rc取值与边界条件和节理几何形貌有关，而位移参数δs＊取值取决于岩石节理初始剪缩阶段的终止。为便于分析，对初始法向应力σn0和法向刚度kn进行量纲归一化处理的公式为式( 3 )和( 4 )，最终拟合结果见表5。

表5 方程拟合参数及决定系数Table 5 Equation fitting parameters and determination coefficient

3.2 极差分析

将表5中的最大剪胀位移ξult和剪胀收敛速率cR作为试验指标，按照每列在同一水平下求和并计算其平均值，分别记为AIi和BIi，其中A代表最大剪胀位移，B代表剪胀收敛速率，I( 取值1～5 )表示影响因素i的5种不同水平。计算出各影响因素的平均值极差，得到最大剪胀位移和剪胀收敛速率影响因素的极差分析结果，见表6。

表6 试验结果的极差分析Table 6 Range analysis of test results

由表6可知，在CNS边界条件下，各因素对最大剪胀位移ξult影响程度从大到小依次为节理的粗糙度系数、初始法向应力、法向刚度；各因素对剪胀收敛速率 cR影响程度从大到小依次为法向刚度、初始法向应力、粗糙度系数。

为便于研究各因素对最大剪胀位移和剪胀收敛速率的影响，以各因素为水平坐标，岩石节理最大剪胀位移ξult和剪胀收敛速率 cR为纵坐标，绘制图5和6。

图5 最大剪胀位移变化趋势Fig. 5 Variation trend of maximum dilation displacement

由图5可知，随着JRC的增大，由于岩石节理最大起伏度不断增大导致节理最大剪胀位移不断增加；随着初始法向应力的增大，岩石节理表面剪切损伤增加导致节理最大剪胀位移不断减小；由于法向刚度的增大，剪切过程中法向应力增大，节理表面损伤变大，最大剪胀位移整体呈逐渐减小的趋势。由图6可知，JRC=8.914时，岩石节理剪胀收敛速率最大；随着初始法向应力和法向刚度的增大，岩石节理剪胀收敛速率不断增大。

图6 剪胀收敛速率变化趋势Fig. 6 Variation trend of dilation convergence rate

3.3 最大剪胀位移量化分析

最大剪胀位移的存在是由于岩石节理表面凸起体尺寸有限，应力作用下产生不同程度的剪切磨蚀，从而削弱了岩石节理的剪胀效应。理想状态下，对于重复规则锯齿节理，在法向应力为0的条件下，预期剪胀位移从0增至2倍的锯齿高度，并按照锯齿的形状周期振荡，这是由于重复锯齿上下表面在不同的剪切位置能完美匹配。在法向应力为非0作用下，锯齿尖端会受到损伤，损伤大小取决于材料强度和法向应力大小，从而削弱岩石节理剪胀。由上述分析可假设最大剪胀位移取决于最大起伏高度ah( 节理剖面线长度范围内最高点与最低点之间的垂直距离 )和应力条件( 由法向刚度和初始法向应力决定 )。理论上，在法向应力为0时，最大剪胀位移ξult等于最大起伏高度ah，当法向应力增加到某一程度时，粗糙节理面被完全破坏，最大剪胀位移接近于0，剪胀效应消失。

以ξultah为应变量，初始法向应力和法向刚度为自变量，整理计算表5中的数据得到表7。

表7 kN ， Nσ ，ah 与ξult ah参数汇总Table 7 Parameter summary of kN ， Nσ ，ah andξult ah

续 表

由表7可知，当法向刚度相同时，随着初始法向应力的增加，岩石节理的最大剪胀位移与最大起伏高度比值ξultah逐渐减小，且ξultah与Nσ间的关系满足式( 5 )。将表7中的数据按照式( 5 )进行曲线拟合，拟合结果如图7所示，拟合所得参数见表8。

表8 不同法向刚度下最大剪胀位移拟合参数Table 8 Fitting parameters of maximum dilation displacement under different normal stiffness

式中，A，B为拟合参数。

分析表7中的数据可知，当初始法向应力相同时，随着法向刚度的增大，岩石节理最大剪胀位移与最大起伏高度比值ξultha逐渐减小。将表7中的数据按照式( 6 )进行曲线拟合，拟合结果如图8所示，拟合所得参数见表9。

表9 不同初始法向应力下最大剪胀位移拟合参数Table 9 Fitting parameters of maximum dilation displacement under different initial normal stresses

图8 不同初始法向应力 kN -ξult ha 散点图及其拟合曲线Fig. 8 kN -ξult ha scatter plots and their fitting curves under different initial normal stresses

式中，C，D为拟合参数。

结合式( 5 )和( 6 )，同时考虑初始法向应力和法向刚度对岩石节理剪胀的影响，ξultah可表示为

式中，a，b，c为最大剪胀位移拟合参数。

采用式( 7 )对表7中的数据进行拟合，拟合所得参数见表10，图9为拟合得到的最大剪胀位移曲面。

表10 最大剪胀位移拟合参数Table 10 Fitting parameters of maximum dilation displacement

图9 最大剪胀位移拟合曲面Fig. 9 Fitting surface of maximum dilation displacement

3.4 剪胀收敛速率量化分析

剪胀位移增长速率与凸起体倾角有关，倾角越大，增长速率越快。而剪胀收敛速率反映了剪胀位移趋于最大值的快慢程度。图10为剪胀收敛速率大小相同但剪胀位移增长速率不同的2条剪胀曲线，剪胀位移越大的曲线，其剪胀位移增长速率越大；但在相同的剪切位移条件下，所对应的剪胀位移与各自最大剪胀位移的百分比相同。因此，只采用凸起体倾角描述 cR欠妥，cR的大小不仅与凸起体倾角有关，还与最大剪胀位移有关。图11为一个理想的规则锯齿节理，将式( 2 )两边进行求导得到

图10 cR 的含义Fig. 10 Meaning ofcR

图11 规则锯齿节理Fig. 11 Regular sawtooth joint

当该锯齿在零法向应力下，刚发生剪胀时，则

结合式( 8 )与( 9 )，可以得到

为了讨论凸起体倾角与最大剪胀位移对收敛速率Rc的影响，结合式( 10 )提出一个新变量taniave表征Rc，其中iave为凸起体平均倾角，计算公式为

式中，N为节理离散点数目；Δx为采样间距；iy为第i个点的高度坐标。

当ξult增加时，Rc逐渐减小；当ξult趋于0时，相当于剪胀位移从0到0，Rc趋于无穷大，意味着未发生剪胀或剪胀位移瞬间达到最大。当iave增加时，岩石节理越粗糙，Rc也随之增加，剪胀速率越大，达到最大剪胀位移所需的剪切位移减小；当iave趋于0时，ξult也趋于0，iave表征的Rc无意义。表11为不同法向刚度条件下Rc与iave的参数汇总，将表中的taniave和Rc进行曲线拟合，如图12所示。

表11 不同法向刚度下 Rc 与 iave参数汇总Table 11 Parameters summary of Rc andtan iave under different normal stiffness

表11 不同法向刚度下 Rc 与 iave参数汇总Table 11 Parameters summary of Rc andtan iave under different normal stiffness

k iave /( ° ) 1 N ξ - i cR /mm-1 tan ult ave 0.25 5.671 0.249 0.374 0.25 8.379 0.562 0.536 0.25 10.973 0.511 0.529 0.25 13.665 0.334 0.397 0.25 16.851 0.286 0.385 0.50 5.671 0.366 0.473 0.50 8.379 0.367 0.469 0.50 10.973 0.327 0.419 0.50 13.665 0.473 0.519 0.50 16.851 0.440 0.487 0.75 5.671 0.487 0.527 0.75 8.379 0.592 0.573 0.75 10.973 0.493 0.532 0.75 13.665 0.330 0.471 0.75 16.851 0.273 0.417 1.00 5.671 0.347 0.508 1.00 8.379 0.352 0.524 1.00 10.973 0.634 0.628 1.00 13.665 0.487 0.569 1.00 16.851 0.448 0.543 1.25 5.671 0.499 0.574 1.25 8.379 0.582 0.623 1.25 10.973 0.476 0.561 1.25 13.665 0.320 0.485 1.25 16.851 0.526 0.605

图12 ia v e - Rc 散点图及拟合曲线Fig. 12 ia v e- Rc scatter plot and its fitting curve

式中，k，d，e，f为岩石节理剪胀收敛速率拟合参数。

用式( 13 )对表11中的数据进行曲面拟合，拟合参数见表12，图13为拟合得到的岩石节理剪胀收敛速率曲面。

表12 岩石节理剪胀收敛速率拟合参数Table 12 Fitting parameters of convergence rate of shear dilatancy of rock joints

图13 岩石节理剪胀收敛速率拟合曲面Fig. 13 Fitting surface of dilation convergence rate of rock joints

3.5 岩石节理剪胀模型

综上可知，CNS条件下岩石节理剪胀的本构模型可表示为

其中，

该剪胀模型的量纲是一致的。图14为采用式( 14 )得到的理论曲线与数值试验曲线的对比，可以看出剪胀理论曲线与数值曲线能较好地吻合。

图14 剪胀理论曲线与数值试验曲线对比Fig. 14 Comparison between theoretical dilation curves and numerical test curves

此外，剪胀角可以表示为

4 结 论

( 1 ) CNS边界条件下岩石节理剪胀段的剪胀位移呈负指数增长，提出了包含最大剪胀位移、剪胀收敛速度和发生剪胀位移为0时的剪切位移的负指数增长剪胀模型。

( 2 ) CNS边界条件下岩石节理的最大剪胀位移随着岩石节理最大起伏高度的增大而增大；随着初始法向应力和法向刚度的增大而减小；岩石节理剪胀收敛速率随着初始法向应力和法向刚度的增大而增大。

( 3 ) 提出的CNS边界条件下岩石节理的剪胀模型同时考虑了节理的最大起伏高度、凸起体倾角、法向刚度和初始法向应力的影响。通过剪胀理论曲线和数值试验曲线的对比，表明该模型具有较好的普适性。

参考文献( References )：

[1] 尹乾，靖洪文，孟波，等. 恒定法向刚度条件下三维粗糙裂隙面剪切力学特性[J]. 岩石力学与工程学报，2020，39( 11 )：2213-2225.YIN Qian，JING Hongwen，MENG Bo，et al. Shear mechanical properties of 3D rough rock fracture surfaces under constant normal stiffness conditions[J]. Chinese Journal of Rock Mechanics and Engineering，2020，39( 11 )：2213-2225.

[2] 崔国建，张传庆，韩华超，等. CNL及CNS条件下结构面剪切特性试验研究[J]. 岩石力学与工程学报，2019，38( S2 )：3384-3392.CUI Guojian，ZHANG Chuanqing，HAN Huachao，et al. Experiment study on shear behavior of artificial joint under CNL and CNS boundary conditions[J]. Chinese Journal of Rock Mechanics and Engineering，2019，38( S2 )：3384-3392.

[3] BAHAADDINI M. Effect of boundary condition on the shear behaviour of rock joints in the direct shear test[J]. Rock Mechanics and Rock Engineering，2017，50( 5 )：1141-1155.

[4] 周辉，程广坦，朱勇，等. 大理岩规则齿形结构面剪切特性试验研究[J]. 岩土力学，2019，40( 3 )：852-860.ZHOU Hui，CHENG Guangtan，ZHU Yong，et al. Experimental study of shear deformation characteristics of marble dentate joints[J]. Rock and Soil Mechanics，2019，40( 3 )：852-860.

根据大坝坝坡稳定分析计算结果表6，死水位及正常水位稳定渗流工况上下游坡抗滑稳定安全系数均大于1，说明坝坡尚未失稳，与工程运行实际相符；正常水位降落工况上游坡抗滑稳定安全系数小于1，说明该工况坝坡失稳，同时滑弧位置与工程实际病害相符。

[5] 金磊磊，魏玉峰. 基于三维形貌和剪胀效应的软-硬节理抗剪强度模型[J]. 工程力学，2020，37( 12 )：180-190.JIN Leilei，WEI Yufeng. Calculation model for the shear strength of soft-hard joints based on three-dimensional morphology and dilatancy effect[J]. Engineering Mechanics，2020，37( 12 )：180-190.

[6] BARTON N. Review of a new shear strength criterion for rock joints[J].Engineering Geology，1973，7( 4 )：287-332.

[7] 赵延林，万文，王卫军，等. 随机形貌岩石节理剪切数值模拟和非线性剪胀模型[J]. 岩石力学与工程学报，2013，32( 8 )：1666-1676.ZHAO Yanlin，WAN Wen，WANG Weijun，et al. Shear numerical simulation of random morphology rock joint and nonlinear shear dilatancy model[J]. Chinese Journal of Rock Mechanics and Engineering，2013，32( 8 )：1666-1676.

[8] 肖维民，余欢，朱占元，等. 薄层充填岩石节理剪胀特性试验研究[J]. 岩土工程学报，2020，42( 8 )：1499-1508.XIAO Weimin，YU Huan，ZHU Zhanyuan，et al. Experimental study on shear dilatancy characteristics of thinly-infilled rock joints[J].Chinese Journal of Geotechnical Engineering，2020，42( 8 )：1499-1508.

[9] GU X F，SEIDEL J P，HABERFIELD C M. Direct shear test of sandstone-concrete joints[J]. International Journal of Geomechanics，2003，3( 1 )：21-33.

[10] 韩观胜. 恒法向刚度条件下岩石多节理结构剪切力学特性研究[D]. 徐州：中国矿业大学，2019.HAN Guansheng. Study on shear mechanical behaviours of multi-joint structures of rock under constant normal stiffness conditions[D].Xuzhou：China University of Mining and Technology，2019.

[11] SHRIVASTAVA A K，RAO K S，RATHOD G W. Numerical simulation of direct shear test on rock joints[C]. Geocongress，2012.

[12] WANG X，KANG H，GAO F. Numerical investigation on the shear behavior of jointed coal mass[J]. Computers and Geotechnics，2019，106：274-285.

[13] PACKULAK T. Laboratory investigation of shear behaviour in rock joints under varying boundary conditions[D]. Kingston，Ontario，Canada：Queen′s University，2018.

[14] 许宏发，李艳茹，刘新宇，等. 节理面分形模拟及JRC与分维的关系[J]. 岩石力学与工程学报，2002，21( 11 )：1663-1666.XU Hongfa，LI Yanru，LIU Xinyu，et al. Fractal simulation of joint surface and the relationship between JRC and fractal dimension[J].Chinese Journal of Rock Mechanics and Engineering，2002，21( 11 )：1663-1666.

[15] 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

[16] 赵延林，万文，王卫军，等. 随机形貌岩石节理剪切-渗流数值模拟和剪胀-渗流模型[J]. 煤炭学报，2013，38( 12 )：2133-2139.ZHAO Yanlin，WAN Wen，WANG Weijun，et al. Shear-seepage numerical simulation and shear dilatancy seepage model of rock joint with random morphology[J]. Journal of China Coal Society，2013，38( 12 )：2133-2139.

[17] MASIVARS D. Water inflow into excavations in fractured rock—A three-dimensional hydro-mechanical numerical study[J]. International Journal of Rock Mechanics and Mining Sciences，2006，43( 5 )：705-725.