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

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

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

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

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

1 線性平行黏結(jié)細觀模型

1.1 常用的線性平行黏結(jié)模型

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

圖1 線性平行黏結(jié)模型組成結(jié)構(gòu)Fig.1 Microstructure of linear parallel bond model

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

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

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

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

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

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

因此改進模型由13個參數(shù)構(gòu)成：黏結(jié)有效模量Eˉ*、線性有效模量/黏結(jié)有效模量E*Eˉ*、黏結(jié)法切向剛度比kˉ*、抗拉強度σˉc、黏聚力/抗拉強度cˉσˉc、摩擦系數(shù)μ、摩擦角φˉ、力矩分配系數(shù)βˉ、填料比例Rf、Weibull分布參數(shù)β、強度比例Rσ、模量比例RE和剛度比比例Rk. 其中Rf表示ball-ball 接觸分組中弱接觸組接觸數(shù)與總接觸數(shù)的比值，抗拉強度、黏聚力和黏結(jié)有效模量按表1賦值后再乘以Weibull分布的隨機變量X，以實現(xiàn)離散效果，Weibull分布中參數(shù)α 設(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è)計是建立在概率論、數(shù)量統(tǒng)計和實踐經(jīng)驗的基礎(chǔ)上，尋求因素水平較優(yōu)組合的一種簡單、高效、快速的試驗設(shè)計方法［13］. 在正交試驗設(shè)計中，各因素的水平根據(jù)先驗知識進行確定，然后進行全部因素水平的組合，再從所有試驗點中依據(jù)正交性選擇部分有代表性試驗點來實施試驗，選岀的這部分試驗點構(gòu)成一張正交表.

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

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

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

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

2.3 PFC3D數(shù)值試驗

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

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

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

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

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

2.4 正交試驗設(shè)計

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

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

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

3 正交試驗結(jié)果分析

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

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

3.1 各因素相關(guān)性分析

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

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

3.2 多元回歸分析

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

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

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

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

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

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

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

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

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

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

基本標(biāo)定步驟如下：

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

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

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

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

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

5 結(jié)論

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

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

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

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