压缩条件下玄武岩柱各向异性及破裂机理数值模拟研究

王永艺，唐春安，夏英杰

(大连理工大学 海岸和近海工程国家重点实验室， 辽宁 大连 116024)

柱状节理岩体是一种由玄武岩浆冷凝收缩而形成的天然岩体，其柱状节理往往将岩体切割成规则或不规则的棱柱体，由于柱状节理构造的存在，其变形和强度具有显著的各向异性特性。此外，部分柱状节理岩体赋存于围压环境下，如金沙江下游的白鹤滩水电站、四川乐山铜街子水电站、云南丽江金安桥水电站，金沙江峡谷段的溪洛渡等水电站均位于柱状节理玄武岩之上，围压环境对其各向异性存在影响。因此，研究不同围压条件下玄武岩柱各向异性及其破裂机理，具有较高的实际和工程应用价值。

针对柱状节理岩体力学各向异性特性，已有相关学者开展相关研究。Xia等[1]提出了一种利用3DP和相似常数精确重建不规则柱状节理岩体结构的合适方法，对重建后的试件开展单轴压缩试验，将其试验结果与现场试验结果进行了比较。Ji等[2]采用水泥、细砂、水、减水剂，制作柱状节理岩体试件，对其开展单轴压缩试验研究，分析不同柱体倾角情况的试件强度变化及破坏特征。肖维民等[3-4]通过单轴压缩试验和三轴压缩试验得到柱状节理岩体在不同柱体倾角下的变形模量和单轴抗压强度，分析柱状节理岩体变形和强度的各向异性特性。柯志强等[5]通过单轴压缩试验，研究柱体倾角和横向节理对岩体各向异性力学特性及破坏机制的影响。郑文棠等[6]采用可变形离散元法建立了柱状节理玄武岩体的三维离散元数值模型，通过数值模拟不同尺寸的承压板试验，探讨了尺寸效应和各向异性对试验成果的影响。崔臻等[7]利用节理网络有限元为工具，研究各结构效应表征参数对柱状节理岩体等效变形模量的影响。倪海江等[8]开展基于离散元的柱状节理岩体等效弹性模量尺寸效应研究。以上研究，在柱状节理岩体的力学特性的认识上，取得了有益的研究成果；但关于不同围压条件下玄武岩柱各向异性及破裂机理，尚需进一步研究。Wang等[9]认为玄武岩柱的围压效应与模型尺寸、柱体倾角、柱体直径、柱体不规则程度、横向节理的错距比、节理力学性质的变化、岩石细观本构的变化、模型边界条件相关，并开展了相应的数值试验研究，但在其平行柱轴方向，有围压作用的玄武岩柱的渐进破裂过程和破坏模式的研究环节中，岩石残余强度系数仅取0.1这一工况。实际工程中，当围压增加，岩石的残余强度系数有可能随之增加。因此，本文将岩石残余强度系数设置为0.5，然后仅改变围压的大小，以研究不同围压条件下不同柱体倾角玄武岩柱的强度、变形特性，以及各自的破裂机理、损伤破裂特征，其研究成果可供对比参考。

单轴及围压条件下，玄武岩柱有强度和变形的各向异性，不同柱体倾角的玄武岩柱的破裂机理和破坏模式存在差异。本文构建不同柱体倾角的玄武岩柱图像，然后基于RFPA3D-CT软件的数字图像处理，将玄武岩柱图像转化为有限元网格模型，并分别赋予节理、岩石的材料力学参数，单轴及围压条件下玄武岩柱数值试验，以研究其强度和变形的各向异性，及其破裂机理与破坏模式。

1 数值模型

1.1 RFPA3D-CT原理简介

RFPA3D-CT软件系统是在RFPA3D-Basic软件[10]的功能基础之上，增加了图片导入、构建数字模型、图像阈值分割、赋予材料参数和构建三维数值模型的功能。本文借助RFPA3D-CT的图像处理功能，将数字图像转化为有限元网格模型，其原理如下。为构建数值模型，需要将图片中的信息转换为建模所需的矢量化数据。数字图像是由正方形像素点组成的，如图1(a)所示，在三维空间中，若认为图像具有一定的厚度，则可将每一个像素点可看作一个有限元网格。将各个像素点的角点坐标转换为相应的矢量空间物理位置(其中每个像素点具有相应的厚度和边长)，并根据像素点灰度值不同，将其归类为不同岩石材料，赋予相应的材料参数。根据上述原理，转化后的有限元网格模型如图1(b)所示。RFPA3D-CT基于细观损伤力学和统计强度理论。细观单元采用最大拉应力准则和Mohr-Coulomb破坏准则。当细观单元的最小主应力超过其单轴抗拉强度时单元产生拉伸损伤；如果细观单元应力状态满足Mohr-Coulomb破坏准则时，细观单元产生剪切损伤，如图2所示。细观单元承载能力随损伤演化过程而降低，在达到破坏准则之后仍保持一定的残余强度。有关RFPA3D-CT详细的原理可参见相关文献[10-15]。

