PFC2D平節理模型細觀參數標定方法

張寶玉， 張昌鎖*， 王晨龍， 郝保欽

(1.太原理工大學 礦業工程學院，太原 030024；2.太原理工大學 應用力學研究所，太原 030024)

1 引 言

顆粒流程序PFC屬于離散元法，可以解決非連續介質力學問題，并能從細觀角度分析巖石的破壞機理。PFC用顆粒集合體來模擬巖石，通過定義顆粒間的黏結接觸模型使得顆粒產生相互作用來重現巖石的力學性質。起初一般采用平行黏結模型(PBM)來模擬巖石，但研究發現PBM模擬巖石有如下缺陷[1-3]。壓拉比過低，不能滿足實際巖石的壓拉比要求；內摩擦角偏小。為此，Potyondy[4]提出了平節理模型(FJM)，將圓形顆粒構造成多邊形，可抑制接觸破壞后顆粒的旋轉，能夠模擬出大壓拉比。因此，本文采用FJM來模擬巖石。

2 平節理模型

平節理模型(FJM)是由晶粒及晶粒間賦予平節理接觸(FJC)構成的，晶粒由圓盤顆粒和名義面(notional surfaces)組成，因此FJM相當于把顆粒構造成了多邊形，可提供足夠的自鎖效應，抑制接觸破壞后顆粒的旋轉，能夠模擬大壓拉比情況。晶粒間的接觸實際上是名義面間的接觸，接觸界面(interface)位于名義面之間。在PFC2D中，FJC是一條線段，該線段平均分成多段，每段代表一個單元，如圖1所示。每個單元的狀態是獨立的，可以是黏結的或是非黏結的。黏結單元作用機理為，當其所受法向拉應力σ+>張拉強度σf時，單元破壞，生成張拉裂紋；當其所受切向應力τ>τf時，單元破壞，生成剪切裂紋。非黏結單元作用機理為，不能承受拉力；當其所受切向應力τ>τf時，沿接觸界面發生滑移。其中，黏結狀態和非黏結狀態下的剪切強度τf分別見式(1,2)，

圖1 平節理接觸模型[4]

(1，2)

3 正交試驗數值模擬

3.1 數值試驗

進行數值試驗包括建立模型和加載兩個部分。

(1) 建立模型

首先，生成上下左右4面墻，墻圍成的區域大小即為模型尺寸(50 mm×100 mm)；然后，按顆粒尺寸及分布形式在墻圍成區域內生成重疊量較大的顆粒，計算彈開顆粒并使模型平衡；接著，通過伺服機制移動墻給模型施加一定的應力，使得顆粒體系均勻；模型中會存在少量懸浮顆粒(顆粒上的接觸數少于3)，通過縮放這些懸浮顆粒的大小使其上的接觸數≥3來達到消除懸浮顆粒的目的；最后，給模型接觸賦予細觀參數并將模型狀態清零(包括顆粒速度和位移等)。建模完成后的數值模型如圖2所示。

圖2 數值模型

(2) 加載

單軸壓縮試驗。刪除模型兩側的墻，給上下兩面墻施加恒定速度來實現加載，通過上下墻監測軸向應力和應變，設置測量圓監測橫向應變，據此可獲取單軸抗壓強度UCS、彈性模量E及泊松比ν。

直接拉伸試驗。刪除模型中所有的墻，在試樣上下端給一定厚度顆粒恒定速度來進行加載，設置測量圓監測軸向應力應變，據此獲取抗拉強度TS。

雙軸壓縮試驗。根據伺服原理控制兩側的墻實現圍壓的施加，給上下兩面墻施加恒定速度實現加載，并通過上下墻監測軸向應力，圍壓σ3和軸向應力σ1呈線性關系，

σ1=kσ3+b

(3)

(4)

(5)

3.2 試驗參數選取

3.2.1 宏觀參數選取

(6)

巖石的起裂強度σc i是判斷圍巖損傷范圍的重要指標。在PFC中可通過裂紋數來確定巖石的σc i,文獻[1，15]將裂紋數達到峰值強度對應裂紋數的10%時所對應的應力視為σc i。

3.2.2 細觀參數選取

綜上，本文研究選取的宏細觀參數列入表1。

表1 選取的宏細觀參數

3.3 正交試驗設計

正交試驗是研究多因素多水平的試驗方法，主要是依據正交性在全面試驗中選擇具有代表性的點(試驗點分布均勻)進行試驗，用正交設計表來安排試驗并進行數據分析，具有高效性，結果可靠簡單易行[16]。

