丹巴石英云母片岩各向异性特征数值模拟研究

周 勇

（中国电建集团华东勘测设计研究院有限公司，浙江杭州，311122）

岩体各向异性是指岩体的性质随方向而异的特性。岩体的各向异性由岩石层理、片理、夹层或定向裂隙系统存在所致。多年来，岩体各向异性问题为力学界和岩土工程界所关注，如何准确描述岩体各向异性特性一直是岩土力学研究的难题。

许多研究者针对岩体各向异性进行了各类理论和试验的深入研究工作。针对层状岩体强度各向异性特征的试验研究规律，Jaeger[1-2]基于摩尔库伦准则发展了单弱面强度理论，为研究固有不连续面对岩体强度的影响提供了有益的起点，Tien Y M[3]、Duveau G[4]等人对此进行了有益的修正。同时，国内外许多学者对不同的层状岩体或节理岩体进行了大量试验研究[5-10]。其中Hoek和Brown[8]通过对片页岩、板岩和砂岩变角度的单轴、三轴室内加载试验，获得各类岩体在不同弱面加载角、不同围岩下的强度各向异性特性试验曲线。刘胜利等[9]开展了武当群绿泥石片岩长方体试样的单轴压缩和圆盘试样的间接拉伸试验，探讨了试样各向异性的力学特性和在不同受力状态下的变形破裂特性，揭示了不同变形破裂的力学机制。张晓平等[10]通过丹巴二云英片岩单轴压缩试验研究，得出片状单轴压缩条件下的裂纹扩展过程存在显著的各向异性，并研究了相应的与片岩裂纹扩展相关的应力门槛值。张春生等[11]采用细观颗粒离散元对石英云母片岩各向异性特征展开了多尺度精细化描述分析。这些理论和试验研究加深了对岩石各向异性特性的认识，并在很大程度上指导了实际工程应用。

笔者基于现场三轴试验成果，采用UDEC[12]方法，对丹巴石英云母片岩的各向异性问题进行数值模拟分析研究，重点考察不同片理方向加载时岩体的强度和变形各向异性特征。之后针对CPD1探洞进行开挖响应模拟，旨在研究考虑岩体各向异性特征后隧洞围岩的破坏机制，并与现场测试结果和工程现象进行对比，用以评价数值模拟方法的适宜性，为丹巴深埋软岩工程的开挖和支护设计提供一定的基础支撑。

1 工程概况及现场三轴试验

1.1 工程概况

丹巴水电站位于四川省丹巴县境内的大渡河干流，电站装机容量1 160 MW，厂坝之间采用长约17.4 km的两条长引水隧洞连接，隧洞直径约14～16 m。目前地质资料显示，引水隧洞尾端的Smx4-2地层为石英云母片岩夹少量云母石英片岩、薄层条带状长英质变粒岩，属于软岩，岩层总体产状为N45°～60°W，NE∠45°～85°，其岩体力学强度较低，具有明显各向异性特征。此段累计长度约2.8 km，约占引水洞总长的17%，分布于引水线路出口段、调压井、压力管道等枢纽部位，最大埋深达1 200 m。

工程面临的深埋软岩问题十分突出，区域地应力水平较高（实测最大主应力值为17～30 MPa），大断面软岩隧洞、调压室施工过程中均可能出现塑性大变形、流变及局部失稳等现象，施工和支护将面临较大困难，因此有必要对丹巴软岩的岩体力学特性进行深入分析研究。在工程前期开展大量的室内物理力学性试验、岩体现场原位试验等试验研究工作[13]，获得的一致性结论为丹巴软岩各向异性突出，因此在开展室内、外各类岩石力学试验时，均分别开展了与片理面呈不同夹角的加载试验研究工作。片理作为石英云母片岩中最普遍的结构面，是导致片岩各向异性宏观特征的主要原因，各类试验均表明，在不同片理方向加载时的岩体强度和变形模量均存在较大差异，这在实际工程中不容忽视。

1.2 现场三轴试验

为深入了解丹巴软岩物理力学特性，在工程现场开展了软岩现场三轴试验[14]。试验支洞选择在CPD1平探洞桩号K0+540 m处，支洞方向与石英云母片岩片理走向近一致，共计13个试验点。坐标系方向规定如下：x向为顺片理面水平方向（支洞洞轴向），y向为近垂直片理面方向（支洞洞径向），z向为近平行片理面方向（铅直向）。试样总体沿石英云母片岩片理面方向制备，但由于片理面不是完全直立，所以试样与片理面之间有一定夹角（约10°）。试样形态及总体布局如图1所示。

