磨礦破碎過程粒度分布的分布式參數(shù)蒙特卡洛動(dòng)力學(xué)模擬及加速方法

盧紹文 蔚潤(rùn)琴 崔玉潔

選礦過程中,顆粒破碎是在外力作用下使物料由大變小的過程,使得原料能夠滿足后續(xù)的選別作業(yè)對(duì)物料的粒度要求,便于物料中有用礦物與脈石解離.顆粒破碎研磨是高耗能工業(yè)過程,運(yùn)行優(yōu)化的主要目標(biāo)是在滿足工藝要求的情況下降低能耗,并使產(chǎn)品粒度穩(wěn)定[1?3].粒度分布是破裂過程的關(guān)鍵工藝指標(biāo),它是不同粒度顆粒在總物料中占比.粒度分布過大不能滿足工藝要求,但是如果物料磨得過細(xì)會(huì)導(dǎo)致能耗過高.粒度分布可以通過激光粒度儀在線檢測(cè),但是這類儀表成本高,而且受到生產(chǎn)環(huán)境干擾影響,需經(jīng)常校準(zhǔn).目前通常采用離線篩析,其主要問題是滯后大,不能滿足在線運(yùn)行優(yōu)化的要求.由于破碎過程的主要擾動(dòng)來自于原料的粒度、硬度的變化,利用破碎過程模擬技術(shù)來預(yù)測(cè)在當(dāng)前的生產(chǎn)條件下最終的粒度分布,可通過運(yùn)行操作優(yōu)化來補(bǔ)償干擾[4?6].此外,經(jīng)過驗(yàn)證的模擬模型還可以用于檢驗(yàn)工藝設(shè)計(jì),對(duì)于優(yōu)化磨礦生產(chǎn)過程控制和設(shè)計(jì)都具有重要應(yīng)用價(jià)值[7?10].

破碎過程粒度分布的模擬主要有三類方法.第一類是穩(wěn)態(tài)模擬方法[11],該方法假設(shè)生產(chǎn)的邊界條件不變,將粒度分布的變化近似為線性過程,主要用于初期的工藝設(shè)計(jì),不能用于在線指導(dǎo)操作.由于磨礦過程機(jī)理極為復(fù)雜,難以建立精確的數(shù)學(xué)模型,一般采用唯像學(xué)建模法[11].20 世紀(jì)40 年代,利用統(tǒng)計(jì)原理研究磨機(jī)內(nèi)物料粒度分布,提出了選擇函數(shù)和破裂函數(shù)兩個(gè)概念[12].其中選擇函數(shù)描述顆粒發(fā)生破裂的概率,破裂函數(shù)刻畫顆粒破碎后得到的子顆粒的粒度分布.Austin 等經(jīng)過實(shí)驗(yàn)給出了選擇函數(shù)和破裂函數(shù)的經(jīng)驗(yàn)公式.文獻(xiàn)[13]提出矩陣動(dòng)力學(xué)模型與總體平衡模型,為磨礦的研究開辟了一個(gè)新的領(lǐng)域.文獻(xiàn)[14]提出磨礦的動(dòng)態(tài)矩陣?yán)硐肽Ｐ?進(jìn)一步推動(dòng)了磨礦數(shù)學(xué)模型的發(fā)展.基于能級(jí)關(guān)系的Bond 模型是比較經(jīng)典的球磨機(jī)模型,但這種模型得不到各粒級(jí)的質(zhì)量分?jǐn)?shù).上述這些模型主要是穩(wěn)態(tài)模型,無法刻畫顆粒破碎過程中粒度分布的演變.

第二類方法為動(dòng)態(tài)模擬方法.磨礦過程的動(dòng)態(tài)模型需要建立粒度分布隨時(shí)間變化的過程模型,輸入物料的性質(zhì)、操作條件,模型輸出每一時(shí)刻粒度分布的模擬值作為輸出.該方法主要采用總量平衡原理建立粒度分布的主方程模型[15].通過磨礦粒度分布的群體平衡動(dòng)態(tài)模型,記錄磨機(jī)中不同位置各個(gè)粒徑的顆粒每時(shí)每刻的狀態(tài),從而描述出磨礦過程中物料在時(shí)間和空間上的粒度分布.但該模型一般是積分微分方程形式,多數(shù)情況下難以獲得解析解[16].

第三類是動(dòng)力學(xué)模擬方法,它是一種隨機(jī)模擬方法[17?18],它用隨機(jī)過程來刻畫顆粒的破碎過程,能夠得到顆粒破碎過程的演化過程.蒙特卡洛動(dòng)力學(xué)模擬方法(以下簡(jiǎn)稱MC 方法)是常用的隨機(jī)模擬方法,它通過模擬每一次破碎事件的發(fā)生,記錄研磨過程中粒度分布的動(dòng)態(tài)變化情況,較真實(shí)模擬了物理過程.文獻(xiàn)[17?18]提出了基于MC 方法的顆粒破碎過程模擬算法.但是在顆粒破碎過程中,大粒徑顆粒會(huì)破碎為多個(gè)小粒級(jí)的顆粒,因此隨著仿真時(shí)間的推進(jìn),需要模擬的顆?？倲?shù)以指數(shù)速度增加,導(dǎo)致計(jì)算量快速增大,效率降低.針對(duì)這一問題,Khalili 等和Smith 等提出了一種定數(shù)量MC(Constant-number Monte Carlo,CNMC)方法[19?20],該算法在處理顆粒的凝并和破碎事件時(shí),通過隨機(jī)從系統(tǒng)中增加或移除顆粒,目的是使系統(tǒng)總的顆粒樣本數(shù)保持不變,從而提高了原有MC方法的效率.趙海波等提出了一種多重MC(Multi Monte Carlo,MMC)方法,通過引入加權(quán)的“虛擬顆?！钡母拍?將實(shí)際顆粒由一顆或者幾顆屬性相近的虛擬顆粒代表并賦予一個(gè)數(shù)目轉(zhuǎn)換權(quán)值,該值就是該虛擬顆粒所代表的當(dāng)?shù)貙?shí)際顆粒的個(gè)數(shù),模擬時(shí)通過調(diào)整數(shù)目轉(zhuǎn)換權(quán)值來保證在發(fā)生凝并和破裂事件之后虛擬顆粒總數(shù)目仍然保持不變.盧紹文在CNMC 基礎(chǔ)上[21]提出一種帶精度反饋和速度控制的動(dòng)力學(xué)模擬算法(ED-KMC-AC),能夠在滿足一個(gè)給定的精度閾值約束下,最大地提高模擬速度.這些方法能夠加快MC 模擬的速度,其本質(zhì)上是通過人為減少系統(tǒng)內(nèi)顆粒數(shù)目實(shí)現(xiàn)的,這將降低模擬的準(zhǔn)確程度.Gillepsie 在文獻(xiàn)[22]中提出了τ-leap近似加速算法,他假設(shè)磨礦過程中多個(gè)顆粒破碎信息可以合并,不需要每一次事件都更新系統(tǒng),通過一個(gè)調(diào)節(jié)因子來控制系統(tǒng)更新的粒度,從而從總體上調(diào)節(jié)模擬的速度和精度.

