李汪洋， 石 崇， 張一平， 陳聞瀟， 張 坤

（河海大學巖土力學與堤壩工程教育部重點實驗室/巖土工程科學研究所，南京 210098）

近年來，基于離散單元法研究巖石與巖體非連續(xù)性質成為巖土類數(shù)值分析的熱點方向. PFC3D是當前應用最廣的顆粒流離散元程序工具，其中線性平行黏結模型經常用來模擬巖石性質，在巖石損傷［1］、圍巖穩(wěn)定性［2］、節(jié)理巖體破壞機制［3-5］等方面已經諸多應用，但該模型在建立宏觀參數(shù)與細觀參數(shù)聯(lián)系時非常困難，導致宏觀試驗與數(shù)值試驗曲線難以吻合.

國內外學者針對土石混合體和巖石等材料，已經嘗試了多種方法來確定細觀參數(shù)和其宏觀力學性質間定性或定量關系. Nardin 等［6］定性地研究了黏性材料細觀參數(shù)及其宏觀特性之間的關系. Potyondy 和Cundall［7］系統(tǒng)總結了平行黏結模型在Lac du Bonnet花崗巖隧道工程模擬中的應用，并給出了模型細觀參數(shù)的選取方法. 劉良軍等［8］基于研究線性平行黏結接觸模型中的細觀參數(shù)與單軸壓縮試驗應力-應變曲線及試樣破壞形式的對應關系，將細觀參數(shù)分為4類. 劉富有等［9］選擇Flatjoint接觸模型作為巖石模擬的基本模型，用單軸壓縮和巴西劈裂數(shù)值試驗測試巖石宏觀參數(shù)，對數(shù)值試驗進行正交設計，采用多因素方差分析和回歸分析研究宏觀和細觀參數(shù)之間的關系. 劉相如等［10］基于平直節(jié)理接觸模型，研究了接觸面單元數(shù)N、平直節(jié)理模量Ec、平直節(jié)理剛度比kr、摩擦系數(shù)μ、平直節(jié)理抗拉強度σc、平直節(jié)理黏聚力c 和平直節(jié)理內摩擦角θ 7種細觀參數(shù)對巖石宏觀力學參數(shù)的影響.

通常采用PFC方法模擬巖石等高黏結材料的單軸抗拉強度都相對較高，導致單軸抗壓強度和單軸抗拉強度比值（UCS UTS）為3～4. 而實際很多巖石的UCS UTS 常常超過10［11］，甚至達到20～30. 目前，多采用抑制顆粒旋轉的簇平行黏結模型［12］來提高UCS UTS 值，但單軸抗壓強度受簇中顆粒數(shù)影響較大. 本研究改進了常規(guī)平行黏結模型，建立了由13個參數(shù)控制的數(shù)值模型開展正交試驗，分析了巖石宏細觀參數(shù)間的定性與定量關系，建立了一種快速達到壓拉比等宏觀參數(shù)的標定方法.

1 線性平行黏結細觀模型

1.1 常用的線性平行黏結模型

在采用數(shù)值方法研究巖石工程的破壞機理，尤其是模擬巖石材料時，線性平行黏結模型［7］使用最為廣泛.

圖1 線性平行黏結模型組成結構Fig.1 Microstructure of linear parallel bond model

線性平行黏結模型更新接觸力和接觸力矩的力-位移法則如下：

1.2 考慮細觀特征分布的顆粒體系改進模型

基于線性平行黏結模型所計算得出的單軸抗拉強度往往較高，導致單軸抗壓強度和單軸抗拉強度之比UCS UTS 為3～4［14］，無法滿足很多巖石實際需求（UCS UTS 常常超過10）的客觀限制. 出于提升壓拉比的需求［15］，對線性平行黏結模型進行改進，將模型產生的ball-ball接觸以不同填充基質配比隨機分為強接觸和弱接觸兩組［16］（圖2）. 不同接觸組的參數(shù)設置見表1. 表1中RE表示弱接觸組黏結有效模量與強接觸組值的比值，Rk表示黏結法切向剛度比與線性法切向剛度比的比值，Rσ表示弱接觸組抗拉強度和黏聚力與強接觸組值的比值.

圖2 模型強弱接觸隨機分布Fig.2 Random distribution of strong and weak contacts in the model

表1 模型不同接觸組參數(shù)Tab.1 Parameters of different contact groups in the model

自然界中的巖石作為一種破壞程度不同的物質，Weibull 分布已被用來描述巖土材料的累積損傷［17-18］.它有兩個優(yōu)點：數(shù)值都為正，并且在距離尺度參數(shù)幾個標準差范圍之外的數(shù)值的概率接近于0，特別是在形狀參數(shù)較大的情況下. 巖土材料的強度參數(shù)均為正，同一礦物在巖土材料中的強度參數(shù)略有差異. 本節(jié)采用Weibull分布來表示礦物間內部膠結強度的隨機分布，威布爾分布的概率密度函數(shù)如下：