图1 现场三轴试验洞位置及试样形态示意图Fig.1 Sketch of tunnel and specimen in the triaxial test

首先，针对现场岩体三轴试样，开展了不同片理加载方向（因试样条件限制，主要为垂直或平行片理面加载）、不同应力水平下的变形试验，获得了垂直或平行片理加载方向上未松弛岩体和松弛岩体的弹形模量。试验中，当试样有较高围压时，可以认为试样处于相对未松弛状态，当围压水平较低时，认为试样处于松弛状态。表1中给出了分析整理后的综合评定结果。

表1 现场试验弹性模量综合评定值[14]Table 1 Elastic modulus in the field triaxial test

此后选取其中7个试验点开展三轴强度试验。此7个试验点在完成变形试验之后，首先开展了卸载（小主应力）破坏试验，用以研究完整岩体的强度特性。然后针对破坏试样，分别开展了重复卸载或重复加载试验，用以研究松弛岩体的强度。通过对现场多组试样进行加卸载三轴试验，获得了各试件破坏时的应力状态曲线，包括图2（a）各试件首次卸载试验的卸载破坏峰值强度点分布、图2（b）单点重复卸载破坏峰值强度点分布。

将由最大、最小主应力表征的强度点进行线性拟合，可得到未松弛岩体和松弛岩体的σ1-σ3关系式：

根据Mohr-Coulomb屈服准则，参照式（2），可计算得各个试样卸载破坏时的峰值强度参数，如表2所示。受现场工程环境和试验条件所限，未开展与片理面呈多角度的现场大型三轴加载试验，上述强度试验可认为是平行于片理方向加载的强度参数。

式中：φ、c分别为完整岩体的内摩擦角和黏聚力。

图2 三轴卸载试验破坏强度及拟合曲线[14]Fig.2 Distributions and fitting curves of the collapse strength in the triaxial cyclic loading tests

表2 现场试验岩体强度综合评定值Table 2 Strength properties in the field triaxial test

2 各向异性数值模拟研究

2.1 UDEC简化模型与参数研究

基于UDEC程序的简化三（单）轴压缩数值模型岩石数值加载试验方案如图3所示，围压（σ3）平行于x方向施加于试件的左右边界，同时在上边界沿y轴负向施加恒定的速度Vy，以产生轴向压缩应力（σ1），使试件逐步受压。结合现场三轴试验模型，基于UDEC模拟的模型尺寸为0.5 m×1 m，加载试件由完整岩体和分布于其中的一组等间距弱面（片理面）组成，图3（b）给出了不同片理方向加载的三种典型方案。岩体和结构面均服从Mohr-Coulomb屈服准则。通过编制相应的FISH命令，控制稳定的加载速度，并监测加载过程中的应力-应变关系变化曲线，由此求得不同片理方向加载时相应的等效变形模量和岩体强度。

图3 简化的二维数值加载模型Fig.3 Simplified 2D numerical simulation loading model

式（3）和（4）分别给出了基于上述数值加载试验获得的试样泊松比和弹性模量求解表达式。式中Lx和Ly为模型尺寸，ux和uy分别为x、y方向上的变形量，νm和Em分别为岩体的泊松比和弹性模量。

通过反复试算，并与石英云母片岩的现场三轴试验结果（见表1和表2）进行对比，反演参数与现场试验结果较吻合。反演的岩体、结构面力学参数见表3和表4。需要指出，UDEC对石英云母片岩的模拟采用岩石基质+片理的方式，因此表3中的岩石基质参数与表1、表2中岩体参数不具备直接的对应关系，岩石基质+片理的综合体的数值试验结果才能与表1和表2中的参数进行对比。

表3 岩石基质参数反演值Table 3 Mechanical parameters of rock

2.2 丹巴石英云母片岩各向异性特征的数值分析

图4为不同片理方向加载时的单轴压缩试验轴向应力-应变关系曲线，可见不同片理方向加载时对应的试件峰值强度存在较大的差异，体现了岩体强度的各向异性特点。当加载方向与片理接近平行或呈大角度（50°～90°）相交时，对应的峰值强度最大，达到15.6 MPa。当加载方向与片理斜交（10°～40°）时，对应的峰值强度减小，β=30°时达到最小值，约为5.5 MPa。

