裴英杰　肖慶飛　石貴明　周強

摘要：球磨機作為礦石磨碎階段的重要設備，對選礦廠提效降耗起著決定性作用。球磨機功耗與能耗分布、技術(shù)參數(shù)與襯板結(jié)構(gòu)參數(shù)設計優(yōu)化是影響球磨機磨礦效果的重要內(nèi)容和技術(shù)難點，離散元法可以準確模擬顆粒運動的形變過程并記錄顆粒的功率及能量分布，在球磨機仿真領(lǐng)域具有良好的應用前景。介紹了離散元法在球磨機仿真過程中應用比較廣泛的顆粒模型、接觸模型及聯(lián)合仿真方法各自的特點及研究進展，以及離散元法在球磨機功率及能量分布研究、技術(shù)參數(shù)優(yōu)化、襯板結(jié)構(gòu)參數(shù)優(yōu)化等方面的應用情況，歸納總結(jié)并指出了球磨機仿真模型精細化、聯(lián)合仿真廣泛應用、球磨機技術(shù)參數(shù)智能調(diào)節(jié)是未來球磨機仿真研究的發(fā)展趨勢。

關(guān)鍵詞：離散元法;球磨機;仿真;顆粒模型;接觸模型;耦合方法

中圖分類號：TD453文獻標志碼：A開放科學（資源服務）標識碼（OSID）：

文章編號：1001-1277（2022）01-0069-09doi：10.11792/hj20220112

基金項目：國家自然科學基金地區(qū)科學基金項目（51964044）;礦冶過程自動控制技術(shù)國家（北京市）重點實驗室開放基金項目（201905）;安徽省重點研究和開發(fā)計劃面上攻關(guān)項目（201904a07020044）;云南省教育廳課題（2019J0738）

引言

球磨機作為礦石磨碎階段的重要設備，具有結(jié)構(gòu)簡單、適應性強等優(yōu)點，現(xiàn)已廣泛應用于國內(nèi)外各類礦山[1]。近年來，隨著礦石“貧細雜”程度不斷加深，選礦廠生產(chǎn)對球磨機設計提出了更高的要求。若球磨機參數(shù)設計不合理，會對磨礦效果產(chǎn)生消極作用，嚴重影響后續(xù)選別指標，并加劇鋼材及電力的損耗[2-3]。為此，學者們進行了一系列的優(yōu)化研究。在對球磨機研究的過程中，采用的研究手段主要分為2類：球磨機試驗與模擬仿真[4-5]。在最初的研究中，基本研究過程是在理論計算后進行球磨機磨礦試驗。近年來，隨著計算機技術(shù)快速發(fā)展，模擬仿真方法越來越多地被應用到球磨機的研究中。目前，球磨機試驗和模擬仿真這2種主要研究手段互相驗證、相輔相成，共同推動球磨機研究領(lǐng)域的快速發(fā)展。

離散元法（Discrete Element Modelling，DEM）可以跟蹤球磨機中單一顆粒的運動并準確模擬顆粒之間的碰撞，已成為球磨機仿真過程最常用的數(shù)值模擬方法[6]。為此，眾多學者對球磨機功率及能量分布、技術(shù)參數(shù)及襯板結(jié)構(gòu)參數(shù)優(yōu)化等方面進行了廣泛的模擬仿真研究，獲得了大量研究成果，對后續(xù)研究人員了解球磨機內(nèi)部的顆粒運動及優(yōu)化各類參數(shù)具有積極的促進作用[7-9]。然而，受計算機能力限制，大多數(shù)仿真研究均經(jīng)過模型簡化處理，導致模擬過程存在一定誤差[10]。適宜的建模仿真可在一定程度上降低仿真誤差，提高球磨機仿真精確性?；诖?，本文針對離散元法在球磨機仿真中的研究進展及應用進行了歸納總結(jié)，并提出了今后離散元法在球磨機仿真研究中的發(fā)展方向，以期為球磨機仿真研究提供一定的研究思路和理論指導。

1離散元法在球磨機仿真中的研究

離散元法于1971年由CUNDALL[11]基于分子動力學理論提出，最初被用來研究單顆粒與離散顆粒群之間的交互作用及其穩(wěn)定性，由于該方法可以模擬巖體非均質(zhì)、不連續(xù)和大變形等特點，逐漸開始應用于巖土力學領(lǐng)域。1992年，離散元法被MISHRA等[12]首次應用到球磨機襯板結(jié)構(gòu)與功耗之間關(guān)系的研究中，但受到當時技術(shù)水平及計算機能力的限制，僅能進行二維建模仿真。1997年，AGRAWALA等[13]采用離散元法對比分析試驗測量與預測所得球磨機功率的差異性，驗證了三維離散元法用于球磨機仿真研究的可行性，為離散元法在球磨機仿真中的應用奠定了良好基礎。

1.1顆粒模型

顆粒模型是指根據(jù)介質(zhì)和礦石顆粒的形狀及特性建立的仿真模型。自19世紀80年代以來，學者們提出包括球體、橢球體[14]、超二次曲面[15]、圓柱體[16]、多面體[17]、多聚球體[18]及多聚橢球體[19]等不同形狀的離散元顆粒模型。為了達到不同的仿真要求，學者們選擇的顆粒模型也各有差異。離散元法模擬仿真應用的部分顆粒模型見表1。