式中：α 和β 分別為控制概率密度函數(shù)均值和峰值的尺度和形狀參數(shù). 威布爾分布的概率分布函數(shù)見下式：

蒙特卡洛方法是一種以概率統(tǒng)計理論為指導的數(shù)值計算方法. 通過蒙特卡洛方法構造（0，1）區(qū)間上均勻分布的隨機數(shù)f，然后生成一個服從威布爾分布的隨機變量X. X的表達式如下：

因此改進模型由13個參數(shù)構成：黏結有效模量Eˉ*、線性有效模量/黏結有效模量E*Eˉ*、黏結法切向剛度比kˉ*、抗拉強度σˉc、黏聚力/抗拉強度cˉσˉc、摩擦系數(shù)μ、摩擦角φˉ、力矩分配系數(shù)βˉ、填料比例Rf、Weibull分布參數(shù)β、強度比例Rσ、模量比例RE和剛度比比例Rk. 其中Rf表示ball-ball 接觸分組中弱接觸組接觸數(shù)與總接觸數(shù)的比值，抗拉強度、黏聚力和黏結有效模量按表1賦值后再乘以Weibull分布的隨機變量X，以實現(xiàn)離散效果，Weibull分布中參數(shù)α 設為1.0. 圖3為不同β 取值下X的概率密度函數(shù)，可見隨著β 的增大，X取值越來越集中于1附近.

圖3 不同參數(shù)值下變量X 的概率密度函數(shù)Fig.3 Probability density function of variable X for different parameter values

2 正交試驗數(shù)值模擬

2.1 正交試驗原理

正交試驗設計是建立在概率論、數(shù)量統(tǒng)計和實踐經驗的基礎上，尋求因素水平較優(yōu)組合的一種簡單、高效、快速的試驗設計方法［13］. 在正交試驗設計中，各因素的水平根據(jù)先驗知識進行確定，然后進行全部因素水平的組合，再從所有試驗點中依據(jù)正交性選擇部分有代表性試驗點來實施試驗，選岀的這部分試驗點構成一張正交表.

2.2 宏觀和細觀試驗參數(shù)

根據(jù)室內常規(guī)三軸試驗曲線，選取彈性模量E、泊松比ν、單軸抗壓強度UCS、單軸抗拉強度UTS、內摩擦角φ、黏聚力c、裂紋損傷應力σcd作為試驗的宏觀參數(shù). 考慮到對巖石高壓拉比的要求，采用宏觀參數(shù)單軸抗壓與抗拉強度之比UCS UTS 來代替單軸抗拉強度UTS .

巖石體積應變壓縮膨脹的拐點被稱為裂紋損傷應力σcd，當應力超過σcd后，即使應力保持不變，也會有越來越多的微裂紋貫通導致巖石破壞，因此裂紋損傷應力σcd可以作為巖石的長期強度［19-20］. σcd對于單軸抗壓強度UCS 來說是一個相對值，本文定義裂紋損傷應力水平指標σcdUCS 作為試驗的宏觀參數(shù).

宏觀和細觀參數(shù)選取如表2. 如果能建立細觀參數(shù)與宏觀參數(shù)間的內在規(guī)律，無疑可以大大簡化參數(shù)標定過程消耗的時間.

2.3 PFC3D數(shù)值試驗

采用PFC3D數(shù)值方法進行壓縮和拉伸數(shù)值研究，三維數(shù)值模型如圖4. 試樣為圓柱體，半徑25 mm，高度100 mm，四周和上下施加墻體進行固定. 試樣顆粒共23 719個，最大粒徑1.5 mm，最小粒徑0.5 mm，孔隙率0.35.

表2 線性平行黏結模型宏細觀參數(shù)Tab.2 Macroscopic and microscopic parameters of linear parallel bond model

圖4 三維數(shù)值模型Fig.4 Three-dimensional numerical model

對模型進行單軸壓縮和單軸拉伸數(shù)值試驗，應力-應變曲線如圖5a和b所示. 圖5c和d為單軸壓縮的體積應變和泊松比與軸向應變的關系曲線. 根據(jù)結果曲線確定宏觀參數(shù)單軸抗壓強度UCS 和單軸抗拉強度UTS，彈性模量E 取軸向應變0.05%和0.15%間的割線斜率，見式（14）；泊松比ν 取軸向應變0.1%和0.2%對應值的平均，見式（15）；裂紋損傷應力σcd取體積壓縮膨脹轉化點對應的軸向應力. 最后進行三軸壓縮試驗，根據(jù)莫爾圓確定黏聚力c 和內摩擦角φ .

圖5 宏觀參數(shù)確定Fig.5 Determination of macroscopic parameters

2.4 正交試驗設計

選取表2中的13個細觀參數(shù)作為正交試驗研究因素，每個因素選取3個水平，各因素的因素水平見表3，7個宏觀參數(shù)作為參數(shù)標定的指標. 選擇正交表L27313安排數(shù)值試驗（表4）.

