劉 歡 任鳳玉 何榮興 李廣輝(東北大學(xué)資源與土木工程學(xué)院，遼寧 沈陽 110819)

PFC(Particle Flow Code)屬于離散單元法，是一種適用于解決非連續(xù)性問題的數(shù)值方法，它僅由顆粒之間的微觀參數(shù)來表現(xiàn)各種材料的宏觀力學(xué)行為[1]。由于PFC可用來模擬礦巖散體的運(yùn)動、旋轉(zhuǎn)以及散體間相互作用等非線性行為，目前已被廣泛地應(yīng)用于研究礦巖散體的各項(xiàng)力學(xué)行為以及采礦工藝中。

王培濤等[2-3]借助PFC研究了無底柱分段崩落法的相關(guān)放礦規(guī)律，得出顆粒間摩擦系數(shù)越大放出橢球體的偏心率越大，并且研究了炮孔邊孔角與端壁所受側(cè)壓力的關(guān)系；程愛平等[4]應(yīng)用PFC模擬了金山店鐵礦低貧化放礦工藝，證明在金山店推行低貧化放礦方法完全可行；孫浩等[5-6]應(yīng)用PFC研究了復(fù)雜邊界和多放礦口條件下崩落礦巖流動特性，表明在復(fù)雜邊界以及多放礦口條件下放出體形態(tài)會產(chǎn)生不同程度的變異，并不是一個規(guī)則的橢球體；鄒曉甜等[7]研究了溜井底部放礦漏斗角對礦石流動性的影響，分析不同放礦漏斗角下的礦石流動性特征，揭示了放礦漏斗角對礦石流動性的影響機(jī)理；井伯祥等[8]研究了覆蓋層中散體顆粒的移動規(guī)律，研究表明覆蓋層內(nèi)小顆粒的平均速度大于大顆粒的速度，且進(jìn)路上方的顆粒速度較快。這些研究均表明，PFC較好地解決了與采礦有關(guān)的各種礦巖散體問題，然而PFC 的原理是通過設(shè)置每個接觸的本構(gòu)模型實(shí)現(xiàn)對材料宏觀性質(zhì)的模擬，如何恰當(dāng)選擇顆粒間的接觸本構(gòu)模型以及模型中的細(xì)觀參數(shù)，使PFC中顆粒的性質(zhì)與礦巖散體的性質(zhì)相匹配，對于研究礦巖散體的各種力學(xué)行為將有重要的意義。

1 PFC接觸本構(gòu)模型

PFC的基本原理是采用介質(zhì)最基本單元(顆粒)、最基本的顆粒運(yùn)動法則(牛頓第二定律)以及顆粒間的接觸本構(gòu)模型來描述介質(zhì)的復(fù)雜力學(xué)行為。不同顆粒間通過接觸點(diǎn)處的力和力矩產(chǎn)生相互作用，在特定條件下不斷更新顆粒間的接觸力和力矩，并通過離散元的方法計(jì)算系統(tǒng)的時間演化，從而為牛頓運(yùn)動方程提供了一種顯式的動態(tài)解決方案。顆粒間接觸本構(gòu)模型體現(xiàn)了接觸顆粒之間的力學(xué)特征并控制著顆粒間相互作用力和力矩的更新。因此，能正確反映礦巖散體的性質(zhì)，需要選擇合理的接觸本構(gòu)模型及其細(xì)觀參數(shù)。目前PFC提供了9種接觸本構(gòu)模型，其中常用來模擬無黏連礦巖散體的接觸本構(gòu)模型有3種，分別為Line Model(線性接觸模型)、Hertz Contact Model(赫茲接觸模型)和Rolling Resistance Linear Model(抗轉(zhuǎn)動線性接觸模型)。當(dāng)然PFC也允許用戶自定義接觸本構(gòu)模型來實(shí)現(xiàn)材料的宏觀本構(gòu)性質(zhì)。

