巖石PFC2D模型細觀參數(shù)確定方法研究

郝保欽，張昌鎖，王晨龍 ，任 超

(1. 太原理工大學 礦業(yè)工程學院，山西 太原 030024；2.太原理工大學 機械與運載工程學院，山西 太原 030024)

0 引 言

常規(guī)巖石室內試驗是研究巖石力學性質的主要手段。受限于試驗設備和測試手段，巖石力學試驗一般只能得到巖石試件的宏觀變形和強度特性，而無法真正揭示巖石細觀層面的力學行為。顆粒流程序(PFC)基于離散單元方法，將巖石類材料視為剛性顆粒的集合體，克服了傳統(tǒng)連續(xù)介質分析方法在研究巖石類非均質材料時的局限性，可以從細觀角度再現(xiàn)裂紋的孕育、擴展以及貫通的時空演化規(guī)律，已成為研究巖石力學問題的有效手段之一[1]。

二維顆粒流程序(PFC2D)通過定義顆粒接觸的本構關系再現(xiàn)巖石宏觀力學性質，接觸黏結模型(CBM)、平行黏結模型(PBM)和平直節(jié)理模型(FJM)是最常用的3種細觀本構模型。不論選取哪種模型，都需要不斷調整細觀參數(shù)來匹配宏觀力學性質。目前大多數(shù)研究普遍采用“試錯法”進行標定[2-3]，調整過程盲目性大、費時費力。國內外學者針對此問題展開了大量研究，進行了一定的改進。劉新榮等[4]以接觸黏結模型中的等效微梁模型為基礎，經(jīng)過詳細的理論推導確定了細觀強度參數(shù)與巖石斷裂韌度之間的關系。由于接觸黏結模型顆粒接觸界面不能抵抗轉矩，平行黏結模型更廣泛地被用于模擬巖石材料。文獻[5-8]將理論分析與數(shù)值模擬相結合，系統(tǒng)地分析了平行黏結模型細觀參數(shù)對宏觀特性的影響，并給出了兩者之間的量化關系。POTYONDY等[1,9-11]研究表明，受顆粒形狀的影響，平行黏結模型不能很好地反映巖石的拉壓強度比和強度包絡線。為此，CHO等[10]分別提出了簇平行黏結模型[10]和平直節(jié)理模型[11]，通過抑制顆粒破壞后的旋轉，使得拉壓強度比顯著減小。鄧樹新等[12]采用PB(Plackett-Burman)試驗設計分析平直節(jié)理模型細觀參數(shù)對宏觀參數(shù)指標的敏感性，PB設計針對每個因素只能選取2個水平，當因素范圍較大時可能會影響分析結果。

上述研究在顆粒流細觀參數(shù)與巖石宏觀力學響應關系的建立方面取得了富有成效的進展，同時也存在一定的不足。以往研究大多集中在單軸壓縮、單軸拉伸和巴西劈裂試驗上，而對于有圍壓的三軸壓縮試驗研究較少。對平行黏結這一最常用模型進行宏細觀參數(shù)研究時，相關學者[1,9-11]做了一定的假設。對于摩擦因數(shù)和黏結摩擦角這些可能對宏觀力學響應有重要影響的細觀參數(shù)并未做深入的研究。鑒于此，針對PFC2D中的平行黏結模型，基于正交設計原理，通過一系列的實驗室試驗(單軸壓縮、單軸拉伸和三軸壓縮)系統(tǒng)全面地分析細觀參數(shù)與宏觀參數(shù)間的關系。

1 正交試驗設計

1.1 正交試驗設計原理

在科學試驗中常會遇到試驗結果受多種因素影響的情況。若采取全面試驗的方法，其工作量將隨因素的個數(shù)呈指數(shù)形式劇增。正交試驗設計是一種優(yōu)化的試驗方法，它利用正交表科學安排因素組合，能夠在明顯減少試驗次數(shù)的前提下得到可靠性較強的試驗結果[13-14]。通過多因素方差分析可以評估眾多因素對試驗指標的影響程度。因此，采用正交試驗設計法對PFC2D平行黏結模型宏細觀參數(shù)關系進行研究，量化參數(shù)對應關系，為參數(shù)標定提供依據(jù)。