上述模型都屬于集總參數(shù)模型,即假設(shè)磨機(jī)是一個(gè)質(zhì)量分布均勻的質(zhì)點(diǎn),物料在磨機(jī)內(nèi)處處均勻.但實(shí)際工業(yè)生產(chǎn)中,工業(yè)磨機(jī)的尺寸日趨大型化,物料在磨機(jī)內(nèi)的駐留時(shí)間更長(zhǎng),物料性質(zhì)的空間分布難以忽略,集總參數(shù)模擬方法已經(jīng)不能準(zhǔn)確刻畫物料在磨機(jī)腔體內(nèi)的駐留和變化[23].針對(duì)這個(gè)問題,Kis 等[23]提出了分布參數(shù)的磨礦過程模型.該模型不僅研究磨機(jī)徑向方向的破裂,同時(shí)也考慮了磨機(jī)軸向方向物料的移動(dòng),考慮了對(duì)流前移和擴(kuò)散后移兩種顆粒消失事件.該方法將磨機(jī)分為若干段虛擬的子磨機(jī),從而得到每一段子磨機(jī)在連續(xù)磨礦過程中物料的變化情況,最后合并得到產(chǎn)品的粒度分布.但是基于該方法的模型需要采用數(shù)值近似方法進(jìn)行求解,由于模型方程非常復(fù)雜,不能保證解的存在.此外,該方法只能解出粒度分布隨時(shí)間變化的近似解,無法描述分布參數(shù)下磨礦模擬過程中顆粒的破裂和移動(dòng)的隨機(jī)特性.針對(duì)這些問題,本文提出一種基于分布式參數(shù)破碎過程模擬的粒度分布終點(diǎn)預(yù)報(bào)方法.我們采用分布參數(shù)模擬方法,將磨機(jī)在軸向分段為若干串級(jí)連接的虛擬的子磨機(jī),由于每個(gè)子磨機(jī)的軸長(zhǎng)足夠小,可以假設(shè)認(rèn)為每段子磨機(jī)內(nèi)物料為充分混合狀態(tài),從而將連續(xù)磨礦過程看作多個(gè)串級(jí)連接的批次子磨過程,但這種方法面臨計(jì)算代價(jià)過大的問題.為了提高模擬效率,在批次磨礦過程的蒙特卡洛動(dòng)力學(xué)模擬方法基礎(chǔ)上,給出一種模擬加速方法.

1 相關(guān)算法及存在的問題

磨機(jī)內(nèi)物料發(fā)生破碎的過程很難精確建模,只能采用概率的方法,即將顆粒的破碎看成一個(gè)隨機(jī)事件,通過人為控制實(shí)驗(yàn)條件下進(jìn)行大量實(shí)驗(yàn)獲得的實(shí)驗(yàn)樣本,采用統(tǒng)計(jì)方法得到不同條件下破碎事件發(fā)生的概率.將顆粒的破碎過程看作是一個(gè)隨機(jī)過程,通過實(shí)驗(yàn)我們可以獲得兩方面重要特性,一個(gè)是破碎的速率,一個(gè)是顆粒破碎后獲得的子顆粒的分布.前者假設(shè)顆粒破碎的發(fā)生符合泊松分布,破碎事件發(fā)生的間隔時(shí)間服從指數(shù)分布,通過選擇函數(shù)的形式給出.后者通過物理實(shí)驗(yàn)獲得,它主要與物料的硬度和材質(zhì)相關(guān),通過破裂函數(shù)的形式給出.但是顆粒破碎速率與當(dāng)前的粒度分布、濃度等多種因素有關(guān).這些因素在生產(chǎn)過程中是不斷變化的,如果再考慮模型的分布參數(shù),模型的形式異常復(fù)雜,難以獲得解析解.

1.1 集總參數(shù)的顆粒破碎蒙特卡洛模擬算法

蒙特卡洛(MC)模擬方法是一種隨機(jī)數(shù)值算法,主要利用隨機(jī)抽樣原理實(shí)現(xiàn)對(duì)任意積分過程的數(shù)值模擬.目前采用蒙特卡洛方法模擬顆粒破碎過程,主要是利用蒙特卡洛抽樣方法來模擬時(shí)變的顆粒破碎概率分布以及破裂函數(shù).蒙特卡洛方法分事件驅(qū)動(dòng)和時(shí)間驅(qū)動(dòng)兩種,它們的數(shù)值模擬能力是等價(jià)的[24].以下給出時(shí)間驅(qū)動(dòng)蒙特卡洛模擬算法.

算法1.時(shí)間驅(qū)動(dòng)蒙特卡洛顆粒破碎算法

1)REPEAT