图4 不同片理方向加载下的应力-应变关系曲线Fig.4 Stress-strain plots of uniaxial tests with different schistosity angle

由于单轴压缩试验应力-应变关系曲线中弹性段的斜率近似等于岩体的弹性模量，图4中显示了随着β角的增大而含弱面试件的等效弹性模量有逐渐降低的趋势，且总体下降幅度较大。图5给出了单轴加载条件下的试件等效弹性模量的具体数值及变化规律。相比平行于片理面加载的弹性模量，垂直于片理加载时下降了近40%，这与表1中现场试验结果相一致。

图6给出了不同围压下试件三轴抗压强度（σ1）与不同片理加载角度间的关系。随着围压的增加，岩体强度的各向异性特点仍然存在，且随加载角β的变化规律与无侧限时相似。图6给出了基于Jeager[1]单一弱面理论获得的单轴压缩下的解析曲线，与UDEC计算结果对照，两者分布规律基本一致。

图5 不同片理方向加载下的弹性模量Fig.5 Elastic modulus of uniaxial tests with different schistosity angle

图6 不同围压、不同片理方向加载下的岩体强度曲线Fig.6 Strength curve of specimen under different confining pressures and different schistosity angle

Jaeger提出的单一弱面强度理论，其主要思想为将岩体假定为含一组定向弱面的各向同性体，基于Mohr-Coulomb准则，可求得沿弱面发生滑动破坏的强度条件：

式中：cw和φw分别为弱面的黏聚力和内摩擦角。

将式（5）对β求偏微分，并令方程式为零，从而可得最小破坏强度发生在：

此外，岩体强度和破坏模式会随弱面倾角的变化而变化，当β趋于0或者β在(90°-φw)至90°范围时，岩体不再沿着弱面滑动，而是产生贯穿弱面的剪切破坏，此时的强度为完整岩体强度。

由于上述基于UDEC的数值模拟采用的假设与Jeager相同，因此计算结果表现出的岩体强度各向异性特点与Jeager提出的理论曲线基本一致。

表5给出了不同片理加载方向下试件的力学指标等效值，包括等效黏聚力、内摩擦角、UCS和弹性模量。其中各片理加载角下的等效黏聚力和内摩擦角根据式（2）获得，其变化规律与单轴压缩强度一致，β约为30°时最小。

3 工程应用

为了勘察丹巴引水隧洞软岩段的工程地质条件，从尾水出口部位往山体内开挖了洞径为3 m的勘探平洞CPD1，并在开挖过程中布置了多个监检测断面，用以评估该软岩围岩的成洞条件。现选择其中一个典型剖面（桩号K0+520），模拟其开挖响应过程。考虑到平洞开挖断面规模相对较小，可直接采用现场三轴试验的相关数据来模拟平洞的开挖响应。结合UDEC程序，将软岩片理结构考虑成一组规则的结构面进行等效模拟，相关参数见表3和表4。具体的计算模型、洞室尺寸及地应力条件如图7所示。

图7 CPD1平洞典型剖面形态及初始条件Fig.7 UDEC model and initial conditions of CPD1 tunnel

图8给出了平洞开挖后的围岩塑性屈服区分布，包括片理面的破坏情况。总体上，塑性屈服在左右边墙及底板分布较广，深度在0.5 m左右，顶拱处塑性屈服深度较小，仅约0.2 m。左侧边墙至左拱肩围岩和右侧边墙下部片理面剪切滑移破坏现象较突出，这与现场观察到的一些破坏现象较一致。平洞规模小、岩体强度低，此洞段现场未发现有剧烈的应力性破坏现象，但左侧边墙、左拱肩及右边墙的零星片状剥落、掉块或鼓胀松弛现象时有发生。从这些现象的位置、破坏模式及程度看，一方面受到了较高的初始地应力场分布特点所影响，同时也与片理和围岩的空间关系密切相关，具体表现为普遍分布的片理面引起的围岩各向异性特征。

