柱状节理岩体原位变形试验力学浅析与模拟

狄圣杰，徐卫亚，王 伟，吴关叶

（1.河海大学 岩土力学与堤坝工程教育部重点实验室，南京 210098；2. 河海大学 岩土工程科学研究所，南京 210098；3.中水顾问集团华东勘测设计研究院，杭州 310014）

1 引 言

白鹤滩水电站是我国拟建的一座千万千瓦级的大型水电工程，也是西电东送的骨干电源点，坝址区玄武岩柱状节理较为发育，如图1所示，并且越来越多的水利工程（如金沙江下游的溪洛渡、乌东德等水电站）均揭露出大量柱状节理玄武岩体，从传统工程地质质量分类来看，属于完整性较差的岩体，其变形性能是否满足高拱坝坝基严格的变形要求是工程界广泛关注的问题。

由于在柱状节理玄武岩地区进行工程建设的需要，柱状节理玄武岩工程特性的研究也得到了开展。石安池等[1]在研究白鹤滩柱状节理玄武岩特点的基础上，对影响柱状节理玄武岩岩体变形特性的因素和变形性能各向异性的原因作了工程介绍和解释说明；徐卫亚[2－3]、孟国涛[4]、郑文棠等[5]对白鹤滩柱状节理岩体本构关系、参数取值、屈服准则及数值方法等各方面进行了详细的研究，并在坝基工程进行了应用[6]；Sharma等[7]建议采用柔性承压板试验来反映节理岩体的变形模量；Unal[8]的研究也表明，在现场测试岩体参数时，必须考虑松动圈对测试结果的影响，强调了中心孔变形试验的重要性；张宜虎[9]、杨春和等[10]对中心孔变形试验资料和循环加卸载岩体变形特性进行了试验解释及应用研究；Palmstrom[11]在分析洞室的承压板试验曲线时进行了论述，同时认为这种初始加载阶段曲线不应当计入整体变形模量的计算中，且对节理岩体的原位测试机变形模量的估算做了较多的理论研究工作；Brady等[12]对大尺寸节理玄武岩体的加卸载规律性进行了详细分析，得出典型性的3个阶段力学行为。本文利用理论分析和数值模拟针对原位变形试验资料对节理刚度取值、单节理岩体加、卸载变形规律、多组贯通节理各向异性特征及柱状节理随机模型模拟等问题进行探讨。

图1 白鹤滩坝址区可视化模型与柱状节理岩体Fig.1 Baihetan visual model dam site and columnar jointed rock

2 基于原位试验资料的节理刚度反分析

结构面参数通常应用室内试验或经验方法进行估算，除此之外还可运用现场中心孔试验法来进行估算[13]。节理刚度可通过节理闭合试验及剪切试验确定，但试样取样、搬运、制备要求严格、繁琐，同时不可避免地存在尺寸效应问题和试件代表性问题。而利用现场中心孔试验反分析节理参数更具说服力。

柔性中心孔原位试验采用4枕柔性承压板中心孔法变形法，试验装置如 2(a)中左图所示，试验时在平硐水平向、铅直向试验面中心部位钻孔，孔径φ75 mm，孔深为6 m，在5 m孔深范围布置多点位移计。采用逐级循环加载并纪录不同深度不同应力水平下的位移W，绘制得到多点位移计的σ-W关系曲线，如2(b)图所示。

图2 白鹤滩柱状节理岩体柔性中心孔法变形试验Fig.2 Center hole deformation test of flexible of columnar jointed rock in Baihetan

根据柱状节理玄武岩中心孔试验成果，9个试验点中分为3个水平和6个竖直试验点，位于平硐PD133和平硐PD36中。平硐PD36中的试验点位于Ⅲ1类地层柱状节理岩层的弱风化带中，为3个水平试验点EC36Z-101、EC36Z-102、EC36Z-103和3个竖直试验点 EC36Z-201、EC36Z-202、EC36Z-203，其中EC36Z-101试验点报废。平硐PD133中的试验点位于为的微风化带中，为3个竖直试验点EC133Z-201、EC133Z-202、EC133Z-203。试验点成果如表1所示。

表1 柔性中心孔法变形试验的试验点资料Table 1 Center hole deformation test of flexible test point data

