水-热耦合下露天矿冻结期靠帮开采边坡稳定性研究

谷浩源 ，韩 流 ，张 健 ，刘世宝 ，高志强 ，吴锋锋

（中国矿业大学 矿业工程学院，江苏 徐州 221116）

煤炭作为我国的主要能源，其开采量和消费量均居能源结构的首位[1]。我国目前拥有露天矿山20 000 余座，超过80%的露天煤矿位于北方季节性冻土地区，冻结期大多在半年以上[2-3]。以内蒙古自治区的露天煤矿为例，因其受季节性因素影响较大，其矿区生产在冬季和夏季对应剥离作业期和剥离停产期，且露天矿边坡处在反复冻融条件下，边坡稳定性也受一定的影响。为了提高采出率、降低剥采比，实现煤炭资源的回收和经济效益的提高，靠帮开采技术得到了广泛的应用[4]。冬季靠帮开采回收端帮压煤，可通过内排压帮提高帮坡角度，帮坡角的变化直接影响煤炭采出率，而边坡稳定性是保证安全开采的前提，因此这两者是靠帮开采的关键参数。

针对边坡稳定性和帮坡角的研究，国内外许多学者进行了研究，取得了丰硕成果。吴榕真等[5]利用FLAC3D软件，结合不同靠帮角度，分析了靠帮开采过程中边坡的稳定性变化，确定了靠帮开采最佳的端帮边坡角；佘长超等[6]对排土场边坡软岩物料进行了冻融循环下剪切破坏试验，得出了冻融循环次数增加将导致边坡稳定性下降的结论；马忠辉等[7]依据冬夏两季边坡岩层的物理力学参数，运用稳定性分析软件，计算两季端帮稳定性系数强化幅度，确定出不同季节下靠帮开采方案；陆翔[8]结合寒区季节性气候变化特征，分析了冻融作用对端帮边坡稳定系数及破坏模式的影响，基于冻结泥岩强度强化特性，提出了寒区露天矿季节性靠帮开采方法。虽然上述学者对露天煤矿端帮靠帮开采进行了深入研究，但针对冻结期冻结层厚度变化对靠帮开采边坡稳定性以及帮坡角的影响，尚未形成统一的认知。

基于此，以宝日希勒露天矿北部端帮为工程背景，基于多孔介质理论和水-热耦合理论，考虑冻结期不同厚度的冻结层对边坡稳定性影响的变化规律，建立软岩边坡稳定性数值计算模型，揭示冻结期靠帮开采边坡稳定性变化规律以及合理的开采时间。

1 水热耦合模型构建

随着露天矿进入冻结期，温度的改变将会引起岩土体内部冰水相变以及水分迁移，而冰水相变所产生的潜热会反作用于温度场，且相变引起的变形也将使得土体力学参数发生变化[9]。针对水-热耦合两相物理场问题，可通过构建二者独立方程，寻求耦合联系方程，为冻结期露天矿边坡水热特性分析提供理论依据。

1.1 温度场控制方程

基于傅里叶热传导定律和能量守恒原理，考虑二维水热耦合问题，并考虑相变热，针对多孔介质冻岩土导热控制微分方程如下[10]：

式中： ρ为岩、土密度，kg/m³；C(θ)为岩、土体积比热容，J/（kg·℃）[2,11]；T为岩、土瞬态温度，℃；t为时间，s； λ(θ)为岩、土导热系数，W/（m·℃）；∇为微分算子，对于二维问题为[∂/∂x,∂/∂y]；L为相变潜热，取值为334.5 kJ/kg；ρi为冰的密度，kg/m³； θ为冻岩、土体积含水率；θi为冻岩、土中冰的体积含量。

1.2 渗流场控制方程

边坡岩土冻融过程中，始终存在自由水流动，其运动过程符合广义达西定律。考虑边坡冻融水分场运动为非饱和多孔介质流动，选用Richards水分迁移方程[12]。考虑冰水相变对渗流场的影响，得冻土中未冻水迁移运动方程：

式中： θu为冻岩、土中未冻水的体积含量；ρw为冻土中水的密度，kg/m³；D(θu)为季冻边坡岩、土体的水力扩散系数；k(θu)为季冻边坡岩、土体在重力加速度方向的土体渗透系数。

为充分描述季冻岩、土体边坡渗流真实情况，TAYLOR[13]在确定水分扩散系数时，引入阻抗因子I表示季冻岩、土体孔隙中冰对未冻水迁移的阻碍，并给出了水分扩散系数的公式：

为简化计算引入相对饱和度S：

