黃 鵬， 顧 峰， 王麗君

（中國電建集團(tuán)成都勘測設(shè)計研究院有限公司，成都 611130）

砂卵石地層主要存在于上更沖洪積、中更冰水堆積而形成的河漫灘、河流階地等. 具有孔隙度大、黏結(jié)性差、顆粒粒度不均一等特點(diǎn)，廣泛分布于平原地區(qū)［1-3］. 當(dāng)隧道穿越砂卵石地層時，由于該類圍巖自穩(wěn)性較差，變形不明顯，開挖后極易發(fā)生掌子面或洞身段的脆性變形及坍塌工程事故［4-6］. 在數(shù)值模擬方法中，顆粒流離散元法因基于分子動力學(xué)的理論及經(jīng)典力學(xué)理論，從細(xì)觀的角度出發(fā)，研究物體的受力、運(yùn)動及能量，從而得到巖土體宏觀的力學(xué)性質(zhì)［7-8］. 因此，顆粒流軟件對于研究砂卵石層隧道圍巖這種散體物質(zhì)、松散堆積體的力學(xué)特性、大變形及斷裂問題提供了新的研究手段和方法［9-11］.

基于離散元的顆粒流數(shù)值仿真技術(shù)在三軸試驗(yàn)的模擬過程中已有一定的運(yùn)用，周健等［9］在室內(nèi)三軸試驗(yàn)成果的基礎(chǔ)上，開展了針對砂土的雙軸數(shù)值模擬試驗(yàn). 羅勇等［12］對砂土的室內(nèi)常規(guī)三軸試驗(yàn)及其剪切帶的形成和發(fā)展進(jìn)行了數(shù)值模擬，并分別對比了不同圍壓下三維顆粒體試樣與室內(nèi)三軸的應(yīng)力應(yīng)變關(guān)系曲線. 張鐸等［13］通過一系列真三軸離散元數(shù)值試驗(yàn)，詳細(xì)地分析了三維應(yīng)力條件下中主應(yīng)力和應(yīng)力路徑對散體材料峰值強(qiáng)度的影響. 劉君等［14］利用二維顆粒流計算機(jī)仿真技術(shù)對堆石料在一定圍壓下的顆粒破碎情況進(jìn)行了數(shù)值模擬，為研究大壩填筑、水庫蓄水情況下的顆粒破碎過程以及顆粒破碎所引起的壩體變形提供了一條途徑. 邵磊等［15］采用三維顆粒流計算程序，對堆石料的大三軸排水剪切試驗(yàn)過程進(jìn)行了數(shù)值模擬，得出顆粒流可較好模擬堆石料的大三軸試驗(yàn).

綜上可以看出，離散元顆粒流模擬技術(shù)在工程中已有應(yīng)用，但是在砂卵石層隧道開挖模擬中應(yīng)用較少，并且對于砂卵石層圍巖在顆粒流模擬中的形態(tài)參數(shù)、接觸關(guān)系等并不明確，從而導(dǎo)致數(shù)值模擬結(jié)果與工程實(shí)際差異較大［16-20］. 本文在前人研究的基礎(chǔ)上，首先研究了顆粒流軟件中二維模型及三維模型細(xì)觀參數(shù)的取值差異，同時分析顆粒形狀及顆粒粒徑對于試驗(yàn)結(jié)果的影響規(guī)律. 基于砂卵石層隧道圍巖分選性差，大小顆粒混雜的特點(diǎn)，采用數(shù)字圖像處理技術(shù)，得到了砂卵石層不同的顆粒形狀，通過三軸剪切試驗(yàn)進(jìn)行標(biāo)定，得到了不同含水率及相對密實(shí)度情況下砂卵石層的細(xì)觀參數(shù)，為砂卵石層隧道開挖的數(shù)值模擬技術(shù)研究奠定基礎(chǔ).

1 細(xì)觀參數(shù)影響規(guī)律

1.1 二維模型與三維模型的影響

在顆粒流軟件中，二維模型和三維模型的計算時間不同，顆粒的數(shù)量也不同，因此，為了研究二維模型和三維模型中細(xì)觀參數(shù)取值的差異及對最終應(yīng)力應(yīng)變關(guān)系的影響情況，分別建立二維和三維情況下的三軸剪切試驗(yàn)?zāi)Ｐ停唧w模型如圖1所示. 細(xì)觀參數(shù)的選取如表1所示.

圖1 二維及三維模型尺寸圖Fig.1 The size diagrams of 2D and 3D model

表1 不同維數(shù)下細(xì)觀參數(shù)取值Tab.1 Value of microscopic parameters in different dimensions

