模擬礦巖散體的PFC細觀參數標定方法

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

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

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

1 PFC接觸本構模型

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

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

Fc=Fl+Fd,

(1)

Mc=Mr,

(2)

式中，Fc為顆粒所受到的接觸力，N；Mc為接觸力矩，N·m；Fl為線性接觸力，N；Fd為阻尼力，N；Mr為抗轉動力矩，N·m。

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

2 礦巖散體物理力學性質標定方法

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

2.1 數值剪切試驗

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

圖1 PFC剪切試驗Fig.1 Shearing test of PFC

2.2 放出體數值試驗

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

圖2 PFC放出體試驗Fig.2 Test of PFC drawn-out ore body

2.3 自然安息角數值試驗

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

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

圖3 PFC圓筒測定試驗Fig.3 Cylinder measurement test of PFC

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

圖4 PFC載壓測定試驗Fig.4 Pressure measurement test of PFC

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

圖5 PFC旋轉測定試驗Fig.5 Rotating measurement test of PFC

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

圖6 PFC塌落式測定試驗Fig.6 Caving measurement test of PFC

綜上所述，在標定接觸本構模型的細觀參數時，推薦采用與測定條件相符合的自然安息角數值試驗。由于PFC接觸模型中接觸的類型包括顆粒與顆粒以及顆粒與墻之間的接觸，故采用數值旋轉式模型來研究PFC接觸本構模型的細觀參數對散體自然安息角的影響。

3 礦巖散體細觀參數的標定

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

3.1 摩擦系數

分別考慮顆粒的摩擦系數與墻的摩擦系數對自然安息角的影響，摩擦系數選取0.1、0.3、0.5、0.7、0.9，抗轉動系數為0.2。結果繪制成摩擦系數與自然安息角關系曲線(如圖7所示)。由圖7可知，隨著顆粒或墻的摩擦系數的增大，自然安息角均呈現出先增大后基本穩定的趨勢，但顆粒的摩擦系數對自然安息角有較顯著的影響。

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

3.2 抗轉動系數

分別考慮顆粒的抗轉動系數與墻的抗轉動系數對自然安息角的影響，抗轉動系數選取0.0、0.2、0.4、0.6、0.8，摩擦系數為0.4。結果繪制成抗轉動系數與自然安息角關系曲線(如圖8所示)。由圖8可知，隨著顆粒或墻的抗轉動系數的增大，自然安息角均呈現出先增大后基本穩定的趨勢，但顆粒的抗轉動系數對自然安息角有較顯著的影響。

3.3 摩擦系數和抗轉動系數

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

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

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

4 結 論

(1)根據PFC接觸本構模型并結合礦巖散體本身的特點，提出采用抗轉動線性接觸模型或將抗轉動線性接觸模型、線性接觸模型以及赫茲接觸模型組合應用分別賦予不同性質的散體顆粒。

(2)闡述了數值剪切試驗、放出體數值試驗以及自然安息角數值試驗的優缺點，推薦采用與測定條件相符合的自然安息角數值試驗來標定PFC細觀參數。

(3)通過自然安息角數值試驗分別研究了抗轉動線性接觸模型中顆粒及墻的摩擦系數和抗轉動系數與散體自然安息角的關系。研究表明：隨著顆粒或墻的摩擦系數和抗轉動系數的增加，自然安息角均呈現出先增大后趨于穩定的趨勢，且顆粒的摩擦系數和抗轉動系數對自然安息角有較顯著的影響。

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

[1] 張社榮，孫 博，王 超，等.雙軸壓縮試驗下巖石裂紋擴展的離散元分析[J].巖石力學與工程學報，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] 王培濤，楊天鴻，柳小波.無底柱分段崩落法放礦規律的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] 王培濤，楊天鴻，柳小波.邊孔角對無底柱分段崩落法放礦影響的顆粒流數值模擬研究[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] 程愛平，許夢國，王 平.金山店鐵礦低貧化放礦數值模擬[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] 孫 浩，金愛兵，高永濤，等.復雜邊界條件下崩落礦巖流動特性[J].中南大學學報:自然科學版，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] 井伯祥，陳 超，張亞賓，等.覆蓋層顆粒移動規律的PFC2D數值模擬[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] 劉興國.放礦理論基礎[M].北京：冶金工業出版社，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.