图1 数字图像转化为有限元网格模型示意图

图2 单轴应力下单元的弹脆性损伤本构关系

1.2 数值模拟验证

数值模拟验证环节，本文采用Ji等[2]、肖维民等[3]的室内物理试验对本文的数值试验进行验证。本文用于数值验证的试件，采用的是宽度50 mm，高度100 mm的矩形试件，平面应变情况；内部的正六棱柱体的直径是10 mm；考虑了平行柱轴方向β=15°的情况；基于RFPA3D-CT，将数字图像转化为有限元计算模型；有限元模型的力学参数取值见表1，其取值参考了玄武岩柱的相关文献资料[1-8]。数值试验采用位移控制加载，加载量为0.005 mm/step，直到试件破坏。数值试验与室内物理试验的试件破坏模式对比，见图3。

表1 数值模拟验证环节，玄武岩柱的岩石及节理的力学参数取值

图3 单轴压缩条件下数值试验与室内物理试验的试件破坏模式对比

1.3 数值模型设置

数值试验中，各个模型的单元尺寸是相同的，其中以3 m×3 m试件为例，该试件的单元数为608 400。图4展示了围压条件下玄武岩柱尺寸效应数值试验的典型设置和边界条件。玄武岩柱的岩石及节理的力学参数取值,见表1。在每一个模型的顶部施加竖向位移荷载；每一步施加的位移量，与模型初始侧向边长的比值，为0.000 017；逐步施加位移荷载，直至试件破坏。

图4 围压条件下玄武岩柱试件的典型设置和边界条件(以β=45°为例)

2 结果和分析

2.1 不同围压条件下不同柱体倾角玄武岩柱的强度及变形特性

由图5可知，在抗压强度方面，对于围压0 MPa～8 MPa的情况，玄武岩柱试件的抗压强度随柱体倾角的增加大致呈U型分布；随着围压的增加，各柱体倾角的试件抗压强度，有较明显的提升；在不同围压条件下，抗压强度的最小值，出现在β=30°的位置。其中，对于围压0 MPa的情况，在β=60°～90°，试件抗压强度变化较为平缓；当围压增加时，在β=60°～75°，试件抗压强度的增加逐渐明显。

图5 不同围压条件下平行柱轴方向玄武岩柱的抗压强度和等效变形模量

在等效变形模量方面，对于围压0 MPa～8 MPa的情况，玄武岩柱试件的等效变形模量随柱体倾角的增加大致呈先减小然后增加再微减的趋势；当围压超过4 MPa，各柱体倾角的试件等效变形模量，变化较不明显；在不同围压条件下，等效变形模量的最小值，出现在β=60°的位置。

由图6可知，在应力应变曲线方面，对于围压0 MPa的情况，存在明显的应力应变曲线分类现象，即β=0°、60°～90°的应力应变曲线，与β=15°～45°的应力应变曲线，区别较明显。其中，β=0°的应力应变曲线的峰前阶段，与β=60°～90°的峰前阶段，存在差异，在跌落阶段大致重叠。β=15°～45°的峰前阶段，大致重叠，而在跌落阶段存在差异。

图6 不同围压条件下平行柱轴方向玄武岩柱的应力应变曲线

对于围压2 MPa的情况，相比其它柱体倾角的试件，β=45°的试件，呈现较明显的延性破坏特征，且其残余强度最大。相比围压0 MPa的情况，围压2 MPa情况下，β=30°与β=60°的残余强度，较为接近。

对于围压4 MPa、6 MPa的情况，β=30°的试件，其延性破坏特征较明显，且其残余强度最大。β=0°与β=90°的应力应变曲线的跌落阶段重叠或接近重叠，而与其它柱体倾角的跌落阶段存在明显差异。

对于围压8 MPa的情况，β=30°延性破坏特征仍较明显，其残余强度仍最大。相比其它围压情况，围压8 MPa的情况下，β=15°、45°、60°、75°，其在应力应变曲线的跌落阶段，发生失稳破坏，而没有呈现之后的残余强度阶段。

2.2 单轴压缩及围压条件下不同柱体倾角玄武岩柱的应力场演化及声发射特征

2.2.1 单轴压缩条件下柱体倾角β=15°玄武岩柱的应力场演化及声发射特征

