层状岩体受压力学特征结构面效应数值分析

鲁海峰，姚多喜，沈 丹，温 亮，王 康

(1.安徽理工大学地球与环境学院，安徽 淮南 232001;2.煤炭资源与安全开采国家重点实验室(中国矿业大学)，北京100083)

煤层底板岩体是沉积岩，其形成过程中的沉积作用和矿物颗粒的择优取向，具有显著的层状结构特点。由于分布有一组占绝对优势的结构面(层面)，层状岩体的变形和强度特性具有明显的各向异性。在承压水上采煤突水预测分析中，底板层状岩体的变形和强度是其重要的输入参数。由于层状岩体强度不仅与完整岩块的性质有关，还与结构面方向和性质有很大的关系，表现出复杂的力学性质，国内外学者作了很多的研究工作。Jaeger 针对节理岩体沿着节理面滑动破裂提出了相应的破坏准则［1］，Tien 等用2 种不同的材料人工预制了3 组不同倾角的层状岩石，研究了横观各向同性体倾角对岩石强度和弹性模量的影响，并针对横观各向同性岩石破坏的2 种不同模式提出了相应的破坏准则［2］。何沛田等从岩石在单轴压应力状态下的性态和岩石工程问题的应用出发，论述了黑灰色钙质页岩岩块层理结构面夹角θ 对其破坏特征、强度和变形的影响［3］。苏志敏等在分析大量试验数据的基础上，探讨了页岩的层状结构面倾角和围压对其强度的影响［4］。张玉军等分析了层状岩体强度异向性对地下洞室稳定性的影响，并提出了相应的经验公式［5］。

以上研究主要通过理论分析和室内试验来进行。近年来，随着计算机技术的不断发展，采用数值模拟方法分析岩石力学问题成为新的发展趋势［6－8］，作为理论分析和试验研究的辅助工具，采用数值分析方法能直观地得到岩体的力学行为特性，从而可为层状结构底板岩体采动破坏分析提供借鉴。为此，本文将层状岩体等效为横观各向同性连续体，运用FLAC3D 数值模拟软件，利用改进的遍布节理模型对层状岩体进行了大量的单轴和三轴压缩试验模拟，并将模拟结果与理论分析结果进行对比，验证数值模拟方法的有效性。

1 数值模型的建立

1.1 横观各向同性弹塑性本构模型

FLAC3D 自身携带的遍布节理模型可作为层状岩体的各向异性屈服准则，但其本构关系仍然是各向同性的，不能够反映层状岩体的变形特征。为此，本文将横观各向同性本构关系引入到遍布节理模型中。在图1所示的横观各向同性体中，以结构面为x'－z'平面，法向为y＇轴，则在局部坐标系x'y'z'中横观各向同性线弹性本构关系采用式(1)表示。

式(1b)中，E1，u1为横观各向同性面(x'oz'平面)内的弹性模量和泊松比;E2，u2为垂直横观各向同性面(y＇轴方向)内的弹性模量和泊松比;G2为垂直横观各向同性面(与x'oz'平面垂直)内的剪切模量。

图1 局部坐标与整体坐标关系

局部坐标系x'y'z'和整体坐标系xyz下的应力关系可由式(2)表示，即

