宋 彧，冷國(guó)杰，楊安玉

（1.浙江國(guó)貿(mào)集團(tuán)東方機(jī)電工程股份有限公司，浙江 杭州 310000；2.南京水利科學(xué)研究院，江蘇 南京 210024；3.水利部農(nóng)村電氣化研究所，浙江 杭州 310012）

1 問題的提出

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

當(dāng)前，常用的顆粒流模擬程序PFC 中，數(shù)值模型的宏觀基本物理特性無法通過直接賦值的形式實(shí)現(xiàn)，只能對(duì)構(gòu)成模型顆粒的幾何特性及顆粒間接觸的細(xì)觀力學(xué)參數(shù)進(jìn)行賦值[7-8]。且必須通過大量試算，標(biāo)定出細(xì)觀參數(shù)，才能使數(shù)值模型表現(xiàn)出所需宏觀物理特征。

文章采用PFC2D 程序建立顆粒流接觸黏結(jié)模型進(jìn)行單軸壓縮模擬，并通過Plackett-Burman 設(shè)計(jì)獲得模型各細(xì)觀參數(shù)的敏感性響應(yīng)，篩選出影響各項(xiàng)宏觀力學(xué)特性的主要細(xì)觀參數(shù)因子，提高顆粒流模型細(xì)觀參數(shù)標(biāo)定的便捷性。

2 Plackett-Burman 試驗(yàn)設(shè)計(jì)

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

P-B 試驗(yàn)設(shè)計(jì)采用線性函數(shù)篩選得出重要因素，不考慮相互影響，其中，函數(shù)方程為：

式中：Y為響應(yīng)函數(shù)值，Xi為設(shè)計(jì)函數(shù)值。β0和βi為回歸系數(shù)，通過回歸系數(shù)的顯著水平P值來反映Xi對(duì)Y的影響程度。

3 基于Plackett-Burman 試驗(yàn)設(shè)計(jì)的模型計(jì)算

3.1 數(shù)值模型建立

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

表1 數(shù)值模型的細(xì)觀參數(shù)表

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

圖1 顆粒流數(shù)值模型圖

3.2 單因子試驗(yàn)確定因子水平

3.2.1 顆粒接觸模量

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

圖2 顆粒接觸模量單因子試驗(yàn)應(yīng)力—應(yīng)變曲線圖

根據(jù)單因子試驗(yàn)，顆粒接觸模量與模型彈性模量的關(guān)系見圖3。

圖3 顆粒接觸模量對(duì)模型彈性模量的影響圖

3.2.2 法向接觸強(qiáng)度/切向接觸強(qiáng)度

對(duì)法向接觸強(qiáng)度/切向接觸強(qiáng)度這一細(xì)觀參數(shù)進(jìn)行單因子試驗(yàn)，法向接觸強(qiáng)度/切向接觸強(qiáng)度分別取10，20，30，40，50，60 MPa，其余參數(shù)與初始模型參數(shù)相同。建立6 個(gè)模型分別進(jìn)行數(shù)值模擬得到應(yīng)力—應(yīng)變曲線（見圖4）。

圖4 法向接觸強(qiáng)度/切向接觸強(qiáng)度單因子試驗(yàn)應(yīng)力—應(yīng)變曲線圖

根據(jù)單因子試驗(yàn)，法向接觸強(qiáng)度/切向接觸強(qiáng)度與模型單軸抗壓強(qiáng)度的關(guān)系見圖5。

圖5 法向/切向接觸強(qiáng)度對(duì)模型單軸抗壓強(qiáng)度的影響圖

3.2.3 顆粒剛度比

對(duì)顆粒剛度比這一細(xì)觀參數(shù)進(jìn)行單因子試驗(yàn)，顆粒剛度比分別取0.5，1.0，1.5，2.0，2.5，3.0，3.5，其余參數(shù)與初始模型參數(shù)相同。建立7 個(gè)模型分別進(jìn)行數(shù)值模擬得到應(yīng)力—應(yīng)變曲線（見圖6）。

圖6 顆粒剛度比單因子試驗(yàn)應(yīng)力—應(yīng)變曲線圖

根據(jù)單因子試驗(yàn)，顆粒剛度比與模型泊松比的關(guān)系見圖7。

圖7 顆粒剛度比對(duì)模型泊松比的影響圖

3.2.4 因子水平確定

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

表2 Plackett-Burman 試驗(yàn)設(shè)計(jì)因子水平表

3.3 Plackett-Burman 試驗(yàn)設(shè)計(jì)與結(jié)果分析

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

本文采用Minitab 軟件生成篩選試驗(yàn)設(shè)計(jì)表，根據(jù)P-B 試驗(yàn)設(shè)計(jì)進(jìn)行12 次計(jì)算分析，設(shè)計(jì)及試驗(yàn)結(jié)果見表3。

表3 Plackett-Burman 試驗(yàn)設(shè)計(jì)及響應(yīng)值表

擬合一階模型為：

對(duì)模型進(jìn)行方差分析（結(jié)果見表4）可知，所擬合的3 個(gè)回歸方程均達(dá)到顯著（P＜0.05），決定系數(shù)R2分別為0.985 0、0.991 0、0.916 0，表明分別用方程擬合7 個(gè)細(xì)觀因子與3 個(gè)宏觀力學(xué)特性之間的關(guān)系是可行的。

表4 回歸模型方差分析表

實(shí)驗(yàn)回歸方程系數(shù)顯著性檢驗(yàn)見表5。

表5 回歸系數(shù)的顯著性檢驗(yàn)表

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

X1(顆粒接觸模量)、X2(顆粒剛度比)對(duì)于彈性模量具有顯著性影響(P＜0.05)，即X1、X2為模型彈性模量的主要影響因素，且影響因素的影響順序?yàn)椋篨1＞X2。由T值可知，X1對(duì)單軸抗壓強(qiáng)度產(chǎn)生的影響為正效應(yīng)，X2對(duì)彈性模量產(chǎn)生的影響為負(fù)效應(yīng)，因此增加X1因子上的水平，減小X2因子上的水平，會(huì)使彈性模量增大。

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

4 結(jié)語

通過P-B 試驗(yàn)設(shè)計(jì)研究模型宏觀力學(xué)特性對(duì)于細(xì)觀參數(shù)的敏感性響應(yīng)，得出影響數(shù)值模型單軸抗壓強(qiáng)度、彈性模量、泊松比的最主要細(xì)觀參數(shù)，即：切向接觸強(qiáng)度、法向接觸強(qiáng)度這2 個(gè)細(xì)觀參數(shù)是數(shù)值模型單軸抗壓強(qiáng)度的主要影響因素，其中，切向接觸強(qiáng)度是首要因素。顆粒接觸模量、顆粒剛度比2 個(gè)細(xì)觀參數(shù)是數(shù)值模型彈性模量的主要影響參數(shù)，其中，顆粒接觸模量是首要因素。最小顆粒半徑、顆粒剛度比2 個(gè)細(xì)觀參數(shù)是數(shù)值模型泊松比的主要影響參數(shù)，其中，最小顆粒半徑是首要因素。

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