目前，在球磨機仿真中最常用的顆粒模型主要是球體、多聚球體這2類。球體顆粒模型由于具有計算成本低的優(yōu)勢，常作為球磨機離散元仿真中的首選顆粒模型，而多聚球體顆粒模型具有可模擬真實顆粒形狀的優(yōu)點，常被用來模擬實際礦物，二者最大的差異表現(xiàn)在仿真時間及仿真精度上。BIBAK等[21]通過粒子填充算法對顆粒建模（步驟見圖1），探究多聚球體顆粒模型和球體顆粒模型動力學之間的差異，研究發(fā)現(xiàn)多聚球體顆粒在高能沖擊作用下的參與度較球體顆粒高5 %，多聚球體顆粒的模擬時間成本較球體顆粒增加35倍。BARRIOS等[27]為鐵礦石球團建立了球體顆粒模型和多聚球體顆粒模型（4個球體），對比分析2種預測模型的休止角及磨礦功率。結(jié)果表明：采用多聚球體顆粒模型仿真測得的休止角角度略低于試驗測得的休止角角度，預測磨礦功率與實測磨礦功率差值小于6 %;當仿真模型為球體顆粒模型時，休止角預測值誤差過高，預測磨礦功率與實測磨礦功率偏差高達30 %。因此，在仿真條件允許的情況下，采用多聚球體顆粒模型仿真得到的試驗數(shù)據(jù)及結(jié)論準確性更高。

在實際磨礦過程中，可能存在數(shù)十億個不同粒度的礦粒在球磨機中同時運動，根據(jù)現(xiàn)階段計算機性能，在離散元仿真中運行數(shù)千萬個粒子頗具難度，模擬實際尺寸及數(shù)量的礦石顆粒幾乎沒有可能。IWASAKI等[28]研究發(fā)現(xiàn)，當球磨機內(nèi)沖擊能量的頻率分布及其在研磨時間內(nèi)的累積量在不同尺寸球磨機中一致時，可以成功實現(xiàn)球磨機及顆粒放大，節(jié)約仿真成本。POWELL等[29]將球磨機筒體長度按比例縮小，球磨機中顆粒數(shù)量由450萬降至11萬，模擬時間大幅降低，有效節(jié)約了仿真成本。因此，在確保仿真精度的前提下，對球磨機及顆粒進行適當縮放，將極大程度地減少離散元仿真所消耗的時間并有效提高計算效率。

1.2接觸模型

離散元仿真過程即顆粒模型運動過程，顆粒運動時必然會引發(fā)顆粒之間的相互碰撞，對顆粒之間接觸作用研究成為離散元仿真過程的重中之重。目前，在離散元仿真中主要采用的2種接觸模型是硬球模型和軟球模型。硬球模型的思想是ALDER等[30]在分子系統(tǒng)的相變研究中提出的，該模型認為顆粒間的碰撞是兩體瞬時碰撞，且顆粒具有絕對剛度，但其常用于氣體輸送等稀疏顆粒系統(tǒng)，不適用于球磨機這類顆粒連續(xù)接觸的致密顆粒系統(tǒng)。目前，球磨機離散元仿真中最常用的是由CUNDALL等[31]提出的軟球模型，該模型通過對受力后發(fā)生的形變及斷裂進行分析來處理顆粒間的接觸力，顆粒間的接觸為柔性接觸，允許存在少量重疊，適用于準靜態(tài)顆粒材料，對顆粒數(shù)量巨大的球磨機離散元仿真系統(tǒng)具有計算效率上的優(yōu)勢。

顆粒接觸時接觸力的計算方法，即顆粒接觸本構(gòu)模型，是離散元仿真顆粒接觸模型的理論核心。目前，應用最廣泛的離散元顆粒接觸本構(gòu)模型是CUNDALL等[31]提出的線性彈簧-阻尼器接觸模型（見圖2），該模型在接觸點假設一組彈簧、阻尼器和滑塊組件，用來表示2個顆粒之間的碰撞，同時考慮了離散顆粒之間接觸力和非接觸力的相互作用。在后續(xù)的離散元仿真研究中，學者們提出了多種不同的本構(gòu)模型。JOHNSON等[32]基于Hertz-Mindlin模型，為顆粒添加一個額外的法向內(nèi)聚力，并考慮接觸物體的表面黏附能對接觸半徑的影響，提出了JKR模型。WALTON等[33]采用部分鎖定彈性形變模型計算法向力，并依據(jù)MINDLIN等[34]的理論在切向加入了摩擦力對彈性體的作用，提出了半鎖彈簧模型。IWASHITA等[35-36]分析了滾動摩擦機制，提出了抗?jié)L動的離散元接觸模型，其便于觀察顆粒滾動的影響。

雖然顆粒接觸本構(gòu)模型眾多，但并非全部都可以契合球磨機仿真研究。為了開發(fā)更適合于球磨機仿真的顆粒接觸本構(gòu)模型，有學者開展了相關(guān)研究。 IWASAKI等[9]忽略球碰撞時滾動摩擦的影響，提出了一種考慮球磨機轉(zhuǎn)速率與球摩擦系數(shù)相關(guān)性的三維濕式球磨機接觸力計算方法，并驗證了該模型的有效性。FENG[37-38]基于接觸體積模型開發(fā)出可以表示任意顆粒形狀的節(jié)能法向接觸模型，該模型可以自動確定任意顆粒形狀的完整法向接觸模型，包括接觸法線方向、接觸點（線）及作用力大小，無需任何其他假設，但使用該模型進行DEM仿真的整體計算效率并沒有實際數(shù)據(jù)可以驗證，需要后續(xù)學者進行相關(guān)工作。

1.3離散元耦合方法

球磨機仿真研究中將離散元法作為最常用的數(shù)值分析手段，但僅使用單一離散元法局限性較大。因此，學者們開發(fā)了不同的離散元耦合方法，其中最常用的是將離散元法與有限元（FEM）、計算流體力學（CFD）及光滑粒子流體力學（SPH）等方法進行耦合。由于以上方法均有各自的特點，它們與離散元法耦合后進行的研究也有所不同。

