基于統計試驗設計下花崗巖兩種本構模型研究

溫燕華，袁周祥，薛少欣

(河北工程大學 水利水電學院，河北 邯鄲 056038)

0 引 言

巖石材料的細觀本構特征，決定了巖石材料的宏觀力學行為特征?；陬w粒流理論的PFC(Particle Flow Code，PFC)軟件可以將細觀力學參數與宏觀力學參數聯系起來，用以模擬巖石的基本力學特性。該方法是P. A. Cundall等在離散元理論基礎上提出的，將顆粒的運動以及顆粒間應力的交互作用作為模擬巖石的基本條件[1]。目前，顆粒流程序中用于表征巖石材料細觀本構的兩種主流模型分別為平行黏結模型(Parallel-Bonded，PB)和平直節理模型(Flat-joint，FJ)。平行黏結模型由POTYONDY等[1-2]創建，針對該模型在模擬巖石力學特性的理論方法以及優勢方面進行過詳細闡述；而蔣明鏡等[3]在研究中發現，平行黏結模型雖然可以較為準確的模擬巖石的脆性破壞特征，但是計算得出的單軸抗拉強度較高，與實際結果存在差異。后續，POTYONDY[4]建立Flat-Joint接觸模型，通過構造顆粒形狀及新的接觸本構關系，使顆粒破壞后的旋轉遭到抑制，顯著提升了材料的壓拉比(USC/TS)，但未指出Flat-Joint接觸模型的適用范圍以及缺陷。劉富有等[5]、陳鵬宇等[6]則針對Flat-Joint模型，研究宏細觀參數間的對應關系，并建立了相應的標定方法，但并未對模型的計算穩定性、破壞特性等方面進行闡述。

盡管上述研究貢獻諸多成果，但仍有待于進行補充與完善。眾所周知，基于顆粒流方法生成巖石材料時，顆粒半徑大小是按均勻分布生成，顆粒間的黏結強度是按正態分布生成。這兩種隨機因素均會對模型的宏觀力學參數與破壞特征存在顯著影響，即影響計算結果的穩定性?；诖?，本文以新疆奇熱哈塔爾過水隧洞片麻狀花崗巖力學參數為基礎，建立二維離散元模型，通過統計學的試驗方法，對比分析PB模型與FJ模型在宏觀力學參數及破壞模式的差異性與穩定性，為后續改進本構模型以及選取合適的本構模型進行巖石數值模擬提供理論依據。

1 研究方法

1.1 物理實驗

以新疆齊熱哈塔爾過水隧洞的片麻狀花崗巖為研究對象，采用TAW-2000微機控制電液伺服巖石三軸試驗機對巖石試塊進行單軸壓縮試驗，加載速率0.002 mm/s，得到巖石的物理力學參數以及破壞模式。見圖1。

圖1 試驗示意圖以及巖石破壞模式Fig.1 Test schematic diagram and the rock failure mode

宏觀物理力學參數見表1，以單軸抗壓強度、彈性模量、泊松比、抗拉強度作為細觀參數反演的主要目標值。

表1 巖石宏觀力學參數Tab.1 Macroscopic mechanical parameters of rock

1.2 模型的建立及計算

本試驗是在顆粒離散元PFC2D軟件平臺上完成。單軸抗壓模型尺寸為長100 mm、寬50 mm。顆粒平均半徑為0.3 mm，最大半徑與最小半徑比值為1.66。通過半徑擴大法來調整模型內的顆粒分布，顆粒半徑在最大半徑與最小半徑之間服從均勻隨機分布。巴西劈裂數值模型半徑為25 mm，顆粒半徑設置與單軸壓縮相同。利用內嵌Fish語言編寫伺服算法控制加載墻體運動，加載速度設置0.02 mm/s，以此實現單軸以及巴西劈裂的數值模擬試驗。

1.3 參數的選取