2)反函數(shù)法求下一破裂事件步長(zhǎng)

3)舍選法選取被破裂顆粒

4)反函數(shù)確定破裂子顆粒的粒級(jí)

5)系統(tǒng)狀態(tài)更新

6)UNTIL 模擬時(shí)間達(dá)到最大

假設(shè)磨機(jī)內(nèi)破碎事件的發(fā)生滿足馬爾科夫特性,粒徑為d的顆粒發(fā)生破碎事件的速率由選擇函數(shù)S(d)定義,事件間隔時(shí)間服從指數(shù)分布.從當(dāng)前時(shí)刻t到下一破碎事件發(fā)生的時(shí)間間隔?t可通過反函數(shù)法得到:

其中,R1為(0,1)區(qū)間隨機(jī)數(shù),Sb為所有計(jì)算區(qū)域內(nèi)顆粒發(fā)生破碎事件的總速率.采用蒙特卡洛舍選法[17]確定被破裂顆粒所在粒級(jí),生成(0,1)區(qū)間隨機(jī)數(shù)R2,在1 到N中隨機(jī)選擇一個(gè)粒級(jí)j,若下式成立,則j為發(fā)生被破裂顆粒所在粒級(jí):

其中,Si為粒級(jí)i的選擇函數(shù).由于破裂函數(shù)一般不隨生產(chǎn)條件變化,仍采用反函數(shù)法確定破裂產(chǎn)生的子顆粒的粒徑,生成(0,1)區(qū)間隨機(jī)數(shù)R3,在粒級(jí)j到N隨機(jī)選擇粒級(jí)k,如果下式成立,則k為破裂產(chǎn)生的子顆粒的粒級(jí).

其中,B(k,j)為破裂函數(shù)的累積量.最后,更新當(dāng)前各個(gè)粒級(jí)顆粒數(shù)量的系統(tǒng)狀態(tài)向量X(t)以及各個(gè)粒級(jí)顆粒的破碎速率S和總破碎速率Sb.

1.2 分布式參數(shù)的顆粒破碎過程離散差分?jǐn)?shù)值算法

算法1 介紹的蒙特卡洛算法屬于集總參數(shù)算法,主要用于批次磨的模擬.對(duì)于連續(xù)磨過程,物料從磨機(jī)入口到出口的傳遞過程中物料顆粒在軸向上不一定分布均勻,所以集總參數(shù)模型不能精確模擬連續(xù)磨過程.針對(duì)這一問題,Kis 等[23]提出一種考慮物料在磨機(jī)中軸向分布的基于離散差分方法的連續(xù)磨粒度分布動(dòng)態(tài)模型.如圖1 所示,首先將磨機(jī)在軸向上等分為J段子磨機(jī),分別為y1,y2,···,yj,···,yJ;徑向上分為N個(gè)粒級(jí),分別為x1,x2,···,xi,···,xN.假設(shè)每段子磨機(jī)內(nèi)物料顆粒分布均勻,對(duì)子磨機(jī)利用批次磨集總參數(shù)模擬方法對(duì)每個(gè)子磨機(jī)同時(shí)進(jìn)行模擬.而與批次磨略有不同的地方是每段子磨機(jī)各個(gè)粒級(jí)的顆粒數(shù)目的變化不僅來自于徑向上顆粒之間的破裂過程,同時(shí)也是相鄰子磨機(jī)之間顆粒相互移動(dòng)的結(jié)果.

圖1 離散徑向分布參數(shù)磨機(jī)模型示意圖Fig.1 Illustration of the distributed parameter grinding mill model

以任意中間段第j段子磨機(jī)的i粒級(jí)為例分析磨機(jī)中間段的物料變化,稱矩形區(qū)域[yj?1,yj][xj?1,xj]為(j,i)單元,u(yj,xi,tn)代表該單元的質(zhì)量,另外定義m(y,x,tn)為物料在tn時(shí)刻的密度,則

該單元物料的質(zhì)量包括四小部分:1)該單元經(jīng)過破裂及先后移動(dòng)后剩余的質(zhì)量;2)上段子磨機(jī)移動(dòng)對(duì)流;3)下段子磨機(jī)移動(dòng)逆流;4)徑向上其他粒徑比較大的粒級(jí)破裂為粒級(jí)顆粒.則該單元在下一時(shí)刻粒度分布可以由以下離散方程表示:

其中,dt表示模擬步進(jìn)時(shí)間,VF表示物料單位時(shí)間向前移動(dòng)比例,VB表示物料單位時(shí)間向后移動(dòng)比例,Pk,i表示單位時(shí)間i粒級(jí)的顆粒破裂到k粒級(jí)的比例.

以(1,i)單元為例分析入料口物料的質(zhì)量變化,與中間段相比該單元物料不存在來自上一段子磨機(jī)的移動(dòng)對(duì)流,但有給料注入物料.因此:

最后,以(J,i)單元為例分析出口處物料的質(zhì)量變化.與中間段相比,出口處不存在下段子磨機(jī)逆流移動(dòng)而來的物料,因此:

如上,通過求解方程得到每一段子磨機(jī)中物料顆粒的動(dòng)態(tài)變化過程,從而模擬整個(gè)磨機(jī)物料的粒度分布的變化情況.

1.3 存在的問題

顆粒破碎蒙特卡洛模擬算法的主要優(yōu)點(diǎn)是它能夠準(zhǔn)確地模擬物料隨時(shí)間變化的過程,但也存在兩方面缺點(diǎn):1)基于集總參數(shù)的蒙特卡洛動(dòng)力學(xué)算法無法反映物料性質(zhì)在磨機(jī)內(nèi)部的差異及動(dòng)態(tài)變化,因此不適于大型磨機(jī)以及連續(xù)作業(yè)過程的模擬;2)這類算法通常需要較大的計(jì)算資源,而且模擬的速度隨著計(jì)算的進(jìn)行不斷衰減,原因見文獻(xiàn)[18].Kis等[23]提出磨機(jī)軸向的離散化近似來解決這一問題.該算法采用時(shí)間差分得到差分方程,然后不斷計(jì)算步進(jìn)時(shí)間內(nèi)物料質(zhì)量的變化,從而預(yù)測(cè)模擬到達(dá)研磨設(shè)定時(shí)間的模擬結(jié)果.Kis 算法的主要問題在于選擇的差分步進(jìn)時(shí)間必須滿足一定的約束條件,而該約束條件是隨著模擬不斷變化的,確定的步進(jìn)時(shí)間可能導(dǎo)致約束條件在模擬過程中被破壞.