(x'，x)为局部坐标轴x＇与广义坐标轴x之间的夹角;［C］T为［C］的转换矩阵。

遍布节理模型的破坏准则包含Mohr－ Coulomb 体内特殊方向上的弱面，根据应力状态、结构面产状以及模型体和结构面的材料特性的不同，屈服可能发生在岩块内，或者发生在结构面上，或者在两个部分都发生。岩块采用Mohr－Coulomb 屈服准则，对于结构面，独立设置局部坐标系，其方位分别用结构面倾角dip 和方位角dd 来表达，同样采用Mohr－Coulomb 屈服准则。

为反映层状岩体的峰后力学性质，本次在横观各向同性弹塑性本构模型中考虑了岩体的硬化－软化特性。根据FLAC3D 模拟应变软化分段线性原则，预先定义软化参数，在塑性应变产生后，部分或所有单元的屈服参数，如粘结力、内摩擦角、剪胀角和抗拉强度都可发生变化，并通过引入塑性参数来实现。对于岩块分别定义塑性参数ks和kt来描述塑性剪切软化行为和塑性拉伸软化行为，它们主要与岩块的塑性剪切应变和塑性拉伸应变密切相关。建议的增量型数学表达式采用如下形式［9］2－46

拉伸硬化增量则由塑性拉伸应变增量

根据ks和kt的具体表达式，则可定义岩块的力学参数为塑性参数ks和kt的函数为

对于结构面，和岩块相同，同样引入塑性参数ksj和ktj来描述结构面的塑性剪切软化行为和塑性拉伸软化行为。结构面的这两个塑性参数也采用与岩块同样的建议方法，它们主要与结构面的塑性剪切应变和塑性拉伸应变密切相关，其表达式为［10］

调用商业软件的源代码，采用C++语言，将横观各向同性的弹性本构关系、Mohr－Coulomb 屈服准则及岩石应变硬化－软化实现，生成.dll 文件后嵌入到FLAC3D软件中，实现过程见文献［9］(4－1，4－18)。

1.2 层状岩体压缩数值模型建立

单、三轴数值模型中存在一组优势软弱结构面，试件采用长方体模型，模型底部长宽为5 cm，高度为10 cm;结构面倾角α 为0°～90°，变化梯度为5°。建立三维试件计算模型，部分倾角模型如图2所示。试件加载方式为位移加载，加载速度为0.002 mm/步，分别记录加载过程中试件应力－应变的关系，数值计算中所用的力学参数见表1。

图2 数值计算模型

表1 数值试验中的层状岩体力学参数

根据分段线性原则，预先定义粘聚力和内摩擦角的软化参数。本文根据所研究岩体特征并参考FLAC3D 对于软化模型参数的建议，当其发生塑性变形后，定义相应的cs、φs、csj、φsj和原始c、φ、cj、φj值与塑性参数的关系如表2所示。

表2 岩块与结构面抗剪强度参数与塑性参数关系

2 结果分析与讨论

2.1 压缩强度与结构面倾角关系

单、三轴压缩计算机模拟试验得到的层状岩体压缩强度与结构面倾角α 的关系如表2 和图3所示。从中可看出，岩样的压缩强度受结构面倾角的变化而变化。岩体的压缩强度随结构面倾角呈现先减小后增大的趋势;当结构面倾角α 为60°时，岩样的压缩强度最小，该角度为层状岩体最不利的结构面倾角;当结构面倾角在α=35°～50°或α=75°～85°时，岩样压缩强度的变化梯度最大，说明岩样的压缩强度受此结构面倾角范围影响较大，该结果与文献［2］1010中的室内试验结果相同。

同时，由Mohr－Coulomb 准则可得，试件沿结构面发生剪切破坏时，结构面破坏(极限平衡)时的主应力应满足的条件

式中:cj、φj为结构面的粘结力和内摩擦角。

从式(8)可看出，cj、φj均为常数，固定σ3不变，当α 趋于90°以及当α 趋于φj时，σ1趋于无限大，当然，σ1不可能无限大，当其达到一定值时岩石材料发生破坏。由此得知，当结构面倾角α 满足φj ＜α＜90°时，岩体沿结构面发生破坏，当α 不满足上述条件时，破坏沿岩石材料内部发生。

对式(8)倾角α 求导，并令导数等于0，可以求得当倾角α=45° +φj/2 时σ1－σ3有最小值，相应的σ1的最小值σ1，min为

式中:Nφj=tan2(45° +φj/2)

同理岩块破坏时满足如下表达式

式中:c、φ 分别为岩块的黏聚力、内摩擦角。

当试样沿结构面发生剪切破坏时，按式(8)可得压缩强度，当试验不沿结构面破坏时，可由式(10)计算得压缩强度，计算结果如表3所示。

表3 层状岩体压缩强度数值计算与理论计算结果对比 MPa

从表3 可以看出，对于不同的结构面倾角，数值计算结果与理论分析的结果基本相同，仅在α=85°时，二者存在一定的差别，最大差别为在围压11 MPa 时，数值解为70.87 MPa，理论解为73.78 MPa，相差2.91 MPa，但两者之间的差别小于4%，其余工况两者之间的差别大都在2%以下。另外从表2 中可看出，层状岩体的力学行为具有显著的围压效应。在围压一定时，压缩强度与倾角α 的关系曲线呈U 型，压缩强度随着α 的增加先减小后增大;对于相同倾角的试件，随着围压的增大，岩石三轴压缩强度不断增大。

计算机模拟试验得到的结论验证了前人的试验及理论推导［10］，同时也说明本文建立的横观各向同性弹塑性模型可以很好的反映层状岩体的强度各向异性。

2.2 压缩变形特征

层状岩体压缩试验得出的等效弹性模量与结构面倾角α 间的关系如图3所示。由图3 可见，层状岩体弹性模量E一般沿平行于结构面方向最大，而沿垂直于结构面方向最小，并随结构面与水平面(或最大主应力)之间夹角的增加而增大(减小)，这一规律与一些试验结果是一致的［11］，表明本文的数值方法在描述层状岩体的变形方面是可行的。

图3 弹性模量与结构面倾角关系

图3 是倾角60°的岩样在1MPa、6MPa、11MPa和16MPa 围压作用下模拟得出的应力－应变全程曲线。由图3 可以看出，随着围压的增加岩样的屈服应力和峰值强度以及残余强度均逐渐增大;另外可看出，随着围压的增大，岩样破裂特性逐渐由脆性破坏向延性破坏转变。为了更清楚地说明围压和倾角的变化对峰值强度和全程曲线峰后段的影响，图5 给出了围压1MPa、不同倾角岩样的应力－应变全程曲线。从图中可以看出，由于岩样结构面倾角不同，压缩曲线的斜率表现不同，岩样可发生沿结构面的剪切破坏，也可发生穿切结构面的剪切破坏，此时，岩体的抗压强度达到最大值，由于结构面倾角不同，导致该最大值也不同。同时从图5 中看出，所有的压缩曲线皆存在较为明显的尖点，说明岩样的脆性特征明显。围压相同时，压缩曲线随倾角的变化，峰后曲线呈现出不同的特性，且对峰值强度以及残余强度都有显著影响。峰值强度相同，残余强度也接近;峰值强度大的，残余强度大。

图4 倾角60°层状岩体压缩应力－应变曲线

图5 围压1MPa 层状岩体压缩应力－应变曲线

2.3 破坏模式

图6 显示的是模拟单轴压缩情况下试件最终破裂时候的位移矢量分布图。从图6 中可以清楚看出，岩石试件各点最终时的位移情况，最大的位移量都出现在破裂面上。另外结合表3 分析，当结构面倾角小于40°和大于85°时，岩体强度基本上由岩块控制，整体较高，岩样在其中央发生鼓胀破坏，且主控破裂面穿越结构面;当倾角在40°和80°之间时，岩体强度整体较低，基本上由结构面控制，呈剪切滑移状，且主控破裂面主要沿结构面展开。

图6 试样破坏位移矢量图

3 结论

通过FLAC3D 数值软件，结合改进的遍布节理本构模型，建立了层状岩体压缩数值模型，分析了单轴、三轴压缩情况下的应力应变响应以及强度特征，通过全文研究，可得出以下结论:

1)数值分析表明，层状岩体压缩强度具有显著的结构面倾角效应。随着倾角的增大，层状岩体的压缩强度呈现先减小后增大。当结构面倾角为40°～80°时，岩体强度整体较低，破裂面主要沿结构面展开，数值试验和理论分析反映的岩体强度变化规律一致，表明本文提出的数值模拟方法是有效的，可用来对层状岩体压缩特征进行分析。同时模拟发现，对于相同倾角的试件，随着围压的增大，岩石三轴压缩强度不断增大。