正交表是正交試驗設計中安排試驗和分析結果的基本工具，一般用Ln(rm)表示。其中，L為正交表代號；n為試驗次數(shù)；r為因素水平；m為因素數(shù)量。正交表具有2個重要性質：表中任一列不同因素水平出現(xiàn)的次數(shù)均相同；表中任意兩列同行水平組成的數(shù)字對出現(xiàn)的次數(shù)相同。

1.2 方差分析

方差分析可以精確量化試驗因素對試驗結果的影響程度，通過計算因素的離差平方和與誤差的離差平方和，并基于此構造F統(tǒng)計量。對于給定的顯著性水平α，可以查F分布表得到對應的臨界值F1-α。當F值大于臨界值F1-α時，認為該因素對試驗結果有顯著影響。

若采用Ln(rm)正交表設計試驗，y1，y2，…，yn為試驗結果，方差分析的基本步驟如下。

1)計算總離差平方和SST。

(1)

其中，

(2)

(3)

2)計算因素j引起的離差平方和SSj。

(4)

其中，

(5)

(6)

式中：Ki為第j列水平號為i的試驗結果之和。

3)計算試驗誤差的離差平方和SSe。

SSe=∑SSk

(7)

式中，SSk為空列離差平方和。

4)計算自由度。SST和SSj的自由度分別為

fT=n-1

(8)

fj=r-1

(9)

誤差的自由度為

fe=∑fk

(10)

式中：fk為空列自由度。

5)計算F值。

(11)

6)顯著性檢驗。對于給定的顯著性水平α，當Fj>F1-α(fj,fe)時，則該因素對試驗結果有顯著性影響。Fj相對F1-α(fj,fe)越大，則該因素對試驗結果的影響越顯著。

2 試驗指標及試驗因素

2.1 試驗指標選取

巖石力學試驗是測試巖石在不同受力條件下(即單軸壓縮、單軸拉伸和三軸壓縮等)的破壞及變形特性參數(shù)，通過巖石力學試驗可以得到巖石材料的力學參數(shù)包括抗壓強度、抗拉強度、彈性模量以及泊松比等。

目前關于參數(shù)標定的大多數(shù)研究都是針對單軸壓縮試驗展開的，通過不斷調整細觀參數(shù)對應力-應變曲線進行擬合。陳鵬宇等[15]研究指出，僅通過單軸壓縮試驗標定得到的顆粒流數(shù)值模型，并不能直接用于研究有圍壓條件下巖石的力學性質。自然界的絕大多數(shù)巖石都處在三向應力狀態(tài)，研究巖石在三向應力狀態(tài)下的強度和變形特征對于指導巖土工程設計具有重要的意義。強度包線可以全面反映三軸試驗的強度指標值，因此選擇它作為參數(shù)標定的依據(jù)。強度包絡線包括黏聚力和內摩擦角2個宏觀力學參數(shù)，對于庫倫強度包線，黏聚力可以由單軸抗壓強度和內摩擦角表示，因此僅選取內摩擦角作為試驗指標。

巖石的抗拉強度是反映巖石力學性質的重要指標之一。巖石的抗拉強度遠小于其抗壓強度，在受荷不大時就可能出現(xiàn)拉伸破壞，因此進行數(shù)值仿真時應將拉伸強度作為重要的標定依據(jù)。POTYONDY等[1,9-10]采用平行黏結模型模擬Lac du Bonnet花崗巖材料時發(fā)現(xiàn)，其抗拉強度與抗壓強度的比值顯著大于實驗室的測量值。但以上的研究對一些參數(shù)如摩擦因數(shù)μ、黏結摩擦角φ以及黏結強度變異系數(shù)Cv等的選取并不充分，它們對于拉壓強度比的影響還不明確。