步進(jìn)時(shí)間的大小直接影響著物料質(zhì)量的變化程度,因此步進(jìn)時(shí)間的大小對(duì)于模擬結(jié)果具有較大的影響.以式(5)為例:等式右邊第一項(xiàng)表示該粒級(jí)內(nèi)剩余物料質(zhì)量,即該粒級(jí)內(nèi)的顆粒發(fā)生移動(dòng)或者破裂事件后剩余的物料質(zhì)量.由于剩余物料質(zhì)量不能為負(fù)值,因此應(yīng)有:

進(jìn)而可得到時(shí)間間隔的約束條件:

該算法時(shí)間步長(zhǎng)的選擇對(duì)于數(shù)值結(jié)果的精度影響至關(guān)重要,當(dāng)分別取不同的步進(jìn)時(shí)間會(huì)得到不同的累積粒度分布曲線,并且差別顯著.在差分方法中,步進(jìn)時(shí)間的選擇需要事先人為決定,導(dǎo)致其并不一定能夠滿足約束條件(10),這時(shí)模擬計(jì)算的結(jié)果是沒有物理意義的.另一方面,離散差分方法雖然能夠得到連續(xù)磨過程內(nèi)物料的質(zhì)量變化,并求得產(chǎn)品的粒度分布,但該方法結(jié)果屬于固定的數(shù)值解,模擬結(jié)果與試驗(yàn)次數(shù)無關(guān),因此它只是反映了連續(xù)磨過程的終點(diǎn)期望結(jié)果.而在連續(xù)磨過程中,顆粒的破裂與移動(dòng)具有隨機(jī)性,固定意義上的數(shù)值解無法描述連續(xù)磨礦的隨機(jī)過程,以及各個(gè)粒級(jí)內(nèi)顆粒的微觀變化過程.

2 分布式參數(shù)磨礦模擬算法

2.1 分布式參數(shù)蒙特卡洛動(dòng)力學(xué)模擬算法

離散差分的方法雖然通過求解方程能夠得到整個(gè)磨機(jī)物料的粒度分布情況,但該種方法不能準(zhǔn)確描述顆粒的動(dòng)態(tài)變化過程.而且磨機(jī)體積趨向于越來越大,里面的物料顆粒數(shù)目比較龐大,不同時(shí)間點(diǎn)各個(gè)位置各個(gè)粒級(jí)顆粒的破碎都是隨機(jī)的,所以本文在離散差分方法的基礎(chǔ)上提出一種基于隨機(jī)算法的分布式參數(shù)磨礦過程模擬方法.

本算法的思路借鑒離散差分算法,將磨機(jī)分為多個(gè)串級(jí)連接的子磨機(jī),如圖1 所示.軸向上對(duì)磨機(jī)按編號(hào)從1 到j(luò)段子磨機(jī).對(duì)于每段子磨機(jī),徑向上劃分為N個(gè)粒級(jí),每個(gè)粒級(jí)用代表粒徑di來表示.并假設(shè)子磨機(jī)內(nèi)物料均勻混合,即子磨機(jī)內(nèi)的粒度分布處處相同;顆粒在子磨機(jī)之間移動(dòng)是在破裂之后發(fā)生的.對(duì)于子磨機(jī)內(nèi)的顆粒破裂過程,我們假設(shè)軸向的粒度差異可以忽略,因此仍采用算法1 給出的時(shí)間驅(qū)動(dòng)蒙特卡洛模擬方法.根據(jù)分段思想,相鄰的子磨機(jī)間有物料的相互傳遞.因此,子磨機(jī)內(nèi)顆粒數(shù)目的增加是由大顆粒破裂或者子磨機(jī)間顆粒的流入而導(dǎo)致的.從微觀角度來觀察,磨機(jī)內(nèi)顆?？赡馨l(fā)生三種事件:向前移動(dòng)事件,是由對(duì)流和物料本身的循環(huán)流動(dòng)引起;向后移動(dòng)事件,是由擴(kuò)散作用引起;破裂事件,是由磨機(jī)研磨作用而導(dǎo)致顆粒破裂至比其粒級(jí)小的粒級(jí)內(nèi).將上述三種事件統(tǒng)稱為顆粒消失事件.

首先,針對(duì)上述離散模型,定義新的系統(tǒng)狀態(tài)矩陣XN×J,其中x(i,j)表示第j個(gè)子磨機(jī)內(nèi)第i個(gè)粒級(jí)的顆粒數(shù)量.由于分布式參數(shù)模型中傾向函數(shù)的計(jì)算需要靠考慮顆粒的對(duì)流前移、擴(kuò)散后移及破裂三種顆粒消失事件,新定義的傾向函數(shù)反映這三種事件發(fā)生的總速率,定義如下:

對(duì)此又分為三種情況:給料段(j=1)存在對(duì)流前移和破裂2 類事件;中間段(j ∈{2,···,J ?1})存在前移、后移和破裂3 類事件;排料段(j=J)存在對(duì)流后移和破裂2 類事件.

其次,根據(jù)傾向函數(shù)計(jì)算步進(jìn)時(shí)間.由于磨機(jī)被等分為J子磨機(jī),不同子磨機(jī)分別采用MC 方法進(jìn)行模擬計(jì)算,由于步進(jìn)時(shí)間的選取只與子磨機(jī)局部狀態(tài)相關(guān),因此導(dǎo)致全局時(shí)間不一致.針對(duì)這一問題,我們每次計(jì)算得到J個(gè)不同的步進(jìn)時(shí)間,并取其中的最小值作為全局系統(tǒng)的推進(jìn)步長(zhǎng).