離散元法與有限元（DEM-FEM）耦合法通常采用DEM分析球磨機內(nèi)礦石與介質(zhì)的運動特性，采用FEM分析球磨機的襯板應力及扭矩。DEM-FEM耦合法分為單向耦合法和雙向耦合法，目前單向耦合法已趨于成熟，而雙向耦合法仍在進一步研究，因此大部分研究常采用單向耦合法分析磨礦介質(zhì)、礦石顆粒和襯板提升條之間的相互作用。李昆塬等[39]采用DEM-FEM耦合法對Φ4.3 m球磨機進行模擬仿真，分析了球磨機中不同面角、高度和頂邊長度的梯形襯板總變形程度與應力的關(guān)系。JONS?N等[40]使用FEM對球磨機中的柔性橡膠襯板及提升條進行建模，使用DEM對磨礦介質(zhì)進行建模，對比分析了仿真預測的撓度值與實際撓度值的關(guān)系。研究發(fā)現(xiàn)，DEM-FEM耦合法仿真獲得的撓度曲線與測量值具有一致性，為后續(xù)學者采用DEM-FEM耦合法模擬球磨機中可形變的橡膠襯板提供了理論基礎。

濕磨作業(yè)仿真需要對礦漿性質(zhì)進行模擬，且由于入料礦石的粒徑較小，無法使用單一離散元法對大量漿料進行建模。離散元法與計算流體力學（DEM-CFD）耦合法則可以通過描述物質(zhì)經(jīng)過網(wǎng)格坐標點時的密度、能量等來表示礦漿運動，且仿真精度高。BEINERT等[41]使用DEM-CFD耦合法對行星式球磨機進行數(shù)值研究，分析了磨礦介質(zhì)的運動、每個時間段的接觸及接觸時的動能。MAYANK等[42]采用CFD對連續(xù)礦漿流建模，采用DEM對離散礦粒進行建模，通過DEM-CFD耦合法對球磨機內(nèi)的礦石顆粒及礦漿進行模擬仿真，發(fā)現(xiàn)仿真結(jié)果與基于PEPT試驗結(jié)果吻合度高。但是，該試驗未考慮礦石顆粒與礦漿之間的相互作用，只考慮了礦石顆粒對礦漿的影響，模型需要進一步完善。

在SPH中粒子之間的關(guān)系根據(jù)當前時刻粒子的位置確定，其隨時間變化而改變，因此該方法亦可用于流體等形變較大的物體，且不會產(chǎn)生類似于拉格朗日網(wǎng)格法（如FEM）嚴重的網(wǎng)格變形導致數(shù)值誤差大的問題。離散元法與光滑粒子流體力學（DEM-SPH）耦合法也常被用于研究礦石顆粒、磨礦介質(zhì)與礦漿之間的交互作用。SINNOTT等[43]使用DEM對礦石顆粒和介質(zhì)建模，使用SPH對礦漿建模，提出一種球磨機中磨礦介質(zhì)和礦漿完全耦合流動的預測模型，為后續(xù)學者研究球磨機中的固相和液相提供了理論基礎。

為了更好地對球磨機運動進行仿真研究，JONSN等[44]將SPH、DEM和FEM 3種方法進行耦合，用DEM和SPH粒子表示介質(zhì)及礦漿，同時使用FEM對球磨機結(jié)構(gòu)進行建模，通過觀察礦漿黏度和密度的差異對扭矩及礦漿壓力分布的影響，考察球磨機模型及礦漿模型的可靠性，結(jié)果表明SPH-DEM-FEM耦合法為預測球磨機中高能區(qū)的體積及優(yōu)化球磨機襯板提升條和運行參數(shù)提供了可能。

2離散元法在球磨機研究中的應用

隨著離散元法研究的不斷深入，其理論體系日漸成熟，該方法在球磨機研究中也得到了更加廣泛的應用，其憑借非連續(xù)性的特點，將理論公式或常規(guī)試驗無法解釋的內(nèi)容更加直觀準確地表現(xiàn)出來。離散元法在球磨機研究中主要應用在功率及能量分布研究、技術(shù)參數(shù)優(yōu)化、襯板結(jié)構(gòu)參數(shù)優(yōu)化等方面，具體研究流程見圖3。

2.1功率及能量分布研究

球磨機運行過程通常需要消耗大量能量，由于磨礦過程中礦石對磨礦介質(zhì)的破碎行為是隨機的，磨礦介質(zhì)可能與礦石發(fā)生碰撞，也可能與其他介質(zhì)或襯板發(fā)生碰撞。當磨礦介質(zhì)與其他介質(zhì)或襯板發(fā)生碰撞的概率較高時，整個磨礦過程中大部分能量都在做無用功，導致能量利用率偏低，造成能源浪費。球磨機的功率及能量分布是由磨礦過程中多種因素決定的，襯板結(jié)構(gòu)及磨礦技術(shù)參數(shù)均會對其產(chǎn)生影響，研究不同因素對磨礦過程中功率及能量分布的影響，不僅可以優(yōu)化產(chǎn)品指標，對選礦廠節(jié)能降耗也至關(guān)重要。

KANO等[45]探究了球磨機磨礦過程中的沖擊能量與磨礦速率常數(shù)之間的相關(guān)性，發(fā)現(xiàn)能量利用率與磨礦速率成正比。CLEARY等[46]對實驗室球磨機建模，對比分析了不同料球比對球磨機功率消耗及能量利用方式的影響。研究發(fā)現(xiàn)，隨著料球比的增加，礦石磨礦消耗能量比例呈線性增加，且隨著礦石填充量的增加，磨礦產(chǎn)品粒度變得更細。IWASAKI等[9，47]基于離散元法探究了球磨機不同轉(zhuǎn)速率與介質(zhì)充填率對球磨機消耗能量的影響。研究發(fā)現(xiàn)，球磨機消耗能量不僅與介質(zhì)充填率成正比，而且隨球磨機速度的增加呈線性增加。