图8中显示了围岩开挖面的破坏现象与片理方向存在较明显的相关性。由于开挖面不同位置与片理方向的交角不一样，表现出了不同的破坏模式、规模，主要可分为三类：①类，位于左右侧墙，受力形式相当于加载方向与片理面呈较小角度（β≈30°）相交，属于弱面滑动破坏类型；②类，位于左侧拱肩处，受力形式相当于加载方向与片理面呈极小角度（β＜10°）相交，岩体在压弯作用下容易出现沿层面劈裂后的折断或屈曲破坏；③类，位于顶拱和右拱肩，受力形式相当于加载方向与片理面呈大角度（β＞50°）相交或垂直片理面，破坏形式以岩体剪切破坏为主。结合图5所示的岩体强度各向异性特点，其中①类加载模式下的岩体强度最低，塑性破坏区深度最大，②类和③类的强度明显要高一些，对应的破坏区深度相对较小。若不计开挖边界形态的影响，UDEC计算获得的围岩塑性破坏区分布特点与此前不同片理方向加载数值试验获得的规律具有一致性，体现了岩体各向异性导致隧洞围岩不同部位失稳型式的不一致。在实际工程中，深入理解此类软岩隧洞中不同部位的变形破坏特征差异性，有助于合理采用针对性的开挖支护方式，保证支护措施的可行性和有效性，提高工程整体安全稳定性。

表5 不同片理方向加载的试件力学指标等效值（未松弛岩体）Table 5 Mechanical equivalent parameters of the specimen with different schistosity angle

图8 CPD1平洞开挖响应结果及实测松动圈轮廓Fig.8 Sketch of plastic failure zone of the simulated and field test results

图8中给出了桩号K0+520断面的两组松动圈（岩体低波速带）测试成果，显示了围岩的松动圈深度约为0.5～0.7 m，其中在顶拱处的深度为0.3～0.5 m。

对比UDEC计算结果，侧墙的塑性破坏区分布范围较好地复核了松动圈测试成果，但计算获得的顶拱塑性屈服区深度与测试结果比，相对偏小。

4 分析与讨论

通过采用UDEC对丹巴石英云母片岩的各向异性特征进行数值模拟分析，并对比现场试验结果、监测结果及相关工程经验，有必要针对以下几个问题进行探讨和说明。

4.1 UDEC模拟片岩各向异性存在的不足

采用UDEC模拟获得的各种加载条件下的岩体强度、变形参数等与现场试验大体上符合，基本体现了片岩的各向异性特点，但也存在不清晰的地方。一方面，现场试验数据量偏少，离散性较大，导致模拟结果的对照性较差；其次，此处模拟中采用了Mohr-Coulomb屈服准则，计算结果也与Jaeger提出的传统单一弱面强度理论解一致，但这并非与室内试验和现有的工程经验完全匹配。

实际上大量试验研究表明，Jaeger破坏准则对沿着弱面的滑动破坏与完整岩石的破坏描述过于简化，特别是从斜交弱面的剪切破坏向沿着弱面滑移破坏的过渡范围内，往往会高估岩体的强度，将此情况下的强度等同于完整岩石介质的强度，忽略了斜穿弱面及沿着弱面滑动的混合破坏模式对整体强度的影响。图9（a）为丹巴室内试验获得的试件各向异性强度曲线，图9（b）为Hoek和Brown[8]获得的页岩、板岩和砂岩在不同弱面加载角、不同围岩下的强度各向异性曲线。经验表明，U型变化的曲线更趋于合理，Jeager破坏准则中平折线的存在使其与实际岩体强度的各向异性规律存在一定的偏差。此外平行于弱面加载和垂直加载时，对应的强度值也不尽相同。

比较CPD1平洞数值模拟塑性破坏区与现场松动圈测试结果存在的差异，除了二维简化计算合理性、监测结果可靠性、爆破对围岩的开挖扰动等影响外，上面阐述的采用UDEC模拟片岩各向异性的方法可能会因高估岩体强度而导致顶拱区域的塑性区深度要小于实际情况。

4.2 岩体各向异性的数值模拟方法

图9 丹巴片岩及其他岩体的强度各向异性试验曲线Fig.9 Strength anisotropy of Danba slates and other rock in the uniaxial or triaxial compression test results