其中對(duì)于發(fā)生消失事件的子磨機(jī),采用算法1計(jì)算當(dāng)前子磨機(jī)的粒度分布.由于增加了前向流動(dòng)和后向流動(dòng)兩種新的事件,因此在用舍選法選擇消失顆粒時(shí),需要根據(jù)式(11)定義的傾向函數(shù)來計(jì)算.

綜上,算法的主要步驟如下:

算法2.顆粒破碎的分布式參數(shù)蒙特卡洛動(dòng)力學(xué)模擬算法

1)REPEAT

2)利用式(11)計(jì)算傾向函數(shù)

3)確定下一顆粒消失事件發(fā)生的子磨機(jī)并根據(jù)式(12)確定全局時(shí)間步長(zhǎng)

4)對(duì)發(fā)生消失事件的子磨機(jī),采用算法1 計(jì)算粒度分布

5)系統(tǒng)狀態(tài)更新

6)UNTIL 模擬時(shí)間達(dá)到最大

2.2 基于τ-leap 的磨礦過程分布式參數(shù)蒙特卡洛模擬加速算法

算法2 實(shí)現(xiàn)的連續(xù)破碎過程模擬能夠精確刻畫每一個(gè)子磨機(jī)的動(dòng)態(tài)變化過程,進(jìn)而得到整個(gè)磨機(jī)中不同時(shí)間不同位置各個(gè)時(shí)間點(diǎn)的粒度分布.但由于該種方法需要記錄每一次事件的發(fā)生,計(jì)算量極大,只能應(yīng)用于小體系的模擬.針對(duì)該問題,本文在算法2 基礎(chǔ)上提出一種基于τ-leap 的模擬加速方法.τ-leap 方法主要廣泛應(yīng)用于化學(xué)反應(yīng)過程的快速模擬.這里將τ-leap 算法引入顆粒破碎過程,通過將模擬時(shí)間劃分為若干時(shí)間片作為步進(jìn)時(shí)間,同時(shí)假設(shè)在每個(gè)步進(jìn)時(shí)間內(nèi)有多個(gè)反應(yīng)事件同時(shí)發(fā)生.即不必每次事件都更新系統(tǒng)狀態(tài),而是根據(jù)步進(jìn)時(shí)間片內(nèi)多個(gè)事件的累積效應(yīng)對(duì)系統(tǒng)狀態(tài)進(jìn)行更新,從而提高模擬速度.

在分布式參數(shù)模擬的框架下采用τ-leap 方法,就必須解決計(jì)算中系統(tǒng)狀態(tài)不一致的問題.在算法2中,子磨機(jī)系統(tǒng)內(nèi)的顆粒破碎和移動(dòng)均被視為該顆粒的消失事件.由于τ-leap 加速方法并不是針對(duì)每個(gè)事件更新系統(tǒng)狀態(tài),將其應(yīng)用于分布式參數(shù)模型時(shí),不同子磨機(jī)的狀態(tài)和時(shí)間推進(jìn)速度并不一致,這就帶來子系統(tǒng)的狀態(tài)同步和時(shí)間的同步問題.此外,MC 模擬過程中需要假設(shè)當(dāng)前子系統(tǒng)的狀態(tài)是獨(dú)立的,但是在時(shí)間段(t,t+τ)內(nèi)各類事件的累積效應(yīng)還需要考慮子磨機(jī)之間存在的物料流動(dòng),所以各個(gè)磨機(jī)之間的狀態(tài)又是相關(guān)的.這就導(dǎo)致τ-leap 算法的思想并不能直接應(yīng)用在分段離散模型中.

針對(duì)上述問題,本文提出一種基于緩沖區(qū)的解決方法.考慮到移動(dòng)事件和破裂事件直接的獨(dú)立性,我們對(duì)每個(gè)子磨機(jī)系統(tǒng)開辟兩個(gè)狀態(tài)同步緩沖區(qū):前向移動(dòng)緩沖區(qū)和后向移動(dòng)緩沖區(qū),如圖2 所示.在計(jì)算過程中,將前向移動(dòng)和后向移動(dòng)的顆粒分別暫存在緩沖區(qū),等到達(dá)下一個(gè)同步時(shí)間點(diǎn)t+τ后,再將緩沖區(qū)的顆粒與相鄰子磨機(jī)進(jìn)行交換.本質(zhì)上,這相當(dāng)于將前/后向的移動(dòng)事件和破裂事件分別進(jìn)行處理.相應(yīng)地,在進(jìn)行破裂事件的處理中,將子磨機(jī)和其附屬的緩沖區(qū)視為一個(gè)獨(dú)立的批次磨過程,并更新狀態(tài)轉(zhuǎn)移矩陣和傾向函數(shù).

圖2 分段子磨機(jī)緩沖區(qū)示意圖Fig.2 The buffer zones for sub-grinding mill

加速算法對(duì)每個(gè)子磨機(jī)的模擬計(jì)算采用基于τ-leap 的批次MC 模擬.子磨機(jī)j的步進(jìn)時(shí)間由以下公式計(jì)算:

其中,M為消失事件的種類,a0為當(dāng)前傾向函數(shù)之和,x為系統(tǒng)狀態(tài),為τ-leap 算法的速度調(diào)節(jié)參數(shù).對(duì)于時(shí)間同步的問題,由于τ-leap 算法的主要作用是加快模擬速度,因此步進(jìn)時(shí)間的同步采用與算法2相反的策略,即選取各個(gè)子磨機(jī)中的最大步長(zhǎng)作為全局系統(tǒng)的步進(jìn)同步時(shí)間點(diǎn).這樣做的目的是盡可能發(fā)揮加速的性能.