對二維模型和三維模型分別進(jìn)行三軸剪切試驗(yàn)，可以得到圍壓為100、200、300及400 kPa下的應(yīng)力應(yīng)變曲線，如圖2所示. 結(jié)果表明，當(dāng)二維和三維模型選取相同的細(xì)觀參數(shù)時，兩者得到的應(yīng)力應(yīng)變曲線并不相同，峰值強(qiáng)度差異也很大，同時，二維模型和三維模型達(dá)到峰值強(qiáng)度所對應(yīng)的軸向應(yīng)變值也不相同. 通過繪制二維模型及三維模型的峰值強(qiáng)度比值隨圍壓的變化關(guān)系圖，如圖3所示. 由圖可知，隨著圍壓的增大，二維模型與三維模型的峰值強(qiáng)度比值逐漸降低，穩(wěn)定在2.0左右.

圖2 不同圍壓下二維和三維模型應(yīng)力應(yīng)變曲線Fig.2 Stress-strain curves of two and three dimensional models under different confining pressures

圖3 峰值強(qiáng)度比值隨圍壓的變化關(guān)系圖Fig.3 Peak strength ratio with confining pressure

通過圍壓和峰值強(qiáng)度可以繪制出二維模型和三維模型的摩爾應(yīng)力圓，采用摩爾庫倫準(zhǔn)則，繪制摩爾庫倫包絡(luò)線，可得到二維和三維模型的黏聚力c和內(nèi)摩擦角φ，如表2所示. 由表2可知，當(dāng)二維模型和三維模型采用相同的細(xì)觀參數(shù)，并進(jìn)行相同尺寸的三軸數(shù)值試驗(yàn)時，兩者的顆粒數(shù)量與計算時間差別非常大，如二維的顆粒數(shù)量有886個，而三維模型的顆粒數(shù)量有19 500個，兩者相差約22倍. 二維模型得到的黏聚力c為207.87 kPa，內(nèi)摩擦角φ為37.02°，而三維模型得到的黏聚力c為598.67 kPa，內(nèi)摩擦角φ為39.95°. 兩個模型間的黏聚力c相差約2.88倍，而內(nèi)摩擦角φ相差卻不大，說明二維和三維模型影響的主要宏觀力學(xué)指標(biāo)為黏聚力c，因此，在之后的參數(shù)標(biāo)定時，應(yīng)當(dāng)選擇合適的模型進(jìn)行標(biāo)定，二維和三維模型的參數(shù)并不能直接套用.

表2 不同維數(shù)下得到的物理力學(xué)參數(shù)指標(biāo)Tab.2 The physical and mechanical parameters obtained under different dimensions

1.2 顆粒粒徑與試驗(yàn)尺寸的影響

在進(jìn)行離散元數(shù)值模擬時，顆粒的大小影響計算的效率，若采用實(shí)際顆粒尺寸，則進(jìn)行大型模型的建立時，需要上百萬個顆粒，計算效率非常低下. 為了研究顆粒尺寸對于試驗(yàn)結(jié)果的影響規(guī)律，以顆粒平均粒徑D 與試驗(yàn)尺寸L 的比值（D/L）作為主要變量，研究不同D/L 情況下，二維雙軸壓縮試驗(yàn)的應(yīng)力應(yīng)變曲線的差異，計算圖示如圖4 所示. 不同計算工況及細(xì)觀參數(shù)取值如表3所示.

對不同D/L 情況下的數(shù)值模型分別進(jìn)行三軸剪切試驗(yàn)，可以得到圍壓為100、200、300 及400 kPa 下的應(yīng)力應(yīng)變曲線，如圖4 所示. 由圖4 可知，當(dāng)D/L 的比值逐漸減小時，所得到的應(yīng)力應(yīng)變曲線逐漸接近，說明顆粒的粒徑對于試驗(yàn)結(jié)果有影響，當(dāng)D/L的比值小于1/60時，所得到的應(yīng)力應(yīng)變曲線差異不大，其峰值強(qiáng)度與對應(yīng)的軸向應(yīng)變值差異均很小. 因此，在進(jìn)行顆粒流軟件進(jìn)行數(shù)值模擬時，應(yīng)當(dāng)控制好顆粒粒徑的大小，過小會使得計算效率變低而對結(jié)果無明顯改善，過大則會影響計算結(jié)果.

表3 顆粒尺寸及平均粒徑細(xì)觀參數(shù)取值Tab.3 Value of mesoscopic parameters of particle size and average particle size

圖4 不同圍壓下應(yīng)力應(yīng)變曲線Fig.4 Stress-strain curves under different confining pressures