结合图7和图8可知，单轴压缩(围压0 MPa)，β=15°，3 m×3 m试件，随着加载的进行，当应力达到应力应变曲线的A点时，玄武岩柱试件内的柱状节理、试件的中上部、下部左右两侧，出现明显的应力集中。当应力加载至B点，试件上部中间的柱状节理，出现开裂，其附近的柱体应力集中较明显。当应力达到峰值点附近的C点时，试件内的柱状节理进一步开裂，此时应力集中主要出现在试件的顶部和底部。当应力降至D点，在试件的顶部和底部，若干柱体边缘，出现应力集中。当应力继续降至E点，在试件的顶部和底部，裂纹萌生发育，在试件的中部，若干柱体边缘，出现应力集中。当应力达到F点，在试件的顶部、中部、底部，若干柱体边缘，裂纹萌生、扩展。在声发射方面，其声发射呈单峰型分布。

图7 单轴压缩(围压0 MPa)，β=15°，3 m×3 m试件的应力应变曲线、声发射特征

图8 单轴压缩(围压0 MPa)，β=15°，3 m×3 m试件的应力场演化

2.2.2 围压2 MPa条件下柱体倾角β=15°玄武岩柱的应力场演化及声发射特征

结合图9和图10可知，围压2 MPa，β=15°，3 m×3 m试件，随着加载的进行，当应力达到应力应变曲线的A点时，试件内的柱状节理，出现较明显的应力集中。当应力达到峰值点附近的B点时，在试件顶部附近的柱状节理，出现开裂。当应力降至C点，试件顶部附近的柱状节理，进一步开裂，且试件顶部附近的柱体，应力集中较明显。当应力继续降至D点，在试件顶部和中部，若干柱体边缘，裂纹萌生发育。当应力达到E点，各柱体裂纹扩展，裂纹尖端处应力集中明显。当应力达到F点，裂纹进一步扩展，破碎区域主要分布在试件的顶部和中部。在声发射方面，其声发射呈单峰型分布。

图9 围压2 MPa，β=15°，3 m×3 m试件的应力应变曲线、声发射特征

图10 围压2 MPa，β=15°，3 m×3 m试件的应力场演化

2.2.3 围压8 MPa条件下柱体倾角β=15°玄武岩柱的应力场演化及声发射特征

结合图11和图12可知，围压8 MPa，β=15°，3 m×3 m试件，随着加载的进行，当应力达到应力应变曲线的A点时，试件内的柱状节理应力集中较微弱。当应力达到峰值点附近的B点时，试件上部和下部的柱状节理应力集中逐渐明显。当应力降至C点，试件顶部出现较明显的应力集中，裂纹萌生发育。当应力继续降至D点，试件顶部的裂纹向下扩展，此外，若干柱体边缘处，出现应力集中。当应力进一步降至E点，试件上部的多个柱体的局部发生破裂，破碎区域相连接，形成一条破碎带。当应力降至F点，原有的破碎带破碎加剧。在声发射方面，其声发射呈单峰型分布。

图11 围压8 MPa，β=15°，3 m×3 m试件的应力应变曲线、声发射特征

图12 围压8 MPa，β=15°，3 m×3 m试件的应力场演化

2.2.4 单轴压缩条件下柱体倾角β=45°玄武岩柱的应力场演化及声发射特征

结合图13和图14可知，单轴压缩(围压0 MPa)，β=45°，3 m×3 m试件，随着加载的进行，当应力达到应力应变曲线的A、B点时，玄武岩柱试件内的柱状节理、试件顶部，应力集中较明显。当应力达到峰值点附近的C点时，试件内的柱状节理压剪、拉伸、滑移。当应力降至D点，在试件内的若干柱体边缘处，有明显的应力集中。当应力继续降至E点，在试件内有两条条带状应力集中区域，在其区域内，若干柱体边缘，裂纹萌生、扩展，裂纹尖端，应力集中。当应力达到F点，各柱体的裂纹进一步扩展。在声发射方面，其声发射大致呈单峰型分布。

图13 单轴压缩(围压0 MPa)，β=45°，3 m×3 m试件的应力应变曲线、声发射特征

2.2.5 围压2 MPa条件下柱体倾角β=45°玄武岩柱的应力场演化及声发射特征

结合图15和图16可知，围压2 MPa，β=45°，3 m×3 m试件，随着加载的进行，当应力达到应力应变曲线的A点时，试件内的柱状节理出现明显的应力集中。当应力达到峰值点附近的B点时，在试件顶部，亦出现应力集中。当应力降至C点，试件内出现明显的条带状应力集中区域。当应力继续降至D点，原有的条带状区域分成两条条带状应力集中区域，其中，在下方的条带状应力集中区域，若干柱体边缘处，裂纹萌生发育。当应力达到E点，各柱体的裂纹扩展，裂纹尖端处，应力集中。当应力达到F点，裂纹进一步扩展。在声发射方面，其声发射呈类似三峰型分布的趋势。