加速算法的主要步驟如下所示:

算法3.顆粒破碎的分布式參數(shù)蒙特卡洛動(dòng)力學(xué)τ-leap 算法

1)REPEAT

2)利用式(11)計(jì)算傾向函數(shù)

3)確定下一顆粒消失事件發(fā)生的子磨機(jī)并根據(jù)式(14)確定全局時(shí)間步長(zhǎng)

4)對(duì)發(fā)生消失事件的子磨機(jī),采用基于τ-leap MC 算法計(jì)算粒度分布,其中τ值由式(13)計(jì)算

5)系統(tǒng)狀態(tài)更新

6)UNTIL 模擬時(shí)間達(dá)到最大

3 算法性能的比較分析

3.1 速度和精度的性能指標(biāo)

動(dòng)力學(xué)模擬算法的性能主要從精度和速度兩個(gè)方面來衡量.本文研究主要采用文獻(xiàn)[21]所提出的兩個(gè)量化指標(biāo)用于以下的分析和討論.首先定義穩(wěn)態(tài)模擬速度(Transient simulation speed,TSS)來度量磨礦過程的模擬速度.TSS 越大說明模擬速度越快.表達(dá)式為:

其中,dti表示第i次與第i ?1 次破裂事件之間模擬時(shí)間間隔,ti表示實(shí)際的計(jì)算時(shí)間.

其次,定義平均粒徑加權(quán)變異系數(shù)(Weighted squared coefficient of variation,WSCV)來度量磨礦過程的模擬精度,WSCV 越小說明模擬精度越好.表達(dá)式為:

其中,v(i)代表第i粒級(jí)礦物顆粒的體積,v(j)代表第j粒級(jí)礦物顆粒的體積.

3.2 模擬計(jì)算實(shí)驗(yàn)

本文分別采用離散差分方法、MC 方法以及τleap 加速方法對(duì)磨礦過程進(jìn)行模擬實(shí)驗(yàn).模型的初始參數(shù)設(shè)置如下:磨機(jī)的長(zhǎng)度L=4.4 m,磨機(jī)分段數(shù)J=10,研磨時(shí)間T=2 min,擴(kuò)散協(xié)同系數(shù)D=0.005,對(duì)流速度u=0.065.為了便于比較,模型中的破裂函數(shù)、選擇函數(shù)參數(shù)以及給料粒度取值與文獻(xiàn)[21]相同,這里不再贅述.

3.2.1 算法精度比較

1)平均粒徑

圖3 所示為磨機(jī)出口物料各粒級(jí)代表粒徑的平均值隨時(shí)間的變化曲線(其中=0.001).從圖中可見磨機(jī)中物料的平均粒徑隨時(shí)間逐漸變小,這反映了磨機(jī)中較大粒徑的顆粒不斷破裂為較小粒徑顆粒的物理過程.從圖中可以看到MC 模擬方法和基于τ-leap 的模擬加速算法計(jì)算得到的平均粒度高度吻合,并不存在系統(tǒng)性的差別.統(tǒng)計(jì)兩種算法的相對(duì)誤差,最大誤差為4.3×10?2,最小誤差為2.13×10?6,平均誤差為8.19×10?4.值得指出的是,在時(shí)刻約0.05 min 之前,可以觀察到模型計(jì)算得到的排料為零.這是因?yàn)榉植际絽?shù)模型考慮了磨機(jī)軸向物料的分布,物流從給料口傳遞到排料口之前磨機(jī)沒有物料排出.這是與集總參數(shù)模型最顯著的不同.

圖3 分布參數(shù)MC 算法(圖中簡(jiǎn)稱MC)和τ-leap 加速算法計(jì)算得到的平均粒徑隨時(shí)間變化趨勢(shì)Fig.3 The averaged particle size over simulation time for distributed parameter MC algorithm and the τ-leap accelerated algorithm

理論和實(shí)驗(yàn)結(jié)果表明基于MC 動(dòng)力學(xué)的模擬算法具有較高的精度[25].我們以分布參數(shù)MC 算法為對(duì)比的基準(zhǔn),對(duì)比τ-leap 加速算法和離散差分?jǐn)?shù)值算法的計(jì)算結(jié)果,見圖4.其中每個(gè)子圖中繪制了該子磨機(jī)模型的輸入–輸出的累積粒度分布曲線.進(jìn)而,將計(jì)算結(jié)果分為兩組比較:第1 組為分布參數(shù)MC 算法和離散差分算法的比較,第2 組為分布參數(shù)MC 算法和τ-leap 加速算法的比較,結(jié)果見表1.從圖4 和表1 中的統(tǒng)計(jì)數(shù)據(jù)可看出,兩組算法的相對(duì)誤差穩(wěn)定在1% 以下,且未表現(xiàn)出明顯的誤差變化趨勢(shì).從第1 組結(jié)果可以看出,分布式參數(shù)MC算法和離散差分?jǐn)?shù)值算法的最終計(jì)算結(jié)果的差距較小.第2 組結(jié)果表明,采用τ-leap 算法對(duì)MC 模擬加速后,仍然確保了較高的最終計(jì)算精度.

2)平均絕對(duì)誤差

圖4 各子磨機(jī)系統(tǒng)的輸入和輸出物料的累積粒度分布曲線對(duì)比Fig.4 The cumulative particle size distribution of the input and output of each sub-mill

表1 離散差分?jǐn)?shù)值算法和τ-leap 方法相對(duì)于分布參數(shù)MC 算法的相對(duì)誤差比較Table 1 The relative errors of the discrete numerical algorithm and the τ-leap algorithm with the MC algorithm

τ-leap 方法確保準(zhǔn)確性的充要條件是傾向函數(shù)在一段時(shí)間內(nèi)不發(fā)生明顯變化.傾向函數(shù)的變化可通過參數(shù)反映.增大相當(dāng)于增大τ-leap 方法充要條件的判斷閾值,進(jìn)而可以更大幅度地加快模擬速度,但同時(shí)也會(huì)帶來更大的誤差;反之,加速幅度降低,但保持較高的精度.