礦巖散體由大量的非均勻巖塊組成，這些巖塊很難被定量地描述其幾何形態(tài)，雖然PFC也可以通過Clump(顆粒簇)將小顆粒組合在一起形成不同形狀的巖塊或?qū)㈩w粒直接黏結(jié)在一起形成巖塊，但去描述這些巖塊的幾何形態(tài)以及形成如此大量的礦巖散體均不太現(xiàn)實(shí)，同時這也會極大地降低計(jì)算機(jī)的計(jì)算速度，尤其當(dāng)顆粒數(shù)目很龐大時模擬極有可能失真。但在離散元方法中可以用滾動摩擦系數(shù)的大小來模擬顆粒的形狀[9]，PFC提供的眾多接觸本構(gòu)模型中，Rolling Resistance Linear Model(抗轉(zhuǎn)動線性接觸模型)[10]增加了抗轉(zhuǎn)動系數(shù)，原理是在接觸點(diǎn)上增加了與接觸顆粒間相對轉(zhuǎn)動時呈線性增加的內(nèi)力矩，這會降低顆粒的轉(zhuǎn)動能力，其與非均勻巖塊的性質(zhì)極其接近。因此，在模擬礦巖散體時推薦采用抗轉(zhuǎn)動線性接觸模型或是將線性接觸模型、赫茲接觸模型和抗轉(zhuǎn)動線性接觸模型組合應(yīng)用分別賦予不同性質(zhì)的散體顆粒?？罐D(zhuǎn)動線性接觸模型的力-位移方程為

Fc=Fl+Fd,

(1)

Mc=Mr,

(2)

式中，F(xiàn)c為顆粒所受到的接觸力，N；Mc為接觸力矩，N·m；Fl為線性接觸力，N；Fd為阻尼力，N；Mr為抗轉(zhuǎn)動力矩，N·m。

在模擬計(jì)算過程中隨著顆粒的運(yùn)動，線性接觸力與阻尼力隨著顆粒間法向與切向接觸位移量的改變而更新，抗轉(zhuǎn)動力矩隨著顆粒間相對轉(zhuǎn)動增量的變化而更新?？罐D(zhuǎn)動線性接觸模型中細(xì)觀參數(shù)主要包括顆粒的有效模量E*，法向和切向剛度比k*，摩擦系數(shù)fric，抗轉(zhuǎn)動系數(shù)rr_fric。

2 礦巖散體物理力學(xué)性質(zhì)標(biāo)定方法

礦巖散體的物理力學(xué)性質(zhì)主要包括散體的密度、塊度及級配、內(nèi)摩擦角和黏聚力，其中散體的密度、塊度及級配可以直接設(shè)定具體值，散體的內(nèi)摩擦角和黏聚力可以借助數(shù)值剪切試驗(yàn)直接獲得，也可以間接通過放出體或自然安息角的大小來反映散體的力學(xué)性質(zhì)。

2.1 數(shù)值剪切試驗(yàn)

剪切試驗(yàn)是測定松散礦巖內(nèi)摩擦角和黏聚力的常用方法，在已知礦巖散體內(nèi)摩擦角、黏聚力以及剪切位移與剪切應(yīng)力變化曲線時，可以借助數(shù)值剪切試驗(yàn)將接觸本構(gòu)模型中細(xì)觀參數(shù)與礦巖散體的物理力學(xué)性質(zhì)參數(shù)相匹配，PFC剪切試驗(yàn)如圖1所示。然而由于礦巖散體的平均塊度較大，目前常規(guī)的剪切試驗(yàn)裝置并不能測定礦巖散體的內(nèi)摩擦角和黏聚力，因此通過數(shù)值剪切試驗(yàn)標(biāo)定接觸本構(gòu)模型中細(xì)觀參數(shù)的方法具有一定的局限性。

圖1 PFC剪切試驗(yàn)Fig.1 Shearing test of PFC

2.2 放出體數(shù)值試驗(yàn)

放出體可以間接反映散體的物理力學(xué)性質(zhì)，放礦過程中礦巖移動的產(chǎn)物就是放出體。放出體是指從出礦口放出的礦石在采場崩落礦巖堆體中原來占有空間位置所構(gòu)成的形體[11]。放出體的形態(tài)與大小反映著礦巖的移動規(guī)律以及散體的流動能力，是散體物理力學(xué)性質(zhì)的一個綜合指標(biāo)。因此可以通過數(shù)值放礦試驗(yàn)，將獲得的放出體形態(tài)和大小與實(shí)際放出體相比較，當(dāng)一致時即可將接觸本構(gòu)模型中細(xì)觀參數(shù)與礦巖散體的物理力學(xué)性質(zhì)參數(shù)相匹配[6]，PFC放出體試驗(yàn)如圖2所示。目前，關(guān)于測定放出體的試驗(yàn)均是在試驗(yàn)室內(nèi)獲得，Castro R等[12]雖然進(jìn)行了較大模型的放礦試驗(yàn)，但與實(shí)際采場放礦還有一定的差距。因此研究放礦的基本規(guī)律時可應(yīng)用該方法，但研究采礦現(xiàn)場問題時該方法并不適用。