FJM中涉及很多細觀參數，正交試驗設計之前，為降低參數標定的難度，確定一些參數的值如下。

表2 因素水平

4 結果分析

表3 正交數值試驗方案及結果Tab.3 Orthogonal numerical test scheme and results

4.1 多因素方差分析

方差分析(F檢驗)用于多個樣本均數差別的顯著性檢驗。根據方差分析計算過程[19]計算出試驗因素及誤差的自由度分別為3和7，查表可知顯著性水平α=0.05時，F0.95(3,7)=4.35,顯著性水平α=0.01時，F0.99(3,7)=8.45；計算各試驗因素的F值并繪成圖3。認為當4.35<因素F值<8.45時，因素對結果有顯著影響；當因素F值> 8.45 時，因素對結果有非常顯著的影響；因素的F值越大，其對結果的影響越大。據圖3分析如下。

圖3 多因素方差分析F值

4.2 回歸分析

根據4.1節的分析，用巖石宏觀參數的主要影響因素對其進行簡單的線性擬合，結果見式(7～12)，可見只有式(7,8,10)的R2大于0.85，擬合效果較好，說明宏細觀參數間的線性關系良好；其余公式的R2則較小，效果差，說明宏細觀參數間的關系復雜，需要進一步分析，因此將試驗結果平均值變化趨勢繪于圖4，分析如下。

圖4 試驗結果平均值變化趨勢

E=0.858Ef-4.677Kf+12.021，R2=0.976

(7)

(8)

UCS=-0.037Ef-10.397Kf+31.312μf+

3.333σf+4.378(cf/σf)-41.644Rs d+

R2=0.819

(9)

Ts=0.007Ef+0.25σf-6.145Rs d+

R2=0.867

(10)

R2=0.607

(11)

R2=0.798

(12)

(1) 圖4(a)顯示，Ef和kf對E的影響趨勢線變化很明顯，且Ef和kf與E均呈良好的線性關系，而E隨其余參數的變化趨勢線基本水平。因此，用E的主要影響因素Ef和kf對其進行簡單線性擬合的效果就很好，見式(7)，對應表4中第1個公式。

表4 宏細觀參數間的擬合公式

5 細觀參數標定流程及實例驗證

5.1 FJM細觀參數標定流程

圖5 標定流程

5.2 實例驗證

劉翌晨[21]從齊熱哈塔爾引水隧洞中鉆取片麻花崗巖試樣進行物理試驗來獲取其宏觀力學參數，試驗結果平均值列入表5，本文以此為依據，進行片麻花崗巖的細觀參數標定。

表5 片麻花崗巖試樣宏觀參數試驗值和模擬值

表6 片麻花崗巖細觀參數標定結果

數值模擬及物理試驗得到的單軸壓縮應力-應變曲線如圖6所示，其中編號為Y-1,Y-2,Y-3和Y-4的曲線分別是物理試驗中4個試樣對應的結果。可以看出，模擬得到的曲線表現的特征與物理試驗得到的曲線表現的特征接近，但數值模擬得到的曲線沒有體現壓密階段特征(應力-應變曲線輕微下凸)，這是由于在建模過程中為得到均勻的顆粒體系進行了伺服，壓密了顆粒體系，因此施加上應力即出現彈性變形特征(應力-應變曲線近似直線)。

圖6 單軸壓縮應力-應變曲線

數值試驗和室內物理試驗試樣破壞結果如圖7所示，可見模擬結果與物理試驗結果吻合較好。單軸壓縮情況下出現貫通試樣上下端的主裂紋，由巖石內部拉伸破壞造成，體現了片麻花崗巖的劈裂張拉破壞特征；三軸壓縮下試樣出現了剪切破裂帶。

圖7 試樣破壞形態

6 結 論

(1) 正交數值試驗結果顯示采用FJM可以模擬出大壓拉比。對試驗結果進行多因素方差分析得到了各細觀參數對宏觀參數的影響程度，確定了影響各宏觀參數的主要細觀參數。

(3) 提出FJM細觀參數標定流程，然后以齊熱哈塔爾片麻花崗巖物理試驗結果平均值為依據進行細觀參數標定，用標定好的模型進行數值試驗，結果顯示模擬得到的宏觀參數值、應力-應變曲線及破壞形式接近物理試驗結果，驗證了本文PFC2D平節理模型細觀參數標定方法的可行性。