宋 彧，冷國杰，楊安玉

（1.浙江國貿集團東方機電工程股份有限公司，浙江 杭州 310000；2.南京水利科學研究院，江蘇 南京 210024；3.水利部農村電氣化研究所，浙江 杭州 310012）

1 問題的提出

天然巖土材料在漫長的地質過程中，形成了異常復雜的結構，成分、大小不一的礦物顆粒以不同的方式結合在一起，聚集為團粒，團粒再堆積起來形成巖土材料，期間由于外部作用，形成節(jié)理裂隙，從而使巖土材料具有結構性。巖土材料的結構性帶來巖土問題的非線性、復雜性[1]。顆粒流法[2-6]是研究巖土體材料機理特征的一種有效的數值分析方法。

當前，常用的顆粒流模擬程序PFC 中，數值模型的宏觀基本物理特性無法通過直接賦值的形式實現，只能對構成模型顆粒的幾何特性及顆粒間接觸的細觀力學參數進行賦值[7-8]。且必須通過大量試算，標定出細觀參數，才能使數值模型表現出所需宏觀物理特征。

文章采用PFC2D 程序建立顆粒流接觸黏結模型進行單軸壓縮模擬，并通過Plackett-Burman 設計獲得模型各細觀參數的敏感性響應，篩選出影響各項宏觀力學特性的主要細觀參數因子，提高顆粒流模型細觀參數標定的便捷性。

2 Plackett-Burman 試驗設計

Plackett-Burman（P-B）試驗設計是由R.L.Plackett和J.P.Burman[9]于1946 年提出。是一種基于非完全平衡塊原理近飽和的2 水平試驗設計方法，通過N次實驗至多可以研究（N-1）個變量，其中每個變量取高、低2 個水平。P-B 試驗設計能夠確定顯著影響因子，從而在眾多的變量因素中快速有效地篩選出最為重要的幾個因素[10]。

P-B 試驗設計采用線性函數篩選得出重要因素，不考慮相互影響，其中，函數方程為：

式中：Y為響應函數值，Xi為設計函數值。β0和βi為回歸系數，通過回歸系數的顯著水平P值來反映Xi對Y的影響程度。

3 基于Plackett-Burman 試驗設計的模型計算

3.1 數值模型建立

研究數值模型宏觀物理特征對于細觀參數的敏感性響應，無需對應于某一真實的巖石試件。數值模型的宏觀物理特征參數為彈性模量、單軸抗壓強度及泊松比。綜合考慮眾多文獻資料，建立一顆粒流接觸黏結模型，模型初始參數取值見表1。

表1 數值模型的細觀參數表

生成的初始模型見圖1，模型尺寸18.58 mm×49.80 mm，共生成顆粒1 079 個。建立上下左右4 面墻（Wall），對模型進行邊界約束。移動上下墻對模型施加軸壓，通過對模型進行數值單軸壓縮試驗可獲取模型的宏觀物理特征參數。

圖1 顆粒流數值模型圖

3.2 單因子試驗確定因子水平

3.2.1 顆粒接觸模量

對顆粒接觸模量這一細觀參數進行單因子試驗，顆粒接觸模量分別取10，20，30，40，50，60，70 GPa，其余參數與初始模型參數相同。建立7 個模型分別進行數值模擬得到應力—應變曲線（見圖2）。

圖2 顆粒接觸模量單因子試驗應力—應變曲線圖

根據單因子試驗，顆粒接觸模量與模型彈性模量的關系見圖3。

圖3 顆粒接觸模量對模型彈性模量的影響圖

3.2.2 法向接觸強度/切向接觸強度

對法向接觸強度/切向接觸強度這一細觀參數進行單因子試驗，法向接觸強度/切向接觸強度分別取10，20，30，40，50，60 MPa，其余參數與初始模型參數相同。建立6 個模型分別進行數值模擬得到應力—應變曲線（見圖4）。

圖4 法向接觸強度/切向接觸強度單因子試驗應力—應變曲線圖

根據單因子試驗，法向接觸強度/切向接觸強度與模型單軸抗壓強度的關系見圖5。

圖5 法向/切向接觸強度對模型單軸抗壓強度的影響圖

