·666* 工程科学学报，第40卷，第6期 KEY WORDS block caving:ore drawing;moisture content;non-Newtonian fluid:isolated movement zone 自然崩落法作为唯一能与露天开采经济效益相 ASTM)提出的岩土颗粒粒级分类标准9，将直径小 媲美的高效地下采矿方法，具有生产能力大，便于组 于19mm的细颗粒流作为似非牛顿流体，大于19 织管理，作业安全及开采成本低等优点，近年来越来 mm的矿岩作为离散体，并提出一种新的自然崩落 越受到采矿工作者的关注-.但自然崩落法在国 法放矿模拟研究方法，即将矿岩细颗粒流与大块离 内有色金属矿山中的应用并不广泛，除与矿山开采 散矿岩分别通过格子波尔兹曼法与离散元法模拟， 的岩体性质相关外，还与我国现阶段自然崩落法放 基于格子波尔兹曼法一离散元法耦合方法对不同含 矿控制理论研究不成熟密切相关，因为放矿效果直 水率下自然崩落法放矿松动体的形态进行细观研 接关系到自然崩落法应用的成败的.放矿过程中， 究.通过将模拟结果与已有研究结论进行对比分 矿岩堆体内产生松动的部分为放矿松动体(isolated 析，验证了基于格子波尔兹曼法一离散元法耦合算 movement zone,IMZ),其形态对于采场结构参数和 法的放矿模型在自然崩落法放矿松动体形态研究中 放矿方式的选择至关重要6 的可靠性与准确性 放矿松动体的形态主要取决于散体崩落矿岩的 流动性及放矿参数，而矿岩含水率是影响散体崩落 1非饱和矿岩颗粒间受力分析 矿岩流动的重要因素⑨.国外Campbellto分析了 崩落矿岩细颗粒间通常含有孔隙水，因此处于 矿岩散体中水的主要存在形式，及其对细小颗粒散 非饱和状态的崩落矿岩散体为固一液一气三相耦合 体摩擦力和粘结力的作用机理；Ketterhagen等nu通 场.此时颗粒表面为润湿状态，相邻湿颗粒会在接 过建立楔形漏斗模型研究了准三维维度下的颗粒流 触点及其附近形成半月状液桥（见图1(a)),进而 态：国内王洪江等☒通过放矿物理模拟试验，得出 将颗粒连接在一起，故矿岩颗粒间的细观水力特性 含水率在10%左右时矿岩抗剪强度较大；张春阳 主要取决于矿岩的含水率.非饱和矿岩颗粒可简化 等)基于相似模型试验研究了含水率对矿岩流动 为湿颗粒的球体模型0，相邻颗粒间的作用力除包 性的影响，得出有效控制矿岩散体的含水量可避免 括离子水化作用力、静电力、范德华力、双电层作用 放矿堵塞与放矿大量贫化等结论.由于在宏观尺度 力等分子尺度的作用力外，还包括颗粒间液桥产生 下无法直接揭示含水矿岩颗粒与液桥的相互作用， 的毛细作用力四.非饱和状态下，颗粒间的毛细作 因此需从细观尺度探究这种相互作用的深层次机 理，进而解释非饱和矿岩宏观流动特性. 用力主要来自于液体的表面的张力及气一液交界处 目前，有限单元法（finite element method, 的静水压力.毛细作用力的计算主要涉及表面张力 FEM)、离散元法(discrete element method,.DEM)、 Y、颗粒半径R、颗粒间距S及固-液接触角0等因素 元胞自动机理论(cellular automaton,CA)及REBOP (见图1(b)),其中，目前毛细作用力通常采用La- (rapid emulator based on PFC)等计算机仿真手段 place--Young方程计算求解u,其计算表达式如下： 在放矿理论研究方面均有不同程度的应用.其中， Ry" R y(tyy(R (1) 离散元法由Cundall于1971年提出并被广泛应 用于岩石力学及颗粒流等领域.但由于当前计算机 式中，△p为气一液交界面的压力差，y=cot(p+), 处理能力的限制，单纯应用离散元法对自然崩落法 y"为y的二阶导数，P为半填充角 放矿进行工程级别的模拟，需生成10量级的颗粒， 相应的总液力计算方程为： 对计算机处理能力要求较高，且模拟过程中可能出 F。=2 myRsin osin(p+0)+TR△psin'p(2) 现高应力区内颗粒接触不稳定的现象，难以在全域 其中，2 myRsin osin(o+8)与πR△o sin分别表 实现离散元模拟a.此外，格子Boltzmann法(Iat- 征表面张力与静水压力. tice boltzmann method,LBM)作为一种新兴的流体 2LBM-DEM耦合算法 数值模拟手段，具有算法简单、并行性高及能处理复 杂边界条件等优点，在多相多组分流、微尺度流动及 2.1LBM方程非牛顿流体化 多孔介质流等领域得到了广泛地发展和应用图 对于二维两相流问题的处理，目前格子波尔兹 因此，基于上述理论与研究，本文根据美国材料与试 曼法通常采用D209(2为空间维数，9为离散速度 验协会(American Society for Testing and Materials, 数)，模型如图2所示.该模型将整个流体域划分为工程科学学报，第 40 卷，第 6 期 KEY WOＲDS block caving; ore drawing; moisture content; non-Newtonian fluid; isolated movement zone 自然崩落法作为唯一能与露天开采经济效益相 媲美的高效地下采矿方法，具有生产能力大，便于组 织管理，作业安全及开采成本低等优点，近年来越来 越受到采矿工作者的关注［1--4］． 但自然崩落法在国 内有色金属矿山中的应用并不广泛，除与矿山开采 的岩体性质相关外，还与我国现阶段自然崩落法放 矿控制理论研究不成熟密切相关，因为放矿效果直 接关系到自然崩落法应用的成败［5］． 放矿过程中， 矿岩堆体内产生松动的部分为放矿松动体( isolated movement zone，IMZ) ，其形态对于采场结构参数和 放矿方式的选择至关重要［6--8］． 放矿松动体的形态主要取决于散体崩落矿岩的 流动性及放矿参数，而矿岩含水率是影响散体崩落 矿岩流动的重要因素［9］． 国外 Campbell［10］分析了 矿岩散体中水的主要存在形式，及其对细小颗粒散 体摩擦力和粘结力的作用机理; Ketterhagen 等［11］通 过建立楔形漏斗模型研究了准三维维度下的颗粒流 态; 国内王洪江等［12］通过放矿物理模拟试验，得出 含水率在 10% 左右时矿岩抗剪强度较大; 张春阳 等［13］基于相似模型试验研究了含水率对矿岩流动 性的影响，得出有效控制矿岩散体的含水量可避免 放矿堵塞与放矿大量贫化等结论． 由于在宏观尺度 下无法直接揭示含水矿岩颗粒与液桥的相互作用， 因此需从细观尺度探究这种相互作用的深层次机 理，进而解释非饱和矿岩宏观流动特性． 目 前，有 限 单 元 法 ( finite element method， FEM) 、离散元法( discrete element method，DEM) 、 元胞自动机理论( cellular automaton，CA) 及 ＲEBOP ( rapid emulator based on PFC) ［14］等计算机仿真手段 在放矿理论研究方面均有不同程度的应用． 其中， 离散元法由 Cundall［15］于 1971 年提出并被广泛应 用于岩石力学及颗粒流等领域． 但由于当前计算机 处理能力的限制，单纯应用离散元法对自然崩落法 放矿进行工程级别的模拟，需生成 108 量级的颗粒， 对计算机处理能力要求较高，且模拟过程中可能出 现高应力区内颗粒接触不稳定的现象，难以在全域 实现离散元模拟［16］． 此外，格子 Boltzmann 法( lat￾tice boltzmann method，LBM) 作为一种新兴的流体 数值模拟手段，具有算法简单、并行性高及能处理复 杂边界条件等优点，在多相多组分流、微尺度流动及 多孔介质流等领域得到了广泛地发展和应用［17--18］． 因此，基于上述理论与研究，本文根据美国材料与试 验协会( American Society for Testing and Materials， ASTM) 提出的岩土颗粒粒级分类标准［19］，将直径小 于 19 mm 的细颗粒流作为似非牛顿流体，大于 19 mm 的矿岩作为离散体，并提出一种新的自然崩落 法放矿模拟研究方法，即将矿岩细颗粒流与大块离 散矿岩分别通过格子波尔兹曼法与离散元法模拟， 基于格子波尔兹曼法--离散元法耦合方法对不同含 水率下自然崩落法放矿松动体的形态进行细观研 究． 通过将模拟结果与已有研究结论进行对比分 析，验证了基于格子波尔兹曼法--离散元法耦合算 法的放矿模型在自然崩落法放矿松动体形态研究中 的可靠性与准确性． 1 非饱和矿岩颗粒间受力分析 崩落矿岩细颗粒间通常含有孔隙水，因此处于 非饱和状态的崩落矿岩散体为固--液--气三相耦合 场． 此时颗粒表面为润湿状态，相邻湿颗粒会在接 触点及其附近形成半月状液桥( 见图 1( a) ) ，进而 将颗粒连接在一起，故矿岩颗粒间的细观水力特性 主要取决于矿岩的含水率． 非饱和矿岩颗粒可简化 为湿颗粒的球体模型［20］，相邻颗粒间的作用力除包 括离子水化作用力、静电力、范德华力、双电层作用 力等分子尺度的作用力外，还包括颗粒间液桥产生 的毛细作用力［21］． 非饱和状态下，颗粒间的毛细作 用力主要来自于液体的表面的张力及气--液交界处 的静水压力． 毛细作用力的计算主要涉及表面张力 γ、颗粒半径 Ｒ、颗粒间距 S 及固--液接触角 θ 等因素 ( 见图 1( b) ) ，其中，目前毛细作用力通常采用 La￾place--Young 方程计算求解［21］，其计算表达式如下: Δp γ = Ｒy″ ( Ｒ2 + y'2 ) 3 /2 + Ｒ y ( Ｒ2 + y'2 ) 1 /2 ( 1) 式中，Δp 为气--液交界面的压力差，y' = cot( φ + θ) ， y″为 y 的二阶导数，φ 为半填充角． 相应的总液力计算方程为: Fn = 2πγＲsin φsin ( φ + θ) + πＲ2 Δp sin2 φ ( 2) 其中，2πγＲsin φsin ( φ + θ) 与 πＲ2 Δp sin2 φ 分别表 征表面张力与静水压力． 2 LBM--DEM 耦合算法 2. 1 LBM 方程非牛顿流体化 对于二维两相流问题的处理，目前格子波尔兹 曼法通常采用 D2Q9( 2 为空间维数，9 为离散速度 数) ，模型如图 2 所示． 该模型将整个流体域划分为 · 666 ·