尾矿库是用于堆存尾砂或其他工业废渣的场所，具有安全风险高、环境污染风险大和管理难度大等特点，是高势能人造泥石流重大危险源与潜在环境污染源（王昆等，2024）。我国现有各类大小尾矿库1万多座，一旦发生尾矿库溃坝和尾矿泄漏等事故，不仅会严重威胁人民群众的生命财产安全，还会对生态环境造成极大的破坏，甚至可能会影响社会安定（张家荣等，2019）。因此，尾矿坝稳定性分析一直是尾矿库安全领域的重要课题。与此同时，推进尾矿资源的安全回采与二次综合利用，不仅能从根源上化解尾矿库的安全隐患与环境问题，还能实现资源循环利用，具有巨大的发展潜力。
尾矿坝稳定性是矿山企业十分关注的问题，国内外学者从静力稳定性、动力稳定性及渗流稳定性等多个方面对其开展研究（杨春和等，2021；许云美等，2023；朱远乐等，2024）。静力稳定性主要研究在尾矿坝自重、堆积荷载等静力条件下，尾矿材料的物理力学性质、坝体结构特征等因素对尾矿坝稳定性的影响（杨春和等，2022）；动力稳定性则研究尾矿坝在面临地震和爆破振动等动力荷载冲击时的稳定性特征（王文松等，2018；李生伟等，2023）。其中，尾矿坝渗流稳定性一直是研究的一个难点，学者们通过室内渗流试验、现场监测以及数值模拟等方法对尾矿坝内渗流场的分布进行分析，研究渗流对坝体稳定性的作用机制，为制定合理的排水减压措施提供科学依据（吴高杰等，2022；Du et al，2024；梁冰等，2025）。在尾矿坝渗流稳定性研究基础上，许多学者开始考虑多场耦合作用对尾矿坝稳定性的影响，并利用数值模拟软件建立尾矿坝多场耦合数值模型，研究尾矿坝在渗流—应力耦合（邓红卫等，2016；Popescu et al，2024）、流—固—化耦合（樊嘉庆，2020）以及热—流—固耦合（李向阳等，2022）作用下的稳定性，分析库水位、化学腐蚀、降雨强度和温度等因素对尾矿坝稳定性的影响规律，为尾矿坝的稳定性分析及后期管理提供理论支持。
近年来，尾矿资源利用受到国内外学者的广泛关注。以往研究表明，尾矿资源潜力巨大，通过推进尾矿综合利用，不仅能够消除尾矿库这一重大危险源，还能够有效促进矿区生态环境保护，对于实现资源循环利用和矿山可持续发展具有重要的战略意义（Araya et al，2021；李鑫等，2024；伍天华等，2024）。随着选矿技术装备的不断进步，国内外学者对尾矿资源安全回采及二次利用进行了研究（Alcalde et al，2018；任珩等，2024）。例如：通过对尾矿库回采方法、回采顺序及回采工艺的系统研究，对比分析了不同回采方案的优缺点，并提出了一种改进的尾矿库回采方案（王新岩，2023；付宪吕等，2024；郭大平等，2023）。在此基础上，通过室内试验（周罕等，2023）和数值模拟（段晓阳等，2023；Pan et al，2023；肖容等，2023）等方法，研究了在尾矿库回采的不同阶段尾矿坝稳定性特征的演化规律，确保尾矿坝在整个回采过程中保持稳定运行。
综上所述，国外学者关于尾矿库资源回采的研究主要集中在回采技术和工艺方法等方面，国内学者则更多关注尾矿库回采过程中的静力稳定性分析，以及不同回采方法对尾矿坝稳定性的影响。然而，在综合考虑渗流—应力耦合作用对尾矿库回采过程中稳定性影响的研究方面，目前国内外学者开展的研究工作相对较少，缺乏系统性成果。在尾矿库回采过程中，干滩长度与浸润线的变化会引发渗流场的改变，从而破坏尾矿库原有的应力平衡状态，导致尾矿坝产生变形破坏。因此，开展尾矿库回采过程中渗流—应力耦合作用下的稳定性研究，对于揭示坝体变形破坏机制和保障尾矿库安全回采具有重要的理论意义和工程参考价值。
1 渗流—应力耦合分析模型研究
在尾矿库回采过程中，尾砂介质内部水头分布发生变化，地下水在水头差驱动下沿着孔隙发生渗流运动，产生渗流压力。渗流压力以体积力的形式作用于尾砂介质，改变其应力场，进而引起位移场的变化。位移场的改变反过来又影响介质的渗透系数和渗流场分布，形成渗流—应力耦合作用。
1.1 渗流场对应力场的作用
渗流场对应力场的作用主要表现为流体在孔隙介质中运动时对固体骨架施加的孔隙水压力，包括静水压力和超静水压力（渗流压力）。
一方面，尾矿库原有地下水及回采工作面基坑中的积水会给尾矿库施加静水压力，它以压力荷载的形式作用于接触面上，其函数表达式为
式中： 为静水压力； 为水的容重； 为位置水头； 为水头高度。
另一方面，尾砂介质中的流体在水头差的作用下发生渗流运动，产生渗流压力，以渗流体积力的形式作用于尾砂介质中，从而影响尾矿库内原有的应力场和位移场，改变其原有的应力平衡状态。根据水力学理论，渗流体积力与水力梯度成正比（李雪梅，2017），其函数表达式为
式中： 为渗流体积力； 分别为渗流体积力沿x、y方向的分力； 分别为水力坡降沿x、y方向的分量。
1.2 应力场对渗流场的作用
随着尾矿库表面尾砂的回采，其内部应力将重新分布，导致某些部分卸荷或荷载减少，形成卸荷效应（如应力释放和体积膨胀），从而使尾矿库原有的应力场和位移场发生变化。应力场和位移场的变化会改变尾砂介质的孔隙结构和分布特征，多孔介质骨架的体积应变会导致孔隙率等孔隙特征发生改变，进而影响尾砂介质的渗透系数，最终导致渗流场的变化。因此，应力场对渗流场的影响主要是通过改变尾砂介质的孔隙特征和渗透系数来实现的。多孔介质的渗透系数可用应力状态函数表示为
此时，受应力场影响的二维渗流微分方程可表示为
式中： 为应力张量； 为渗透系数。
1.3 渗流—应力耦合模型在COMSOL中的实现
COMSOL多物理场仿真软件可通过添加多个物理场接口实现多物理场耦合分析，具有易操作和可视化强的特点，被广泛应用于渗流—应力耦合分析模型研究中。本文研究中假设地下水服从达西定律，尾矿库尾砂为各向同性等效连续介质，因此渗流场选择达西定律接口，应力场选择固体力学接口，以此搭建尾矿库渗流—应力耦合分析模型。
在COMSOL软件进行渗流—应力耦合计算中，考虑孔隙水压力作用的达西定律表达式为
式中： 为达西速度； 为渗透率； 为流体的黏度； 为孔隙水压力； 为流体密度； 为重力加速度。
固体力学和达西定律接口进行耦合计算时，通过达西定律接口计算出模型中的渗流场分布情况，并将其作为初始条件导入到固体力学接口中，以此实现渗流场与应力场的耦合分析。由于尾矿库浸润线以上部分尾砂处于非饱和状态，因此用理查兹方程来描述流体在尾砂介质中的非饱和渗流过程，考虑有效应力的控制方程为
式中： 为含水量； 为有效饱和度； 为储水系数； 为介质饱和渗透率； 为相对渗透率； 为渗流源汇项； 为Biot系数； 为体积应变量。
2 尾矿库数值模型的构建
2.1 工程概况
河南省某钼矿尾矿库总体上为山谷型尾矿库，局部属傍山型。该尾矿库目前堆积坝顶标高约为1 420.76 m，现状总坝高为60.50 m，全库容约为85×104 m3，有效库容约为76.5×104 m3。初期坝为碾压式土石坝，坝高为14.76 m，坝顶长度为69.6 m，顶宽为3.17 m，外坡比为1∶2.48，内坡比为1∶1.75。堆积坝采用上游式筑坝，堆积坝高约为45.74 m，水平长度为260.62 m，堆积坝整体外坡比约为1∶5.62，共堆有18级子坝。子坝坝面覆土植草，植被生长茂密，护坡良好。该尾矿库现状等别为四等库，根据《尾矿库安全规程》（GB39496-2020），正常运行时尾矿坝稳定性计算的荷载包括运行期正常库水位时的稳定渗透压力、坝体自重以及坝体和坝基中的孔隙水压力，坝体正常运行的最小安全系数为1.25。
本次设计回采的尾砂按照成分可划分为2个部分进行利用，通过取样化验，具备再选价值的尾砂输送至选矿厂，通过二次造浆后进行再选；不具备再选价值的尾砂则用作建筑材料，用于建设或修整矿区内道路。
2.2 尾矿库回采方式
尾矿库采用坝前干式与库内湿式的联合回采方式，按照先内后外、先库后坝、先上后下和分层开采的顺序进行回采，如图1所示。库内湿式回采区域选择水枪与漂浮式砂泵船联合工艺进行尾砂回采，由库内向坝前方向推进，设计回采分层高度为2 m，回采坑坡比按1∶3放坡以防止坡面垮塌，回采采坑水面以上的边坡高度应小于等于3 m。坝前干滩部分选择干式回采，采用推土机将干滩区堆积的尾砂推送至湿采区再回采，在回采过程中需保持滩面坡度向库内倾斜大于等于1%，干滩长度大于等于50 m。在回采初始阶段，库内湿式回采区域先开采1~2层，形成必要滩面形态，然后按库内湿采区与坝前干采区有序循环进行回采，回采过程中干采区域始终高于湿采区域2 m以上。
2.3 模型建立及参数设定
结合矿山企业提供的工程地质剖面图及尾砂沉积规律，建立尾矿库坝体稳定性二维分析模型，如图2所示。构建模型时，固体力学模块选择德鲁克—普拉格本构模型匹配莫尔—库伦准则进行计算，对模型底部施加固定约束，左右两侧施加辊支撑约束，同时对整个模型施加y方向的重力。渗流场则选择达西定律接口开展计算研究，材料选择非饱和多孔介质，持水模型选择V-G模型，根据工程地质剖面图中浸润线的位置确定压力水头边界，模型底部为无流动边界。尾矿库坝体材料物理力学参数根据室内试验成果和现有的工勘报告等相关资料综合确定，计算所需材料的物理力学参数见表1。V-G模型计算需要用到模型参数α、n、l和残余含水率θr等参数，相关参数根据室内试验并参考经验数值选取，如表2所示。在渗流—应力耦合计算时，首先通过达西定律模块计算出尾矿库内部孔隙水压力分布，然后将渗流场计算得到的孔隙水压力作为荷载施加到应力场。通过固体力学和达西定律接口进行耦合计算，基于有效应力原理和强度折减原理分析尾矿库在回采过程中的渗流场、应力场及稳定性变化情况。
表1 尾矿库坝体材料物理力学参数Table 1 Physical and mechanical parameters of the tailings dam materials |
| 材料名称 | 容重/（kN·m-3） | 孔隙率 | 杨氏模量/MPa | 泊松比 | 内聚力/kPa | 内摩擦角/（°） | 渗透系数/（cm·s-1） |
|---|---|---|---|---|---|---|---|
| 尾中（细）砂 | 21.6 | 0.54 | 320 | 0.335 | 30 | 27 | 1.5×10-3 |
| 尾粉砂 | 22.3 | 0.53 | 210 | 0.326 | 10 | 20 | 3.75×10-4 |
| 尾粉土 | 21.0 | 0.51 | 160 | 0.320 | 15 | 12 | 2.5×10-5 |
| 碾压土石坝 | 20.6 | 0.02 | 780 | 0.350 | 180 | 34 | 5.0×10-1 |
| 基岩 | 26.2 | 0.015 | 5 000 | 0.350 | 500 | 40 | 7.5×10-5 |
表2 V-G模型计算所需基本参数Table 2 Basic parameters required for V-G model calculation |
| 材料名称 | α/（1·m-1） | n | l | θr |
|---|---|---|---|---|
| 尾中（细）砂 | 0.91 | 2.01 | 0.5 | 0.034 |
| 尾粉砂 | 0.54 | 2.32 | 0.5 | 0.037 |
| 尾粉土 | 0.15 | 2.53 | 0.5 | 0.040 |
| 碾压土石坝 | 0.95 | 1.81 | 0.5 | 0.023 |
| 基岩 | 0.37 | 1.66 | 0.5 | 0.034 |
3 尾矿库回采过程渗流与应力耦合分析
根据尾矿库回采方案，尾矿库回采过程中干滩长度应大于等于50 m，因此取回采过程推进0，50，100，150，200，230 m共6种工况，利用COMSOL软件分析回采过程中尾矿库坝体的渗流规律、应力变化、变形特征和坝体稳定性。在尾矿库回采过程中，不同推进距离工况下尾矿库坝体的最大孔隙水压力、最大Mises等效应力、最大剪应力、最大位移、最大等效塑性应变和安全系数见表3。由于尾矿库回采尚处于可行性研究阶段，目前矿山企业已根据相关研究成果开展尾矿库回采先期准备工作，后续将采用尾矿库回采过程中的现场监测数据与模拟结果进行对比，验证数值模拟的可靠性。
表3 回采过程中尾矿库参数数值模拟结果Table 3 Numerical simulation results of tailings pond parameters during the process of mining |
| 推进距离/m | 最大孔隙水 压力/kPa | 最大Mises 等效应力/kPa | 最大剪应力/kPa | 最大位移/mm | 最大等效塑性应变 | 安全系数 |
|---|---|---|---|---|---|---|
| 0 | 225.402 | 311.622 | 341.332 | 123.960 | 0.0407 | 1.881 |
| 50 | 234.382 | 324.087 | 356.647 | 125.397 | 0.0412 | 1.861 |
| 100 | 241.328 | 322.609 | 354.931 | 130.005 | 0.0427 | 1.845 |
| 150 | 250.383 | 323.721 | 356.413 | 142.189 | 0.0468 | 1.809 |
| 200 | 265.433 | 324.785 | 357.551 | 148.963 | 0.0487 | 1.783 |
| 230 | 277.866 | 334.612 | 369.501 | 153.616 | 0.0500 | 1.743 |
3.1 尾矿库回采过程渗流规律分析
在尾矿库回采过程中，随着回采区域由库内向坝前推进，基坑中的水将向尾砂介质中渗流，导致尾矿库的渗流场发生变化。图3所示为尾矿库回采过程中尾矿库坝体的饱和度分布情况。图中红色区域为饱和区域，蓝色区域为非饱和区域。由图3可以看出，随着回采区域向前推进，尾矿库原本非饱和区域变为饱和区域，饱和区域范围逐渐增大，且有向堆积坝表层蔓延的趋势；非饱和区域范围逐渐减小，主要集中在堆积坝表层和初期坝区域。随着尾矿库从0 m回采至230 m，尾矿库干滩长度由280 m减小至50 m，饱和区域面积由10 811 m2增大至12 455 m2（不考虑基岩），增幅为15.21%。图4所示为尾矿库回采推进0 m和230 m这2种工况下尾矿库坝体的孔隙水压力分布云图。为了查看尾矿库回采过程中坝体应力情况，对模型中的基岩部分进行隐藏，该操作不影响相关数值计算结果。由图4可知，尾矿库坝体的孔隙水压力随着地层深度的增加而增大，最大孔隙水压力位于尾粉土层底部与基岩结合处。
这是由于尾粉土的渗透性较低，阻碍了地下水的快速排出，导致孔隙水压力在该区域积聚。图5中灰色曲线为尾矿库回采过程中尾矿库坝体的最大孔隙水压力变化情况。随着尾矿库从0 m回采至230 m，尾矿库坝体的最大孔隙水压力逐渐增大，压力值从225.402 kPa增大至277.866 kPa，增幅为23.28%。孔隙水压力增大的原因主要是随着回采工作的推进，基坑中的水增多并向深部渗流，导致下部尾砂中的水无法及时排出，从而形成超静水压力。孔隙水压力的增大会降低尾砂介质的有效应力，导致土体的抗剪强度降低，甚至引发液化现象，进而增加尾矿库变形和破坏的风险，严重影响尾矿坝的整体稳定性。
3.2 尾矿库回采过程应力变化分析
von Mises应力是一种基于剪切应变能的等效应力，可用于快速确定尾矿坝的最危险区域。图6所示为尾矿库回采推进0 m和230 m这2种工况下坝体的von Mises应力分布云图。由图6可以看出，von Mises应力随土层深度的增加而增大，其最大值出现在初期坝内侧及堆积坝最大埋深处。这表明该区域为应力集中区，是尾矿坝的最危险区域。由于初期坝的强度远高于堆积坝，因此堆积坝最大埋深处更容易发生变形和破坏。图5中红色曲线为尾矿坝的最大von Mises应力变化曲线。由图5可以看出，尾矿坝的最大von Mises应力随回采推进距离的增加而增大。其中，0~50 m和200~230 m回采段的斜率较大，主要原因是：回采0~50 m区域时，回采工作面中的积水开始渗流到尾砂介质中，尾矿库原本非饱和区的饱和度提高，孔隙水压力增大，导致von Mises应力显著增大；回采200~230 m区域时，孔隙水压力的进一步增大导致尾砂介质的有效应力降低，从而加速了von Mises应力的增大。von Mises应力的增大会降低尾矿坝的抗剪强度，增大变形和破坏的风险，进而影响其整体稳定性。
Tresca剪应力则主要用来表征尾矿坝的最大剪应力，从而判断其发生剪切破坏的最危险区域。图7所示为尾矿库回采推进0 m和230 m这2种工况下坝体的Tresca剪应力分布云图。与von Mises应力分布特征类似，Tresca剪应力也随着土层深度的增加而增大，最大值出现在初期坝内侧及堆积坝最大埋深处。图5中蓝色曲线为尾矿坝的最大Tresca剪应力变化曲线。由图5可知，尾矿坝的最大Tresca剪应力随着回采推进距离的增加而增大，其变化趋势与von Mises应力基本一致。其中，0~50 m和200~230 m回采段的斜率较大，主要原因是：回采0~50 m时，回采工作面中的积水导致孔隙水压力增大，土体的饱和度提高，剪应力显著增大；回采200~230 m区域时，浸润线的上升进一步加剧了孔隙水压力的增大，导致剪应力加速上升。Tresca剪应力的增大会增加尾矿坝发生剪切破坏的风险，特别是在堆积坝最大埋深处，由于土体强度较低，更容易发生破坏。
3.3 尾矿库回采过程变形特征分析
在实际工程中，尾矿坝的位移是判断其稳定性的重要指标。图8所示为尾矿库回采推进0 m和230 m这2种工况下坝体的位移分布云图。由图8可以看出，尾矿坝最大位移出现在堆积坝底部与初期坝交界处。这是由于该区域的孔隙水压力、von Mises应力和Tresca剪应力均较大，对土体的变形和破坏产生了显著影响。随着回采区域向前推进，尾矿坝最大位移所在区域从堆积坝底部向堆积坝顶部扩展。图9中灰色曲线展示了回采过程中尾矿坝的最大位移变化情况。由图9可以看出，尾矿坝的最大位移随着回采推进距离的增加而增大，当回采区域从0 m推进至230 m时，尾矿坝的最大位移从123.960 mm增大至153.616 mm，增幅为23.92%。位移的显著增大会增加尾矿坝发生失稳破坏的风险。在尾矿库回采过程中，由于渗流场与应力场的耦合作用，坝体的潜在滑动面表现为典型的圆弧形滑动面，并有向坝体顶部蔓延的趋势。这种滑动面的形成主要是由于坝体内部发生剪切破坏，特别是在堆积坝底部与初期坝交界处，由于应力集中和孔隙水压力增大，土体的抗剪强度显著降低，因此更容易发生滑动破坏。
3.4 尾矿库回采过程坝体稳定性分析
图10所示为尾矿库回采0 m和230 m这2种工况下坝体的等效塑性应变云图。由图10可以看出，塑性区主要位于尾矿库堆积坝与基岩及初期坝结合处，塑性应变最大值出现在尾矿库堆积坝底部与初期坝顶部交界处。这是由于该区域的应力集中和孔隙水压力增大，导致土体的抗剪强度显著降低，从而容易发生塑性变形。图9中红色曲线展示了尾矿库回采过程中坝体的最大等效塑性应变变化情况。由图9可以看出，尾矿库坝体的最大等效塑性应变随着回采推进距离的增加而增大。当回采区域从0 m推进至230 m时，最大等效塑性应变由0.0407增大至0.0500，增幅为22.85%。最大等效塑性应变的变化趋势与最大位移的变化趋势一致，表明塑性变形是导致坝体位移增大的主要原因。
运用强度折减理论，借助COMSOL软件的参数辅助扫描功能，计算得到尾矿库回采过程中尾矿库坝体的安全系数，如表3所示。图9中蓝色曲线为尾矿库回采过程中尾矿库坝体的安全系数变化情况。由表3和图9可知，随着回采区域由0 m推进至230 m，尾矿库坝体的安全系数从1.881下降至1.743，降幅为7.34%。安全系数的下降主要是由于回采过程中干滩长度减小，导致孔隙水压力增大和潜在滑动面的摩擦系数降低。此外，安全系数的下降趋势在回采后期呈现加速现象。这是由于随着回采区域的推进，尾矿库湿采区域的水渗流到尾砂介质中，导致坝体内部侵蚀加剧和浸润线抬高，有效应力进一步降低，进而引起安全系数加速下降。根据数值模拟结果，尾矿库回采过程中的安全系数均超过《尾矿库安全规程》（国家市场监督管理总局，2020）规定的最小值（1.25），且具有一定富余。这说明尾矿坝在回采过程中处于稳定状态，研究结果可为尾矿库回采工程提供重要参考。
4 结论
以河南省某钼矿尾矿库作为研究对象，基于渗流—应力耦合理论，运用COMSOL软件建立了尾矿库的渗流—应力耦合分析模型，系统研究了尾矿库回采过程中坝体的渗流规律、应力变化、变形特征及稳定性演化规律，得出如下结论：
（1）在渗流—应力耦合作用下，随着尾矿库回采向前推进，尾矿库的干滩长度减小，饱和区域逐渐扩大，孔隙水压力、von Mises应力和最大剪应力均逐渐增大，导致尾矿坝整体位移和等效塑性应变逐渐增大，尾矿坝的安全系数逐渐减小，尾矿库的稳定性下降。
（2）在尾矿库回采过程中，尾矿坝的安全系数由1.881下降至1.743，但仍远大于《尾矿库安全规程》（国家市场监督管理总局，2020）规定的四等坝安全系数（1.25），说明尾矿库回采过程中的稳定性满足规范要求，在满足安全生产管理条件下，坝体可正常运行。本文研究成果能够为尾矿库回采稳定性及可行性分析提供理论支持。
甘公网安备 62010202000672号
