摘要: 在目前的城市公路建设过程中，越来越多的公路的建设方式采用高架桥的方式进行，这种建设方式可以更好地利用地面以上的空间，在资源型城市中建设高架桥往往会遇到采空区问题。本文以某资源型城市在某公路快速化改造过程中遇到的采空区问题为背景，以FLAC3D模拟软件为工具，研究该项目在采空区变形过程中的规律以及原治理方案即注浆方案的合理性，对比注浆前后采空区区域内的地表变形、应力、以及塑性区的变化规律，为类似的采空区区域内建设高架桥的工程提供建设和治理的建议。
Abstract:
In the current process of urban highway construction, more and more highways are constructed in the form of viaducts, which can make better use of the space above the ground. The construction of viaducts in resource-based cities often encounters goaf problems. This paper takes the goaf problems encountered in a resource-based city during the rapid transformation of a highway as the background, and uses FLAC3D simulation software as a tool, studies the regularity of the project in the deformation process of the goaf and the rationality of the original treatment scheme, namely the grouting scheme, and compares the variation regularity of the surface deformation, stress, and plastic zone in the goaf area before and after the grouting, so as to provid
1. 引言
随着现代社会的发展，我们每年都要从地球内部开采出数量巨大的矿产品，采空区就是这类开采活动带来的破坏之一。采空区是指地下矿产被人为挖掘而遗留的空洞，由于开采过程中，土体初始平衡应力遭到破坏，致使采空区产生逐渐变形，其上方土体缓慢沉降，很可能导致采空区上方建筑物由于地基失稳而破坏，进而威胁人们的生命财产安全。我国不仅幅员辽阔，矿产资源储量丰富，并且分布广泛。以产煤大省山西省为例，该省全省面积的七分之一(约2万多公里)已采空，部分地面塌陷严重，并且沉陷区面积正以每年94 km2的速度较快增长 [1] [2] 。
一般来讲，工程建设在勘察设计、选址环节中遇到采空区区域时，都会遵循避让的原则，但有时碍于线路的要求和建设成本等因素的考虑，工程项目有时会无法避免地选址在采空区区域 [3] 。
早在19世纪，国外就已经有人对矿产开采带来的地面影响进行了研究。1838年，比利时的工程师哥诺特 [4] 提出了关于开采沉陷的第一个理论——“垂线理论”，随后有了“法线理论”与“拱形理论”。再后来，在1947年，前苏联学者阿维尔申将塑性理论引入到开采沉陷的研究中，1954年，波兰学者李特维尼申提出了开采沉陷的随机介质理论，1961年，Brauner提出了水平移动影响函数 [5] 。Berry和Sales分别在1961年和1963年对开采沉陷区地表位移进行二维研究，获得了二维的地表的下沉和水平变形规律 [6] [7] 。到了70年代，Jone S等人 [8] 研究分析了因采矿造成的塌落对公路的影响；紧接着80年代，Jone S [9] 、Drum [10] 、Wang M.C. [11] 等人对采矿过程及因开采矿产资源而形成的采空空洞进行了研究，并分析了其对上部建筑物地基的影响和可能造成的危害。上个世纪90年代以后，由于科学技术水平的不断发展和进步，开采沉陷学的研究进入到突飞猛进的阶段。一方面，人们运用计算机技术操作复杂的算法，使得过去无法计算的问题成为可能；另一方面，一些新的数值分析计算方法，包括有限元法、离散元法和边界元法等，均应用到了开采沉陷的计算中。张杰 [12] 利用相似模拟实验、理论分析、数值模拟等手段对采空区覆岩结构及失稳特征进行了研究，得到了煤柱和上覆地层的失稳规律以及采动裂隙的失稳规律。刘垚鑫 [13] 等通过理论推导得出了倾斜采空区的内应力场表达式，并且给出了内应力场宽度的表达式。并通过数值模拟，研究了倾斜特厚煤层条件下，沿空掘巷时煤柱宽度对塑性区分布与垂直应力分布的影响，综合给出了最佳煤柱宽度。近几年，模型试验搭配耦合的数值计算方法成为研究工程问题的热门方法 [14] 。
2. 工程概况
2.1. 地形地貌
项目区地貌单元为古黄河泛滥形成的冲洪积平原，地势较平坦。地面标高约35~44.0 m，项目区附近由于煤炭开采活动，使地表形成大面积塌陷并积水成塘，塌陷区水深可达5~6 m。
2.2. 地质条件
项目区地处冲(洪)积平原地带，地表被第四系土层覆盖，第四系下伏基岩主要由石炭系地层构成。改造工程下伏基岩地层由老到新为奥陶系马家沟组～阁庄组、石炭系本溪组～太原组、二叠系下统山西组～下石盒子组、二叠系上统上石盒子组，其中奥陶系以灰岩、白云岩为主，石炭系本溪组～太原组地层岩性以灰岩、碎屑岩夹煤层为主，二叠系下统山西组～上石盒子组以碎屑岩夹煤层为主，均被第四系所覆盖。项目区7煤上覆基岩地层主要为二叠系下统山西组~下石盒子组以碎屑岩夹煤层为主，各组地层工程地质性质如下：二叠系下统山西组~下石盒子组属碎屑岩软质岩石亚类，岩性主要为砂岩、页岩夹煤层，表层风化呈粘土状。砂岩岩芯呈灰白色，泥质、钙质胶结，层状构造，层理发育，泥质充填，岩芯呈长柱状、短长柱状，完整砂岩单轴抗压强度平均值41.77 MPa；砂质页岩岩芯呈灰色、灰黑色，强～中风化，泥质胶结，含砂质，层状构造，局部为砂岩页岩互层，层理发育，微解理裂隙发育，泥质充填，岩芯呈碎块状、短柱状、长柱状，完整页岩单轴抗压强度平均值3.69 MPa；页岩岩芯呈灰色、灰黑色，中风化，泥质胶结，层状构造，层理发育，微解理裂隙发育，泥质充填，岩芯呈碎块状、短柱状、长柱状，砂岩页岩互层平均值为8.36 MPa。
2.3. 水文条件
项目区及周边下部含水层有第四系松散岩类孔隙含水层、碎屑岩类裂隙含水层和碳酸岩类裂隙岩溶含水层。
第四系厚70~72 m，与下伏的上石盒子组地层呈不整合接触。自上而下可划分为上部松散砂层、孔隙潜水含水层组、中部砂层孔隙承压含水层组及底部粘土砾石孔隙承压含水层组，上、中部砂层含水层富水性相对较好，底部粘土砾石层富水性较弱。其中上部粉砂土层为黄泛区的冲、淤积物，厚度0~16.40 m，平均8.20 m。该层松软、空隙大、富水性中等，水位埋深一般在3 m左右，地下水位随季节变迁略有变化，变化幅度约为1.0~2.0 m，为区内民用主要水源；底部粘土砾石孔隙承压含水层是区内各基岩含水层的主要补给水源。
碎屑岩类裂隙含水层自上而下可划分为上石盒子组底部奎山砂岩裂隙承压含水层、下石盒子组砂岩裂隙承压含水层组、下石盒子组底部分界砂岩裂隙承压含水层及山西组砂岩裂隙承压含水层组。
上石盒子组底部奎山砂岩裂隙承压含水层：砂岩厚4.6~41.7 m，平均为13.9 m，其为灰~灰白色的中~粗粒含砾砂岩，整合接触于下伏地层，分布稳定，富水性中等，透水性较强，属富水性中等的含水层。下石盒子组砂岩裂隙承压含水层组：砂岩平均厚度为63.10 m，单层砂岩最厚为35.04 m，砂岩富水性弱~中等，渗透性差且不均一，地下水赋存受构造控制且以静储量为主，属富水性弱~中等的含水层。下石盒子组底部分界砂岩裂隙承压含水层：砂岩厚13 m左右，为浅绿、浅灰色中~粗粒含砾砂岩，属富水性弱~中等的含水层。山西组砂岩裂隙承压含水层组：砂岩1~6层，砂岩平均厚度为35.64 m，单层砂岩最小仅为0.79 m，多为灰色细~中粒砂岩，是开采7、9煤层的直接充水含水层，本组砂岩富水性弱，渗透性差，砂岩裂隙水以静储量为主，易于疏干。
3. 数值模拟
3.1. 模型建立
首先根据研究区域岩土层分布高度，岩土层参数，巷道分布信息，确定好模型的大小。分别建立桩、土层、岩层、巷道模型的计算模型。但是由于FLAC3D在前处理方面稍有不足，其在建立复杂的非对称的三维立体模型时会变得异常困难，所以在前处理阶段时，采用同为ITASCA公司开发的Rhinoceros建模软件进行模型的建立，并初步划分网格，再使用Griddle(通用网格处理器插件)对其进行网格细化并通过其内置的转换语言，将其转换为FLAC3D能够读取的网格文件形式再将所得到的网格文件导入到FLAC3D中，得到了我们所要的计算模型。图1为建立的模型示意图。
选取的本构模型为摩尔库伦本构模型，物理力学参数的参考资料为《工程地质手册(第五版)》，摩尔–库伦弹塑性模型共有8个参数，分别是体积模量、粘聚力、剪胀角、内摩擦角、泊松比、剪切模量、抗拉强度、弹性模量。表1为模型中用到的材料力学性能参数。
Table 1. Table of rock and soil material parameters
表1. 岩土体材料参数表
3.2. 模拟结果
3.2.1. 位移分析
截取全部土层三个方向的位移与云图进行分析：
首先从数值大小上来讲，竖向位移的数值明显比其他两个反向的变形要大，竖向位移从图中可以得知大多数在0~6 mm附近左右，而地表在X方向上的位移均在2 mm以内，最小值接近于0，Y方向上的的最大值为1.2 mm左右，只是在桥墩两侧由于桥墩的影响有轻微变形，其他地方的Y向变形接近于0。
从位移规律上来看，三个方向的位移分布均有明显的规律，其中，X方向的位移具有中间大，两边小的特征，而且从图2和图3中可以看出，在Y的正方向上，远离线路中心桥墩的X向位移大，靠近桥墩处的X向位移较小，由于注浆范围是以道路桥墩的中心线为基准确定的，所以在不注浆区域其变形量会稍大。图4中云图的不连续部分，可能也是由于桥墩所在地方导致的，Z向位移的规律变现为从左往右，数值依次减小，这可能是由于在同时跨落时，和位移与Z向的变形规律类似，都是从左往右，位移量依次减小，因为三个方向中，Z向的位移分量最大，其他两个方向的位移量相比Z方向很小，所以和位移的规律基本与Z方向的规律相似。
Figure 2. X-direction deformation of soil surface
图2. 土层地表X向变形
Figure 3. Y-direction deformation of soil surface
图3. 土层地表Y向变形
Figure 4. Z-direction deformation of soil surface
图4. 土层地表Z向变形
3.2.2 应力分析
1) 地表竖向应力分析
由图5中的应力分析云图可知，采空区垮落对地表的应力竖向应力分布的规律影响不大，无论是竖向应力还是水平向应力均按深度梯度增加。
2) 倾斜采空区局部最大主应力分析
图6为倾斜采空区局部的最大主应力在Y = 0剖面上的云图，在FLAC3D中，规定，正应力为拉应力，负应力为压应力，由图中云图的信息可知，最大正应力为2.7 × 106为正应力也就是拉应力，可知采空区上覆岩层垮落后有部分围岩处于受拉状态。
Figure 6. Distribution nephogram of local maximum principal stress in inclined goaf
图6. 倾斜采空区局部最大主应力分布云图
图7为塑性区分布图，由图可知，采空区垮落后一定范围内的岩体进入塑性屈服状态，图例中n表示现状态依然为屈服状态，p表示过去曾经屈服过但是现在应力已经降低到屈服应力以下，shear表示剪应力达到抗剪屈服强度，tension表示拉应力达到抗拉屈服强度。
Figure 7. Distribution diagram of plastic zone on section Y = 0
图7. Y = 0剖面上的塑性区分布图
4. 结论
通过FLAC3D的数值模拟，对采空区区域垮落后的位移和应力分析，可知垮落后的土层位移均有明显的规律性，竖向位移大于水平向位移，应力方面，采空区的垮落对土层的应力影响不大，但是对采空区区域的影响较大，在采空区上方会有部分区域的受压状态转变为受拉状态，采空区垮落后对应区域都有较大的塑性区域。