啟裂強度是裂紋萌生和穩(wěn)定擴展的標志，它對估算巖體隧洞松動圈范圍、指導圍巖支護有重要的意義。室內試驗常通過計算裂紋體積應變[16]，或結合聲發(fā)射試驗結果確定巖石試件的啟裂強度[17]。采用PFC模擬巖石試件可以記錄試驗加載過程中的裂紋數(shù)目及擴展情況，可以結合裂紋數(shù)目判斷試件的啟裂強度。POTYONDY等[1]通過裂紋數(shù)目的變化來確定啟裂強度，認為裂紋數(shù)目達到峰值強度裂紋數(shù)目1%時所對應的強度即為啟裂強度。數(shù)值模型相對半徑范圍較大，為兼顧所有數(shù)值模型，選取10%為啟裂強度的判斷條件。

2.2 試驗因素確定

采用平行黏結模型(PBM)模擬巖石室內試驗。平行黏結模型將巖石表示為具有有限厚度圓盤的組合，圓盤之間通過黏結鍵相連，圓盤之間的接觸以及黏結區(qū)域具有有限剛度。當施加的局部應力(即張力、剪切力或旋轉力矩)超過規(guī)定的鍵的強度時，這些鍵斷裂形成一個微裂紋。在這些過程中，PFC不需要任何描述屈服后響應流動規(guī)律來控制斷裂行為，只需要顆粒的運動規(guī)律、顆粒變形的簡單本構規(guī)律和顆粒鍵斷裂的屈服規(guī)律。平行黏結模型參數(shù)主要包括顆粒細觀參數(shù)和黏結細觀參數(shù)，具體見表1。

表1 PFC2D模型主要的宏細觀力學參數(shù)

2.3 試驗因素范圍

在顆粒離散元模型中，黏結破壞后顆粒間的滑移與摩擦因數(shù)μ有關，POTYONDY等[1]認為摩擦因數(shù)可能會影響材料的峰后行為，建議對于巖石材料μ取0.5。CHO、程愛平等[10,18]認為顆粒間的摩擦力會抑制顆粒的運動，從而提高材料的強度。為了研究μ對宏觀參數(shù)的影響，設定μ的取值范圍為0.2～0.8，這與文獻[10,12,15]的取值相近。

文獻[19]表明能夠體現(xiàn)材料統(tǒng)計平均性質代表性體積單元的尺寸至少應該是顆粒尺寸的7～8倍。周喻等[20]研究發(fā)現(xiàn)，當RES=(L/Rmin)[1/(1+Rmax/Rmin)]≥10，即L/Rmin≥26.6(Rmax/Rmin=1.66)時，計算結果離散性不大。尹小濤等[21]研究表明，材料內尺度比非常小時，它對模型特性幾乎無影響，但會造成計算時間急劇升高。兼顧計算精度和計算效率，設定L/Rmin取值范圍為30～120。

陳鵬宇等[15,23]利用單軸或三軸數(shù)值試驗展開了黏結內摩擦角φ與宏觀力學參數(shù)關系的研究，結合他們的研究成果，將φ的取值設定為0°～60°。

變異系數(shù)Cv是黏結強度標準差與強度均值的比值，反映了模型細觀強度的不均勻性。阿比爾的等[6]研究發(fā)現(xiàn)Cv<0.2時，模型強度參數(shù)趨于穩(wěn)定，Cv取值過大則會造成強度參數(shù)偏小，因此將Cv的取值設定為0.2～0.8。

3 細觀參數(shù)敏感性分析

3.1 正交試驗設計方案及計算結果

根據(jù)因素個數(shù)與因素水平，采用L32(49)正交表設計試驗。依據(jù)因素水平可設置32組參數(shù)組合?；趨?shù)組合分別進行單軸壓縮、直接拉伸和三軸壓縮試驗(圍壓分別為5、10、15、20 MPa)，共計192組數(shù)值試驗，具體結果見表3。數(shù)值試樣長度H=100 mm，寬度L=50 mm，密度ρ=2 700 kg/m3。