通過圍壓和峰值強(qiáng)度可以繪制出不同D/L情況下的摩爾應(yīng)力圓，采用摩爾庫倫準(zhǔn)則，繪制摩爾庫倫包絡(luò)線，可得到不同D/L下的黏聚力c和內(nèi)摩擦角φ，如表4所示. 由表4可知，隨著顆粒尺寸D/L的減小，顆粒的數(shù)量增加，但是峰值強(qiáng)度有著減小的趨勢，當(dāng)D/L的值小于1/60時，顆粒的尺寸對于峰值強(qiáng)度的影響減弱，當(dāng)D/L 的值大于1/30 時，顆粒的尺寸效應(yīng)明顯，如D/L 為20 時所得到的黏聚力c 大概是D/L 為20 時的1.48 倍.但是，隨著D/L的增大，顆粒的數(shù)量及計算所需時間也在呈指數(shù)增加，為了保證后續(xù)模型計算的有效性，顆粒的尺寸不宜太小或太大，應(yīng)當(dāng)保持在D/L為60～100之間為宜，顆粒數(shù)量在10 000左右效果較好.

表4 不同比值下峰值強(qiáng)度及物理力學(xué)指標(biāo)表Tab.4 Peak strength and physical mechanics indexes at different ratios

1.3 顆粒形狀的影響規(guī)律

在離散元軟件中，基本的顆粒形狀有ball 和clump 兩種，其中，ball 顆粒為圓盤或球體，clump 為采用Bubblepack算法，用不同半徑的球形單元對不規(guī)則幾何邊界進(jìn)行填充得到的基本顆粒. 為了研究顆粒形狀對于最終結(jié)果的影響規(guī)律，本文建立了三種不同的顆粒形狀，分別為圓形、方形和橢圓形為基本顆粒形狀的數(shù)值模型. 顆粒的形狀及填充效果如圖5所示，模型中的細(xì)觀參數(shù)取值如表5所示.

圖5 不同顆粒形狀及填充效果局部放大圖Fig.5 Different particle shapes and filling effects

對不同顆粒形狀下的數(shù)值模型分別進(jìn)行三軸剪切試驗(yàn)，可以得到圍壓為100、200、300、400 kPa下的應(yīng)力應(yīng)變曲線，如圖6所示.

由圖6可知，當(dāng)顆粒形狀為橢圓形時，得到的峰值強(qiáng)度最大，其對應(yīng)的軸向應(yīng)變也最大，當(dāng)顆粒形狀為圓形時，所得到的峰值強(qiáng)度最小，對應(yīng)的軸向應(yīng)變也最小. 因此，當(dāng)顆粒的形狀越扁平時，顆粒與顆粒之間的咬合力增加，接觸面積也增大. 所得到的峰值強(qiáng)度也越大. 所以在用顆粒流軟件進(jìn)行數(shù)值模擬時，若所模擬的對象形狀不均勻，則應(yīng)考慮形狀對于計算結(jié)果的影響.

圖6 不同圍壓下應(yīng)力應(yīng)變曲線Fig.6 Stress-strain curves under different confining pressures

通過圍壓和峰值強(qiáng)度可以繪制出不同顆粒形狀下的摩爾應(yīng)力圓，采用摩爾庫倫準(zhǔn)則，繪制摩爾庫倫包絡(luò)線，可得到不同顆粒形狀下的黏聚力c 和內(nèi)摩擦角φ，如表6 所示. 由表6 可知，在不同顆粒形狀但相同顆粒面積的情況下，橢圓形顆粒和方形顆粒所得到的黏聚力c 和內(nèi)摩擦角φ值相近，但兩者得到的峰值強(qiáng)度卻不相同，這可能是因?yàn)閮烧叩男螤畈煌?，影響了試?yàn)過程中破壞的過程，但是兩者都具有一定棱角，且其黏結(jié)強(qiáng)度和抗拉強(qiáng)度相同，導(dǎo)致了最終得到了相近的黏聚力和內(nèi)摩擦角. 這也側(cè)面反映了單純地考慮顆粒形狀對黏聚力c和內(nèi)摩擦角φ值的影響并不全面，顆粒情況應(yīng)當(dāng)有一指標(biāo)，在本例中方形和橢圓形該指標(biāo)相同，所以才會導(dǎo)致相近c(diǎn)、φ值的出現(xiàn). 而方形與橢圓形和圓形的對比，則說明了顆粒形狀對于試驗(yàn)結(jié)果的影響.

表6 不同比值下峰值強(qiáng)度及物理力學(xué)指標(biāo)表Tab.6 Peak strength and physical mechanics indexes at different ratios

