摘要：以延安市宝塔区杨家岭毛泽东旧居旁的靠崖式单孔窑洞为研究对象，依据当地黄土的力学特性，按弹塑性理论，使用MidasGTSNX岩土有限元分析软件，对靠崖式单孔黄土窑洞的静力性能进行分析，以揭示窑洞的破坏形式和规律并提出加固窑洞的相应方法和措施。
关键词：黄土窑洞；静力性能；有限元分析；加固方法
延安位于陕西北部，处在典型的黄土分布区域。当地黄土覆盖深厚稳固，经过100多万年的沉积，局部黄土层厚度达到二三百米以上，有的地方呈现出陡峭的黄土直壁，加之气候干旱少雨，当地居民便利用这些自然资源在土崖中开挖洞穴作为居室[1-2],“窑洞”这种独特的建筑民居形式由此而生。传统的黄土窑洞具有冬暖夏凉、施工简单、节约材料、成本低廉等特点，在中国经济困难时期是当地民众主要的住宅形式。随着经济社会的发展和人民生活水平的提高，窑洞所存在的采光差、通风不畅、交通不便、基础设施不完善等缺点，已无法满足现代生活的需要。2013年7月，延安经历了一场百年不遇的持续强降雨，致使十几万间窑洞垮塌或损毁，在干旱条件下良好运转多年的依山而建的窑洞，在连续降雨中暴露出严重的安全缺陷。为确保民众安全，“弃窑进楼”在政府的引导下正有序进行，窑洞面临着被遗弃的困境[3].
然而，延安窑洞是地域文化的体现和生土建筑的典型代表，也是中国革命的见证。从传承窑洞建筑文化和政治经济角度来看，保护现存窑洞建筑具有重要的现实意义。要提高现有窑洞的抵御灾害能力，必须研究窑洞在建造施工阶段的静力性能，为进行合理的加固保护提供可靠的理论依据。
陕北地区窑洞按照修筑方式不同可分为靠崖式土窑（直接从黄土崖面向内挖掘的土窑）、明箍窑（直接在平地上箍造的砖窑和石窑，上部覆盖一定厚度的土）、接口窑（在土窑口接上一段砖石拱圈）[4].目前对于黄土窑洞静力方面研究已取得不少成果[5~17],然而研究者通常将实际模型简化为二维平面应变问题或正方体三维模型进行分析，计算模型与实际有所出入。为更逼近实际情况，获得更精确的解答和更深入的研究，本文以延安市宝塔区杨家岭毛泽东旧居旁的靠崖式单孔窑洞为研究对象，依据当地黄土的力学特性，按弹塑性理论，使用MidasGTSNX岩土有限元分析软件，建立单孔窑洞三维立体模型（包含窑口前方岩体），模拟黄土窑洞的开挖过程，分析窑洞的位移大小和塑性区的变化趋势，探讨窑洞的破坏形式和规律，依据分析结果，提出加固窑洞的相应方法和措施。
1Midas软件简介
MidasGTSNX是一款利用尖端计算机图形处理和分析技术研发而成的针对岩土领域的通用有限元分析软件。它可以进行静力分析、施工阶段分析、动力分析、边坡稳定、渗流和固结等类型分析，提供给用户多种单元类型和材料本构模型，能真实地模拟施工阶段开挖过程、地基应力状态（自应力和构造应力）、围岩状态下各种结构作用以及岩土的非线性特点，以便得到比较真实的结果。
使用MidasGTSNX软件分析问题的步骤与常用有限元软件基本类似，即：建立模型→定义材料属性→网格划分→施加荷载→边界约束→定义分析类型（线性静态、施工阶段、时程等）→计算分析→查看结果。MidasGTSNX软件与常用软件ANSYS、FLAC3D等相比，有简单易操作、计算周期短等优点。
2靠崖式生土窑洞静力性能
2.1窑洞模型的建立
靠崖式单孔黄土窑洞剖面尺寸示意图见图1.对延安市宝塔区杨家岭毛泽东旧居旁的单孔窑洞（图2）进行实地测量，测量的数据列于表1.
靠崖式窑洞一般是将自然的黄土崖面铲削成75°~85°的坡度，然后在崖壁上凿挖成一定高度、宽度和进深的土洞。为了保持窑体的稳定性，窑洞的顶部拱形一般为半圆或尖圆。这里用MidasGTSNX软件建模时，崖面坡度采用80°，窑洞拱形为单心圆弧。有限元法处理无限域或半无限域问题通常是在有限的区域内进行离散。为了使有限元的离散误差和计算结果满足一定精度要求，又要保证边界条件尽可能接近实际情况，因此必须划定周围岩土合理的计算范围，一般选取的计算区域是沿洞径各个方向均不小于3~4倍洞径[16~17].本文计算模型中的周边土体范围取相应方向跨径的3倍，故确定的土体外围尺寸为宽、高各20m,窑口前部土体深度为20m,后部深度为27m.计算模型见图3.
2.2材料参数的确定和材料属性的定义
因靠崖式黄土窑洞的稳定性取决于黄土的强度，一般均建造在力学性质较为稳定的Q2（离石黄土）和Q3（马兰黄土）黄土坡体中。由于所能收集到的陕北地区黄土的基本物理力学性质资料数据很多，而且同一窑洞模型在不同土质参数下的位移区和塑性区的变化规律和趋势大致相同，为减少计算工作量，本文参考文献[18-19]中的延安市黄土物理力学参数，在窑洞模拟分析中，选取这两种黄土物理性质和力学性质的平均值，确定的黄土材料参数取值列于表2.
这里的土体材料属性采用理想弹塑性本构关系，土体塑性屈服后服从Mohr-Coulomb破坏准则，该准则通常用于描述土体和岩石的剪切破坏，也是MidasGTSNX软件对一般岩体推荐使用的准则。
2.3网格的划分
为更细致地查看窑洞周围的受力情况，在划分网格时，先将窑洞周围的网格尺寸设置的小些，然后采用自动划分实体网格法，使得其余部分网格尺寸可适当大一些，以减少计算量。本文在窑口周围将线尺寸控制设置为0.2m,其余部分设置为1.0m,划分的网格如图4所示。
2.4模拟所用的单元类型、约束条件和荷载
在用MidasGTSNX进行有限元分析时，采用三维实体单元（SolidElement）模拟岩土地层。实体单元可利用四节点、六节点或八节点构成，本文采用四节点单元，即三角锥单元。三角锥单元（四节点单元）计算的位移结果比较准确，但是应力结果精确度较其他实体单元差。使用二维板单元来模拟窑洞衬砌，板单元可以考虑平面受拉、受压、平面受剪、平面外受弯及厚度方向的剪切。
由于计算模型中的土体边界距离窑洞较远，根据圣维南原理，可认为边界处无位移发生。因而模型中土体的两个侧面的水平位移即x向位移为0,正面和背面土体y向位移为0,土体底面z向位移为0,窑口所在平面以及土体顶面为自由平面，这样的边界条件与实际相符。窑洞除承受周围土体自重外，并无其他外荷载，故只考虑自重荷载。
2.5开挖过程模拟
窑洞在开挖过程中是逐步向前掘进的，而且边挖边用木板、木桩加以支护。为了模拟窑洞开挖的动态过程，本文采用设置不同的施工阶段来模拟。对窑洞模型的开挖分为6步进行，每次掘进1.1m,即将开挖进深分别设为0~1.1、1.1~2.2、2.2~3.3、3.3~4.4、4.4~5.5、5.5~6.6m.在静力分析时，将分析工况设置为施工阶段分析，这是因为在MidasGTSNX中，施工阶段分析采用的是累加模型，即每个施工阶段都继承了上个施工阶段的分析结果，并累加了本施工阶段的分析结果。因此，先要确定土体在窑洞没有开挖前的初始应力及位移状态，然后计算窑洞在开挖之后的应力场及位移场。本文在定义施工阶段时，将初始阶段的应力状态和位移设置为零，开挖的每阶段中将增加相当于增加的单元自重荷载。
2.6结果
分析在运行分析前，首先定义施工阶段相关的分析工况，利用分析控制功能进行施工阶段分析的细部设定。这里不勾选应力分析初始阶段，确认最后计算阶段指定为“最后阶段”,然后点击运行分析。计算完成后，可以查看模型在开挖各个阶段的各个方向的位移、应力及应变云图，也可以提取关键位置的具体数值进行分析。
2.6.1位移分析 针对窑洞结构安全稳定问题，本文只考察水平x方向位移和竖直z方向位移。最终阶段的水平向位移云图和竖直向位移云图结果显示：水平方向的最大位移出现在窑洞左、右直壁偏下方处的位置，数值相近，方向内向；而最大竖直向下位移发生在拱顶处，窑洞底面中心部位有隆起位移。现将拱顶部位、洞室底部的最大竖向位移（z向）和窑洞左、右直壁部位的最大水平位移（x向）在不同开挖深度时的数值列于表3中。可见拱顶竖向位移和水平位移随着开挖进深的扩大有逐渐增大的趋势。
由于在不同开挖深度时拱顶的最大竖向位移均出现在同一节点处，而此节点即为窑口拱顶点，将此节点定义为1点，绘制出1点的竖向位移随开挖深度增大的变化规律（图5），可以看出窑洞在开挖初期，由于自重作用，拱顶下沉位移增加很快，故应在开挖初期做好支护工作。在窑口左右直臂距底部三分之一处的节点定义2点和3点，绘制出2、3点的水平位移随开挖深度增大的变化规律（图6）。可以看出，左直臂处发生向右的水平位移，右直臂处发生向左的水平位移，而且左右数值接近对称，开挖初始阶段位移增量很大，到后期趋于平稳。这是由于土体被挖开后，原有土体失去最初的平衡状态，使得两侧土体向窑洞内空界面发生移动。因此，开挖过程中应用木板、木桩等进行支护，以防止窑洞坍塌或土层掉块现象。
2.6.2塑性区分析 为观察窑洞在开挖过程中可能出现破坏的位置，应对塑性区的分布情况进行研究，因为破坏先是从出现塑性区的地方开始的。当最大等效塑性应变值为0时，认为没有进入塑性区；不0为时，则认为这一区域进入了塑性区，且数值越大，说明越不安全。通过对不同阶段等效塑性应变云图的观察和对比发现，最先出现塑性区的是洞室两侧直臂靠近底脚处，塑性区主要分布在窑洞直臂处及直臂与拱圈交接处附近区域内，有逐渐向拱顶延伸的趋势。最大等效塑性应变值随开挖深度的变化规律如图7所示，可以看出，随着开挖深度的加大，塑性区分布范围在不断扩大，最大等效塑性应变值也在逐渐增大，在前3个阶段增大趋势较快，后面增加逐渐趋于平缓。说明随着开挖的进行，窑洞周围土体从弹性阶段进入塑性阶段，但仍有一定的承载力而继续变形。
【图7】
3窑洞加固方法
从以上对靠崖式窑洞静力性能的分析可以看出，窑洞的薄弱部位是窑口、拱顶、直臂与拱圈交接及底脚处。在修建时应做好前期的支护工作，最好是边挖边支，防止窑口塌陷、底脚破坏；同时在使用中加强排水、防水和防潮工作，以免土体遇水或潮气使自重变大，窑顶出现裂缝或坍塌。对于现存靠崖式窑洞，因窑址的选择已不能改变，应增强山体的稳定，如加砌护坡或采用防护技术，防止土体下滑、崖壁崩塌等一系列地质灾害。同时还要做窑洞衬砌对原有窑洞进行加固[20],如做土坯衬砌、普通烧结砖衬砌、毛石衬砌等，还可以应用高强度复合材料进行抹面处理。
本文通过Midas软件中的二维板单元模拟衬砌，对比了做30cm厚的土坯衬砌、12cm厚的烧结砖衬砌、24cm厚的烧结砖衬砌及30cm厚的毛石衬砌时窑洞拱顶最大位移（表4）。有限元分析结果显示加设衬砌后，窑洞周围土体受到约束，拱顶竖向位移和塑性区范围均有减小，加固效果明显。若从经济角度考虑，用土坯做衬砌费用最少，用普通砖做衬砌费用次之，砖衬砌越厚，加固效果越好，而毛石衬砌费用则最贵。综合考虑可以优先选用24cm厚的普通砖作为衬砌材料。
4结语
本文利用岩土有限元数值分析软件，对陕北延安宝塔区杨家岭革命旧址中毛泽东旧居旁的靠崖式单孔窑洞进行三维立体模拟，依照弹塑性理论，分析了靠崖式生土窑洞的静力性能。结果显示，最大竖直向下位移发生在拱顶处，水平方向的最大位移出现在窑洞左、右直壁偏下方处的位置；竖向位移和水平位移随着开挖进深的扩大有逐渐增大的趋势，而且开挖初期竖向位移和水平位移均增长较快；塑性区主要分布在窑洞直臂及直臂与拱圈交接处附近区域内，有逐渐向拱顶延伸的趋势，而最先出现塑性区的是洞室两侧直臂靠近底脚处；最大等效塑性应变值随开挖深度增加而逐渐增大，并且前期增大趋势较快，后期趋于平缓。根据以上分析结果，判断出窑洞的破坏形态和薄弱部位，提出对现有窑洞采用保护窑洞崖体、窑内加设衬砌的加固方法。综合考虑加固效果及经济条件，认为加设24cm普通砖衬砌是加固窑洞的优选措施。
参考文献：
[1]王继唐。窑洞民居防灾[J].灾害学，1993,8（1）：86-90.
[2]杨志威。浅议我国窑洞建筑的现状与未来[J].西北建筑工程学院学报，1990（3）：130-140.
[3]郝艳娥，张扬，杨红霞。延安市革命旧址中窑洞的保护问题分析[J].建筑设计管理，2015（12）：56-58.
[4]张厶文，郝艳娥，李奇，等。陕北窑洞结构有限元分析及加固方法探讨[J].建材发展导向，2014（24）：65-67.
[5]任侠。黄土窑洞弹塑性有限元分析[J].兰州铁道学院学报，1993,12（4）：1-7.
[6]闫月梅。石砌窑洞拱圈的受力分析和截面计算[J].工程力学，1996（增刊）：663-666.
[7]张玉香，田兴运。窑洞破坏的原因及其防治措施[J].西北农林科技大学学报：自然科学版（增刊），2004,32:145-149.
[8]童丽萍，张晓萍。生土窑居的存在价值探讨[J].建筑科学，2007,23（12）：7-9.
[9]童丽萍，韩翠萍。传统生土窑洞的土拱结构体系[J].施工技术，2008,37（6）：113-116.
[10]童丽萍。传统生土窑居的灾害及民间防灾营造[J].建筑科学，2008,24（12）：7-9.
[11]童丽萍，韩翠萍。黄土窑居自支撑结构体系的研究[J].四川建筑科学研究，2009,35（2）：71-73.
[12]黄才华，王泽军。窑洞建筑的结构分析[J].长春工程学院学报：自然科学版，2009,10（1）：46-48.
[13]马琳瑜。中国民居的拱券结构研究[D].西安：西安建筑科技大学，2007.
[14]姬栋宇，谢湘赞。弹塑性力学理论在靠崖窑结构中的应用[J].中国建设教育，2012,4（4）：35-37.
[15]王玺。考虑节理影响的黄土窑洞稳定分析[D].西安：西安建筑科技大学，2009.
[16]胡晓坤。甘肃中东部黄土窑洞的静力特性研究[D].兰州：兰州理工大学，2011.
[17]陈莉粉。黄土地区窑洞建筑中结构稳定性的研究[D].西安：西安科技大学，2012.
[18]张乘千。延安市降雨与黄土滑坡相关性分析[D].西安：长安大学，2014.
[19]孟龙。黄土窑洞灾害形成机理数值模拟研究[D].西安：长安大学，2011.
[20]许方义。用ADINA软件对庆阳地区黄土窑洞的模拟计算与分析[D].兰州：兰州理工大学，2010.