圖2 PFC放出體試驗(yàn)Fig.2 Test of PFC drawn-out ore body

2.3 自然安息角數(shù)值試驗(yàn)

自然安息角是指自然濕度下的松散礦巖，在某一特定條件下堆積，其自然坡面和水平面所形成的最大傾角[13]。自然安息角是松散礦巖顆粒間摩擦角、黏聚力和重力平衡的一個重要標(biāo)志，同時自然安息角也是試驗(yàn)室模擬礦巖散體與現(xiàn)場礦巖散體力學(xué)相似的準(zhǔn)則[13]。因此通過自然安息角數(shù)值試驗(yàn)，當(dāng)測得的值與實(shí)際自然安息角相同時，即可將接觸本構(gòu)模型中細(xì)觀參數(shù)與礦巖散體的物理力學(xué)性質(zhì)參數(shù)相匹配。然而對于同一種礦巖散體使用不同的測定方法時，測得的自然安息角也存在差異，需要依據(jù)實(shí)際的測定條件、測定方法或采礦方法來建立與之匹配的數(shù)值模型。根據(jù)不同的測定條件，可建立如下幾種數(shù)值模型。

(1)有靜壓力，無邊壁的影響。該條件下常采用圓筒式測定裝置測定散體的自然安息角，因此建立的數(shù)值模型如圖3所示。測得的自然安息角可用來標(biāo)定卸礦場、廢石場中散體的物理力學(xué)性質(zhì)。

圖3 PFC圓筒測定試驗(yàn)Fig.3 Cylinder measurement test of PFC

(2)有靜壓力，有邊壁的影響。該條件下常采用載壓式測定裝置測定散體的自然安息角，建立的模型如圖4所示。其中關(guān)于壓力的模擬可借助PFC伺服控制系統(tǒng)控制頂部墻體與顆粒間的接觸力，以此來模擬施加的壓力。測得的自然安息角可用來標(biāo)定有底柱采礦法、無底柱采礦法以及留礦采礦法(大量放礦的初期)中散體的物理力學(xué)性質(zhì)。

圖4 PFC載壓測定試驗(yàn)Fig.4 Pressure measurement test of PFC

(3)無靜壓力，有邊壁的影響。該測定條件下常采用旋轉(zhuǎn)式測定裝置測定散體的自然安息角，建立的數(shù)值模型如圖5所示。測得的自然安息角可用來標(biāo)定空場采礦法、留礦采礦法(大量放礦的后期)中礦巖散體的物理力學(xué)性質(zhì)。

圖5 PFC旋轉(zhuǎn)測定試驗(yàn)Fig.5 Rotating measurement test of PFC

(4)有滾動和滑動摩擦力的影響。該測定條件下采用塌落式測定裝置，建立的數(shù)值模型如圖6所示。測得的自然安息角可用來標(biāo)定礦倉、溜井中散體的物理力學(xué)性質(zhì)。

圖6 PFC塌落式測定試驗(yàn)Fig.6 Caving measurement test of PFC

綜上所述，在標(biāo)定接觸本構(gòu)模型的細(xì)觀參數(shù)時，推薦采用與測定條件相符合的自然安息角數(shù)值試驗(yàn)。由于PFC接觸模型中接觸的類型包括顆粒與顆粒以及顆粒與墻之間的接觸，故采用數(shù)值旋轉(zhuǎn)式模型來研究PFC接觸本構(gòu)模型的細(xì)觀參數(shù)對散體自然安息角的影響。

3 礦巖散體細(xì)觀參數(shù)的標(biāo)定

抗轉(zhuǎn)動線性接觸模型中有效模量E*與巖石顆粒的變形指標(biāo)彈性模量有關(guān)，并且隨著彈性模量的增加而增加，法向和切向剛度比k*與巖石的泊松比有關(guān)，這2個參數(shù)可通過相關(guān)的力學(xué)試驗(yàn)獲得。在其他參數(shù)一定時，通過調(diào)整摩擦系數(shù)fric和抗轉(zhuǎn)動系數(shù) rr_fric完成模擬礦巖散體的PFC細(xì)觀參數(shù)標(biāo)定。數(shù)值旋轉(zhuǎn)式模型如圖5，模型尺寸為200 mm×400 mm，顆粒的半徑為1～3 mm、3～5 mm、5～7 mm、7～10 mm且所占的體積份數(shù)分別為55%、25%、15%、5%，密度為2 700 kg/m3，有效模量為1.0×107Pa，法向和切向剛度比為1.0。