分析白鹤滩柱状节理玄武岩的岩体强度变化规律可知：同一岩层内，弱风化带的变形模量低于微风化带的变形模量；同一风化层的相同深度范围内，受倾斜柱状节理面影响，竖直方向的变形模量低于水平方向的变形模量；同一风化层的不同深度范围内，受卸荷影响，浅部岩体的变形模量低于深部岩体的变形模量。符合以上规律的柔性中心孔试验点为：EC36Z-103、EC36Z-201、EC36Z-202、EC36Z-203、EC133Z-201、EC133Z-202，考虑到松弛圈的厚度约为70～120 cm，概化原位试验地层为两层：试验点以孔深86 cm以内为松弛柱状节理岩体，孔深86 cm以下为未松弛柱状节理岩体。

由以上试验点资料可以估算不同风化程度、松弛程度柱状节理面刚度系数的有效个数为：

（1）弱风化层松弛岩体，可估计节理面参数，水平方向1组，竖直方向3组；

（2）弱风化层未松弛岩体，可估计节理面参数，竖直方向1组；

（3）微风化层松弛岩体，可估计节理面参数，竖直方向2组；

（4）微风化层可估计节理面参数，竖直方向 2组。

由于结构面的变形是非线性且不可恢复的，而岩石的弹性变形是可恢复的，因此，利用某级循环荷载下的应力差和相应的变形差可估算结构面的刚度[3]。如图2所示，S1为柱状节理面的平均间距，可以用柱体直径代替；S2为横节理的平均间距，可以用柱体乘以纵横比得出。利用中心孔试验成果，对于水平孔和铅直孔分别有如下关系式：

式中：Δhi为不同的深度测点距离差；、′为不同深度段的平均应力；Δ Wi+1-ΔW 为每段岩体整体变形量减去岩石受压而产生的压缩量，即结构面的变形量；、根据Boussinesq弹性理论近似求得，承压面积中心以下深为zi处的铅直应力为

式中：q为受荷表面的均布压力；a为中心孔试验板半径；zi为深度。

假定1、2组结构面有相同的刚度，通过每两级加载间各测点的应力差和位移差，根据式（1）～（4）可估算出单一节理面的刚度系数，反分析成果见表2。

表2 柔性中心孔试验点原位试验资料反分析成果Table 2 Flexible center hole test points, the back analysis of in situ test results

3 节理岩体加载、卸载力学变形特性

初始加载阶段，表现为对应斜率较低，但斜率增长速率很快的上凹曲线。对于这种上凹曲线的物理机制，试验点总是会受到开挖爆破的扰动而产生卸荷以致松动张开，因此，在这个阶段节理均表现为压密闭合，不能代表原岩无扰动下的状态。岩体压密阶段斜率为斜率Ⅰ，之后进入加载阶段，斜率为斜率Ⅱ，为较平直的一条准直线，表现了岩块的压缩变形及节理的压缩和剪切变形共同作用的结果；初始卸载阶段表现为斜率为斜率Ⅲ的曲线，斜率相对较高，后继卸载阶段对应较缓的曲线，斜率为斜率Ⅳ[4,14]。

取含一条封闭节理的岩体进行力学分析，节理闭合，建立概念模型如图 3(a)所示，岩体尺寸长为B、宽为 W、高为 H，将岩块视为各向同性的弹性体，岩块的压剪刚度分别是km、kn，弹性模量为E，岩体等效刚度为k，节理面法向与加载方向成倾角β，节理的法向刚度和切向刚度分别为Kjn和Kjs，摩擦角为φ，封闭节理长l，沿其迹线总长为L，模型上施加作用力F。单节理岩体加卸载曲线如图3(b)所示，PD36典型承压板试验曲线如图3(c)所示，由于节理的相互作用，初始卸载曲线不为直线回弹，表现为斜率渐小的曲线回弹。

假设岩石是理想弹性体，节理刚度值不随应力而改变，节理岩体中各点的应力状态与完整岩石应力状态相同，这对于较致密、刚度较大的玄武岩体是合适的。先研究贯通节理岩体的变形，沿加载方向岩体的轴向变形量为两部分产生的变形之和u=urn+ujs+ujn，刚度满足以下关系式：

式中：urn、ujs、ujn分别为岩块、节理切向和节理法向产生的沿加载方向的变形；kjn和kjs为由节理法向刚度Kjn和切向刚度Kjs对沿加载方向刚度的贡献，kj为作用力与变形的关系，单位为N/m，Kj为应力与变形的关系，单位为N/m3。节理和岩块产生的沿加载方向的变形为

图3 单节理岩体概念模型及典型荷载变形曲线示意图Fig.3 Conceptual model of a single rock mass and the typical load deformation curve diagram

对应的刚度分别为