2)数值分析表明，层状岩体的等效弹性模量沿平行于结构面方向最大，而垂直于结构面方向最小，并随结构面与水平面之间夹角的增加而增大。

3)数值分析表明，随着围压的增加，岩样的屈服应力和峰值强度以及残余强度均逐渐增大，同时岩样破裂特性逐渐由脆性破坏向延性破坏转变。同时随着结构面倾角α 的增大，压缩曲线斜率逐渐增大。

［1］JAAGER J C.Shear failure of anisotropic rocks［J］.Geology Magazine，1960(97):65－72.

［2］TIEN Y M，TSAO P F.Preparation and mechanical properties of artificial transversely isotropic rock［J］.International Journals of Rock Mechanics＆Mining Science，2000，37(6):1 001－1 012.

［3］何沛田，黄志鹏.层状岩石的强度和变形特性研究［J］.岩土力学，2003，24(S1):1－5.

［4］苏志敏，江春雷，Ghafoori M.页岩强度准则的一种模式［J］.岩土工程学报，1999，21(3):311－314.

［5］张玉军，刘谊平.层状岩体的三维弹塑性有限元分析［J］.岩石力学与工程学报，2002，21(11):1 615－1 619.

［6］梁正召，唐春安，李厚祥，等.单轴压缩下横观各向同性岩石破裂过程的数值模拟［J］.岩土力学，2005，26(1):57－62.

［7］尹立明，张培森.双结构面对岩体强度影响的模拟分析［J］.采矿与安全工程学报，2010，4(27):600－603.

［8］周科峰，李宇峙，柳群义.层状岩体强度结构面特征的数值分析［J］.中南大学学报:自然科学版，2012，43(4):1 424－1 428.

［9］ITASCA CONSULTING GROUP Inc.FLAC3D(Version 2.1)users manual［M］.USA.Itasca Consulting Group Inc.，2003:1－18.

［10］黄书岭，徐劲松，丁秀丽，等.考虑结构面特性的层状岩体复合材料模型与应用研究［J］.岩石力学与工程学报，2010，29(4):743－756.

［11］席道瑛，陈林，张涛.砂岩的变形各向异性［J］.岩石力学与工程学报，1995，14(l):49～58.