Study on Surrounding Rock Deformation and Slope Stability in AiGong Cave Based on FLAC3D
刘勇兵① LIU Yong-bing；杨志法② YANG Zhi-fa；
魏雪云① WEI Xue-yun；王婷静① WANG Ting-jing
（①衢州学院建筑工程学院，衢州 324000；②中国科学院地质与地球物理研究所，北京 100029）
（①College of Civil Engineering and Architecture，Quzhou University，Quzhou 324000，China；
②Institute of Geology and Geophysics，Chinese Academy of Sciences，Beijing 100029，China）
摘要：本文主要借助FLAC3D软件建立了龙游县艾公洞古洞室三维模型，针对洞室模型下的应力和位移两方面进行围岩变形及其边坡的稳定性数值分析。通过研究发现，洞室顶板底中部产生的拉应力远小于岩石的抗拉强度，证实了顶板的拉弯受力满足承载力的安全要求，洞室围岩的破坏主要是在重力作用下的剪切错动破坏，且具有长期缓慢发展的趋势，本研究可为该洞室边坡的加固提供重要的理论和现实指导意义。
Abstract: In this paper, the 3D model of the ancient AiGong Cave is established with FLAC3D software, so numerical analysis of surrounding rock deformation and its slope stability is carried out in terms of stress and displacement. Through the research, it is found that the tensile stress generated in the middle part of cave roof is much smaller than the tensile strength of the rock, which confirms that the tensile strength of the cave roof meets the safety requirements of the bearing capacity. But deformation failure of the surrounding rock is possible and is mainly the shearing fault under gravity and has a long-term slow development trend. This study provides important theoretical and practical guide for the reinforcement of this cave.
关键词：古洞室；围岩；变形；数值模拟；边坡稳定性
Key words: ancient cave；surrounding rock；deformation；numerical simulation；slope stability
中图分类号：TD32 文献标识码：A 文章编号：1006-4311（2019）25-0226-06
0 引言
龙游艾公洞毗邻龙游县小南海镇定埠村附近的古衢江边，与濲波岩、翠光岩摩刻、龙游石窟等，一起被称为小南海石室群，均被列为国家级重点文物保护单位。艾公洞与濲波岩之间还有一个较小的无名洞。濲波岩、无名洞和艾公洞三者之间大致相对位置如图1所示。据龙游县文化广电新闻出版局、博物馆副馆长朱土生的考证，艾公洞是由康熙年间（1662~1722年）知县艾朝东主持修建，在当时主要是文人欣赏衢江美景和吟诗作画之所[1]。因艾公洞上方不远处存在一幢四层兵工厂建筑，其围岩的工程地质条件也较差、水文地质条件的不断变化以及岩体的风化剥蚀等，使艾公洞围岩及其所在边坡发生了多处开裂，局部脱落、掉块、局部崩塌等现象。但是其边坡的受力比较复杂，影响其稳定的关键因素未知，为了古艾公洞的安全和长期稳定，作为国家重点文物保护对象，现急需对艾公洞洞室边坡稳定性进行研究。
1 工程概况
经实地查勘，艾公洞洞室上方不远处于1980年初建有一幢四层（层高均为3.3m）兵工厂平屋面建筑，该建筑距离洞室最近一边北端墙约7m。洞室南端洞口至衢江主河槽间已淤积了大量的泥沙，测得洞口距现今衢江边水的距离约为47.4m。洞室所处地层主要为白垩系上白垩统衢县组（K2q）的砖红色泥质粉砂岩[2]，且洞室内部西面出露有一层粉砂质泥岩（图2），厚约0.3~0.45m，位于离洞室底面1.8~2.45m，并向东面延伸了约3.0m距离。洞室前方（南面）为一陡坡地形，坡角大致在60°~70°，该陡坡地形自洞室两边沿纵向各延伸约300多米以上。洞室顶部（坡顶面）高出淤积的河滩表面近9.0m，洞室内上缘至陡坡面最短距离约1.0m，顶板厚度约为1.3~2.0m。
现洞室内部两侧面（东、西面）均出现了多条平行的斜向剪切裂缝，两侧面发展非常明显的各有6条，并且两侧裂缝近似平行地向坡脚开展延伸，并且东边墙距离洞室附近发展分布有一条较为主要的、长度最长的剪切-开裂面。而对于西边墙来说（图2），5条裂缝均只延伸至二层方梁孔位置，这或许是由于泥质夹层起到了缓冲分散应力集中的作用，且夹层的倾角约为21°，对于边坡来说是比较有利的。洞室各平面图和剖面详图见图3。
2 数值模型、计算参数与边界条件
2.1 数值模型
结合龙游艾公洞的实际破坏情况和变形特点，现采用FLAC3D软件对艾公洞古洞室围岩及其所在边坡进行数值分析，对其稳定性作出评价。经实地测得洞室和兵工厂建筑所处在边坡的尺寸和形状、地层及结构面等分布情况，在建模时，考虑洞室上方不远处的兵工厂建筑物荷载对洞室及其所在边坡的影响。三维数值计算模型坐标采用沿洞室平面南北方向为y轴，朝北为正，其南北方向总长度为28.25m；沿洞室平面东西方向为x轴，朝东为正，其东西方向总宽度为22.1m，垂直地表水平面方向取作z轴，向上为正，沿z轴向下取至平地表以下16.98m，各方向尺寸也均保证大于相应洞室方向尺寸的3倍。将洞室南北长度范围及向建筑物方向2.4m长范围内的单元全部进行加密分布[3]。计算模型共划分有468399个单元，491523个节点，单元多采用六面体单元（如图4）。
2.2 计算参数
经试验分析，模拟计算时采用的岩土体参数如表1所示。
2.3 边界条件
由于洞室前方淤积有一层厚厚的粉砂层，沿洞室前方延伸达四十多米至衢江河槽，因此模型的南面以河滩表面为界，河滩表面以上为自由边界，以下为单向水平方向固定。东西两面及北面（兵工厂建筑物方向）边界均为单向水平方向固定，底边界（zz方向）为洞室底面以下8.4m处为固定边界。上边界为地表，取作自由边界，将兵工厂建筑物荷载转化为均布压力，施加在建筑物作用位置处。
3 计算结果与实际破坏情况比对分析
通过数值模拟计算分析后，主要对模型下的应力和位移两个方面进行分析。关于应力方面，主要针对模型下顶板四角、中部及两边各裂缝上端点的应力大小等作分析研究；关于位移方面，则主要针对模型下的顶板四角、顶板中部、两边各裂缝上端点的位移等进行分析。
3.1 洞室应力场分析
据现场实地分析，龙游石窟古洞室破坏多数是因剪切滑移、受力不足导致的，且以沿着破裂面大面积坍塌为主。
3.1.1 洞室顶板板底四角点与中部点单元主应力分析
考虑到艾公洞洞室顶板底部的中央部分承受较大拉应力，属于较易威胁洞室稳定的重点部位之一，加上顶板为水平分布，上面生长有一定的小灌木，其根系可能在其洞室顶板表面产生沿着洞室边墙的垂直裂缝，因此也可能会导致顶板沿着洞室四周边墙的塌落破坏[4]。故在洞室顶板板底的四角、中部部位附近设置相应的监测点，求得到这些部位的单元最小主应力σmin和最大主应力σmax进行分析，其结果如表2及图5~图6所示。
②所建立的整体坐标系取洞室北端墙底部中点为坐标原点，自西向东为x轴正向，自南朝北为y轴正向，自洞室底至地表面为z轴正向。
从表2的计算结果表明：洞室模型顶板板底四角单元主应力为压应力，且靠北一端的值均较靠南一端的值大，其原因主要是洞口部位岩层厚度较薄，而靠北一端除了岩层较厚外，其距离兵工厂建筑物较近也有一定关系，因此墙角处不容易导致破坏，现场也没有发现任何裂缝。而在洞室模型顶板板底中部单元单元最小主应力σmin和最大主应力σmax均表现为拉应力，主要顶板设计成了水平方向，其东西横向跨度为5.3m，南北纵向为7.25m，在表层黏土和顶板岩石自重作用下，顶板底部承受较大的弯矩，在顶板底中部产生拉应力。但最大拉应力8.63×104Pa也要远小于岩石的抗拉强度1.54×106Pa，这也正是洞室顶板板底至今没有发现裂缝的原因。模型在设定求解5000时间步长后的单元体区域最终得到的主应力等值线图分别如图7~图8所示。从各图中可以明显看出模型可能出现破坏区域均是指向洞口方向斜向增加，且兵工厂建筑对洞室的受力还是有一定的影响。
3.1.2 洞室东、西边墙裂缝上端点处的单元主应力分析
考虑到艾公洞洞室两边均已发育有多条斜向裂缝，各裂缝的分布长度、开裂宽度都不一样，而且东边墙裂缝分布长度明显大于西边墙，但两边裂缝大致均近似平行斜向边坡临空面，且均形成了三条主要裂缝。裂缝⑥离北端墙最近且发展最小，在此不作为分析对象。现分析各主要裂缝后期发展情况，以便作出相关预测与加固处理，其结果如表3。
从表3分析结果表明，洞室模型东、西边墙裂缝上端的单元应力基本上相当，且均表现为以压为主，这说明边坡的破坏主要是在重力作用下的剪切错动破坏[5]。但由于西边墙存在泥质夹层，起到了过度缓冲作用[6]，致使裂缝未能向洞底面延伸。
3.2 洞室位移场分析
3.2.1 洞室顶板板底四角点与中部点位移分析
为说明顶板的沉降情况，以及探明顶板是否会发生薄层剥落，考虑到顶板可能发生塌顶的可能，忽略x向与y向变形对z向的影响，现对洞室顶板板底四角处与中间点位置进行z方向的位移观测，以便对求得相互之间的沉降差值作出定性判断。洞室顶板板底四角点与中部点处的位移?驻z计算结果如表4及图9、图10所示。
从表4及图9、图10计算结果表明：洞室模型北端墙的沉降量要大于南端墙的，其原因位兵工厂建筑作用的结果。而对于洞室西南墙角的沉降量要大于东南墙角的，这可能是西边存在一层较薄的泥岩夹层，其本身变形压缩量要大一点，抵消了总沉降量的一部分，边坡发生剪切错动时，导致裂缝发展长度较东边要短。且洞室模型顶板中部的沉降量最大，顶板中部点相对于四角点平均产生的扰度为1.45×10-4/5.3=1/36552，由此可见实际产生的挠度非常小，对洞室顶板不会造成多大影响。
3.2.2 洞室东、西边墙裂缝上端点处的位移分析
表5及图11、图12给出了东、西边墙裂缝上端点处的z方向、y方向位移计算结果。图13和图14所示为洞室模型z方向和y方向位移等值线图。
表5及图11~图12可以看出：洞室东、西边的各裂缝上端点处z方向位移值基本接近，沉降位移基本稳定，由于兵工厂建筑影响靠北端墙位置稍微大些。而各裂缝y方向（朝洞室南面）的水平位移值靠洞口越近位移越大，且各y方向位移值均具有不断发展的趋势，说明了洞室的破坏主要是沿着陡坡的破裂面长期缓慢滑动为主。在长时间发展变化中，洞室将有可能沿着陡坡的破裂面最终产生剪切错动破坏。
4 结论
借助FLAC3D软件三维建模，对经历几百年变化之后，考虑洞室后方不远处的兵工厂建筑的艾公洞古洞室进行了数值模拟分析。通过分析计算，对洞室顶板四角、中部及两边各裂缝上端点的应力和位移进行分析与对比发现：
①应力方面，洞室模型顶板板底四角单元主应力以压应力为主，拉应力相对很小，且靠北一端的值均较靠南一端的值大，说明兵工厂建筑对洞室的受力是有一定的影响；模型顶板底中部产生的拉应力远小于岩石的抗拉强度，说明洞室顶板板底至今没有发现裂缝，这就说明边坡的破坏主要是在重力作用下的剪切错动破坏。
②位移方面，顶板中部点相对于四角点产生的实际挠度非常小，对洞室顶板不会造成多大影响；各裂缝上端朝洞室南面方向的位移值靠洞口越近位移越大，且具有不断增大的趋势，说明洞室的破坏主要是沿着陡坡的破裂面长期缓慢剪切错动为主。
参考文献：
[1]杨志法，岳中琦，李丽慧.龙游石窟大型古地下工程洞室群科学技术问题研究[M].北京：科学出版社，2010.
[2]浙江省地质矿产局.浙江省岩石地层[M].武汉：中国地质大学出版社，1996.
[3]孙书伟，林杭，任连伟.FLAC3D在岩土工程中的应用[M]. 北京：中国水利水电出版社，2011.
[4]刘传正.崩塌滑坡灾害风险识别方法初步研究[J].工程地质学报，2019，27（1）：88-89.
[5]兰恒星，仉义星，伍宇明.岩体结构效应与长远程滑坡动力学[J].工程地质学报，2019，27（1）：108-122.
[6]周晓飞，孙金山，刘贵应.基于幂函数的边坡岩体泥质夹层长期强度直剪蠕变试验研究[J].安全与环境工程，2018，25（1）：6-11. |