式中：c(θu) 为 比水容量，1/m；kS为饱和岩土渗透系数，cm/s；a0、m、l为与岩、土本构相关的参数，取值参考Van Genuehten（VG）滞水模型[14]；θr为 残余含水率； θs为饱和含水率。

固液比BI相关概念的提出[11]，使得季冻边坡岩、土体由冰水相变引起的温度、体积含水量和体积含冰量三者变化产生关联。其计算式如下：

式中：Tf为初始土体冻结温度，℃；B为土质含盐量变化经验系数[15]。

1.3 模型及验证

基于COMSOL Multi-physics 有限元数值模拟软件，具有的多物理场耦合能力以及非线性微分方程组求解能力，通过对其进行二次开发，实现边坡岩、土体水热力场耦合模拟，软件中内置系数型偏微分方程（PDE）形式如下式：

式中：ea为方程质量系数；u为方程自变量；da为 阻尼或质量系数；c为扩散系数； α为守恒通量对流系数； γ为守恒通量源； β为对流系数；a为吸收系数；f为源项。

将式（7）、式（8）代入式（1）中，整理后得温度场控制方程为：

为验证模型开发的合理性，参考魏道凯[16]对一维土柱冻胀试验的数据进行测试。土柱宽0.14 m，高0.4 m，采用与原文相同的四边形单元进行网格划分，水分场温度场边界条件与原文一致。模型计算所需参数为[17]：①材料：黏性土；②Tf：-1.33 ℃；③B：0.56；④a0：2.65；⑤m：0.25；⑥l：0.5；⑦S：0.76；⑧θs：0.432；⑨θu：0.33；⑩ks：10-8cm/s。

不同时刻温度和未冻水含量的试验值与模拟值结果对比如图1。

图1 模型试验值与模拟值对比Fig.1 Comparison of model test values and simulation values

从结果来看，模拟值得出的温度和未冻水含量变化规律与试验值一致，都呈现随时间增加而下降的趋势，而且无论是温度还是未冻水含量变化，在同一时刻其试验值均大于模拟值，这是由于数值模拟相较于实际试验而言，过于理想化，实际试验无法做到保温隔热，避免能量损失。综上所述，所采用的水热耦合计算模型与实际试验数据整体上贴切，证实了所开发的水热耦合模型描述冻土水热特性的合理性。

2 冻结期边坡稳定性

北方露天矿多以土质边坡为主，每年10 月份开始进入冻结期，岩土孔隙水受温度影响冷结为冰。随着低温时间延长，冻结层向深部扩展形成季节性冻土。此过程主要分为初冬冻融交替、冻结发展以及融冻3 个阶段[18]。研究冻结期（冻结发展阶段）冻结层厚度对边坡稳定性的影响，可为冻结期露天矿靠帮开采边坡稳定性分析提供理论依据。

2.1 冻结期岩土体强度参数变化

土和岩石都属于多孔介质。以土为例，土是由固体颗粒、水和气体组成。非冻结期和冻结期土体结构如图2。

图2 非冻结期和冻结期土体结构Fig.2 Soil structure in non-freezing period and freezing period

天然状态下，土是三相体系，而在形成冻土后就变为“四相”体系。冻土形成过程中，温度的改变将会引起冰水相变以及水分迁移，冰水相变所产生的潜热会反作用于温度场，相变引起的变形也将影响土体力学参数，而且冻结后土的强度增加往往是大于其天然状态下强度和冰强度之和[19]。对于土质边坡稳定性而言，抗剪强度是主要的影响因素，其中黏聚力和内摩擦角是关键参数[20]。宝日希勒露天矿区冻结期和非冻结期地层物理力学参数见表1。

表1 非冻结期和冻结期边坡地层物理力学参数Table 1 Physical and mechanical parameters of slope strata in non-freezing period and freezing period

2.2 冻结期温度对冻结层厚度的影响

根据宝日希勒露天矿采掘平面图选择某一北部端帮剖面，通过在Auto-CAD 中画出宝矿北帮冻结期模型，模型长600 m，高253 m，进入冻结期后，边坡台阶会因为外部温度降低的影响形成冻土层，冻结期北帮模型示意图如图3。

图3 冻结期北帮模型示意图Fig.3 North slope model of freezing period

研究冻岩土试样的试验不少，但涉及边坡冻岩土体的研究不多，研究参数相对而言较为稀缺，因此只能假设、简化或参见其它学者的研究成果。假设边坡岩体热力学和水力学参数同岩石试样，部分岩体热、水力学参数可用物理力学参数近似的岩样代替。各参数的取值参考相关文献中进行了详细讨论[12,21-24]，北帮边坡水热学参数见表2，北帮边坡水力学参数见表3。

