張中衛(wèi), 程 曉， 李 潮

(長江大學(xué)城市建設(shè)學(xué)院，湖北荊州 434000)

?

PFC2D自定義模型使用研究

張中衛(wèi), 程 曉， 李 潮

(長江大學(xué)城市建設(shè)學(xué)院，湖北荊州 434000)

論文研究了基于visual stdio 2013自定義接觸模型的實現(xiàn)方法，并嘗試將自定義接觸模型用在模擬中，分析內(nèi)置接觸模型與自定義接觸模型之間替代關(guān)系。得出結(jié)論：自定義接觸模型無法與常規(guī)的顆粒生成方法一起運算，常規(guī)的顆粒生成方法指半徑擴大法、explosive生成法；自定義接觸模型只能與固定坐標(biāo)生成的顆粒一同使用，巖土模擬中受限；自定義接觸模型在顆粒接觸剛度很小的情況下可作為內(nèi)置接觸模型的一種替代。

PFC2D； 接觸模型； 細觀參數(shù)； 數(shù)值模擬； c++

1 PFC2D自定義模型介紹

PFC2D(particle flow code for 2dimention)是一種對顆粒散體進行模擬的有效數(shù)值手段。核心算法是離散元法，其中離散單元指的是每一顆小顆粒，顆粒和顆粒之間是按照粘結(jié)接觸聯(lián)系在一起的。一個完整的接觸模型包括接觸強度模塊、滑動模塊和粘結(jié)模塊三部份。其中，接觸強度模塊用于給定顆粒和顆粒接觸的法向強度和切向強度，滑動模塊是用來限定滑動強度，粘結(jié)模塊用來規(guī)定顆粒和顆粒之間的粘結(jié)是點粘結(jié)還是平行粘結(jié)?，F(xiàn)在大多數(shù)研究使用的都是PFC內(nèi)置的接觸模型。

PFC2D同時也支持用戶自己定義接觸模型，具體的做法是運用c++將接觸模型相關(guān)的函數(shù)和參量編譯為dll文件，然后在PFC2D內(nèi)置命令流中調(diào)用dll文件。自定義接觸模型相對于PFC命令流來說是相對獨立的，它有自己的名稱、數(shù)據(jù)結(jié)構(gòu)、定義類型。自定義模型與內(nèi)置接觸模型功用體現(xiàn)在被調(diào)用后可將用戶顆粒間的特征數(shù)據(jù)傳遞到顆粒粘結(jié)中并自己定義顆粒粘結(jié)破壞強度，以及粘結(jié)破壞后顆粒的軌跡。

2 自定義接觸模型國內(nèi)外研究現(xiàn)狀

PFC由于能從細觀方面研究宏觀的表現(xiàn)，從而被國內(nèi)外學(xué)者廣泛采用。但是幾乎沒有研究者用到自定義接觸模型。王剛[1]等利用PFC2D程序從宏觀和細觀角度探討節(jié)理在不同邊界條件下剪切過程中的力學(xué)演化規(guī)律和破壞機制；Axel Ruiken[2]等在指出前人關(guān)于由實驗孔隙率推導(dǎo)2D情況下的孔隙率方法的不通用性的基礎(chǔ)上提出了自己由3D孔隙率得到2D孔隙率的方法；徐文杰[3]等提出了一種基于數(shù)字圖像處理的非均質(zhì)巖土材料細觀結(jié)構(gòu)PFC2D數(shù)值計算模型自動生成方法。

莊德林[4]等首次采用了顆粒流程序PFC2D模擬垮落區(qū)充填中煤層采出、頂板垮落、充填的整個動態(tài)過程；孔祥臣[5]等針對土石混合料振動擊實過程中顆粒運動規(guī)律和密度形成機理問題引入PFC2D顆粒離散元方法,模擬了一定級配條件下不同土石比的土石混合料的振動擊實過程。

朱煥春[6]在PFC 及其在礦山崩落開采研究中的應(yīng)用一文中對PFC2D進行了功能介紹。其中提到PFC2D是幾乎用命令流操作的，所以上手難度很大，研究一般在高校中。另外也有提到“構(gòu)建 PFC 模型和進行相應(yīng)的運算準(zhǔn)備工作，必須使用 PFC 的二次開發(fā)功能，可通過編程實現(xiàn)上述目的?！比欢诮鼉赡甑难芯恐?，PFC模擬都是使用內(nèi)置的接觸模型。論文對當(dāng)下學(xué)者對PFC研究中沒有使用自定義接觸模型的原因進行研究。

3 構(gòu)建自定義接觸模型

