含水率对放矿松动体形态的细观影响

工程科学学报，第40卷，第6期：665672,2018年6月 Chinese Journal of Engineering,Vol.40,No.6:665-672,June 2018 DOI:10.13374/j.issn2095-9389.2018.06.003:http://journals.ustb.edu.cn 含水率对放矿松动体形态的细观影响 李 涛”，吴爱祥，王洪江”，尹升华”，冯云田2》 1)北京科技大学土木与资源工程学院，北京1000832)斯旺西大学工程学院，斯旺西SA28PP ☒通信作者，E-mail:ustblt(@163.com 摘要为从细观尺度研究矿岩含水率对自然崩落法放矿松动体形态的影响，对非饱和矿岩颗粒间的受力进行了分析，并分 别将放矿场内细颗粒流与大块离散矿岩利用格子波尔兹曼法与离散元法处理，基于格子波尔兹曼法一离散元法耦合算法建立 自然崩落法放矿模型，得出含水率与放矿松动体形态间的关系，并通过将模拟结果与已有研究结论进行对比分析，验证了基 于格子波尔兹曼法一离散元法耦合算法的放矿模型准确性及可靠性.研究表明：矿岩含水率对放矿松动体形态影响显著，在同 等刊矿岩放出质量分数情况下，随着含水率的增大，放矿松动体高度呈先增大后减小的趋势，放矿松动体形态先逐渐变为细长 型再逐渐恢复，放矿松动体形态变化的含水率临界值在10%左右. 关键词自然崩落法：放矿：含水率：非牛顿流体：松动体 分类号TD853 Influence of moisture content on the shape of isolated movement zone in mesoscale LI Tao,WU Ai-xiang,WANG Hong jiang",YIN Sheng-hua,FENG Yun-tian 1)School of Civil and Resource Engineering,University of Science and Technology Beijing,Beijing 100083,China 2)School of Engineering,Swansea University,Swansea SA2 8PP,UK Corresponding author,E-mail:ustblt@163.com ABSTRACT Block caving mining has the advantages of large production capacity,safe operation,low mining cost,and easy organi- zation and management.Therefore,it is considered as one of the natural alternatives to the current open cut mining.The migration law of ore-rock and the design of ore drawing parameters during ore drawing are the key to block caving mining.The caving rocks usually contain pore water between the fine particles,therefore,the caving rock in an unsaturated state is a solid-liquid-gas three-phase cou- pling field.The electric double layer force as well as the ionic hydration force,electrostatic force,and Van der Waals force result in the interaction between adjacent particles on a molecular scale.However,the interaction is also resulted in the force of the liquid bridge.To explore the influence of moisture content on the shape of an isolated movement zone in block caving method under ore draw- ing in mesoscale,the stress between unsaturated ore particles was analyzed.The fine particles flow and discrete bulks of ore-rocks were investigated using lattice Boltzmann method (LBM)and discrete element method (DEM),respectively.Based on the LBM-DEM cou- pling algorithm,the ore-drawing models of the block caving method were constructed,and the relationship between the moisture content and the shape of the isolated movement zone was obtained.The accuracy and reliability of the ore-drawing models based on the LBM- -DEM coupling algorithm were verified by comparing the simulation results with the existing research results.The results indicate that the moisture content of ore-rock has significant effects on the shape of the isolated movement zone.When the same ore-rock mass frac- tion is discharged,with an increase of moisture content,the height of isolated movement zone increases and then decreases,and the shape of isolated movement zone becomes slender and then is restored gradually.The critical value of moisture content during the change of isolated movement zone shape is about 10%. 收稿日期：2017-0805 基金项目：国家重点研发计划资助项目(2016YFC0600704:2016YFC0600709)