傳統(tǒng)的能量評估指標均假定球沿軸向均勻分布，然而實際情況并非如此，若假定球沿軸向均勻分布無法準確地反映球磨機的磨礦效果。HUANG等[48]將試驗研究和離散元法相結(jié)合，通過在球磨機中加入不同尺寸的鋼球，發(fā)現(xiàn)球磨機中大小不同的鋼球往往會彼此分開，而相同類型的鋼球會沿軸聚集（見圖4）。因此，有必要通過區(qū)域總能量來研究球磨機的能量分布情況，以分析球磨機中不同區(qū)域的磨礦效果。

學者們常將碰撞能量作為縮放仿真模型的參照，以達到提高仿真效率的目的。依據(jù)能量對球磨機筒體按比例縮放的方法，最早是由BOND[49]提出的基于能量的經(jīng)驗性縮放球磨機尺寸的方法。但是，該方法缺乏理論基礎，最終磨礦產(chǎn)品粒度分布預測值與實際值存在約20 %的誤差。學者們基于碰撞能量對縮放模型的研究均是在該模型的基礎上不斷完善，隨著研究的不斷深入，球磨機仿真模型的縮放模型研究已基本成熟。JIANG等[10]建立了不同半徑球磨機中顆粒運動的DEM模型，對單個粒子的運動軌跡和能量變化進行分析，發(fā)現(xiàn)采用減小模型尺寸并改變顆粒粒度的方法，可以減小模擬計算的規(guī)模并提高計算效率。SHI等[50]開發(fā)了一種基于能量的球磨機尺寸減小模型，該模型采用基于能量的特定方程確定磨礦產(chǎn)品粒度分布，且充分考慮球磨機的幾何形狀、裝料量及球磨機的生產(chǎn)能力等運行條件的影響。IWASAKI等[28]運用DEM對鋼球的沖擊能量進行數(shù)值分析，探究球磨機直徑與轉(zhuǎn)速率對沖擊能量的影響。研究發(fā)現(xiàn)，在不同尺寸的球磨機中，當沖擊能量的頻率分布與能量累積量一致時，可以實現(xiàn)球磨機及顆粒的放大。

2.2技術(shù)參數(shù)優(yōu)化

在球磨機運行過程中，影響磨礦效果的因素眾多，但技術(shù)參數(shù)最易調(diào)整，通過優(yōu)化技術(shù)參數(shù)改善磨礦效果是成本最低的方式。分析技術(shù)參數(shù)對球磨機中介質(zhì)運動狀態(tài)及磨礦產(chǎn)品粒度組成的影響，對降低選礦廠能耗及提高產(chǎn)品質(zhì)量具有重要意義。離散元法對球磨機的不同技術(shù)參數(shù)進行微觀模擬，可以更直觀地了解球磨機內(nèi)部運動情況，有助于對球磨機技術(shù)參數(shù)進行優(yōu)化，具體優(yōu)化研究流程見圖5。

OROZCO等[51]使用鍵合單元法開展離散元仿真研究，以礦石平均粒度和比表面積為評價指標，對比分析了鋼球數(shù)量與球磨機轉(zhuǎn)速率對磨礦過程的影響（見圖6）。研究發(fā)現(xiàn)，磨礦速度隨相同尺寸鋼球數(shù)量的增加而增大，但鋼球數(shù)量過多會引起鋼球之間的非彈性碰撞，增加能耗、降低磨礦效率。李騰飛等[8]通過DEM對球磨機建模，研究不同轉(zhuǎn)速率下球磨機內(nèi)介質(zhì)的運動狀態(tài)及不同粒度礦石與介質(zhì)之間的碰撞情況，發(fā)現(xiàn)球磨機內(nèi)介質(zhì)的拋落角、下落高度、落回瞬間的動能及顆粒之間的碰撞能均與球磨機轉(zhuǎn)速率有關(guān)，且轉(zhuǎn)速率越高，越有利于破碎大尺寸及較硬礦石。d/d0—產(chǎn)品與給礦平均粒徑比值S/S0—產(chǎn)品與給礦平均比表面積比值n—球磨機轉(zhuǎn)數(shù)Nb—鋼球數(shù)量

球磨機磨礦技術(shù)參數(shù)之間大多有交互作用，僅對一種參數(shù)進行分析通常很難達到優(yōu)化的目的，因此需要多方面考慮。CLEARY[52]采用DEM探究了充填率、介質(zhì)及礦石粒度分布對球磨機功耗及礦石運動狀態(tài)的影響。孫鵬[53]采用DEM研究了充填率、轉(zhuǎn)速率、鋼球尺寸及料球比對球磨機單位功耗和介質(zhì)與物料間沖擊能量的影響。劉剛等[54]針對轉(zhuǎn)速率、充填率、球料比等參數(shù)，以磨礦工作效率、破碎能力為性能指標，采用均勻試驗方法，對球磨機內(nèi)的物料進行了仿真運動分析及技術(shù)參數(shù)優(yōu)化。