表2 北帮边坡水热学参数Table 2 Hydrothermal parameters of north slope

表3 北帮边坡水力学参数Table 3 Hydraulic parameters of north slope

模型整体采用三角形单元进行网格划分，对局部加密处理，保证计算结果的精确性。根据地质资料，设置初始相对饱和度为0.6，模型水分场边界条件设为0 通量，两侧温度边界条件为绝热边界，下边界为温度与当地年平均气温保持一致的恒温带，模型上边界采用Dirichlet 温度边界条件，参考宝日希勒露天矿当地年平均气温，选择10 月中旬作为冻结期起点，得出上边界温度T（t）表达式：

式中：t为冻结期时间，d。

计算结果表明，坡面各处冻结厚度大致相同，以最低冻结温度-0.54 ℃所在位置作为冻结层厚度判断标准，最大冻结厚度不超过5 m。为此，选择高程+252 m 坡顶面中点为特征点，沿坡面内法线方向做向下延伸5 m 的截线，提取不同时间截线的温度作为冻结层厚度的判断依据。不同时期的边坡冻结层厚度如图4。

图4 不同时期的边坡冻结层厚度Fig.4 Thickness of frozen layer of slope at different periods

由图4 结果可知：初始冻结期为0 d 时，即每年10 月中旬，此时边坡整体温度为5 ℃，未产生冻结；随着时间增加，在冻结期60 d，即12 月中旬时，边坡冻结层厚度从0 m 增加到1.82 m；在前145 d 内，冻结层厚度随时间不断增加属于冻结发展阶段，最大可达3.45 m，这与宝日希勒露天矿地质勘探报告中最大冻结厚度3.5 m 接近，印证了模型开发的合理性；从145 d 至160 d，冻结层厚度维持在3.4 m 左右，出现稳定冻结；随着气温的逐渐回升，180 d 时边坡整体冻结厚度出现不连贯现象，在厚度为2.38 m 到3.13 m 处温度达到了冻结温度，冻结层包含在边坡内部。这是因为此时边坡处于融冻阶段，初期冻土层停止向下发展，冻土层厚度达到稳定，后期冻土层由底部和表层向冻土层消融解冻，在坡体内部形成不连贯冻结层，因此可以确定160 d 以后边坡处于融冻阶段，直至190 d 后，即第2 年4 月中下旬边坡冻层完全融化。

2.3 不同冻结层厚度下边坡稳定性

目前，关于边坡稳定性的计算方法，主要有极限平衡法、概率分析法、数值计算法[25]。其中，极限平衡法是工程领域中最基本的方法。GEOslope 是工程领域用于分析边坡稳定性十分常用的软件[26]，其主要是依据极限平衡法的分析原理，对边坡稳定性系数进行求解。基于稳定性系数Fs的极限平衡法有多种，常用于边坡稳定性分析的 方 法 有Morgenstern-Price 法[27]、Spencer 法[28]、JanbuGeneralized 法[29]、Bishop 法[30]。因矿山实际生产情况复杂，选择上述4 种方法进行计算，根据不同方法结果进行对比，选取较为合适且相对保守的稳定系数作为边坡稳定系数。

利用水热耦合模型得出的不同时期冻结层厚度，结合GEO-slope 软件进行冻结期边坡稳定系分析。利用不同方法对北帮边坡进行不同时期不同厚度冻结层下边坡稳定性分析，不同冻结层厚度下北帮稳定性系数见表4。

表4 不同冻结层厚度下北帮稳定性系数Table 4 Stability coefficients of north slope with different thickness of frozen layer

由表4 可知：随着冻结时间的增加，冻结层的发育，无论是哪种计算方法，其边坡稳定性都是在增加的，说明在冻结期岩体强度的提升，冻结层厚度的变化将对边坡稳定性有一定影响。依照选取较为合适且相对保守的稳定系数为选取原则，可以得出冻结期靠帮开采最大稳定系数为1.336，相较于非冻结期1.245 提升了7.3%。从表4 中可以看出边坡稳定性系数在2 月中旬到3 月中旬增幅为0.07%，可认为在此期间的稳定性系数最大且变化不大，靠帮开采的时机最佳。

3 冻结期靠帮开采帮坡角度