自定義接觸是靠PFC2D主程序調(diào)用用戶編譯的dll文件來實現(xiàn)的。構(gòu)建自定義接觸的步驟有兩步：(1)取得dll文件關(guān)聯(lián)文件；(2)dll文件編譯。

3.1 取得dll文件關(guān)聯(lián)文件

編譯dll文件需要.cpp文件以及.h頭文件、外部關(guān)聯(lián)項。其中.cpp文件和.h頭文件可以從user manual中得到，直接嵌套運行不會通過，需要進行調(diào)試；外部關(guān)聯(lián)項是PFC2D程序自帶的cmodlewrv.dll和同名的object file library。

3.2 dll文件編譯

編譯dll文件的平臺需要為c++6.0以上版本，以vs2013為例，創(chuàng)建win32 dll工作項目，關(guān)聯(lián)以上提到的關(guān)聯(lián)文件，調(diào)試報錯為0后，下拉菜單生成解決方案生成dll文件。

4 不同顆粒生成法與自定義接觸模型的嵌套

PFC2D生成模型的方法有：直接用gen命令生成少量顆粒、半徑擴大法、explosive顆粒生成法、坐標(biāo)生成法。但是不同生成方法采用不同的參數(shù)值，因此可能導(dǎo)致無法在試樣中使用自定義模型。接下來對上述三種顆粒生成法與自定義接觸模型的耦合性進行比較研究。

顆粒流建模首先要確定微觀參數(shù)，通過試錯法，論文中的微觀參數(shù)的選取如表1所示[7]。

表1 微觀參數(shù)的選取

4.1 半徑擴大法

半徑擴大法的思路是先生成規(guī)定數(shù)目的顆粒填充限定的范圍，此時顆粒的半徑大小是在用戶設(shè)定半徑范圍內(nèi)的一個隨機值，顆粒與顆粒之間有重疊，然后定義半徑擴大系數(shù)，也叫乘子。隨著cycle的進行，顆粒的半徑逐漸增加，最終形成相互點接觸的顆粒試樣。

按照自定義模型的插入步驟調(diào)用編寫好的dll文件，檢查命令流可以看出需要賦予顆粒的性狀已經(jīng)完成，但是在加載階段cycle計算開始時，系統(tǒng)報錯：property 0 not found in usercm1 model。證明半徑擴大法生成的試樣缺失自定義模型需要的參數(shù)。

4.2 explosive 顆粒生成法

explosive顆粒生成法的思路是先生成規(guī)定數(shù)目的顆粒，此時生成的顆粒相對于半徑擴大法來說更加不均勻，顆粒間的重合更加明顯。explosive生成法的特點就是在cycle運算的開始，互相重疊的顆粒能迅速形成點接觸，并且這種生成方法可以控制生成試樣的孔隙率。重復(fù)半徑擴大法中的自定義模型嵌套，同樣報錯參數(shù)缺失。

4.3 坐標(biāo)生成法

坐標(biāo)生成方法是通過指定顆粒坐標(biāo)然后使用顆粒生成命令直接生成的。填充一份10×20試樣一般需要2 000左右顆粒，如果每個顆粒都要用戶自己定義坐標(biāo)會很繁瑣，所以這種顆粒生成方法適用范圍很小。

4.3.1 模型方案的建立

模型通過定義顆粒坐標(biāo)的方式生成3個顆粒，并且用4個墻體將顆粒圍起來，然后在顆粒接觸部分分別導(dǎo)入內(nèi)置接觸模型和自定義接觸模型。微觀參數(shù)的選取按照表2選取[8]，生成模型見圖1。

表2 微觀參數(shù)的選取

圖1 坐標(biāo)法建模

4.3.2 加壓與記錄

加壓采用賦予上部墻體固定速度0.001 m/s的方式，記錄顆粒間的平均非均衡力和平均接觸應(yīng)力，并繪出圖像(圖2)。

圖2 顆粒間應(yīng)力(kn=ks=1)

5 自定義接觸模型中參數(shù)影響

從圖2得出，內(nèi)置接觸模型和自定義接觸模型的均不平衡力-時步線和均接觸應(yīng)力-時步線近乎一致。當(dāng)kn=ks=1時，自定義接觸模型可以取代內(nèi)置接觸模型。

接下來參數(shù)進行修改，來確定顆粒連接參數(shù)的選取是否會影響自定義接觸模型對內(nèi)置接觸模型的取代。

5.1 顆粒法向剛度和切向剛度的影響

保持摩擦系數(shù)為0.1的情況下，保持法向剛度和切向剛度相等，并分別取為0.3 Pa、0.5 Pa、1 Pa、3 Pa,繪出均不平衡力-時步曲線和均接觸應(yīng)力-時步曲線(圖3～圖5)。