球磨機技術(shù)參數(shù)優(yōu)化的另一個研究熱點是球磨機預測模型。造成磨礦效果差的部分原因是缺乏監(jiān)控及預測球磨機磨礦效果的適當方法，使得磨礦過程很難被控制。而離散元法能夠有效預測球磨機磨礦效果及能耗[18]。PEDRAVES等[55]為有效控制球磨機磨礦過程，充分考慮球磨機轉(zhuǎn)速率、待磨顆粒性質(zhì)及充填率等因素，對DEM仿真生成的扭矩數(shù)據(jù)頻譜進行分析，研究發(fā)現(xiàn)，球磨機的負載扭矩信號足以體現(xiàn)球磨機的負載水平，并提出了一種基于扭矩數(shù)據(jù)頻譜分析球磨機充填率的方法。DE CARVALHO等[56]以礦石粒度分布、轉(zhuǎn)速率、充填率、球磨機直徑及介質(zhì)尺寸為自變量，建立預測礦石破碎率的函數(shù)模型，研究發(fā)現(xiàn)，使用該模型的預測結(jié)果與HERBST等[57]提出的按比例關(guān)系放大程序得到的結(jié)果一致，預測效果明顯。

2.3襯板結(jié)構(gòu)參數(shù)優(yōu)化

球磨機襯板是礦石及介質(zhì)提升的動力來源，球磨機的襯板結(jié)構(gòu)對礦石顆粒的運動狀態(tài)及磨礦效果有很大的影響。襯板的結(jié)構(gòu)種類眾多，根據(jù)材質(zhì)和形狀不同主要分為8類（見圖7）：波形襯板、梯形襯板、階梯形襯板、楔形襯板、平襯板、波紋形襯板、橡膠襯板、磁性襯板。在實際生產(chǎn)中，球磨機襯板提升條形狀及數(shù)量等參數(shù)設計不合理時，會導致球磨機出現(xiàn)磨礦效率低、鋼球消耗高及磨礦產(chǎn)品質(zhì)量差等問題，對產(chǎn)品指標產(chǎn)生惡劣影響。

ROSALESMAR?N等[58]對4種不同面角襯板提升條進行離散元仿真，分析球磨機襯板幾何形狀對礦石破碎率和襯板提升條磨損率的影響。CHIMWANI等[59]采用離散元法對球磨機襯板形狀進行了優(yōu)化，發(fā)現(xiàn)沿球磨機軸向逐漸改變襯板提升條形狀可以實現(xiàn)不同的破碎及磨礦要求。WU等[60]分別對不同尺寸球磨機的襯板磨損情況開展仿真試驗研究，發(fā)現(xiàn)大型球磨機中鋼球?qū)σr板的沖擊作用與小型球磨機相差無幾。張學東等[61]基于仿生耦合理論和磨粒磨損相關(guān)理論，依據(jù)生物耐磨特征的尺寸參數(shù)和比例關(guān)系，分別對不同的非光滑襯板提升條形狀和現(xiàn)有襯板提升條形狀開展仿真對比試驗，發(fā)現(xiàn)建立的生物非光滑襯板提升條較現(xiàn)有襯板提升條耐磨效果更好。

當現(xiàn)場襯板形狀設計不合理或者磨損嚴重導致球磨機磨礦效果差時，若無法對球磨機內(nèi)的產(chǎn)品進行實時觀測，易造成嚴重的能源浪費，影響磨礦產(chǎn)品指標。因此，針對球磨機內(nèi)襯板對介質(zhì)及礦石顆粒的影響，實時精確控制球磨機內(nèi)的磨礦過程及對襯板磨損進行預測具有重要意義。POWELL等[29]將磨礦速度作為預測球磨機襯板提升條磨損的影響因素，提出一種基于磨礦速度預測球磨機襯板磨損演變的方法。BOEMER等[62]對球磨機的軸向切片進行三維DEM仿真，開發(fā)了一種基于DEM預測襯板表面磨損演變過程的程序，使球磨機襯板仿真研究變得更直觀，但襯板形狀變化的精度受到計算機能力的限制。REZAEIZADEH等[63]基于影響襯板磨損的主要參數(shù)，開發(fā)了一種預測襯板磨損率的模型，該模型能根據(jù)球磨機物料的密度、幾何特性和球磨機轉(zhuǎn)速率的變化，預測不同襯板提升條的使用壽命。

綜上可見，國內(nèi)外學者應用離散元法在功率及能量分布研究、技術(shù)參數(shù)優(yōu)化、襯板結(jié)構(gòu)參數(shù)優(yōu)化等方面積累了大量研究經(jīng)驗。在中國選礦廠提效降耗的背景下，該方法在球磨機仿真領(lǐng)域具有良好推廣應用前景。

3離散元法在球磨機仿真中的發(fā)展趨勢

目前，全球?qū)ΦV石的需求量仍居高不下，礦產(chǎn)資源的開發(fā)利用逐漸趨向于集約化、規(guī)?；椭悄芑?，顛覆傳統(tǒng)的生產(chǎn)技術(shù)與研究方法不斷革新。如前文所述，離散元法在球磨機設計優(yōu)化中具有不可替代的作用，其加快了球磨機的發(fā)展進程及新產(chǎn)品研制的速度。隨著計算機行業(yè)的快速發(fā)展及各種先進技術(shù)的應用，該方法必將展示出更強大的設計指導能力及試驗替代能力。但是，目前該方法在球磨機仿真領(lǐng)域仍處于探索階段，數(shù)據(jù)質(zhì)量、測量精度等技術(shù)方面的短板依然飽受業(yè)界質(zhì)疑，本文參照眾多學者的研究結(jié)果，結(jié)合當前離散元法在球磨機仿真中的研究難題，梳理出在球磨機研究領(lǐng)域離散元法的模型精度、仿真方式及技術(shù)參數(shù)控制3個方面的發(fā)展趨勢：