巖石宏觀力學參數與細觀力學參數之間存在高度非線性關系，兩者之間的函數關系無法用精確的數學表達式表示。因此需通過數值力學試驗，進行宏細觀匹配演算。本文針對上述兩種細觀本構，開展單軸壓縮試驗和巴西劈裂試驗，單軸試驗獲取的宏觀力學參數包括彈性模量、抗壓強度、泊松比，巴西劈裂試驗獲取的宏觀力學參數為抗拉強度。當該細觀本構的數值模擬結果與目標值的相對誤差在5%以內，則以該組細觀參數為基準研究該細觀本構數值計算結果的穩定性。需要說明的是，由于PB模型無法準確計算出抗拉強度，因此在宏細觀匹配演算過程中不考慮該力學參數。根據上述方案，細觀參數的計算結果見表2，由表2中的細觀參數列表可知，FJ本構比PB本構多出一個接觸片段參數。

表2 兩種本構模型細觀參數Tab.2 Mesoscopic Model Mesoscopic Parameters

1.4 試驗方案

依據研究目的，針對PB本構與FJ本構各隨機生成10組用于單軸壓縮試驗的試樣以及10組用于巴西劈裂試驗的試樣。隨機因子分別設置為1001，2001，3001，4001，5001，6001，7001，8001，9001，10001。

由此可知，數值試驗結果數據為成組設計定量資料，區組因素為細觀本構。本文將采用描述性統計方法，針對兩種本構模型，計算分析其宏觀力學參數的穩定性以及破壞特征。采用定量統計檢驗的方法，針對兩種本構模型，檢驗其宏觀力學參數的統計學差異性。選用統計檢驗方法的原則如下：若各組數據資料滿足獨立性、正態性與方差齊性，則優先選用檢驗功效較高的參數檢驗——T檢驗，否則采用相應的非參檢驗——Wilcoxon秩和檢驗[7]。根據文獻[8-11]的建議，利用Shapiro-Wilk檢驗和Levene檢驗對資料的正態性和方差齊性進行檢驗。

2 試驗結果及分析

2.1 宏觀力學參數的穩定性分析

如試驗方案所述，通過力學試驗，獲取兩種細觀本構的宏觀力學參數，采用描述性的統計方法進行分析。圖2展示的是抗壓強度、彈性模量、泊松比以及壓拉比的數值試驗結果。由圖2中數據可知，兩種細觀本構的抗壓強度、彈性模量、泊松比的均值數據相差較小，并未體現出較大的數值差異；而壓拉比的均值則表現出較大的數值差異，即FJ模型的壓拉比均值要遠大于PB均值。同時，由各樣本數據分布特征、標準差(Standard Deviation，SD)以及變異系數(Coefficient of variation, CV)可看出，FJ模型的各力學參數的離散程度均大于PB模型。FJ模型除了彈性模量的變異系數小于5%以外，其它力學參數的變異系數均高于5%，以壓拉比最為明顯。PB模型除了壓拉比的變異系數略大于5%以外，其它力學參數的變異系數均低于5%。

以上研究結果說明，兩種本構模型對力學參數的均值影響很小，即均能體現出花崗巖的宏觀力學特性。二者主要區別在于力學計算結果的穩定性方面，PB模型的力學計算穩定性質明顯優于FJ模型；但是由于無法準確標定抗拉強度，在模擬巖石的真實壓拉比方面，FJ模型要優于PB模型。

圖2 花崗巖數值模型力學參數結果Fig.2 Results of mechanical parameters of granite numerical model

2.2 宏觀力學參數的統計檢驗分析

由表3可知，兩種細觀本構模型的抗壓強度、彈性模量、泊松比的樣本數據均滿足正態性與方差齊性，因此可以采用T檢驗來比較兩組樣本均值是否存在顯著的統計學差異。本文將顯著性檢驗水平定為0.05，假設檢驗為兩組樣本的總體均值相等。由統計檢驗結果可知，這3類宏觀力學參數的均值，在兩種本構模型之間并未存在顯著的統計學差異，即兩種細觀本構在描述同一種巖石材料時，在抗壓強度、彈性模量、泊松比上具有較高的一致性。但抗拉強度均值及壓拉比均值在兩種本構模型之間存在顯著的統計學差異，即FJ模型的壓拉比遠高于PB模型，前者更能客觀體現出巖石材料的抗拉強度。