3.1 摩擦系數(shù)

分別考慮顆粒的摩擦系數(shù)與墻的摩擦系數(shù)對自然安息角的影響，摩擦系數(shù)選取0.1、0.3、0.5、0.7、0.9，抗轉(zhuǎn)動系數(shù)為0.2。結(jié)果繪制成摩擦系數(shù)與自然安息角關(guān)系曲線(如圖7所示)。由圖7可知，隨著顆?；驂Φ哪Σ料禂?shù)的增大，自然安息角均呈現(xiàn)出先增大后基本穩(wěn)定的趨勢，但顆粒的摩擦系數(shù)對自然安息角有較顯著的影響。

圖7 摩擦系數(shù)與自然安息角的變化曲線Fig.7 Curve of the friction coefficient with the repose angle

3.2 抗轉(zhuǎn)動系數(shù)

分別考慮顆粒的抗轉(zhuǎn)動系數(shù)與墻的抗轉(zhuǎn)動系數(shù)對自然安息角的影響，抗轉(zhuǎn)動系數(shù)選取0.0、0.2、0.4、0.6、0.8，摩擦系數(shù)為0.4。結(jié)果繪制成抗轉(zhuǎn)動系數(shù)與自然安息角關(guān)系曲線(如圖8所示)。由圖8可知，隨著顆?；驂Φ目罐D(zhuǎn)動系數(shù)的增大，自然安息角均呈現(xiàn)出先增大后基本穩(wěn)定的趨勢，但顆粒的抗轉(zhuǎn)動系數(shù)對自然安息角有較顯著的影響。

3.3 摩擦系數(shù)和抗轉(zhuǎn)動系數(shù)

為了進(jìn)一步研究摩擦系數(shù)與抗轉(zhuǎn)動系數(shù)對自然安息角的影響，分別選取顆粒的摩擦系數(shù)為0.1、0.3、0.5、0.7、0.9，抗轉(zhuǎn)動系數(shù)為0.2、0.4、0.6，保證墻的摩擦系數(shù)和抗轉(zhuǎn)動系數(shù)不變，模擬結(jié)果詳見圖9。由圖9可知，隨著顆粒的摩擦系數(shù)和抗轉(zhuǎn)動系數(shù)的增大，自然安息角呈增長的趨勢，且不同的摩擦系數(shù)與抗轉(zhuǎn)動系數(shù)組合可以獲得相同的自然安息角。圖9中自然安息角的范圍從22°～50°，與實(shí)際松散礦巖自然安息角相符。

圖8 抗轉(zhuǎn)動系數(shù)與自然安息角的變化曲線Fig.8 Curve of the rolling resistance coefficient with the repose angle

圖9 抗轉(zhuǎn)動系數(shù)和摩擦系數(shù)對自然安息角的影響Fig.9 Influence of rolling resistance coefficient and friction coefficient on the repose angle

4 結(jié) 論

(1)根據(jù)PFC接觸本構(gòu)模型并結(jié)合礦巖散體本身的特點(diǎn)，提出采用抗轉(zhuǎn)動線性接觸模型或?qū)⒖罐D(zhuǎn)動線性接觸模型、線性接觸模型以及赫茲接觸模型組合應(yīng)用分別賦予不同性質(zhì)的散體顆粒。

(2)闡述了數(shù)值剪切試驗(yàn)、放出體數(shù)值試驗(yàn)以及自然安息角數(shù)值試驗(yàn)的優(yōu)缺點(diǎn)，推薦采用與測定條件相符合的自然安息角數(shù)值試驗(yàn)來標(biāo)定PFC細(xì)觀參數(shù)。

(3)通過自然安息角數(shù)值試驗(yàn)分別研究了抗轉(zhuǎn)動線性接觸模型中顆粒及墻的摩擦系數(shù)和抗轉(zhuǎn)動系數(shù)與散體自然安息角的關(guān)系。研究表明：隨著顆?；驂Φ哪Σ料禂?shù)和抗轉(zhuǎn)動系數(shù)的增加，自然安息角均呈現(xiàn)出先增大后趨于穩(wěn)定的趨勢，且顆粒的摩擦系數(shù)和抗轉(zhuǎn)動系數(shù)對自然安息角有較顯著的影響。