工程科学学报，第 40 卷，第 6 期: 665--672，2018 年 6 月 Chinese Journal of Engineering，Vol． 40，No． 6: 665--672，June 2018 DOI: 10． 13374 /j． issn2095--9389． 2018． 06． 003; http: / /journals． ustb． edu． cn 含水率对放矿松动体形态的细观影响 李 涛1) ，吴爱祥1) ，王洪江1) ，尹升华1) ，冯云田2) 1) 北京科技大学土木与资源工程学院，北京 100083 2) 斯旺西大学工程学院，斯旺西 SA2 8PP 通信作者，E-mail: ustblt@ 163． com 摘 要 为从细观尺度研究矿岩含水率对自然崩落法放矿松动体形态的影响，对非饱和矿岩颗粒间的受力进行了分析，并分 别将放矿场内细颗粒流与大块离散矿岩利用格子波尔兹曼法与离散元法处理，基于格子波尔兹曼法--离散元法耦合算法建立 自然崩落法放矿模型，得出含水率与放矿松动体形态间的关系，并通过将模拟结果与已有研究结论进行对比分析，验证了基 于格子波尔兹曼法--离散元法耦合算法的放矿模型准确性及可靠性． 研究表明: 矿岩含水率对放矿松动体形态影响显著，在同 等矿岩放出质量分数情况下，随着含水率的增大，放矿松动体高度呈先增大后减小的趋势，放矿松动体形态先逐渐变为细长 型再逐渐恢复，放矿松动体形态变化的含水率临界值在 10% 左右． 关键词 自然崩落法; 放矿; 含水率; 非牛顿流体; 松动体 分类号 TD853 收稿日期: 2017--08--05 基金项目: 国家重点研发计划资助项目( 2016YFC0600704; 2016YFC0600709) Influence of moisture content on the shape of isolated movement zone in mesoscale LI Tao1) ，WU Ai-xiang1)  ，WANG Hong-jiang1) ，YIN Sheng-hua1) ，FENG Yun-tian2) 1) School of Civil and Ｒesource Engineering，University of Science and Technology Beijing，Beijing 100083，China 2) School of Engineering，Swansea University，Swansea SA2 8PP，UK  Corresponding author，E-mail: ustblt@ 163． com ABSTＲACT Block caving mining has the advantages of large production capacity，safe operation，low mining cost，and easy organi￾zation and management． Therefore，it is considered as one of the natural alternatives to the current open cut mining． The migration law of ore-rock and the design of ore drawing parameters during ore drawing are the key to block caving mining． The caving rocks usually contain pore water between the fine particles，therefore，the caving rock in an unsaturated state is a solid-liquid-gas three-phase cou￾pling field． The electric double layer force as well as the ionic hydration force，electrostatic force，and Van der Waals force result in the interaction between adjacent particles on a molecular scale． However，the interaction is also resulted in the force of the liquid bridge． To explore the influence of moisture content on the shape of an isolated movement zone in block caving method under ore draw￾ing in mesoscale，the stress between unsaturated ore particles was analyzed． The fine particles flow and discrete bulks of ore-rocks were investigated using lattice Boltzmann method ( LBM) and discrete element method ( DEM) ，respectively． Based on the LBM--DEM cou￾pling algorithm，the ore-drawing models of the block caving method were constructed，and the relationship between the moisture content and the shape of the isolated movement zone was obtained． The accuracy and reliability of the ore-drawing models based on the LBM- -DEM coupling algorithm were verified by comparing the simulation results with the existing research results． The results indicate that the moisture content of ore-rock has significant effects on the shape of the isolated movement zone． When the same ore-rock mass frac￾tion is discharged，with an increase of moisture content，the height of isolated movement zone increases and then decreases，and the shape of isolated movement zone becomes slender and then is restored gradually． The critical value of moisture content during the change of isolated movement zone shape is about 10% ．

·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 ·

