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

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

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

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

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

1 LBM-DEM耦合框架

1.1 LBM非牛頓流體處理

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

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

為將碎石細(xì)顆粒擬作非牛頓流體，需對LBM方程進(jìn)行非牛頓流體化處理。在處理流固耦合的問題時，LBM方程標(biāo)準(zhǔn)形式為［15］：

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

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

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

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

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

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

針對IMB法，通過引入一個與格子固含率有關(guān)的加權(quán)函數(shù)β來修正LBM，其表達(dá)式為：

對應(yīng)式（7），修正后的LBM方程為：

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

2 流體與離散顆粒耦合力

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

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

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

式中，β為加權(quán)函數(shù)。

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

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

式中，?=Fp,Tf,up.

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

3.1 計算參數(shù)及數(shù)值模型

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

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

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

3.2 模擬結(jié)果及分析

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

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

根據(jù)圖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點。上述壓力分布特點原因在于放礦場內(nèi)顆粒介質(zhì)的離散性與隨機(jī)性主要通過顆粒間的摩擦與碰撞產(chǎn)生動量傳遞與能量耗散。臨近放礦口處，礦巖顆粒運動速度快，礦巖間的碰撞與摩擦作用明顯，而遠(yuǎn)離放礦口位置的顆粒運動速度較慢，礦巖間碰撞阻尼作用較弱。放礦場內(nèi)礦巖動力學(xué)響應(yīng)機(jī)理研究結(jié)果與文獻(xiàn)［5］、［18］等研究結(jié)論一致，并對研究放礦成拱原因、放礦場邊界壓力及控制底部結(jié)構(gòu)壓力等具有指導(dǎo)作用。

4 結(jié)論

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

（2）從放礦動力學(xué)角度對放礦場內(nèi)單個散體礦巖進(jìn)行受力分析，得出放礦過程中其主要受重力及顆粒間的碰撞與摩擦阻力作用。按照動量定理，推導(dǎo)出顆?？偰Σ亮傲赜嬎愎剑⑼ㄟ^對顆粒運動方程進(jìn)行離散求解，得出相應(yīng)的顆粒運動方程。

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

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