圖3 顆粒間應(yīng)力(kn=ks=0.3 Pa)

圖4 顆粒間應(yīng)力(kn=ks=0.5 Pa)

圖5 顆粒間應(yīng)力(kn=ks=3 Pa)

從圖2～圖5可以得到以下結(jié)論：顆粒間的連接剛度kn與ks相等，并且維持在一個較小的數(shù)值范圍內(nèi)(0 Pa，1 Pa)，此時可以認(rèn)為自定義接觸模型可以完全代替內(nèi)置接觸模型。當(dāng)ks=ks≥1 Pa時，內(nèi)置接觸模型和自定義接觸模型的顆粒間應(yīng)力差距較大，此時不能用自定義接觸模型代替內(nèi)置接觸模型。

另外，數(shù)值模擬表明當(dāng)kn=ks≥10 Pa時，顆粒間的應(yīng)力大小會差5個數(shù)量級，所以未將圖像未放在一起比較。

5.2 顆粒間的摩擦系數(shù)的影響

接下來通過設(shè)定摩擦系數(shù)為0.1、1.2，控制kn=ks=1 Pa來確定顆粒間的摩擦系數(shù)是否會影響自定義接觸模型對內(nèi)置接觸模型的代替(圖6、圖7)。

圖6 顆粒間應(yīng)力(fric=0.1)

圖7 顆粒間應(yīng)力(fric=1.2)

圖6和圖7進行對比，摩擦系數(shù)從0.1上升成1.2后，接觸應(yīng)力-時步曲線的峰值和線型都發(fā)生了明顯變化，證明摩擦系數(shù)的變動也影響自定義接觸模型對內(nèi)置接觸模型的替代。

6 結(jié) 論

在PFC2D用戶手冊中，自定義接觸模型作為二次開發(fā)手段被提出來，它的建立需要用戶掌握一定的c++語言，給PFC2D的建模增加了難度，并且經(jīng)過數(shù)據(jù)分析，自定義接觸取代內(nèi)置接觸模型的條件也較苛刻，所以在近幾年的PFC2D的研究中沒有關(guān)于自定義接觸模型的使用和論述。另外通過參數(shù)控制得出以下結(jié)論：

(1)使用自定義接觸模型替代內(nèi)置接觸模型的條件是：顆粒位置是通過坐標(biāo)定義生成；顆粒間的連接剛度kn與ks相等，并且維持在一個較小的數(shù)值范圍內(nèi)(0 Pa，1 Pa)。

(2)自定義接觸模型無法與半徑擴大法生成的顆粒一起運算也無法與explosive法生成的顆粒一同運算。

(3)顆粒間摩擦系數(shù)的變化影響自定義接觸模型對內(nèi)置接觸模型的替代。

[1] 王剛，張學(xué)朋，賈宇靜.基于顆粒離散元法的巖石節(jié)理面剪切破壞細觀機理[J].中南大學(xué)學(xué)報， 2015， 46(4)：1-3.

[2] Zhijie Wang,Axel Ruiken ,Felix Jacobs and Martin Ziegler. A new suggestion for determining 2D porosities in DEM studies[J]. Geomechanics and Engineering, 2014,7(6): 665-678.

[3] 徐文杰,胡瑞林,王艷萍. 基于數(shù)字圖像的非均質(zhì)巖土材料細觀結(jié)構(gòu)PFC2D模型[J].煤炭學(xué)報， 2007， 32(4)：1-3.

[4] 莊德林,李興尚,許家林. 垮落區(qū)注漿充填壓實特性的PFC2D模擬試驗[J].采礦與安全工程學(xué)報， 2008， 25(1)：1-3.

[5] 孔祥臣, 陳謙應(yīng),賈學(xué)明. 土石混合料振動擊實試驗的PFC2D模擬研究 [J].重慶交通學(xué)院學(xué)報, 2005, 24(1)：1-4.

[6] 朱煥春. PFC及其在礦山崩落開采研究中的應(yīng)用[J].巖石力學(xué)與工程學(xué)報， 2006， 25(9)：1-3.

[7] 彭博，陳玉明，袁利偉. 基于PFC2D馬鞍山鐵礦2號露天采場邊坡穩(wěn)定性分析[J].采礦工程， 2014， 35(12)：2-3.

[8] Cundall P A. PFC2D User's Manual(version3.0)[Z].Minnesota:Itasca Consulting Group Inc.

TU441+.3

A

[定稿日期]2016-04-17

[作者信息]張中衛(wèi)(1992～)， 男， 碩士研究生，研究方向為巖土數(shù)值模擬。