表3 因素水平表Tab.3 Factor level table

表4 正交數(shù)值試驗方案Tab.4 Orthogonal numerical test scheme

3 正交試驗結果分析

按照表4進行PFC3D單軸壓縮、單軸拉伸及三軸壓縮（圍壓1 MPa）的數(shù)值試驗，得到試驗結果（表5）.宏觀參數(shù)中壓拉比UCS UTS 范圍為3～19，可以很好地滿足實際巖石的要求.

表5 正交數(shù)值試驗結果Tab.5 Orthogonal numerical test results

3.1 各因素相關性分析

采用皮爾遜積矩相關系數(shù)計算各細觀參數(shù)和宏觀參數(shù)之間的相關性，如式（16），計算結果見表6.

表6 皮爾遜積矩相關系數(shù)Tab.6 Pearson product-moment correlation coefficient

3.2 多元回歸分析

根據(jù)圖6可以得到各宏觀參數(shù)的主要影響因素，分析如下：

根據(jù)上述分析，用巖石宏觀參數(shù)的主要細觀影響因素對其進行多元線性回歸分析，擬合公式見式（17）～（23）.

圖6 各宏觀參數(shù)對應相關系數(shù)Fig.6 Correlation coefficients corresponding to each macroscopic parameter

從上述擬合公式可以看出，除式（20）外，其他公式的R2大于0.7，擬合效果較好，說明宏細觀參數(shù)間的線性關系比較明顯；式（20）的R2較小，擬合效果一般，說明宏細觀參數(shù)間關系復雜，不能簡單地用線性關系來描述，需進一步分析.

圖7為宏觀參數(shù)φ 與對應各細觀參數(shù)的試驗平均值的關系曲線，散點表示試驗值，點線表示各組宏觀參數(shù)試驗結果的平均值，橢圓表示置信度為95的置信區(qū)間.

從圖7可以看出，μ 和Rf與內摩擦角φ 呈現(xiàn)出一定的對數(shù)關系，對其進行非線性曲線擬合，再結合其他變量進行擬合后的結果如式（24）.

綜上，經過修正后的宏細觀參數(shù)間的擬合公式見式（25）～（31），E 單位為GPa，UCS、c 單位為MPa，φ 單位為（°）：

圖7 宏觀參數(shù)φ 與各細觀參數(shù)試驗平均值關系曲線Fig.7 Curves between macroscopic parameter φ and mean values of test results of each microscopic parameter

4 細觀參數(shù)快速標定實例

依據(jù)文獻［21］對灰砂巖常規(guī)三軸壓縮試驗結果得到的宏觀力學參數(shù)進行模型細觀參數(shù)標定，UCS UTS值取為12.

基本標定步驟如下：

3）進行PFC初步數(shù)值模擬，通過單軸壓縮、單軸拉伸以及三軸壓縮試驗得到宏觀參數(shù). 根據(jù)初步模擬值和試驗值的差值以及宏細觀參數(shù)間的趨勢關系調整細觀參數(shù)，減小誤差，最終達到理想效果. 細觀參數(shù)最終取值見表7.

表7 灰砂巖模型細觀參數(shù)最終取值Tab.7 Calibration values of microscopic parameters of grey sandstone model

數(shù)值模擬與室內物理試驗得到的宏觀參數(shù)對比如表8，誤差除σcdUCS 外均在7%以內，表明快速標定方法是可行的.

圖8為單軸壓縮情況下，數(shù)值模擬與室內物理試驗得到的應力-應變曲線，模擬結果與室內試驗結果基本相符.

圖8 單軸壓縮應力應變曲線Fig.8 Stress-strain curves under uniaxial compression

5 結論

以改進的線性平行黏結模型作為顆粒接觸本構模型，對PFC3D模型進行正交數(shù)值試驗，研究了線性平行黏結模型13個細觀參數(shù)與巖石7個宏觀參數(shù)間的關系，并探討了PFC3D線性平行黏結模型細觀參數(shù)的標定方法，得出以下主要結論：

1）通過對線性平行黏結模型的改進措施，可以有效提高數(shù)值模擬結果的壓拉比UCS UTS，滿足實際工程中巖石的需要.

2）關聯(lián)性研究表明，宏觀參數(shù)與模型細觀參數(shù)關系較復雜，個別參量無主導因素. 此時可先假定一部分參數(shù)取值，如E*Eˉ*、kˉ*、βˉ、Rf等，再利用標定出的規(guī)律確定主要細觀參數(shù)，可以大大簡化標定過程.

3）利用所得出宏觀、細觀參數(shù)之間的關系，依據(jù)灰砂巖的室內物理試驗結果對線性平行黏結模型的細觀參數(shù)進行顆粒流標定，模擬結果與室內試驗結果基本相符，表明可用該方法標定巖石力學參數(shù).