放礦場動力學響應機理研究的LBM-DEM耦合方法

常寶孟 杜翠鳳 韓 斌 李 濤 初曉峰

（北京科技大學土木與環(huán)境工程學院，北京100083）

自然崩落法具有生產能力大、便于組織管理、作業(yè)安全及開采成本低等優(yōu)點，是唯一能與露天開采經濟效益相媲美的高效地下采礦方法［1-4］。自然崩落法放礦是在崩落圍巖覆蓋條件下進行的，覆巖中的廢石極易摻雜于礦石中，并與礦石一起經放礦口放出，從而造成放礦貧化。此外，放礦控制不當可能造成礦巖崩落不均勻，易造成放礦口堵塞及井下空氣沖擊波等危害，嚴重影響井下人員與設備安全及礦山生產效益，因此有必要對放礦過程進行理論研究。

自然崩落法放礦過程的本質是崩落礦巖在放礦漏斗內運動的動力學過程［5］，國內外學者提出的橢球體放礦、類橢球體放礦及隨機介質放礦等放礦理論主要是基于放礦運動學方面的研究，其研究前提是將崩落礦巖視為理想散體，與礦山生產實際存在出入，因而具有一定的局限性和不可靠性。目前，國內外關于放礦動力學的研究主要有：賴森華等［5］在分析放礦過程中礦巖移動規(guī)律的基礎上，建立了礦巖塊體力學模型和相應的動力學方程，從理論上對放礦動力學過程進行了探討；王昌漢［6］根據實驗研究結果，建立了松散礦巖從漏斗口均勻放出時的力學方程，從散體力學角度對放礦力學體系進行了完善；F Melo等［7-8］基于Bergmark-Roos理論分析了放礦過程中礦巖塊體的受力狀態(tài)，建立了放出體極坐標運動方程。由于放礦過程本身的復雜性及研究手段的局限性，放礦動力學相關研究成果較少且發(fā)展緩慢。

針對自然崩落法放礦研究，離散元法（Discrete Element Method，DEM）是巖土力學中最常用的有效模擬手段。但由于受模擬顆粒數量及計算機性能的限制［9］，難以在全域上采用離散顆粒進行模擬。按照美國ASTM（Unified Soil Classification System）巖土顆粒粒級分類標準［10］，相較于粒徑單向尺寸大于1 m的礦巖，顆粒直徑小于19 mm的碎石可視為細顆粒。劉傳平等［11-12］基于“擬流體”思想，通過顆粒斜槽流實驗獲得了顆粒流黏性的本構關系，其數學模型預測值與實驗結果的相對誤差在13%以內，表明顆粒流擬非牛頓流體研究是可行的。而格子Boltzmann法（Lattice Boltzmann Method，LBM）作為一種新興的流體數值模擬方法，具有高并行性及能夠處理復雜邊界條件等顯著優(yōu)勢，近年來在多相多組分流、粒子懸浮流及多孔介質流等方面得到了廣泛的發(fā)展和應用［13-15］。根據上述理論與研究，本研究將自然崩落法放礦過程中大塊散體礦巖及碎石細顆粒分別視為離散體與非牛頓流體，并基于LBM-DEM耦合方法實現自然崩落法放礦過程模擬，從動力學角度對放礦場內崩落礦巖的力學響應機理進行分析，為放礦理論研究提供新的研究思路與技術支持。

1 LBM-DEM耦合框架

1.1 LBM非牛頓流體處理

在LBM粒子速度離散模型中，典型的二維D2Q9模型（如圖1所示）的離散速度ei為：

式中，C=dx/Δt為格子速度，Δt為時間步長。

為將碎石細顆粒擬作非牛頓流體，需對LBM方程進行非牛頓流體化處理。在處理流固耦合的問題時，LBM方程標準形式為［15］：

式中，（fix,t）為t時刻位于x處以離散速度ei運動的粒子分布函數；τ為無量綱松弛時間；fieq為局部平衡態(tài)方程。D2Q9模型對應的平衡態(tài)分布函數具有如下形式：

式中，ωi為數值積分公式的權系數，為格子聲速；u為速度。宏觀層面的流體密度ρ、速度u和壓力P可由粒子分布函數的零階和一階速度矩定義：