考虑到岩体普遍存在的不均一性、各向异性、不连续性等特点，要通过数值模拟方式表征岩体各向异性的室内或现场试验曲线特性，则数值模型中必须能够反映岩体的真实结构特点。若在计算中仅简单等效为横观各向同性的连续介质，或将岩体等效为均质岩体和弱面的联合体的方式，则显然忽略了岩体中赋存的其他构造或缺陷对岩体各向异性的影响。而这很难在有限差分、有限元等方法中加以考虑，当前可尝试采用基于细观颗粒离散元法的PFC进行描述［11］。颗粒离散元法中，因细观颗粒层次上的不连续性诱发的复杂宏观非连续行为并不需要复杂的本构方程来实现，这意味着颗粒离散元模型可能比其他数值模拟技术更切合实际，可适当地表征岩体的破坏特性，借此就可以直接对真实的岩体结构开展各向异性分析工作。

4.3 尺寸效应对岩体各向异性的影响

现场三轴试验参数不能完全地解释CPD1平洞开挖响应现象的另一个重要原因可能是岩体尺寸效应的影响。在采用现场试验数据进行CPD1平洞开挖模拟时，其尺度是否与开挖尺度相匹配值得思考。从小尺度上看，如室内加载实验（0.05 m×0.1 m）和现场三轴试验（0.5 m×0.5 m×1 m）等，片岩可等效为横观各向同性材料，这样有助于简化分析工作。但当工程研究尺度增加后，如CPD1平洞开挖尺寸3 m×3 m、引水隧洞的直径14～16 m、软岩圆筒型调压室的直径30 m时，岩体中包含的结构面将大幅增加，不仅包含了片理面，还可能有优势结构面或其他无规律性结构面等存在，这将导致不同尺度岩体的各向异性特性也表现出不同的规律。尽管随着岩体内结构面数量的增加，岩体强度各向异性可能向等效各向同性转换，但单一或少数优势结构面导致的各向异性影响却不能忽视。

对CPD1平洞桩号K0～920 m洞段的构造裂隙发育情况进行了统计分析，共计426条结构面，区域节理裂隙中等发育，节理密度见图10所示。丹巴软岩段发育有两组优势结构面，其中一组的产状为N40°～70°W，NE∠45°～80°，另一组的产状为N45°～75°E，ES∠40°～60°。可见，随着研究尺度的增加，这两组优势结构面带来的岩体各向异性问题将逐步显现，需要在大尺度工程问题中加以考虑。

图10 丹巴软岩段节理密度图Fig.10 Joint density map of Danba schist

一般而言，进行较大规模的现场试验可得到相似尺度的岩体强度和变形特性，并有助于对大尺度工程问题的认识。但此类试验难度大、成本高，不利于开展大量的试验工作。通常情况下工程人员多通过经验的方法评估岩体物理力学参数，如Hoek等人提出的GSI系统等，但其并未体现大尺度岩体的各向异性特征。当前，可采用数值模拟试验分析此类问题，如通过建立多组结构面组成的节理裂隙网络模型（DFN），可以较真实地反映此类岩体的非连续各向异性特点，但对复杂结构面的表现方式、本构关系、空间分布形态以及力学参数等的表征均还存在难度，有待进一步深入研究。

5 结语

通过对丹巴石英云母片岩各向异性特征开展深入分析研究工作，有如下认识：

（1）现场试验结果显示了丹巴石英云母片岩具有明显的各向异性性质。基于UDEC对此类岩体各向异性特性进行了模拟分析，得到了不同片理方向下加载试验的力学指标等效值。当加载方向与片理接近平行或呈大角度（50°～90°）相交时，对应的峰值强度最大可达到15.6 MPa；当加载方向与片理斜交（10°～40°）时，峰值强度减小，其中β=30°时约5.5 MPa；垂直于片理加载时，试件的等效弹性模量相比平行于片理面加载下降了近40%。数值计算获得的片岩各向异性规律与试验结果以及Jaeger提出的单一弱面强度理论基本一致。

（2）CPD1平洞开挖响应特征表明，隧洞围岩不同部位的破坏形式与片理和开挖边界的几何切割关系密切，考虑片岩各向异性特性可以较好地解释现场围岩开挖现象。深入理解此类软岩隧洞中不同部位处的变形破坏特征差异性有助于合理采用针对性的开挖支护方式，保证支护措施的可行性和有效性。

（3）UDEC模拟结果能一定程度上反映丹巴片岩的各向异性特点，可反映岩体变形各向异性特征，也可较好地反映沿弱面滑移破坏时的低强度特征，但并不能很好地模拟含弱面岩体的穿层破坏现象，后续有待更深入研究。