其中,pmc(t)為t時(shí)刻下分布式MC 方法的粒度分布輸出,ptau(t)是基于τ-leap 的加速方法得到的t時(shí)刻的輸出,為1 范數(shù).

3)加權(quán)變異系數(shù)

式(16)定義的加權(quán)變異系數(shù)是對(duì)過程隨機(jī)性的衡量指標(biāo),它能夠反映模擬算法對(duì)于微觀動(dòng)力學(xué)過程的模擬精度.表3 所示為出口段子磨機(jī)輸出在采用τ-leap 的加速方法后,其加權(quán)變異系數(shù)WSCV與分布式MC 方法之誤差在不同時(shí)間段的值,其中列出了τ-leap 方法下不同取值得到的結(jié)果.可以看出兩種方法的WSCV 值均會(huì)隨著模擬過程減小,這主要是由于破裂過程導(dǎo)致系統(tǒng)內(nèi)整體顆粒個(gè)數(shù)逐漸增加,即式(16)右邊括號(hào)內(nèi)減去的S(i)x(i)y一項(xiàng).但是,WSCV 的減小速度會(huì)越來越慢,這是因?yàn)橄到y(tǒng)內(nèi)的顆粒的數(shù)量增加到一定數(shù)目以后,顆粒粒度的分布范圍收窄,從而部分抵消了增加顆粒的數(shù)目對(duì)于WSCV 減小的貢獻(xiàn),即式(16)右邊括號(hào)內(nèi)的求和項(xiàng).具體分析過程參見[21].

圖5 不同值所對(duì)應(yīng)的Mp值隨時(shí)間變化圖Fig.5 The values of Mpduring the simulation under different

表2 基于τ-leap 的加速算法在不同取值下與分布式參數(shù)MC 算法的Mp的相對(duì)誤差(%)Table 2 The relative difference of Mp(%)of the τ-leap acceleration algorithm comparing with the MC algorithm

表2 基于τ-leap 的加速算法在不同取值下與分布式參數(shù)MC 算法的Mp的相對(duì)誤差(%)Table 2 The relative difference of Mp(%)of the τ-leap acceleration algorithm comparing with the MC algorithm

表3 基于τ-leap 的加速算法在不同取值下與分布式參數(shù)MC 算法的Mp的相對(duì)誤差(%)Table 3 The relative difference of Mp(%)of the τ-leap acceleration algorithm comparing with the MC algorithm

表3 基于τ-leap 的加速算法在不同取值下與分布式參數(shù)MC 算法的Mp的相對(duì)誤差(%)Table 3 The relative difference of Mp(%)of the τ-leap acceleration algorithm comparing with the MC algorithm

圖6 基于τ-leap 的加速算法在不同取值下的WSCV 與分布式參數(shù)MC 算法相對(duì)誤差Fig.6 The relative difference of WSCV of the τ-leap acceleration algorithm comparing with the MC algorithm

此外,由于本文所描述的對(duì)象為“冷車”啟動(dòng),初始時(shí)磨機(jī)中并沒有待研磨物料,因此在模擬初期顆粒數(shù)目比較少,使得傾向函數(shù)較小,當(dāng)顆粒發(fā)生消失事件時(shí),傾向函數(shù)產(chǎn)生很大的波動(dòng),破壞了τ-leap條件.所以,可以從圖中看到,τ-leap 方法與MC 方法相比的WSCV 結(jié)果在模擬初始時(shí)存在較大的差距,但隨著模擬時(shí)間的進(jìn)行,由于系統(tǒng)內(nèi)的顆?？倲?shù)目的增多,兩種方法的差距越來越小.這說明:1)基于τ-leap 的加速模擬算法的結(jié)果會(huì)隨著模擬時(shí)間不斷逼近相應(yīng)的MC 模擬結(jié)果;2)在相對(duì)較短的模擬周期內(nèi),二者的相對(duì)誤差較大.以上兩點(diǎn)趨勢(shì)可以從圖6 中更明顯地看出.

3.2.2 τ-leap 算法的加速性能

本節(jié)主要考察基于τ-leap 算法的加速性能.主要根據(jù)式(15)定義的穩(wěn)態(tài)模擬速度指標(biāo)TSS 來分析.由于磨機(jī)中顆粒的移動(dòng)和破裂的隨機(jī)性導(dǎo)致TSS 曲線波動(dòng)較大.所以對(duì)每個(gè)0.1 分鐘模擬時(shí)間段內(nèi)的瞬時(shí)TSS 求平均值,結(jié)果見表4.可以看到采用τ-leap 加速方法能夠明顯提高模擬計(jì)算的速度,如圖7 所示.其中,調(diào)節(jié)因子可以用于調(diào)節(jié)加速的性能.例如,當(dāng)取0.01 時(shí),基于TSS 計(jì)算的加速比最大可達(dá)227 倍,最小也可以取得30 倍的加速.減小會(huì)降低加速比,極端情況下,如當(dāng)取0.0001時(shí),τ-leap 幾乎退化成MC 算法.此外,可以看到τ-leap 的加速比隨著模擬過程越來越大.這是由于MC 和τ-leap 算法本身的特點(diǎn)導(dǎo)致的.由文獻(xiàn)[21]的分析可知,基于MC 的顆粒破裂動(dòng)力學(xué)模擬的速度隨著系統(tǒng)中顆粒數(shù)目的增加而衰減;而τ-leap 算法的效率幾乎不受系統(tǒng)顆粒數(shù)目的影響.這解釋了我們看到的加速比越來越大的現(xiàn)象.

圖7 基于τ-leap 的加速算法相對(duì)于分布式參數(shù)MC 算法的TSS 加速倍數(shù)Fig.7 The magnitude of speedup achieved by τ-leap over the MC algorithm

4 結(jié)論