图16 围压2 MPa，β=45°，3 m×3 m试件的应力场演化

2.2.6 围压8 MPa条件下柱体倾角β=45°玄武岩柱的应力场演化及声发射特征

结合图17和图18可知，平行柱轴方向，围压8 MPa，β=45°，3 m×3 m试件，随着加载的进行，当应力达到应力应变曲线的A点时，试件内柱状节理呈现微弱的应力集中。当应力达到峰值点附近的B点时，柱状节理的应力集中较明显。当应力降至C点，试件顶部亦出现应力集中情况。当应力继续降至D点，在试件内，形成一条明显的条带状应力集中区域，此外，试件顶部，裂纹萌生发育。当应力进一步降至E点，原有的条带状应力集中区域，有分化成两条条带状应力集中区域的趋势，在其区域内，若干柱体边缘处，裂纹萌生扩展。当应力降至F点，裂纹进一步扩展，试件内形成两条条带状破碎区域，破碎加剧。在声发射方面，其声发射呈双峰型分布。

图17 围压8 MPa，β=45°，3 m×3 m试件的应力应变曲线、声发射特征

图18 围压8 MPa，β=45°，3 m×3 m试件的应力场演化

2.3 围压8 MPa条件下不同柱体倾角玄武岩柱的损伤破裂特征分析

由图19可知，围压8 MPa，柱体倾角β=0°，试件上部的柱状节理开裂，此外，试件上部有若干条条带状损伤破裂区域。当β=15°，柱状节理发生损伤，试件顶部若干个柱体损伤破裂，试件中上部有一条损伤破碎区域。当β=30°，柱状节理损伤较明显，试件顶部若干个柱体损伤破裂，试件中上部有一条损伤破碎区域，相比β=15°的情况，β=30°的条带状损伤破碎区域较细。当β=45°，试件内形成两条明显的条带状损伤破碎区域。当β=60°，试件内形成两条明显的条带状损伤破碎区域，其中，下方的条带状损伤破碎区域较宽。当β=75°、90°，试件顶部有若干条条带状损伤破碎区域。

图19 围压8 MPa，β=0°～90°，3 m×3 m试件的损伤破裂特征

3 结 论

(1) 不同围压条件下，不同柱体倾角的玄武岩柱强度的各向异性。对于围压0 MPa～8 MPa的情况，玄武岩柱试件的抗压强度随柱体倾角的增加大致呈U型分布； 随着围压的增加， 各柱体倾角的试件抗压强度，有较明显的提升；在不同围压条件下，抗压强度的最小值，出现在β=30°的位置。其中，对于围压0 MPa的情况，在β=60°～90°，试件抗压强度变化较为平缓；当围压增加时，在β=60°～75°，试件抗压强度的增加逐渐明显。在等效变形模量方面，在不同围压条件下，等效变形模量的最小值，出现在β=60°的位置。

(2) 以柱体倾角β=15°的玄武岩柱为例，研究其在不同围压条件下的破裂机理及破坏模式：① 单轴压缩条件下(围压为0 MPa)，加载的初始阶段，玄武岩柱试件内的柱状节理、试件的中上部、下部左右两侧，出现明显的应力集中；随着加载的进行，柱状节理逐渐开裂，若干柱体边缘应力集中，裂纹萌生发育；继续加载，在试件的顶部、中部、底部，若干柱体边缘，裂纹萌生、扩展；② 当围压为8 MPa，加载的初始阶段，试件内的柱状节理应力集中由较微弱发展为较明显；随着加载的进行，试件顶部出现较明显的应力集中，裂纹萌生发育；继续加载，试件上部的多个柱体的局部发生破裂，破碎区域相连接，形成一条破碎带。

(3) 以柱体倾角β=45°的玄武岩柱为例，研究其在不同围压条件下的破裂机理及破坏模式：① 单轴压缩条件下(围压为0 MPa)，加载的初始阶段，玄武岩柱试件内的柱状节理、试件顶部，应力集中较明显；随着加载的进行，试件内的柱状节理压剪、拉伸、滑移，若干柱体边缘处，有明显的应力集中；继续加载，试件内有两条条带状应力集中区域，在其区域内，若干柱体边缘，裂纹萌生、扩展；② 当围压为8 MPa，加载的初始阶段，柱状节理的应力集中由微弱发展为较明显；随着加载的进行，在试件内，形成一条明显的条带状应力集中区域；继续加载，各柱体裂纹萌生扩展，试件内形成两条条带状破碎区域，破碎加剧。