表2 各參數(shù)因素水平

對表3的正交試驗設計結果進行初步的分析，數(shù)值模型試件的拉壓強度比σt/σc是關注的重點。巖石的抗拉強度一般為抗壓強度的1/25～1/4，平均為1/10[24]。對表3中的結果進行統(tǒng)計可以得到拉壓強度比在0.09～0.46，這基本滿足了大部分巖石的標定需求。同時也表明采用平行黏結模型模擬巖石試件時，全面地對細觀參數(shù)進行合理取值就可以得到較小的拉壓強度比。

表3 基于正交試驗設計的數(shù)值計算方案及結果

3.2 正交設計結果分析

對數(shù)值計算結果進行方差分析見表4。根據(jù)式(5)和式(8)分別計算因素和誤差的離差平方和。根據(jù)公式(10)可計算得各因素自由度f=r-1=4-1=3，誤差自由度fe=3。計算選定顯著性水平α=0.01，F(xiàn)0.99(3,3)=29.46。當F≥F0.99(3,3)時，認為因素j對試驗結果有非常顯著影響，記作“**”。對于顯著性水平α=0.05，F(xiàn)0.95(3,3)=9.28。當Fj≥F0.95(3,3)時，認為因素j對試驗結果有顯著性影響，記作“*”。為更直觀地分析試驗因素對試驗指標的影響程度，將方差分析結果繪制成柱狀圖，如圖1所示。

表4 方差分析

圖1 多因素方差分析F統(tǒng)計量

結合表4和圖1的方差分析統(tǒng)計結果可以看出，試驗因素對試驗指標影響程度各不相同，具體分析如下：

3)圖1c顯示，對試塊泊松比ν有非常顯著影響的只有剛度比kn/ks這1個細觀參數(shù)。雖然Ec和L/Rmin達到了著影響水平，但他們對ν的影響可以忽略。

6)圖1f有且僅有一個試驗因素Cv對啟裂強度與單軸抗壓強度比值σci/σc有非常顯著的影響。通過多因素方差分析，可以從眾多的細觀參數(shù)中篩選出對試驗指標有顯著性影響的因素。對顯著性影響因素與試驗指標進行回歸分析，可以得到宏細觀參數(shù)間的對應關系：

(12)

ν=0.144ln(kn/ks)+0.080 (R2=0.70)

(13)

(14)

(15)

(16)

(17)

4 標定方法及應用

4.1 標定方法

圖2 細觀參數(shù)標定流程

4.2 實例驗證