对于未贯通节理，不考虑岩桥与节理接触部位的尖端应力集中等因素，对上述贯通节理的式子进行修正，岩桥与节理的联合刚度采用加权平均等效考虑，节理的贯通率定义为0 ＜η=l/ L ＜1，如图4所示，则：

式中：Krn、Krs、Kjn、Kjs分别为岩石和节理的压剪刚度，分别由其压剪试验确定；和为考虑节理与岩桥单元的等效刚度。对于未贯通节理单元，需用式（12）和（13）替代式（6）～（11）中的Kjn和Kjs。

故加载阶段满足以下关系：

式中：T为节理端部产生的拉力；N为力学等效连续节理上法向力；f为沿节理面的摩擦阻力；us′为端部产生的沿节理面的阻变形；u′为岩体变形总量。

图4 非贯通节理介质Fig.4 Non-jointed media

将式（16）代入式（17）得加载阶段曲线的斜率为

当竖向荷载按一定速率被移除时，则节理岩体将产生一个反方向的回弹变形，节理面上的摩擦阻力作用相反，当摩擦阻力足够大，能阻止节理面向上滑移时，反弹值只有岩体的变形值而无节理的滑移值，则 斜率I II=k/ BW 。同理卸载满足：

如果取模型为单位长度尺寸的正方体，则3者斜率蜕化为

柱状节理倾角一般在 70°以上，则斜率Ⅳ＜斜率III ，斜率II＜斜率III 。一个完整的加载循环包括加载阶段岩体变形，节理面存在非线性力学行为；初始卸载变形阶段只存在岩体变形回弹；后继卸载阶段岩体变形回弹节理面存在非线性力学行为，如图3(b)所示，表明概念模型是能够反映承压板试验成果的。

对于多节理岩体的循环加载试验，表现为初始卸载阶段为斜率递减的卸载曲线，各级卸载曲线不沿加载曲线路径。在初始卸载阶段才能抑制节理的非弹性变形，可代表节理岩体真实的弹性力学行为，在计算节理岩体等效弹性模量时，应采用承压板原位试验中的初始卸载阶段曲线。

一般荷载低于岩石的屈服强度，认为其滞后完全是由于节理的非线性行为导致，如图 3(c)所示，原位典型上凹型试验曲线具有层理裂隙等结构面的非均质岩体特征，反映了随着压应力的增加结构面逐渐被压密，加载模量随之增大的趋势，另一方面岩体各条卸载曲线在卸载到应力水平较低时，曲线斜率明显减小，可以说明卸载后裂隙发生了松弛，发生张开和扩展。

4 基于变形解析式的岩体各向异性特征

对于一定尺度下的岩体，在某应力作用下，由式（6）～（8）可得多组节理情况下加载方向产生的总变形满足以下关系式：

含一条L长度节理的岩体变形量为

对于随机分布贯通节理的岩体变形量普遍的形式为

式中：m为节理的条数，对于有厚度节理单元，用节理单元的厚度代表其间距；βi为第i条节理法向与加载方向的夹角；H为加载方向的岩体尺寸；W为垂直加载方向岩体尺寸； Licosβi/W反映对加载方向的等效变形量的贡献，即有效长度比，不同节理贯长度的影响其实是由于节理面接触面积的影响。求得岩体平行于加载方向的等效变形模量为

由于引入了加载方向与节理法向的夹角，针对不同的方向加载，可以反映变形模量的变化规律和岩体的各向异性特征。当节理为两组等密正交节理时，如图5所示，正方形岩块尺寸10 m，节理为正交，间距为0.5 m，岩块弹性模量为50 GPa，节理法向刚度为50 GPa/m，k为切向刚度与法向刚度比值，根据式（26）～（27）得到不同节理倾角时的等效弹性模量极坐标。

图6(a)为两组正交节理法向刚度相等时，岩体弹性模量随节理切向刚度与法向刚度比值k的变化曲线。图6(b)为第1组节理法向刚度为其正交的第2组节理2倍时的曲线，图中 k=Ks/Kn，曲线代表弹性模量随节理切向刚度与法向刚度比值从0.1变化至 5.0时对应的曲线。各向异性曲线特征与 Ki-Bok[15]揭示的规律一致。