表3 宏觀力學參數的統計分析結果Tab.3 Statistical analysis results of macroscopic mechanical parameters

2.3 破壞模式分析

圖3與圖4分別為FJ模型與PB模型在不同壓拉比下對應的破壞模式圖。由圖3可知，PB模型所表現的破壞模式以單斜面剪切滑移為主，同時在剪切帶附近形成多條次生裂紋帶。次生裂紋帶方向趨向于加載應力方向，這說明次生裂紋帶多為試件橫向膨脹所形成的張拉破壞。由圖4可知，FJ模型表現的破壞模式主要表現為兩種：①單斜面剪切破壞；②隨著壓拉比的增加，剪切帶與試件軸向夾角逐漸減小，呈現為軸向劈裂以及表層剝離的破壞模式，即這些破壞主要是由于橫向拉應力超過巖石抗拉極限所引起的。

圖3 不同壓拉比下，PB模型的主要破壞模式Fig.3 Main failure modes of PB model under different draw ratios

圖4 不同壓拉比下，FJ模型的主要破壞模式Fig.4 Main failure modes of FJ model under different draw ratios

將黏結鍵斷裂傾角按10°間隔進行分組，統計每組傾角范圍內的斷裂數占總裂紋數的比值，繪制相應的裂紋傾角玫瑰圖，見圖5。由圖5可看出，兩種細觀本構模型的破壞裂紋傾角集中在45°～145°之間，均表現出以斜剪破壞為主，但亦有拉伸破壞。傾角在90°附近的裂紋數，FJ模型要高于PB模型，由上述分析可知，這是由壓拉比較高的條件下所導致。

圖5 不同本構模型硬巖的裂紋傾角玫瑰圖Fig.5 Crack dip rose of hard rock with different constitutive model

由以上分析可知，PB模型在表現花崗巖宏觀破裂模式方面較為穩定，以斜剪破壞為主，但會伴生次生的張拉裂紋，這主要是由于PB模型壓拉比變化范圍較小所致。FJ模型在宏觀破裂模式上，會隨著壓拉比的變化，呈現斜剪向拉伸破壞模式的轉變，即更能表現出巖石材料的非均質性。

3 結 論

以新疆奇熱哈塔爾過水隧洞片麻狀花崗巖物理實驗力學參數為基礎，建立顆粒流花崗巖數值模型，在考慮顆粒半徑及細觀黏結強度隨機分布的條件下，通過單因素成組的試驗設計，對比分析FJ本構模型與PB本構模型在宏觀力學參數及破壞模式的差異性與穩定性，共得出以下結論：

1) 兩種本構模型均能做到較為準確描述片麻狀花崗巖材料宏觀力學參數，如彈性模量、抗壓強度、泊松比，但PB模型力學參數的變異系數均低于FJ模型，即從計算穩定性角度來講，PB模型計算結果更具有確定性。

2) PB模型由于自身理論的缺陷，相對于FJ模型未能表現出花崗巖巖石材料真實的壓拉比這一力學性質。因此，利用PB模型研究均質性較強且抗拉強度較高的硬巖材料的力學行為時，更容易發現普遍性規律，具有一定優勢；而對于均質性較弱且壓拉比較高的硬巖材料，FJ模型更為適合。

3) 兩種細觀本構模型的破壞裂紋傾角集中在45°～145°之間。PB模型在單軸壓縮狀態下，其花崗巖材料的宏觀破壞模式以斜剪破壞為主，在宏觀裂紋帶兩側會伴生次生的張拉裂紋，破壞模式相對穩定，這主要由壓拉比及PB模型計算穩定性所決定。FJ模型在宏觀破裂模式上，會隨著壓拉比的變化，呈現斜剪向拉伸破壞模式的轉變，張拉裂紋數要高于PB模型，其宏觀破壞模式較為不確定，更能表現出巖石材料的非均質性。

4) 針對以上研究結果，建議今后在進行花崗巖等硬性巖石材料的數值試驗時，可根據研究內容與巖石特性，選擇合適的本構模型，為后續研究巖石破壞機理提供參考價值。