1）球磨機仿真模型的精細化發(fā)展。目前，球磨機仿真受到計算機能力限制，大多仿真研究對礦石顆粒的運動及襯板結(jié)構(gòu)建立了簡化模型。然而，隨著計算機處理能力的提高、軟件算法的優(yōu)化，球磨機及磨礦系統(tǒng)的仿真模型將更加符合實際，具體表現(xiàn)為礦石顆粒形狀及顆粒間的接觸更符合實際、礦石顆粒粒度分布更確切、球磨機結(jié)構(gòu)部件更真實、仿真效果更適合指導現(xiàn)場。

2）多種仿真方法聯(lián)合仿真的廣泛應用。根據(jù)現(xiàn)階段研究，離散元法與其他數(shù)值仿真方法耦合的目的在于提高仿真精度、對球磨機進行綜合分析。多種仿真方法聯(lián)合仿真可以有效解決單一方法的局限性，對礦漿及球磨機結(jié)構(gòu)進行綜合研究，獲得更加準確的仿真結(jié)果，為球磨機發(fā)展開拓新的視角。

3）球磨機技術(shù)參數(shù)的智能調(diào)節(jié)。目前，離散元法在球磨機中的仿真研究主要集中于能耗、技術(shù)參數(shù)及結(jié)構(gòu)性能的優(yōu)化研究。仿真模擬局限于球磨機的性能預測、襯板磨損失效及磨礦機理的研究，但對球磨機運行時的參數(shù)自適應調(diào)節(jié)研究較少。工業(yè)大數(shù)據(jù)、人工智能、數(shù)字孿生等智能技術(shù)能夠為球磨機仿真研究提供新思想、新思路。通過建立球磨機磨礦參數(shù)的智能調(diào)節(jié)系統(tǒng)，可實時監(jiān)測球磨機的運行狀態(tài)，并動態(tài)調(diào)整技術(shù)參數(shù)，提高球磨機的磨礦效果及襯板和介質(zhì)的使用壽命。

4結(jié)論

自離散元法首次在球磨機仿真中應用至今已有近30年的時間，在此期間研究人員對球磨機的離散元仿真進行了各種研究，加快了球磨機的發(fā)展進程及新產(chǎn)品研制的速度，使球磨機的理論研究及應用均取得了長足進展，但目前離散元法在球磨機仿真中仍然存在著模型精度低及仿真時間長等問題。因此，本文重點論述了各離散元顆粒模型、接觸模型和聯(lián)合仿真方法的特點，并從功率及能量分布研究、技術(shù)參數(shù)優(yōu)化、襯板結(jié)構(gòu)參數(shù)優(yōu)化3個方面介紹了離散元法在球磨機仿真中的應用情況，提出球磨機仿真模型的精細化發(fā)展、球磨機聯(lián)合仿真方法的廣泛應用、球磨機技術(shù)參數(shù)的智能調(diào)節(jié)是離散元法在球磨機仿真應用中未來的發(fā)展趨勢。隨著工業(yè)大數(shù)據(jù)、人工智能、數(shù)字孿生等智能技術(shù)交叉融合，離散元法在球磨機仿真領(lǐng)域的智能創(chuàng)新應用，勢必成為智慧礦山體系建設中不可或缺的重要技術(shù)手段。

[參 考 文 獻]

[1]JAYASUNDARA C T，ZHU H P.Impact energy of particles in ball mills based on DEM simulations and data-driven approach[EB/OL].（2021-09-24）[2021-10-05].https：∥www.sciencedirect.com/science/article/abs/pii/S0032591021008561.

[2]WANG T，ZOU W J，XU R J，et al.Assessing load in ball mill using instrumented grinding media[J].Minerals Engineering，2021，173：107198.

[3]YU J W，QIN Y H，GAO P，et al.An innovative approach for determining the grinding media system of ball mill based on grinding kinetics and linear superposition principle[J].Powder Technology，2021，378（Part A）：172-181.

[4]GAO P，ZHOU W T，HAN Y X，et al.Enhancing the capacity of large-scale ball mill through process and equipment optimization：an industrial test verification[J].Advanced Powder Technology，2020，31（5）：2 079-2 091.

[5]MARIJNISSEN M J，GRACZYKOWSKI C，ROJEK J.Simulation of the comminution process in a high-speed rotor mill based on the feeds macroscopic material data[J].Minerals Engineering，2021，163：106746.

[6]耿凡，袁竹林，孟德才，等.球磨機中顆粒混合運動的數(shù)值模擬[J].熱能動力工程，2009，24（5）：623-629，682.

[7]BIAN X L，WANG G Q，WANG H D，et al.Effect of lifters and mill speed on particle behaviour，torque，and power consumption of a tumbling ball mill：experimental study and DEM simulation[J].Minerals Engineering，2017，105：22-35.

[8]李騰飛，林蜀勇，張博，等.不同轉(zhuǎn)速率下球磨機內(nèi)鋼球的碰撞研究[J].中南大學學報（自然科學版），2019，50（2）：251-256.

[9]IWASAKI T，YAMANOUCHI H.Ball-impact energy analysis of wet tumbling mill using a modified discrete element method considering the velocity dependence of friction coefficient[J].Chemical Engineering Research and Design，2020，163：241-247.

[10]JIANG S Q，YE Y X，TAN Y Q，et al.Discrete element simulation of particle motion in ball mills based on similarity[J].Powder Technology，2018，335：91-102.

[11]CUNDALL P A.A computer model for simulating progressive large-scale movements in blocky rock systems[C]∥ISRM.Proceedings of the symposium of the international society for rock mechanics.Lisbon：ISRM，1971：128-132.

[12]MISHRA B K，RAJAMANI R K.The discrete element method for the simulation of ball mills[J].Applied Mathematical Modelling，1992，16（11）：598-604.

[13]AGRAWALA S，RAJAMANI R K，SONGFACK P，et al.Mechanics of media motion in tumbling mills with 3d discrete element method[J].Minerals Engineering，1997，10（2）：215-227.