李涛等：含水率对放矿松动体形态的细观影响 ·667· b 图1非饱和状态下颗粒间液桥力模型.()半月状液桥：(6)毛细作用力影响因素 Fig.1 Model of liquid bridge forces between particles in unsaturated state:(a)semilunar liquid bridge:(b)capillary force influencing factors 图2LBM标准格子及D2Q9模型 Fig.2 Standard LBM lattice and D209 model 互相正交的网格，位于网格节点处的流体粒子具有 离散速度e:(i=0~8),其速度矢量表达式表征 o=wop(1-3u 2c24u} 如下： 3 2c(ew)2-3 9 f=ωP1+3 e;'u e:= (i=1,…,8) 0 (i=0) (5) (m(）a-) (6=1,2,3,4) 式中，p为流体密度，u为平衡态质点速度，w:为数 (m(-+）(-) 值积分公式的权系数： (i=5,6,7,8) 4 1 1 (6) (3) w0=9w1234=9w5.61.8=36 式中，c为流体粒子的迁移速度，且c=l/△t,l为格 崩落矿岩中的细颗粒流具有似非牛顿流体特 子间距，△t为时间步长 性四，在格子波尔兹曼法中，非牛顿流体的松弛时 每次时步迭代，流体质点按以下密度分布函数 间x,与牛顿流体的松弛时间0存在如下关系a: 分布： .=3+(-3）E (7) f(x+e:△t,t+△t)= 式中，E为非牛顿流体发生剪切流动时的应变率，n f(x,)-上(x,）-(x,] To 为幂律指数. (i=0,…,8) (4) 由式(4)与式(7)联立可得适用于非牛顿流体 式中：x为计算域中的一个格点；t为时间：f为速度 的密度分布函数方程： 分布函数；T。为松弛时间，量纲为一："(x,t)为局部 f(x+e:△i,t+△)= 平衡态分布函数，其计算表达式为： x,)-1(x,)-(x)] (8) T

李 涛等: 含水率对放矿松动体形态的细观影响 图 1 非饱和状态下颗粒间液桥力模型． ( a) 半月状液桥; ( b) 毛细作用力影响因素 Fig． 1 Model of liquid bridge forces between particles in unsaturated state: ( a) semilunar liquid bridge; ( b) capillary force influencing factors 图 2 LBM 标准格子及 D2Q9 模型 Fig． 2 Standard LBM lattice and D2Q9 model 互相正交的网格，位于网格节点处的流体粒子具有 离散速度 ei ( i = 0 ～ 8) ，其速度矢量表达式表征 如下: ei = 0 ( i = 0) (c ( cos π( i － 1) ) 2 ， ( sin π( i － 1) ) ) 2 ( i = 1，2，3，4) 槡2 (c ( cos π( i － 5) 2 + π ) 4 ， ( sin π( i － 5) 2 + π ) ) 4 ( i = 5，6，7，8        ) ( 3) 式中，c 为流体粒子的迁移速度，且 c = l /Δt，l 为格 子间距，Δt 为时间步长． 每次时步迭代，流体质点按以下密度分布函数 分布: fi ( x + eiΔt，t + Δt) = fi ( x，t) － 1 τ0 ［fi ( x，t) － f eq i ( x，t) ］ ( i = 0，…，8) ( 4) 式中: x 为计算域中的一个格点; t 为时间; fi 为速度 分布函数; τ0 为松弛时间，量纲为一; f eq i ( x，t) 为局部 平衡态分布函数，其计算表达式为: f eq 0 = ω0 ρ ( 1 － 3 2c 2u·u ) f eq i = ωiρ ( 1 + 3 c 2 ei ·u + 9 2c 2 ( ei ·u) 2 － 3 2c 2u·u ) ( i = 1，…，8        ) ( 5) 式中，ρ 为流体密度，u 为平衡态质点速度，ωi 为数 值积分公式的权系数: ω0 = 4 9 ; ω1，2，3，4 = 1 9 ; ω5，6，7，8 = 1 36 ( 6) 崩落矿岩中的细颗粒流具有似非牛顿流体特 性［22］，在格子波尔兹曼法中，非牛顿流体的松弛时 间 τ* 与牛顿流体的松弛时间 τ0 存在如下关系［16］: τ* = 1 2 + ( τ0 － ) 1 2 E' n － 1 ( 7) 式中，E'为非牛顿流体发生剪切流动时的应变率，n 为幂律指数． 由式( 4) 与式( 7) 联立可得适用于非牛顿流体 的密度分布函数方程: fi ( x + eiΔt，t + Δt) = fi ( x，t) － 1 τ* ［fi ( x，t) － f eq i ( x) ］ ( 8) · 766 ·