随着 Ks/ Kn比值减小而增强，并在β接近45°时达到极值。当两组节理刚度特性不同时，即水平向节理刚度为竖直向节理的2倍时，呈扁平化特征；由于两组节理间距相同，节理刚度相同，则表现为各向同性特征，对应为球形。且 Ks/Kn=1时，虽然两组节理间距相同，但节理刚度不同，仍表现为各向异性特征。由于式（26）～（27）是通式，对于任意倾角贯通节理同样可以求解，如图6(b)所示。对于未贯通节理岩体，根据等效刚度来确定，将式（12）～（13）中等效刚度代入通式（26）～（27）中即可得到。坝址区柱状节理玄武岩体的柱体与其主轴夹角约为15°，柱元直径S1约为0.2 m，柱体长度S3约为1 m，岩块弹性模量取为65.1 GPa，节理法向刚度与切向刚度分别取 284.48 GPa/m和99.31 GPa/m。概化模型如图7所示，通过式（27）计算得到主轴坐标系下弹性模量E3和E1分别为53 GPa和30.4 GPa，主轴偏转15°后xy坐标系下Ey与Ex分别为31.6 GPa和27.1 GPa。

图5 正交节理岩体与随机节理岩体Fig.5 Orthogonal jointed rock with random jointed rock

图6 二维正交节理岩体各向异性特征曲线Fig.6 Two-dimensional orthogonal anisotropy curves of jointed rock

图7 柱状节理岩体柱体偏转示意图Fig.7 Deflection diagram of columnar jointed rock column

5 承压板变形试验数值模拟

Hart等[16]指出，BWIP规则六边形柱状节理岩体原位变形试验揭示出的滞后效应、应变分布的不均匀性及刚度的围压效应都可以归结为节理网络的非线性行为造成。Hart等采用数值试验分析了规则六边形柱状节理岩体原位变形试验中揭示的滞后行为，认为柱状节理的转动和滑移是柱状节理岩体非线性行为的根本原因，如图8所示。

图8 BWIP规则六边形柱状节理数值模拟结果[16]Fig.8 Numerical simulation results of BWIP rules hexagonal columnar joints[16]

此处基于不规则柱元形状的离散元数值构建方法，采用Voronoi构建语法，在6 m×3 m的区域内，采用平均边长为0.2 m，迭代次数小于100来控制柱面形状以四边形和五边形为主，生成与现场勘查较一致的柱状节理数值模型，如图9所示。

承压板宽度为2 m，参数选取如表3所示，岩体节理刚度初始加载阶段采用松弛节理刚度Kn0和Ks0，为初始加载压密阶段的参数，压密后采用未松弛节理刚度Kn和Ks。分4级加、卸载，每级以2 MPa为单位加、卸载，反映的是白鹤滩柱状节理岩体垂直于柱体方向应力-变形关系。建立多组随机柱状节理模型，得到多组加卸载曲线，可得到多组数据的平均变形模量与弹性模量值。

竖向变形与应力云图如图10所示，模拟的加、卸载曲线如图11所示。

图9 原位横截面素描图与生成的Voronoi模型Fig.9 Sketches of the cross-section of in-situ generation diagram and Voronoi model

表3 不规则柱状节理参数Table 3 Parameters of irregular columnar joints

图10 竖向变形与应力云图Fig.10 Vertical deformations and stress nephograms

图11 原位加卸载模拟曲线Fig.11 Simulation curves of in-situ loading and unloading

卸载阶段结束后，节理岩体存在不可恢复的永久变形，与原位试验初始卸载曲线吻合较好，单节理力学概念模型的初始卸载曲线是直线，复杂节理情况下为多折线型曲线，斜率逐渐减小，表明模拟的柱状节理加、卸载力学行为及滞后特性效果较好。竖向应力从表面向深部随机延伸，这有别于各向同性材料,显现出各向异性特性。在应力较大时可能形成节理滑移的应力临界条件，可理解为块体在应力作用下局部排列方式的改变，由平动和偏转导致。根据资料统计[4]，Ⅲ1类地层微风化带左岸柱状节理岩体垂直于柱体方向的变形模量测试值为 18.86～36.25 GPa，弹性模量测试值为28.83～41.20 GPa，计算变形模量为20.69 GPa，弹性模量为33.3 GPa，其值在试验值范围内。

6 结 论

（1）由于柱状节理刚度值较难确定，室内试验取样难度也较大，利用原位测试成果反算其值是一种比较合适的方法；

（2）对节理岩体加、卸载力学性能分析，从概念模型研究入手，解释其非线性力学行为，同时基于概念模型得到的变形公式是一通式，可以用于等效变形参数的计算和节理岩体各向异性的分析；

（3）应用由现场测试结果反算的节理刚度值和离散元方法进行了柱状节理的模拟计算，规律性与现场试验结果一致，所得弹性模量和变形模量均在实测值范围内，文中所采用方法可为涉及柱状节理的工程计算及参数取值提供参考。

