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