2 砂卵石層細(xì)觀參數(shù)標(biāo)定

2.1 砂卵石層顆粒粒徑的確定

本次試樣取自川西南某水電站庫區(qū)岸灘公路隧道，為灰白色卵礫石土，顆粒大小混雜，分選性差. 對試樣進(jìn)行顆粒篩分試驗(yàn)，可以得到顆粒粒徑級配曲線圖，如圖7所示. 因此，在顆粒流軟件中進(jìn)行相應(yīng)的數(shù)值模擬時，采用相同的級配關(guān)系，其顆粒粒徑分布及形狀如表7所示.

圖7 顆粒粒徑級配曲線Fig.7 Particle size grade fitting line

表7 離散元顆粒形狀、分布及粒徑參數(shù)取值Tab.7 The parameters of shape，distribution and particle size of discrete element particles

2.2 砂卵石層形狀的確定

砂卵石層具有分選性差，顆粒分布極不均勻的特點(diǎn)，若采用傳統(tǒng)的有限元方法進(jìn)行計算，則相對誤差會很大，無法體現(xiàn)顆粒的離散性. 為了更好體現(xiàn)砂卵石層的顆粒離散的特點(diǎn)，采用離散元顆粒流軟件進(jìn)行相應(yīng)的數(shù)值模擬. 在本次數(shù)值模擬中，基于clump 來模擬實(shí)際砂卵石層，采用Matlab 中的Candy 算法求取圖像的邊界坐標(biāo)，從而可以獲得砂卵石層圖像的顆粒形狀，得到的數(shù)字圖像處理結(jié)果.

2.3 砂卵石層細(xì)觀參數(shù)標(biāo)定

通過對不同相對密實(shí)度及不同含水率的試樣進(jìn)行大型三軸剪切試驗(yàn)，可以得到不同工況條件下的峰值強(qiáng)度及對應(yīng)的軸向應(yīng)變，具體如表8所示. 在進(jìn)行細(xì)觀參數(shù)的標(biāo)定時，峰值強(qiáng)度完全相同的情況幾乎只是理想狀態(tài)，所以需要有一定的誤差允許，在本文中，選擇峰值強(qiáng)度與試驗(yàn)數(shù)據(jù)相差不超過5%作為標(biāo)準(zhǔn)進(jìn)行標(biāo)定，因此可以得到不同工況條件下砂卵石層的峰值強(qiáng)度限值. 根據(jù)上述的標(biāo)定原則，不斷地改變顆粒流軟件的細(xì)觀參數(shù)，使之滿足條件，得到的細(xì)觀參數(shù)標(biāo)定結(jié)果如表9所示.

表8 砂卵石層大型三軸試驗(yàn)結(jié)果Tab.8 Results of large-scale triaxial test on sand-pebble bed

表9 砂卵石層細(xì)觀參數(shù)標(biāo)定結(jié)果Tab.9 Calibration results of sand and pebble bed microparameters

3 結(jié)論

1）相同細(xì)觀參數(shù)下，二維模型和三維模型所得到的應(yīng)力應(yīng)變曲線不同，且隨著圍壓的增大，二維模型和三維模型所得到的峰值強(qiáng)度的比值逐漸趨于2.0. 二維模型和三維模型所影響的主要宏觀力學(xué)性質(zhì)為黏聚力c，對于內(nèi)摩擦角φ 的影響很小.

2）顆粒粒徑D與試樣尺寸L的比值對于得到的應(yīng)力應(yīng)變曲線有影響，且當(dāng)顆粒粒徑與試樣尺寸比值大于1/60時，其影響減弱，但是計算效率降低，因此在實(shí)際計算時，D/L的比值應(yīng)當(dāng)在1/60附近，不宜過大或過小.

3）顆粒形狀對于得到的應(yīng)力應(yīng)變曲線有很大影響，顆粒的越扁平，其峰值強(qiáng)度越高，顆粒越接近圓形，其峰值強(qiáng)度越低；在不同顆粒形狀但相同顆粒面積的情況下，橢圓形顆粒和方形顆粒所得到的黏聚力c和內(nèi)摩擦角φ 值相近，但兩者得到的峰值強(qiáng)度卻不相同.

4）通過與不同相對密實(shí)度和含水率情況下砂卵石層大型三軸剪切試驗(yàn)所得到的峰值強(qiáng)度進(jìn)行對比，得到了不同相對密實(shí)度及不同含水率情況下的砂卵石層細(xì)觀參數(shù)，為砂卵石層隧道開挖模擬技術(shù)奠定了基礎(chǔ).