(4)模擬結(jié)果中自然安息角的變化范圍為22°～50°，與實(shí)際松散礦巖自然安息角相符。

[1] 張社榮，孫 博，王 超，等.雙軸壓縮試驗(yàn)下巖石裂紋擴(kuò)展的離散元分析[J].巖石力學(xué)與工程學(xué)報(bào)，2013，32(S2)：3083-3091.

Zhang Sherong，Sun Bo，Wang Chao，et al.Discrete element analysis of crack propagation in rocks under biaxial compression[J].Chinese Journal of Rock Mechanics and Engineering，2013，32(S2)：3083-3091.

[2] 王培濤，楊天鴻，柳小波.無底柱分段崩落法放礦規(guī)律的PFC2D模擬仿真[J].金屬礦山，2010(8)：123-127.

Wang Peitao，Yang Tianhong，Liu Xiaobo.PFC2D numerical simulation of ore drawing rule with pillarless sublevel caving[J].Metal Mine，2010(8)：123-127.

[3] 王培濤，楊天鴻，柳小波.邊孔角對無底柱分段崩落法放礦影響的顆粒流數(shù)值模擬研究[J].金屬礦山，2010(3)：12-16.

Wang Peitao，Yang Tianhong，Liu Xiaobo.Particle flow numerical simulation investigation on influence of lateral opening angle on ore drawing with sublevel pillarless caving[J].Metal Mine，2010(3)：12-16.

[4] 程愛平，許夢國，王 平.金山店鐵礦低貧化放礦數(shù)值模擬[J].金屬礦山，2011(1)：31-34.

Cheng Aiping，Xu Mengguo，Wang Ping.A numerical simulation of low dilution ore drawing in Jinshandian Iron Mine[J].Metal Mine，2011(1)：31-34.

[5] 孫 浩，金愛兵，高永濤，等.復(fù)雜邊界條件下崩落礦巖流動特性[J].中南大學(xué)學(xué)報(bào):自然科學(xué)版，2015，46(10)：3782-3788.

Sun Hao，Jin Aibing，Gao Yongtao，et al.Flow characteristics of caved ore and rock under complex boundary conditions[J].Journal of Central South University:Science and Technology，2015，46(10)：3782-3788.

[6] Jin A B，Sun H，Ma G W,et al.A study on the draw laws of caved ore and rock using the discrete element method[J].Computers and Geotechnics，2016，80：59-70.

[7] 鄒曉甜，劉艷章，張丙濤，等.溜井底部放礦漏斗角對礦石流動性的影響研究[J].金屬礦山，2016(12)：160-164.

Zou Xiaotian，Liu Yanzhang，Zhang Bingtao，et al.Research on effect of drawing funnel angle at the bottom of ore-pass on ore fluidity[J].Metal Mine，2016(12)：160-164.

[8] 井伯祥，陳 超，張亞賓，等.覆蓋層顆粒移動規(guī)律的PFC2D數(shù)值模擬[J].金屬礦山，2016(6)：43-48.

Jing Boxiang，Chen Chao，Zhang Yabin，et al.PFC2D numerical simulation in velocity law of covering layer particles[J].Metal Mine，2016(6)：43-48.

[9] Wensrich C M，Katterfeld A.Rolling friction as a technique for modelling particle shape in DEM[J].Powder Technology，2012，217：409-417.

[10] Kazuyoshi I，Masanobu O.Rolling resistance at contacts in simulation of shear band development by DEM[J].Journal of Engineering Mechanics，1998，124(3)：285-292.

[11] 劉興國.放礦理論基礎(chǔ)[M].北京：冶金工業(yè)出版社，1995.

Liu Xingguo.Base of Ore-drawing[M].Beijing：Metallurgical Industry Press，1995.

[12] Castro R，Trueman R，Halim A.A study of isolated draw zones in block caving mines by means of a large 3D physical model[J].International Journal of Rock Mechanics and Mining Sciences，2007，44(6)：860-870.

[13] 周 驥.對松散礦石自然安息角的研究[J].有色金屬，1983，35(2)：24-29.

Zhou Ji.An investigation on angles of natural repose of bulk ores[J].Nonferrous Metals，1983，35(2)：24-29.