磨礦是選礦過程中最重要的工藝過程之一,對(duì)于選礦過程的運(yùn)行品質(zhì)至關(guān)重要.而物料的粒度分布是反映磨礦過程中顆粒破碎情況的直接工藝指標(biāo).本文以磨礦過程為背景,介紹了兩種顆粒研磨過程的動(dòng)力學(xué)模擬算法.首先提出一種分布式參數(shù)的連續(xù)磨蒙特卡洛動(dòng)力學(xué)模擬算法,進(jìn)而針對(duì)該算法效率較低的問題,提出了基于τ-leap 的模擬加速算法.本文的主要貢獻(xiàn)是:

表4 基于τ-leap 的加速算法相對(duì)于分布式參數(shù)MC 算法的TSS 加速倍數(shù)Table 4 The magnitude of speedup achieved by τ-leap over the MC algorithm

1)提出分布式參數(shù)蒙特卡洛動(dòng)力學(xué)模擬算法.將磨機(jī)分為多個(gè)串級(jí)連接的子磨機(jī).對(duì)于每段子磨機(jī),徑向上劃分為若干粒級(jí),假設(shè)子磨機(jī)內(nèi)物料均勻混合而且顆粒在子磨機(jī)之間移動(dòng)是在破裂之后發(fā)生的.那么子磨機(jī)內(nèi)的破裂過程的模擬采用批次磨的時(shí)間驅(qū)動(dòng)蒙特卡洛模擬方法實(shí)現(xiàn).根據(jù)分段思想,相鄰的子磨機(jī)間有物料的相互傳遞.因此,子磨機(jī)內(nèi)顆粒數(shù)目的增加是由大顆粒破裂或者子磨機(jī)間顆粒的流入而導(dǎo)致的.從微觀角度來觀察,磨機(jī)內(nèi)顆?？赡馨l(fā)生三種事件:向前移動(dòng)事件,是由對(duì)流和物料本身的循環(huán)流動(dòng)引起;向后移動(dòng)事件,是由擴(kuò)散作用引起;破裂事件,是由磨機(jī)研磨作用而導(dǎo)致顆粒破裂至比其粒級(jí)小的粒級(jí)內(nèi).將上述三種事件統(tǒng)稱為顆粒消失事件.以此為基礎(chǔ)定義系統(tǒng)狀態(tài)矩陣、傾向函數(shù)和基于蒙特卡洛的模擬計(jì)算方法.

2)基于τ-leap 的磨礦過程分布式參數(shù)蒙特卡洛模擬加速算法.分布式參數(shù)蒙特卡洛動(dòng)力學(xué)模擬算法能夠精確刻畫每一個(gè)子磨機(jī)的動(dòng)態(tài)變化過程,進(jìn)而得到整個(gè)磨機(jī)中不同時(shí)間不同位置各個(gè)時(shí)間點(diǎn)的粒度分布.但由于該種方法需要記錄每一次事件的發(fā)生,計(jì)算量極大,只能應(yīng)用于小體系的模擬.針對(duì)這個(gè)問題,我們提出基于τ-leap 的模擬加速方法.將模擬時(shí)間劃分為數(shù)個(gè)時(shí)間片作為步進(jìn)時(shí)間,同時(shí)假設(shè)在每個(gè)步進(jìn)時(shí)間內(nèi)有多個(gè)反應(yīng)事件同時(shí)發(fā)生.即不必每次事件都更新系統(tǒng)狀態(tài),而是根據(jù)步進(jìn)時(shí)間片內(nèi)多個(gè)事件的累積效應(yīng)對(duì)系統(tǒng)狀態(tài)進(jìn)行更新,從而提高模擬速度.為了解決狀態(tài)不一致的問題,我們創(chuàng)新性地提出了一種基于緩沖區(qū)的解決方法.考慮到移動(dòng)事件和破裂事件直接的獨(dú)立性,對(duì)每個(gè)子磨機(jī)系統(tǒng)開辟兩個(gè)狀態(tài)同步緩沖區(qū):前向移動(dòng)緩沖區(qū)和后向移動(dòng)緩沖區(qū).在計(jì)算過程中,將前向移動(dòng)和后向移動(dòng)的顆粒分別暫存在緩沖區(qū),等到達(dá)下一個(gè)同步時(shí)間點(diǎn)后,再將緩沖區(qū)的顆粒與相鄰子磨機(jī)進(jìn)行交換.本質(zhì)上,這相當(dāng)于將前/后向的移動(dòng)事件和破裂事件分別進(jìn)行處理.相應(yīng)地,在進(jìn)行破裂事件的處理中,將子磨機(jī)和其附屬的緩沖區(qū)視為一個(gè)獨(dú)立的批次磨過程,并更新狀態(tài)轉(zhuǎn)移矩陣和傾向函數(shù).

對(duì)本文工作加以總結(jié)如下:

1)本文提出徑向離散化的模型框架來解決磨機(jī)連續(xù)磨過程的分布式參數(shù)蒙特卡洛動(dòng)力學(xué)模擬問題,從仿真結(jié)果來看,算法能夠精確地刻畫顆粒的隨機(jī)破碎和徑向方向的流動(dòng);但是蒙特卡洛動(dòng)力學(xué)模擬的計(jì)算速度遠(yuǎn)大于差分?jǐn)?shù)值算法.

2)針對(duì)分布式參數(shù)蒙特卡洛動(dòng)力學(xué)算法的速度問題,提出基于τ-leap 加速算法來解決這一問題.通過仿真實(shí)驗(yàn)以及算法精度與速度的比較分析可知,基于τ-leap 的方法在本質(zhì)上是一種對(duì)蒙特卡洛動(dòng)力學(xué)算法的精度和速度的平衡機(jī)制,通過調(diào)節(jié)參數(shù)可以實(shí)現(xiàn)速度和精度的取舍.

3)本文所提出的方法,可以用于球磨機(jī)或棒磨機(jī)的分布參數(shù)動(dòng)力學(xué)模擬.