以錦屏二級水電站深埋白山組大理巖[25-26]為研究對象，采用上述方法進行細觀參數(shù)標定。對比大理巖力學試驗及聲發(fā)射試驗結果，驗證標定方法的可行性。表5再現(xiàn)了細觀參數(shù)標定過程，共進行了5次調整過程。式(12)—式(17)反算可以得到步驟1的初始細觀參數(shù)，此時對應的宏觀參數(shù)與目標值接近，但仍有一定的差距。這是因為細觀參數(shù)數(shù)目眾多，回歸分析時只考慮了對宏觀參數(shù)有顯著性影響的細觀參數(shù)，而影響不顯著的參數(shù)可能對試驗結果有一定的影響。 步驟2～6參考圖2的標定流程，對細觀參數(shù)進行了適當?shù)恼{整。將試驗平均值與模擬值匯總在表6中，除內摩擦角以外，其他模擬結果與目標值之間的誤差均在7%以內，可認為標定結果符合要求，結束標定。

表5 大理巖宏觀參數(shù)試驗值與模擬值

表6 大理巖細觀參數(shù)標定過程

將數(shù)值模擬試驗和物理試驗獲得的單軸壓縮應力-應變曲線繪制在圖3中。可以看出兩者曲線基本重合，但存在一定的區(qū)別。大理巖試件中存在的微裂隙會在壓力作用下產生閉合，應力-應變曲線會出現(xiàn)明顯的壓密階段。而數(shù)值模擬在生成顆粒時采用了一定伺服壓力，使產生的顆粒集合體接觸致密、均勻，因此其應力應變關系峰前階段近似為一條直線。

圖3 物理試驗和數(shù)值試驗的單軸壓縮應力-應變曲線對比[26]

莫爾-庫侖強度準則因其簡單實用而在巖土工程領域得到廣泛應用，可以用來表征強度包絡線。該準則在主應力平面上呈線性關系，即

(18)

因此，若已知σ1與σ3間的關系就可以求得黏聚力C和內摩擦角θ。

(19)

(20)

式中：k為極限主應力(σ1-σ3,σ1為最大主應力，σ3為最小主應力)關系曲線斜率[24]；σc為極限主應力關系曲線在縱坐標軸上的截距。

圖4為數(shù)值試驗與室內試驗極限主應力關系曲線對比，試驗曲線參考文獻[26]中σ1-σ3擬合曲線。根據(jù)擬合曲線，由式(18)—式(20)可計算得到巖石的內聚力和摩擦角，將計算結果標注在圖4中?？梢钥闯?，當圍壓較小時，數(shù)值模型預測的極限軸向應力與試驗結果較接近，而隨著圍壓的增大，他們之間的差值也逐漸增大。這可能是由于PFC基于離散單元方法，將巖石類材料視為圓形/球形顆粒的集合體，而實際巖石礦物顆粒形狀則為不規(guī)則的多面體。形狀的不同可能會導致顆粒間的咬合作用產生差異，而圍壓的增加則放大了這種差異。POTYONDY 等[1,12]研究也指出顆粒離散元法在極限主應力關系曲線的擬合上存在一定困難。

圖4 極限主應力關系曲線[26]

將數(shù)值試驗得到的試件破壞形式實際物理試驗進行對比，如圖5所示。對比發(fā)現(xiàn)不同加載條件下的數(shù)值試驗與物理試驗結果吻合性較好。圖5a單軸壓縮試驗試件的破壞形式表現(xiàn)為沿剪切帶的整體滑移破壞，同時與主要滑移面斜交方向產生1個次生破壞面。圖5b圍壓18 MPa下的三軸壓縮試驗試塊均出現(xiàn)了宏觀剪切破壞帶(其他圍壓條件下試件破壞形式也一致)。

圖5 破壞形式對比[26]

5 結 論

1)通過對正交試驗結果進行多因素方差分析，精確估計了各細觀參數(shù)對宏觀參數(shù)的影響程度。巖石單軸抗壓強度主要受黏結切向強度和黏結強度比的影響；巖石泊松比主要受剛度比的影響，它們之間呈對數(shù)關系，泊松比隨著剛度比增大而增大；對巖石的宏觀彈性模量有顯著影響的細觀參數(shù)，影響程度由高到低為：細觀彈性模量、相對半徑及黏結強度比；巖石拉壓強度比主要受黏結內摩擦角和黏結強度變異系數(shù)的影響，合理地調整細觀參數(shù)可以使拉壓強度比降至0.1左右，滿足大部分巖石的要求；對巖石內摩擦角有顯著影響的細觀參數(shù)類型最多，影響程度由高到低依次為：黏結內摩擦角、摩擦因數(shù)、黏結切向強度以及細觀彈性模量。

2)根據(jù)方差分析結果，綜合考慮多個細觀參數(shù)對宏觀參數(shù)的影響，進一步建立了顯著影響因素與試驗指標間的多元函數(shù)關系。

3)根據(jù)宏細觀參數(shù)間的關系，確定了細觀參數(shù)的標定順序，提出了具體的標定方法。并采用該方法對錦屏白山組大理巖力學參數(shù)進行了標定。數(shù)值試驗與物理試驗的宏觀力學參數(shù)、破壞模式以及應力-應變曲線基本對應，驗證了細觀參數(shù)標定方法的可行性。

4)由于采用圓形顆粒，平行黏結模型對極限主應力關系曲線的擬合效果較差。