3.2.3 顆粒剛度比

對顆粒剛度比這一細觀參數進行單因子試驗，顆粒剛度比分別取0.5，1.0，1.5，2.0，2.5，3.0，3.5，其余參數與初始模型參數相同。建立7 個模型分別進行數值模擬得到應力—應變曲線（見圖6）。

圖6 顆粒剛度比單因子試驗應力—應變曲線圖

根據單因子試驗，顆粒剛度比與模型泊松比的關系見圖7。

圖7 顆粒剛度比對模型泊松比的影響圖

3.2.4 因子水平確定

根據單因子試驗結果，接觸強度的-1 和+1水平分別設為20 MPa 和50 MPa，顆粒接觸模量的-1 和+1 水平分別設為20 GPa 和60 GPa，顆粒剛度比的-1 和+1 水平分別設為1 和3，其余參數水平根據常規(guī)設置選取。最終，P-B 試驗設計因子水平范圍見表2。

表2 Plackett-Burman 試驗設計因子水平表

3.3 Plackett-Burman 試驗設計與結果分析

Minitab 是一種具有實驗設計功能的統(tǒng)計軟件，是一種研究與處理多因素試驗的科學方法。Minitab 具有體積小、功能強大、操作簡單、易于掌握、計算精確、兼容性好等特點，在影響因素眾多且主因素不確定的情況下，通過篩選試驗設計，經過少量試驗即可找出主因素[11]。

本文采用Minitab 軟件生成篩選試驗設計表，根據P-B 試驗設計進行12 次計算分析，設計及試驗結果見表3。

表3 Plackett-Burman 試驗設計及響應值表

擬合一階模型為：

對模型進行方差分析（結果見表4）可知，所擬合的3 個回歸方程均達到顯著（P＜0.05），決定系數R2分別為0.985 0、0.991 0、0.916 0，表明分別用方程擬合7 個細觀因子與3 個宏觀力學特性之間的關系是可行的。

表4 回歸模型方差分析表

實驗回歸方程系數顯著性檢驗見表5。

表5 回歸系數的顯著性檢驗表

由表5 可知，X4(法向接觸強度)、X5(切向接觸強度)對于單軸抗壓強度具有顯著性影響(P＜0.05)，即X4、X5為模型單軸抗壓強度的主要影響因素，且影響因素的影響順序為：X5＞X4。由T值可知，X4、X5對單軸抗壓強度產生的影響均為正效應，因此增加X4、X5因子上的水平，會使單軸抗壓強度增大。

X1(顆粒接觸模量)、X2(顆粒剛度比)對于彈性模量具有顯著性影響(P＜0.05)，即X1、X2為模型彈性模量的主要影響因素，且影響因素的影響順序為：X1＞X2。由T值可知，X1對單軸抗壓強度產生的影響為正效應，X2對彈性模量產生的影響為負效應，因此增加X1因子上的水平，減小X2因子上的水平，會使彈性模量增大。

X2(顆粒剛度比)、X7(最小顆粒半徑)對于泊松比具有顯著性影響(P＜0.05)，即X2、X7為模型泊松比的主要影響因素，且影響因素的影響順序為：X7＞X2。由T值可知，X2、X7對泊松比產生的影響均為正效應，因此增加X2、X7因子上的水平，會使泊松比增大。

4 結語

通過P-B 試驗設計研究模型宏觀力學特性對于細觀參數的敏感性響應，得出影響數值模型單軸抗壓強度、彈性模量、泊松比的最主要細觀參數，即：切向接觸強度、法向接觸強度這2 個細觀參數是數值模型單軸抗壓強度的主要影響因素，其中，切向接觸強度是首要因素。顆粒接觸模量、顆粒剛度比2 個細觀參數是數值模型彈性模量的主要影響參數，其中，顆粒接觸模量是首要因素。最小顆粒半徑、顆粒剛度比2 個細觀參數是數值模型泊松比的主要影響參數，其中，最小顆粒半徑是首要因素。

文中闡述的影響各宏觀力學特性的主要細觀參數及影響順序，能夠便捷地調整試算方向及參數，有效降低細觀參數選取的復雜性，為后續(xù)類似顆粒流數值計算提供參考。