韩流等[31]、尚涛等[32]针对冻结期靠帮开采提出了条分式靠帮开采方法，在保证边坡稳定的前提下，提高了煤炭回收量的同时降低了剥采比，结合宝日希勒露天矿北帮实际情况冻结期选择不剥离覆岩的靠帮开采[20,33]。

宝日希勒矿区每年的10 月初进入冻结期，于第2 年4 月中旬解冻，最大冻土厚度可达3.45 m。选择在冻结2 月中旬至3 月中旬期间进行靠帮开采，此时的边坡稳定系数最高，靠帮开采角度相对而言较大，端帮煤炭回收量较多。通常选取2倍工作帮台阶高度作为内排压帮的台阶高度，采用条分式靠帮开采，选择30 m 作为压端帮底脚台阶的宽度，并用与端帮底脚高度相同的内排台阶及时跟进压帮，台阶高度为25 m。构建的端帮稳定性分析模型如图5。

图5 冻结期北帮稳定性分析模型Fig.5 Analysis model of north slope stability during freezing period

宝日希勒露天矿原北帮帮坡角为26°，当冻结层厚度为3.4 m 时稳定系数为1.336。分别对帮坡角为27°、28°、29°以及30°的端帮进行稳定性计算并得出稳定系数，得出的北帮在不同帮坡角时边坡稳定性系数见表5。由表5 可知：随着帮坡角不断增大，其对应的稳定系数会随之减小。

表5 不同帮坡角的稳定系数Table 5 Stability coefficients of different slope angles

根据边坡稳定性的判定依据可知，当滑体上的抗滑力（矩）与下滑力（矩）比值大于1，即抗滑力大于下滑力，此时滑体处于安全状态，边坡稳定；反之，边坡失稳。在满足规范端帮边坡稳定系数1.10 前提下，确定靠帮开采的最佳帮坡角为29°。

4 工程实践

根据上述研究可知，冻结期内由于边坡表面形成了强度较大的冻结层，边坡整体稳定性提升，稳定系系数得到增大。在此基础上，对宝日希勒露天矿冻结期采用内排压帮宽度30 m，内排台阶高度25 m 的不剥离靠帮开采方式，进行冻结期靠帮开采。靠帮开采合理帮坡角为29°，通过CAD计算得可回收的煤层横截面面积为180 m2，根据地质资料描述，采用350 m 的可采长度计算，可回收煤炭8.4 万t，以 150 元/t 的市场价格计算，可创造1 260 万元的经济价值，由此可见基于冻结期软岩强度提升的靠帮开采产生的经济效益是十分显著的。

融冻阶段坡体含水量如图6。

图6 融冻阶段坡体含水量Fig.6 Water content of slope in thawing and freezing stage

冬春交替时节，气温逐渐上升，当最低气温达到0 ℃及以上，露天矿边坡岩土进入解冻期，如图6（a）；随着时间增加，气温回暖，冻结层开始融化使得边坡岩体含水量增大，坡体水压力增大，滑坡危险系数增加，如图6（b）。

解冻层岩石黏聚力为20.5 kPa、内摩擦角为32.3°，密度为2.20 t/m³。将解冻层岩石力学参数代入到Geo-slope 中，得出解冻期帮坡角为29°时端帮稳定系数为1.084，无法达到规范端帮边坡稳定系数1.10 的要求。为满足安全要求，可通过增加内排台阶宽度的方法，使其内排压帮体积增大，以保证解冻后端帮稳定系数的增加。解冻期端帮保障措施如图7，图中△V为增大内排台阶宽度所需的废弃物量。

图7 解冻期端帮保障措施Fig.7 Guarantee measures for north slope in thawing period

台阶宽度以5 m 为1 个步距递增，经多次稳定性计算后，确定内盘宽度在原有基础上增加10 m，即内排宽度为40 m 时，此时的稳定系数为1.122，满足规范端帮边坡稳定系数1.1 的要求。

5 结 语

1）针对冻结期冻结层厚度受温度、水分变化影响，利用水热耦合控制方程构建了1 个分析冻结期边坡冻结层厚度变化的模型。

2）边坡稳定性系数与冻结层厚度呈现正相关，冻结期内稳定性系数随着冻结厚度增加而增大，降低而减小。靠帮开采最佳时机为2 月中旬至3月中旬，冻结层厚度达到3 m 以上，最大为3.45 m，稳定性系数1.336。

3）得出了最佳时机进行靠帮开采，合理帮坡角为29°，与原设计26°相比，多回收煤炭8.4 万t，创造1 260 万元经济价值，考虑解冻期边坡稳定性降低，通过增大内排台阶宽度40 m 的方式，保障了边坡过渡期安全。