[14]HE S Y，GAN J Q，PINSON D，et al.Particle shape-induced axial segregation of binary mixtures of spheres and ellipsoids in a rotating drum[J].Chemical Engineering Science，2021，235：116491.

[15]YANG S L，WANG H，WEI Y G，et al.Flow dynamics of binary mixtures of non-spherical particles in the rolling-regime rotating drum[J].Powder Technology，2019，361：930-942.

[16]FENG Y T，HAN K，OWEN D R J.A generic contact detection framework for cylindrical particles in discrete element modelling[J].Computer Methods in Applied Mechanics and Engineering，2017，315：632-651.

[17]XIE C H，SONG T，ZHAO Y Z.Discrete element modeling and simulation of non-spherical particles using polyhedrons and super-ellipsoids[J].Powder Technology，2020，368：253-267.

[18]WEERASEKARA N S，LIU L X，POWELL M S.Estimating energy in grinding using DEM modelling[J].Minerals Engineering，2016，85：23-33.

[19]ZHANG B，REGUEIRO R，DRUCKREY A，et al.Construction of poly-ellipsoidal grain shapes from SMT imaging on sand，and the development of a new DEM contact detection algorithm[J].Engineering Computations，2018，35（2）：733-771.

[20]METZGER M J，GLASSER B J.Simulation of the breakage of bonded agglomerates in a ball mill[J].Powder Technology，2013，237：286-302.

[21]BIBAK Z，BANISI S.A combined physical and DEM modelling approach to investigate particle shape effects on load movement in tumbling mills[J].Advanced Powder Technology，2021，32（2）：916-930.

[22]李同清.球磨機礦磨介質(zhì)動力學行為研究[D].徐州：中國礦業(yè)大學，2018.

[23]王肖江.基于離散元法的武山銅礦半自磨介質(zhì)優(yōu)化研究[D].昆明：昆明理工大學，2017.

[24]李鴻程，董為民，趙德化，等.礦物顆粒形狀對球磨機磨礦性能的影響[J].中國粉體技術(shù)，2013，19（2）：70-73.

[25]GOVENDER N，RAJAMANI R，WILKE D N，et al.Effect of particle shape in grinding mills using a GPU based DEM code[J].Mine-rals Engineering，2018，129：71-84.

[26]CLEARY P W，MORRISON R D.Comminution mechanisms，particle shape evolution and collision energy partitioning in tumbling mills[J].Minerals Engineering，2016，86：75-95.

[27]BARRIOS G K P，DE CARVALHO R M，KWADE A，et al.Contact parameter estimation for DEM simulation of iron ore pellet handling[J].Powder Technology，2013，248：84-93.

[28]IWASAKI T，YABUUCHI T，NAKAGAWA H，et al.Scale-up methodology for tumbling ball mill based on impact energy of grinding balls using discrete element analysis[J].Advanced Powder Technology，2010，21（6）：623-629.

[29]POWELL M S，WEERASEKARA N S，COLE S，et al.DEM modelling of liner evolution and its influence on grinding rate in ball mills[J].Minerals Engineering，2011，24（3/4）：341-351.

[30]ALDER B J，WAINWRIGHT T E.Phase transition for a hard sphere system[J].Journal of Chemical Physics，1957，27（5）：1 208-1 209.

[31]CUNDALL P A，STRACK O D L.A discrete numerical mode for granular assemblies[J].Géotechnique，1979，29（1）：47-65.

[32]JOHNSON K L，KENDALL K，ROBERTS A D.Surface energy and the contact of elastic solids[J].Proceedings of the Royal Society A：mathematical，Physical and Engineering Sciences，1971，324：301-313.

[33]WALTON O R，BRAUN R L.Viscosity，granular-temperature，and stress calculations for shearing assemblies of inelastic，frictional disks[J].Journal of Rheology，1998，30（5）：949-980.

[34]MINDLIN R D，DERESIEWICZ H.Elastic spheres in contact under varying oblique forces[J].Journal of Applied Mechanics—Transactions of the Asme，1953，20（3）：327-344.

[35]IWASHITA K，ODA M.Rolling resistance at contacts in simulation of shear band development by DEM[J].Journal of Engineering Mechanics，1998，124（3）：285-292.

[36]ODA M，KONISHI J，NEMAT-NASSER S.Experimental micromechanical evaluation of the strength of granular materials：effects of particle rolling[J].Mechanics of Materials，1982，1（4）：269-283.

[37]FENG Y T.An energy-conserving contact theory for discrete element modelling of arbitrarily shaped particles：basic framework and general contact model[J].Computer Methods in Applied Mecha-nics and Engineering，2021，373：113454.

[38]FENG Y T.An energy-conserving contact theory for discrete element modelling of arbitrarily shaped particles：contact bolume based model and computational issues[J].Computer Methods in Applied Mechanics and Engineering，2021，373：113493.

[39]李昆塬，童昕，李占福，等.基于DEM-FEM耦合的球磨機梯形襯板參數(shù)研究[J].福建工程學院學報，2017，15（6）：523-527.

[40]JONSN P，PLSSON B I，TANO K，et al.Prediction of mill structure behaviour in a tumbling mill[J].Minerals Engineering，2010，24（3/4）：236-244.

[41]BEINERT S，F(xiàn)RAGNIRE G，SCHILDE C，et al.Analysis and modelling of bead contacts in wet-operating stirred media and planetary ball mills with CFD-DEM simulations[J].Chemical Engineering Science，2015，134：648-662.

[42]MAYANK K，MALAHE M，GOVENDER I，et al.Coupled DEM-CFD model to predict the tumbling mill dynamics[J].Procedia IUTAM，2015，15：139-149.