（4）值得指出的是，基于解析式的节理刚度反算及岩体各向异性分析，不能考虑围压和节理交错影响，同时岩体在构造应力和二次应力场的作用下也会引起局部应力集中现象，并且个别柱状节理会在加载过程中产生屈服强化，这些都会影响到岩体应变不均和滞后效应，在实际工程中需要考虑到。

[1] 石安池, 唐鸣发, 周其健. 金沙江白鹤滩水电站柱状节理玄武岩岩体变形特性研究[J]. 岩石力学与工程学报,2008, 27(10): 2079－2086.SHI An-chi, TANG Ming-fa, ZHOU Qi-jian. Research of deformation characteristics of columnar jointed basalt at Baihetan hydropower station on Jinsha River[J]. Chinese Journal of Rock Mechanics and Engineering, 2008,27(10): 2079－2086.

[2] 徐卫亚, 郑文棠, 宁宇, 等. 柱状节理坝基岩体三维各向异性数值计算[J]. 岩土力学, 2010, 31(3): 949－955.XU Wei-ya, ZHENG Wen-tang, NING Yu, et al. 3D anisotropic numerical analysis of rock mass with columnar joints for dam foundation[J]. Rock and Soil Mechanics, 2010, 31(3): 949－955.

[3] 徐卫亚, 狄圣杰, 郑文棠, 等. 白鹤滩水电站上坝线左拱座楔形体安全性能分析[J]. 岩石力学与工程学报,2011, 30(5): 910－916.XU Wei-ya, DI Sheng-jie, ZHENG Wen-tang, et al.Safety performance analysis of rock wedges under left skewback of upstream dam line in Baihetan Hydropower Station[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(5): 910－916.

[4] 孟国涛. 柱状节理岩体各向异性力学分析及其工程应用[D]. 南京: 河海大学, 2007.

[5] 郑文棠. 不规则柱状节理岩石力学及在高边坡坝基岩石工程中的应用[D]. 南京: 河海大学, 2008.

[6] DI Sheng-jie, XU Wei-ya, NING Yu, et al.Macro-mechanical properties of columnar jointed basaltic rock masses[J]. Journal of Central South University of Technology, 2011, 18(6): 2143-2149.

[7] SHARMA V M, SINGH R B, CHAUDHARY R K.Comparison of different techniques and interpretation of the deformation modulus of rock mass[J]. Proc. Indian Geotechnical Conference IGC. Visakhapatnam, 1989,I: 439.

[8] UNAL E. Determination of in situ deformation modulus:New approaches for plate-loading tests[J]. International Journal of Rock Mechanics and Mining Sciences, 1997,34(6): 897－915.

[9] 张宜虎, 石安池, 周火明, 等. 中心孔变形试验资料的解释与应用[J]. 岩石力学与工程学报, 2008, 27(3): 489－595.ZHANG Yi-hu, SHI An-chi, ZHOU Huo-ming, et al.Explanation and application of data obtained in deformation test of central hole[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(3): 489－595.

[10] 杨春和, 马洪岭, 刘建锋, 等. 循环加、卸载下盐岩变形特性试验研究[J]. 岩土力学, 2009, 30(12): 3562－3568.YANG Chun-he, MA Hong-ling, LIU Jian-feng, et al.Study of deformation of rock salt under cycling loading and unloading[J]. Rock and Soil Mechanics, 2009,30(12): 3562－3568.

[11] PALMSTROM A, SINGH R. The deformation modulus of rock masses—comparisons between in situ tests and indirect estimates[J]. Tunnelling and Underground Space Technology, 2001, 16(2):115－131.

[12] BRADY B H G., CRAMER M L, HART R D.Preliminary analysis of a loading test on a large basalt block[J]. Int. J. Rock Mech. Min. Sci. & Geomech.Abstr., 22: 345－348.

[13] 刘佑荣, 唐辉明. 岩体力学[M]. 武汉: 中国地质大学出版社, 1999.

[14] Itasca Consulting Group Inc. Verification problems and example applications, version 4.0[M]. Minneapolis: [s. n.],2003.

[15] KI-BOK M. Fractured rock masses as equivalent continua——A numerical study[D]. Stockholm: KTH University,2004.

[16] HART R D, CUNDALL P A, CRAMER M L. Analysis of a loading test on a large basalt block[J]. Research &Engineering Applications in Rock Masses, Proceedings of the 26th US Symposium on Rock Mechanics, 1985,2: 759－768.