在低馬赫數假設前提下，通過Chapman-Enskog展開由式（2）格子Boltzmann方程可正確恢復出Navier-Stokes方程，并得到非牛頓流體黏性與模型參數的關系［16］為：

式中，v為非牛頓流體黏度；e為應變率；τ*與τ0分別為非牛頓流體與牛頓流體松弛時間。

1.2 流體—移動顆粒邊界耦合控制策略

要實現放礦場內流體與離散顆粒的耦合框架，重點在于解決細顆粒流與移動散體礦巖的邊界耦合問題。Noble等［17］提出針對LBM-DEM耦合的浸入運動邊界法（Immersed Moving Boundary，IMB）。為解決流固邊界耦合動量不連續(xù)性問題，并使移動顆粒受力計算更平順，引入新物理參數格子固含率ε，其取值取決于固體節(jié)點所覆蓋的節(jié)點單元比例，圖2中固體節(jié)點的格子固含率ε為填充區(qū)域與其邊框的面積比值，fi（xf）為流體節(jié)點，f-i（xs）為固體節(jié)點。

針對IMB法，通過引入一個與格子固含率有關的加權函數β來修正LBM，其表達式為：

對應式（7），修正后的LBM方程為：

通過在LBM修正方程中引用非牛頓流體松弛時間τ*，可解決碎石細顆粒流擬非牛頓流體及流體與移動顆粒邊界耦合的問題。

2 流體與離散顆粒耦合力

自然崩落法放礦過程中，散體礦巖顆粒的運動存在隨機性與突發(fā)性，難以對每一塊體每一時刻的運動狀態(tài)和力學作用過程進行精確描述和分析。根據崩落礦巖的散體流動性及放礦場的邊界條件，可對放礦場內礦巖顆粒在放礦過程中的受力狀態(tài)進行分析，如圖3所示。

放礦場內單個礦巖顆粒運動遵循牛頓第二運動定律，礦巖顆粒在重力mg作用下，與細顆粒流及其他相鄰礦巖顆粒間存在摩擦力Ff，且Ff為阻尼作用，阻礙礦巖塊體向放礦口運動。除此之外還存在碰撞作用，由于塊體間的碰撞作用存在隨機性，因而將其簡化為垂直于Ff方向的合力N。根據圖3的受力分析，放礦場內單個礦巖顆粒的動力學方程可表達如下：

式中，m和J分別為礦巖顆粒質量與慣性矩；θ?為角加速度；Ff和Tf分別為作用于礦巖顆粒上的總摩擦力與總力矩。流體粒子與固體顆粒交互作用時，兩者間將產生動量交換，根據流體粒子與固體顆粒之間的動量交換可計算作用力。按照動量定理，顆粒受到的總摩擦力和力矩分別為：

式中，β為加權函數。

獲得作用于顆粒上的力和力矩后，分別對顆粒的運動方程進行離散求解，散體顆粒運動方程為：

式中，Mp和Fp分別為顆粒質量和作用于顆粒上的總力；up為顆粒運動速度；xc為顆粒質心。其中，t時刻顆粒受到的作用力、力矩和顆粒速度三者平均值為：

式中，?=Fp,Tf,up.

3 算例：自然崩落法單漏斗放礦試驗

3.1 計算參數及數值模型

國內某礦山已探獲銅金屬量為480萬t，礦化帶長 2 300 m，寬 600～800 m，面積約 1.09 km2，呈穹窿狀。由于該礦具有巖層破碎、礦石品位較低、礦體相對較為厚大等特點，故采用自然崩落法全地下開采。開采區(qū)礦體及圍巖主要由石英二長斑巖與閃長玢巖組成，根據開采技術條件研究成果整理，礦體及圍巖力學參數見表1。

在生成自然崩落法放礦數值模型時，需對放礦場內的不同塊度礦巖進行分類處理：根據顆粒斜槽流實驗結果［12-13］將細顆粒流視為非牛頓流體，在時間與空間統(tǒng)計平均意義下，其細觀結構不需考慮，并通過LBM方法處理；大塊散體礦巖在放礦漏斗內的運動視為離散體運移，并通過DEM方法模擬其運動狀態(tài)。由于在實際放礦場中，崩落礦巖顆粒是非球形的，且塊度大小不均，此類復雜介質通常通過定義當量直徑Dp來表示其中顆粒的大小，Dp可根據不同的統(tǒng)計平均方法確定。在此統(tǒng)計平均意義下，崩落散體礦巖可表示為按高斯分布的直徑為Dp的二維圓盤介質（顆粒粒級服從高斯分布）?？紤]模型計算量和計算精度的要求，散體礦巖及流體域參數分別按表2、表3選取。