[43]SINNOTT M D，CLEARY P W，MORRISON R D.Combined DEM and SPH simulation of overflow ball mill discharge and trommel Flow[J].Minerals Engineering，2017，108：93-108.

[44]JONSN P，STENER J F，PLSSON B I，et al.Validation of a model for physical interactions between pulp，charge and mill structure in tumbling mills[J].Minerals Engineering，2015，73：77-84.

[45]KANO J，MIO H，SAITO F.Correlation of grinding rate of gibbsite with impact energy of balls[J].AIChE Journal，2000，46（8）：1 694-1 697.

[46]CLEARY P W，MORRISON R D.Understanding fine ore breakage in a laboratory scale ball mill using DEM[J].Minerals Engineering，2011，24（3/4）：352-366.

[47]穆海芳，何康，韓君，等.基于離散元法的球磨機介質(zhì)運動的仿真研究[J].濟寧學院學報，2019，40（2）：12-16.

[48]HUANG P，DING Y F，WU L，et al.A novel approach of evaluating crushing energy in ball mills using regional total energy[J].Powder Technology，2019，355（6）：289-299.

[49]BOND F C.The third theory of comminution[J].Transactions of AIME Mining Engineering，1952，193：484-494.

[50]SHI F N，XIE W G.A specific energy-based size reduction model for batch grinding ball mill[J].Minerals Engineering，2015，70：130-140.

[51]OROZCO L F，NGUYEN D H，DELENNE J Y，et al.Discrete-element simulations of comminution in rotating drums：effects of grinding media[J].Powder Technology，2020，362：157-167.

[52]CLEARY P W.Charge behaviour and power consumption in ball mills：sensitivity to mill operating conditions，liner geometry and charge composition[J].International Journal Mineral Processing，2001，63（2）：79-114.

[53]孫鵬.基于離散元素法的巖土球磨機仿真與分析[J].軟件，2019，40（8）：145-148，189.

[54]劉剛，張冠宇，李婧錫，等.基于離散元的球磨機工作參數(shù)優(yōu)化[J].礦山機械，2014，42（7）：79-82.

[55]PEDRAVES F，NORNIELLA J G，MELERO M G，et al.Frequency domain characterization of torque in tumbling ball mills using DEM modelling：application to filling level monitoring[J].Powder Technology，2018，323：433-444.

[56]DE CARVALHO R M，TAVARES L M.Predicting the effect of operating and design variables on breakage rates using the mechanistic ball mill model[J].Minerals Engineering，2013，43/44：91-101.

[57]HERBST J A，F(xiàn)UERSTENAU D W.Scale-up procedure for continuous grinding mill design using population balance models[J].International Journal of Mineral Processing，1980，7（1）：1-31.

[58]ROSALES-MAR?N G，ANDRADE J，ALVARADO G，et al.Study of lifter wear and breakage rates for different lifter geometries in tumbling mill：experimental and simulation analysis using population balance model[J].Minerals Engineering，2019，141：105857.

[59]CHIMWANI N，BWALYA M M.Using DEM to investigate how shell liner can induce ball segregation in a ball mill[J].Minerals Engineering，2020，151：106311.

[60]WU M，WANG V.Modeling ball impact on the wet mill liners and its application in predicting mill magnetic liner performance[J].Minerals Engineering，2014，61：126-132.

[61]張學東，董為民，周海燕，等.具有仿生特征的球磨機提升條耐磨性數(shù)值模擬[J].有色金屬（選礦部分），2017（6）：56-62.

[62]BOEMER D，PONTHOT J P.A generic wear prediction procedure based on the discrete element method for ball mill liners in the cement industry[J].Minerals Engineering，2017，109：55-79.

[63]REZAEIZADEH M，F(xiàn)OOLADI M，POWELL M S，et al.A new predictive model of lifter bar wear in mills[J].Minerals Engineering，2010，23（15）：1 174-1 181.

作者簡介：裴英杰（1997—），男，山東平度人，碩士研究生，研究方向為碎磨理論與工藝;昆明市五華區(qū)一二一大街文昌路68號，昆明理工大學國土資源工程學院，650093;E-mail：z179180802@163.com

通信作者，E-mail：13515877@qq.com，13888532894

裴英杰，肖慶飛，石貴明，周強（1.昆明理工大學國土資源工程學院; 2.礦冶過程自動控制技術(shù)國家（北京市）重點實驗室;3.玉溪師范學院化學生物與環(huán)境學院

Research progress and application of discrete element modeling in ball mill simulationPei Yingjie，Xiao Qingfei，Shi Guiming，Zhou Qiang

（1.School of Land Resource Engineering，Kunming University of Science and Technology;

2.State（Beijing） Key Laboratory of Process Automation in Mining & Metallurgy Process;

3.School of Chemical Biology and Environment，Yuxi Normal University）

Abstract：As important equipment in the ore grinding and crushing stage，the ball mill plays a decisive role in improving efficiency and reducing consumption of the concentrator.Ball mill power consumption and energy consumption distribution，the design optimization of technical parameters and liner structural parameter are important content and technical difficulties that affect the grinding effect of ball mills.Discrete element modeling can accurately simulate the deformation process of particle motion and record the power and energy distribution of particles，so that in the field of ball mill simulation，it has good application prospects.The paper introduces the characteristics of the widely used particle model，contact model，and joint simulation method of the discrete element modeling in the ball mill simulation process and their research progress，as well as? the application of ball mill power and energy distribution research，technical parameter optimization，liner structural parameter optimization，summarizes and points out that the refinement of the ball mill simulation model，the wide application of joint simulation，and the intelligent adjustment of ball mill technical parameters are the future ball mill simulation research development trend.

Keywords：discrete element modeling;ball mill;simulation;particle model;contact model;coupling method