·668 工程科学学报，第40卷，第6期 2.2LBM-DEM耦合方法 权函数B来修正LBM,其表达式为： 根据离散元理论，放矿场中的大块离散颗粒可 ε(T-1/2) 采用考虑了接触力大小的软球模型，颗粒间的碰撞 B=1-E+(r-1D) (9) 及运动通过树形搜索算法计算P四.在LBM一DEM 式中，x为松弛时间，将式(9)代入式(8)，即可得适 耦合场中，处理流体与静止颗粒及流体与固定边界 用于非牛顿流体的LBM-DEM耦合修正方程： 情况时，通常假设流固交界点处流体与固体节点 f(x+e:△t,t+△t)=f(x,t）- (或固体边界)速度相等，利用反弹格式处理. (1-B®(x,)-f(x,0]+(10) T LBM-DEM耦合方法的关键在于对流体与移动颗粒 的耦合处理.为解决流固边界耦合动量不连续性问 式中，为修正后的密度分布函数.根据不同含水 率所对应的细颗粒流运动黏度，通过式(10)即可实 题，并使移动颗粒受力计算更平顺，Noble与Torc- 现自然崩落法放矿过程中矿岩细颗粒流与大块离散 zynski提出针对LBM-DEM耦合的浸入运动边界 矿岩的耦合计算 法(immersed moving boundary,IMB),并引入新物理 参数格子固含率ε（见图3），其取值取决于固体节 3数值计算及分析 点所覆盖的节点单元比例.图3中单个固体节点的 3.1数值计算模型 格子固含率ε即为粉色填充区域与绿色边框的面积 国内某矿山具有岩层破碎、矿石品位较低、矿体 比值，T与V。分别为颗粒角速度及平移速度， 相对较为厚大等特点，对采矿方法进行综合经济技 f(x)、f:(x)分别为相邻流体节点与固体节点的 术比较后，采用自然崩落法全地下开采.开采区矿 合力. 体及围岩主要由石英二长斑岩与闪长玢岩组成，其 颗粒表面 物理力学参数见表1.目前已探获铜金属量为480 万t,矿化带长2300m,宽600~800m,面积约1.09 km2,呈穹窿状，铜品位0.20%~3.74%，平均 0.57%.放矿口尺寸为5m×4.2m,最大放矿高度 为110m.本文将以该矿矿岩特性选取数值模型物 ●固体节点 理力学参数，研究结论主要是针对该类矿山的放矿 ●流体节点 松动体形态变化特性，但在含水率对放矿松动体形 态影响方面仍具有普适性.根据该矿放矿参数，分 T 别在不同含水率下生成单漏斗放矿模型，计算模型 图3IMB耦合模型及格子固含率定义 如图4所示.其中，图4(a)为单漏斗放矿模型网格 Fig.3 IMB scheme and definition of local solid ratio s 图，根据网格间距，模型内流体域划分为550×250 针对MB法，通过引入与格子固含率相关的加 个正方形网格 表1某铜矿岩石物理力学参数 Table 1 Rock physical and mechanical parameters in one copper mine 块体密度/ 弹性模量/ 单轴抗压/ 抗拉强度/ 内聚力/ 内摩擦角/ 岩石名称 (g-cm-3) 泊松比 GPa MPa MPa MPa () 石英二长斑岩 2.70 54.58 0.27 127.96 7.07 22.06 47.31 闪长玢岩 2.76 58.68 0.25 185.67 12.29 22.70 41.18 为研究矿岩含水率对放矿松动体形态的影响， 块离散矿岩通过LBM与DEM模拟，利用Fortran95 分别考虑矿岩含水率w为0、5%、10%及15%的情 程序设计语言编写算例代码，并在Microsoft Visual 况，与之对应的细颗粒流运动黏度v分别为0、 Studio平台对代码进行编译，通过Tecplot软件对模 0.18、0.37与0.20m2·s.模型中矿岩颗粒半径为 型计算结果进行后处理 0.5~1.25m,且颗粒粒径服从高斯分布.其他具体 3.2模型可靠性分析 模型参数见表2及表3.计算模型基于LBM一DEM Castro等通过大型物理放矿实验研究得出， 耦合算法生成，分别将放矿场内矿岩细颗粒流与大 在矿岩不含水的情况下（ω=0），放矿松动体高度

工程科学学报，第 40 卷，第 6 期 2. 2 LBM--DEM 耦合方法 根据离散元理论，放矿场中的大块离散颗粒可 采用考虑了接触力大小的软球模型，颗粒间的碰撞 及运动通过树形搜索算法计算［23］． 在 LBM--DEM 耦合场中，处理流体与静止颗粒及流体与固定边界 情况时，通常假设流固交界点处流体与固体节点 ( 或固体边界) 速度相等，利用反弹格式处理［17］． LBM--DEM 耦合方法的关键在于对流体与移动颗粒 的耦合处理． 为解决流固边界耦合动量不连续性问 题，并使移动颗粒受力计算更平顺，Noble 与 Torc￾zynski［24］提出针对 LBM--DEM 耦合的浸入运动边界 法( immersed moving boundary，IMB) ，并引入新物理 参数格子固含率 ε( 见图 3) ，其取值取决于固体节 点所覆盖的节点单元比例． 图 3 中单个固体节点的 格子固含率 ε 即为粉色填充区域与绿色边框的面积 比值，T 与 Vc 分 别为颗粒角速度及平移速度， fi ( xf ) 、f － i ( xs) 分别为相邻流体节点与固体节点的 合力． 图 3 IMB 耦合模型及格子固含率定义 Fig． 3 IMB scheme and definition of local solid ratio ε 针对 IMB 法，通过引入与格子固含率相关的加 权函数 β 来修正 LBM，其表达式为: β = ε( τ － 1 /2) 1 － ε + ( τ － 1 /2) ( 9) 式中，τ 为松弛时间，将式( 9) 代入式( 8) ，即可得适 用于非牛顿流体的 LBM--DEM 耦合修正方程: fi ( x + eiΔt，t + Δt) = fi ( x，t) － 1 τ* ( 1 － β) ［fi ( x，t) － f eq i ( x，t) ］+ βf m i ( 10) 式中，f m i 为修正后的密度分布函数． 根据不同含水 率所对应的细颗粒流运动黏度，通过式( 10) 即可实 现自然崩落法放矿过程中矿岩细颗粒流与大块离散 矿岩的耦合计算． 3 数值计算及分析 3. 1 数值计算模型 国内某矿山具有岩层破碎、矿石品位较低、矿体 相对较为厚大等特点，对采矿方法进行综合经济技 术比较后，采用自然崩落法全地下开采． 开采区矿 体及围岩主要由石英二长斑岩与闪长玢岩组成，其 物理力学参数见表 1． 目前已探获铜金属量为 480 万 t，矿化带长 2300 m，宽 600 ～ 800 m，面积约 1. 09 km2 ，呈 穹 窿 状，铜 品 位 0. 20% ～ 3. 74% ，平 均 0. 57% ． 放矿口尺寸为 5 m × 4. 2 m，最大放矿高度 为 110 m． 本文将以该矿矿岩特性选取数值模型物 理力学参数，研究结论主要是针对该类矿山的放矿 松动体形态变化特性，但在含水率对放矿松动体形 态影响方面仍具有普适性． 根据该矿放矿参数，分 别在不同含水率下生成单漏斗放矿模型，计算模型 如图 4 所示． 其中，图 4( a) 为单漏斗放矿模型网格 图，根据网格间距，模型内流体域划分为 550 × 250 个正方形网格． 表 1 某铜矿岩石物理力学参数 Table 1 Ｒock physical and mechanical parameters in one copper mine 岩石名称 块体密度/ ( g·cm － 3 ) 弹性模量/ GPa 泊松比 单轴抗压/ MPa 抗拉强度/ MPa 内聚力/ MPa 内摩擦角/ ( °) 石英二长斑岩 2. 70 54. 58 0. 27 127. 96 7. 07 22. 06 47. 31 闪长玢岩 2. 76 58. 68 0. 25 185. 67 12. 29 22. 70 41. 18 为研究矿岩含水率对放矿松动体形态的影响， 分别考虑矿岩含水率 ω 为 0、5% 、10% 及 15% 的情 况，与之对应的细颗粒流运动黏度 ν 分 别 为 0、 0. 18、0. 37 与 0. 20 m － 2·s． 模型中矿岩颗粒半径为 0. 5 ～ 1. 25 m，且颗粒粒径服从高斯分布． 其他具体 模型参数见表 2 及表 3． 计算模型基于 LBM--DEM 耦合算法生成，分别将放矿场内矿岩细颗粒流与大 块离散矿岩通过 LBM 与 DEM 模拟，利用 Fortran 95 程序设计语言编写算例代码，并在 Microsoft Visual Studio 平台对代码进行编译，通过 Tecplot 软件对模 型计算结果进行后处理． 3. 2 模型可靠性分析 Castro 等［3］通过大型物理放矿实验研究得出， 在矿岩不含水的情况下( ω = 0) ，放矿松动体高度 · 866 ·

李涛等：含水率对放矿松动体形态的细观影响 ·669· (a) (b) 图4单漏斗放矿模型.()网格模型：(b)几何模型 Fig.4 Isolated drawpoint model:(a)mesh model:(b)geometric model 表2流体域基本参数 Table 2 Basic parameters of the fluid domain 漏斗高， 漏斗宽， 放矿口宽， 漏斗倾角， 流体密度， 时间步长， 网格间距， him w/m wo/m 81() p/(kg'm-3) △ls lm 110 50 5 58 2760 1.0×10-2 0.2 表3矿岩颗粒基本参数 Table 3 Basic parameters of ore-ock 颗粒数量 R/m 摩擦系数，∫ 法向刚度，k/(Nm1) 颗粒密度，p/(kgm3) 1000 0.5-1.25 0.2 1×108 2700 (Hz)与放矿量(m)满足以下方程： 80 Hiz (m)=Ho(1-e-m/mu)+gm (11) ■模拟结果 一理论曲线 式中，H。mH及专均为常数系数，其中，H。及mH分 ◆ 别表示Hzm关系曲线呈指数形式增长阶段结束 50 时的松动体高度和矿岩放出质量，专为Hzm关系 40 曲线呈线性增长阶段的线性增长率. 30 为验证数值计算模型的可靠性及准确性，基于 10 Levenberg--Marquardt算法对w=0时的Hazm关系 0 曲线进行非线性拟合，其拟合方程为式(11)，具体 -10 参数拟合结果见表4（表中拟合参数右侧括号内的 0 150300450600750900 放矿量，mt 数据为相应的拟合误差值). 图5放矿松动体高度理论曲线与模拟数据对比（ω=0） 表4中ω=0时模拟方案的拟合优度2值接近 Fig.5 Comparison between the theoretic curve and simulation data 于1，且对应的放矿松动体高度理论曲线与模拟数 for the height of isolated movement zone (@=0) 据高度拟合（见图5），说明计算模型与放矿实际情 3.3不同含水率下放矿松动体形态对比分析 况比较一致，即验证了基于LBM一DEM耦合算法在 以矿岩放出质量分数8=18%时刻为例，不同 放矿模拟中的可靠性. 含水率下放矿漏斗内不同矿岩界面形态如图6所 表4式(11)参数拟合结果 示.根据图6可知，不同含水率下放矿场内细颗粒 Table 4 Fitted parameters in Eq.(11) 流的流动速度均存在由放矿口轴线向两侧逐渐减小 Holm m日作 5/(mt-1) 的现象，且放矿口处最大流速可达4.5mms.含 19.847(±1.64)147.6（±7.1）10（±0.59） 0.97 水率在0~15%范围内，随着含水率的不断增加，放

李 涛等: 含水率对放矿松动体形态的细观影响 图 4 单漏斗放矿模型． ( a) 网格模型; ( b) 几何模型 Fig． 4 Isolated drawpoint model: ( a) mesh model; ( b) geometric model 表 2 流体域基本参数 Table 2 Basic parameters of the fluid domain 漏斗高， h /m 漏斗宽， w/m 放矿口宽， w0 /m 漏斗倾角， θ /( °) 流体密度， ρ /( kg·m － 3 ) 时间步长， Δt / s 网格间距， l /m 110 50 5 58 2760 1. 0 × 10 － 2 0. 2 表 3 矿岩颗粒基本参数 Table 3 Basic parameters of ore-rock 颗粒数量 Ｒ/m 摩擦系数，f 法向刚度，kn /( N·m － 1 ) 颗粒密度，ρp /( kg·m － 3 ) 1000 0. 5 ～ 1. 25 0. 2 1 × 108 2700 ( HIMZ ) 与放矿量( m) 满足以下方程: HIMZ ( m) = H0 ( 1 － e － m/mH ) + ξm ( 11) 式中，H0、mH 及 ξ 均为常数系数，其中，H0 及 mH 分 别表示 HIMZ--m 关系曲线呈指数形式增长阶段结束 时的松动体高度和矿岩放出质量，ξ 为 HIMZ--m 关系 曲线呈线性增长阶段的线性增长率． 为验证数值计算模型的可靠性及准确性，基于 Levenberg--Marquardt 算法对 ω = 0 时的 HIMZ--m 关系 曲线进行非线性拟合，其拟合方程为式( 11) ，具体 参数拟合结果见表 4( 表中拟合参数右侧括号内的 数据为相应的拟合误差值) ． 表 4 中 ω = 0 时模拟方案的拟合优度 r 2 值接近 于 1，且对应的放矿松动体高度理论曲线与模拟数 据高度拟合 ( 见图 5) ，说明计算模型与放矿实际情 况比较一致，即验证了基于 LBM--DEM 耦合算法在 放矿模拟中的可靠性． 表 4 式( 11) 参数拟合结果 Table 4 Fitted parameters in Eq． ( 11) H0 /m mH /t ξ /( m·t － 1 ) r 2 19. 847( ± 1. 64) 147. 6( ± 7. 1) 10( ± 0. 59) 0. 97 图 5 放矿松动体高度理论曲线与模拟数据对比( ω = 0) Fig． 5 Comparison between the theoretic curve and simulation data for the height of isolated movement zone ( ω = 0) 3. 3 不同含水率下放矿松动体形态对比分析 以矿岩放出质量分数 δ = 18% 时刻为例，不同 含水率下放矿漏斗内不同矿岩界面形态如图 6 所 示． 根据图 6 可知，不同含水率下放矿场内细颗粒 流的流动速度均存在由放矿口轴线向两侧逐渐减小 的现象，且放矿口处最大流速可达 4. 5 mm·s － 1 ． 含 水率在 0 ～ 15% 范围内，随着含水率的不断增加，放 · 966 ·