根據第1節(jié)中LBM-DEM耦合框架，分別將放礦場內細顆粒流與大塊崩落礦巖通過LBM與DEM模擬，利用Fortran語言編寫程序代碼，按礦山放礦場參數設置模型參數，并在Visual Studio平臺對算例代碼進行編譯。生成的單漏斗放礦耦合模型如圖4所示。

3.2 模擬結果及分析

圖5為編譯程序經Tecplot軟件處理后獲得的放礦場內不同時刻壓力分布云圖。從圖5中可以看到，在時間t分別為5 s、10 s及15 s時，放礦口周邊壓力較大，一直保持在3.002 5×107～3.007×107N水平，而放礦口上部壓力較小，在2.991×107N左右，且放礦場內壓力隨崩落礦巖的放出表現為脈動分布形式。其主要原因在于崩落礦巖從放礦口放出時存在松散和壓縮2種變形，放礦口上部為松散變形區(qū)，礦巖壓力較小，周邊為壓縮變形區(qū)，礦巖壓力較大。此外，液—固兩相流中顆粒相的細觀結構與宏觀力學特性均十分復雜，流場中離散顆粒是其運動的基本單元，在放礦過程中表現為隨機脈動的運動形式，而細顆粒流擬為非牛頓流體，其動力學特性在統(tǒng)計平均意義下具有宏觀上的連續(xù)性與規(guī)律性，因而會出現放礦場內壓力脈動分布現象。

為進一步研究放礦場內不同位置的壓力分布特征，分別以放礦口水平為x軸，以放礦口軸線為y軸建立直角坐標系，分別取A（-5，1）、B（0，10）、C（0，50）3點進行監(jiān)測，壓力監(jiān)測結果如圖6所示。

根據圖6可知，放礦口處監(jiān)測點A壓力在3.002 5×107～3.007×107N間波動，放礦口軸線上B、C2點壓力值波動區(qū)間比較接近，區(qū)間為2.993 5×107～2.998×107N，A點壓力值明顯大于B、C2點。上述壓力分布特點原因在于放礦場內顆粒介質的離散性與隨機性主要通過顆粒間的摩擦與碰撞產生動量傳遞與能量耗散。臨近放礦口處，礦巖顆粒運動速度快，礦巖間的碰撞與摩擦作用明顯，而遠離放礦口位置的顆粒運動速度較慢，礦巖間碰撞阻尼作用較弱。放礦場內礦巖動力學響應機理研究結果與文獻［5］、［18］等研究結論一致，并對研究放礦成拱原因、放礦場邊界壓力及控制底部結構壓力等具有指導作用。

4 結論

（1）將碎石細顆粒流與大塊崩落礦巖分別視為非牛頓流體與離散體，并分別通過LBM與DEM模擬，基于LBM-DEM耦合框架可以很好地實現自然崩落法單漏斗放礦模擬，同時有效降低離散元模擬規(guī)模，為放礦模擬提供了新的研究手段。

（2）從放礦動力學角度對放礦場內單個散體礦巖進行受力分析，得出放礦過程中其主要受重力及顆粒間的碰撞與摩擦阻力作用。按照動量定理，推導出顆?？偰Σ亮傲赜嬎愎?，并通過對顆粒運動方程進行離散求解，得出相應的顆粒運動方程。

（3）根據本研究提出的計算分析方法，得出礦巖放出過程中放礦口上部存在松散與壓縮2種變形，放礦口周邊壓力在3.002 5×107～3.007×107N水平，放礦口上部壓力在2.991×107N左右，且放礦場內壓力因散體礦巖的隨機脈動形式隨時間做脈動分布。

（4）通過設置不同監(jiān)測點，得出A、B、C3點處壓力值均在一定區(qū)間內波動，且A點壓力要明顯大于B、C2點，與放礦場內顆粒介質的離散性與隨機性主要通過顆粒間的摩擦與碰撞產生動量傳遞與能量耗散有關。