土工袋雙軸壓縮的離散元數(shù)值模擬

高軍軍,湯 雷,劉斯宏,李 建,王恩準(zhǔn)

(1.南京水利科學(xué)研究院材料結(jié)構(gòu)研究所,江蘇 南京 210029;2.河海大學(xué)水利水電學(xué)院,江蘇 南京 210098; 3.南京瑞迪高新技術(shù)有限公司,江蘇 南京 210029)

土工袋技術(shù)是近年來通過系列試驗(yàn)研究、理論分析及工程應(yīng)用而發(fā)展起來的一項(xiàng)地基處理新技術(shù)。Matsuoka等[1-2]通過一系列的室內(nèi)試驗(yàn)研究并驗(yàn)證了土工袋的加筋機(jī)理,即將散狀土體裝入土工編織袋,利用壓縮過程中袋子產(chǎn)生的張力約束袋內(nèi)土體,達(dá)到提高袋內(nèi)土體強(qiáng)度的目的。土工袋已在房屋基礎(chǔ)、公路與鐵路路基、溝槽回填及擋土墻等工程中得到廣泛應(yīng)用。目前對土工袋的研究主要以室內(nèi)試驗(yàn)為主：劉斯宏等[3]通過現(xiàn)場平板載荷試驗(yàn)驗(yàn)證了土工袋具有比一般土體高得多的承載能力;高軍軍等[4-5]通過室內(nèi)試驗(yàn)和數(shù)值模擬驗(yàn)證土工袋可顯著提高淤泥土的強(qiáng)度,其極限抗壓強(qiáng)度約為普通混凝土的1/10;李玲君等[6]探討了不同袋內(nèi)充填材料對土工袋動力特性的影響;相關(guān)學(xué)者還通過鋁棒摩擦試驗(yàn)研究了不同條件下柔性土工袋層間摩擦與剛性體間摩擦特征的區(qū)別[7]。

對土工袋的數(shù)值模擬研究相對較少,劉斯宏等[8]將袋子張力的作用等效為附加應(yīng)力,利用彈塑性有限元法對加固后的軟基承載力進(jìn)行了數(shù)值模擬;Ansari等[9]將土和袋子分別采用不同的模型利用有限元法研究單個土工袋在豎向荷載下的力學(xué)反應(yīng)。但土工袋是一種復(fù)雜的復(fù)合型材料,存在非均質(zhì)性和介質(zhì)不連續(xù)性,尤其是不同材料介質(zhì)間的接觸問題,采用基于連續(xù)介質(zhì)理論的有限元法反映其真實(shí)力學(xué)特性的難度較大,土工袋固有的間隙也直接制約有限元法計算的有效進(jìn)行。故而采用基于非連續(xù)介質(zhì)理論的顆粒離散元法對土工袋的力學(xué)特性進(jìn)行數(shù)值模擬來驗(yàn)證土工袋的加筋機(jī)理有一定的現(xiàn)實(shí)意義。本文對土工袋在雙軸壓縮條件下的力學(xué)特性進(jìn)行離散元數(shù)值模擬,從細(xì)觀角度分析土工袋內(nèi)部應(yīng)力分布及發(fā)展規(guī)律。

1 土工袋模型的生成

考慮真實(shí)土體顆粒分布的空間隨機(jī)性,本文采用蒙特卡洛隨機(jī)算法生成一定顆粒粒徑及顆粒級配的袋內(nèi)土體樣本,其中土顆粒間采用彈簧-阻尼器-滑塊接觸模型,袋子顆粒間采用皮筋-阻尼器接觸模型。

1.1 蒙特卡洛隨機(jī)算法

蒙特卡洛隨機(jī)算法又稱計算機(jī)隨機(jī)模擬方法,它是一種基于隨機(jī)數(shù)的計算方法[10-12],在大數(shù)定理和中心極限定理等概率論的基礎(chǔ)上,對隨機(jī)變量進(jìn)行方差判定來分析隨機(jī)問題的不確定程度。

設(shè)X1、X2、…、Xn(n為隨機(jī)變量的總數(shù))為獨(dú)立且均勻分布的隨機(jī)變量序列,其期望值為E,標(biāo)準(zhǔn)差為S,則正態(tài)隨機(jī)數(shù)可表示為

(1)

顆粒生成程序中,根據(jù)邊界條件和顆粒數(shù)可以求出顆粒的平均半徑,輸入比平均半徑較小的一個值R0和預(yù)設(shè)標(biāo)準(zhǔn)差S0,通過式(2)就可以得到第i個顆粒粒徑Ri,如果Ri小于0,則回到隨機(jī)算法重新生成100個處在區(qū)間[0,1]內(nèi)的隨機(jī)變量并重新得到Q100,根據(jù)式(2)得到新的粒徑值,重復(fù)判斷,直到Ri大于0,進(jìn)入下一個循環(huán)生成第i+1個顆粒的粒徑Ri+1。

Ri=R0+S0Q100

(2)

根據(jù)顆粒單元的交叉(重合)情況通過分級衰減粒徑大小的方法再重新定位顆粒單元,如此往復(fù)直至生成較滿意的顆粒系統(tǒng)。此顆粒系統(tǒng)需滿足兩個條件:①顆粒在給定空間內(nèi)具有各向均勻性;②顆粒體的初始重疊量不宜過大。

1.2 顆粒體生成

首先確定矩形試樣尺寸為長度40 cm、高度10 cm(土工袋的合理尺寸[13])。如果顆粒最大粒徑為9.0 mm,生成初始顆粒數(shù)為1 250,再根據(jù)衰減系數(shù)對發(fā)生重疊的顆粒進(jìn)行分級粒徑衰減,經(jīng)過5級衰減后,得到粒徑分別為9.00 mm、7.20 mm、6.12 mm、5.20 mm、4.16 mm、3.33 mm。此時生成的顆粒樣本雖然不存在重疊現(xiàn)象,但是部分顆粒間隙較大,存在顆粒騰空現(xiàn)象,需要再進(jìn)行自重堆積得到新的顆粒位置狀態(tài)。此時顆粒樣本在頂部并不平整,會影響豎向均布荷載的施加,可在其頂部先施加10 kPa的預(yù)荷載,同時限制側(cè)向變形為0,待穩(wěn)定后得到試驗(yàn)最終樣本。

最終生成的樣本包含1 250個顆粒,其中顆粒數(shù)按粒徑從大到小分別為:180、148、156、189、184、393,顆粒級配曲線如圖1所示,最終的顆粒孔隙比為0.20,不均勻系數(shù)為2.1。

圖1 顆粒級配曲線

1.3 顆粒接觸模型

在顆粒離散單元中,宏觀系統(tǒng)的力學(xué)行為通過細(xì)觀顆粒間的接觸模型體現(xiàn)出來。基于牛頓第二定律,通過顆粒間的作用力求得各顆粒的加速度,加速度對時間進(jìn)行一次積分即可求得顆粒的速度,加速度對時間進(jìn)行二次積分則可求得顆粒的位移。

在土工袋的離散元模擬計算中,袋內(nèi)土顆粒間的接觸簡化為彈簧-阻尼器-滑塊模型,如圖2(a)所示,彈簧主要模擬顆粒間的壓縮和切向所用,滑塊主要模擬顆粒間的滑移,而阻尼器主要是為了消耗系統(tǒng)中多余的能量使系統(tǒng)快速趨向平衡狀態(tài);土工袋通過在土顆粒的外面增加一層具有張力的小顆粒模擬,袋子顆粒間的接觸簡化為皮筋-阻尼器模型,如圖2(b)所示,其間的接觸模型與土顆粒相同,但只有法向接觸,無切向接觸,且只受拉不受壓。

圖2 袋子DEM數(shù)值模擬顆粒接觸模型

2 無袋子情況下的土體強(qiáng)度

為了反映土工袋對袋內(nèi)土體的加筋作用(即提高袋內(nèi)土體強(qiáng)度),先對無袋子情況下近似同等體積的土體進(jìn)行數(shù)值模擬,分析其在雙軸壓縮條件下的力學(xué)強(qiáng)度。

表1 土體試樣DEM計算參數(shù)

模擬試驗(yàn)在4塊相對獨(dú)立、互不影響的剛性加載板的范圍內(nèi)進(jìn)行,如圖3所示。顆粒密度為2.7 g/cm3,顆粒間的摩擦角為20°;假定顆粒與剛性加載板之間無摩擦,避免加載板對試樣端部變形及試樣內(nèi)部應(yīng)力的影響。具體計算參數(shù)如表1所示。

圖3 無袋子情況下的土體試樣雙軸壓縮模擬

加載過程中對試樣施加以50 kPa的穩(wěn)定圍壓,再分別對頂部和底部的加載板施加以大小為0.01 cm/s、方向相反的豎向等應(yīng)變率加載速度。

圖4為無袋子包圍的土顆粒試樣偏應(yīng)力及體應(yīng)變與軸向應(yīng)變之間的關(guān)系曲線。從圖中可以看出土體在雙軸壓縮過程中應(yīng)力隨著應(yīng)變非線性增長。而土體的體應(yīng)變在壓縮過程中先發(fā)生微小的剪縮現(xiàn)象,隨即開始剪脹,且體應(yīng)變隨軸向應(yīng)變的增大而增大,當(dāng)土體發(fā)生破壞后,剪脹值有所降低。

圖4 無袋子土體的應(yīng)力-應(yīng)變關(guān)系曲線

圖5 土顆粒試樣內(nèi)部力鏈分布

圖5為土顆粒試樣達(dá)到極限強(qiáng)度時的力鏈結(jié)構(gòu)圖,線條的粗細(xì)代表顆粒間接觸力的大小。力鏈?zhǔn)穷w粒物質(zhì)力學(xué)的基本概念,孫其誠等[14]認(rèn)為顆粒在外荷載作用下相互擠壓形成接觸網(wǎng)絡(luò),該網(wǎng)絡(luò)是支撐外荷載的物質(zhì)基礎(chǔ),就如同土力學(xué)中土骨架的概念一樣[15-16]。從圖5可看出,顆粒間的強(qiáng)力鏈(即較粗力鏈)分布較均勻,基本上與大主應(yīng)力方向平行,而兩側(cè)顆粒間的力鏈較弱,說明土顆粒試樣的破壞由豎向荷載主導(dǎo),這與事實(shí)相符。

為了通過摩爾圓來分析土體的強(qiáng)度特性,又對該試樣作了不同圍壓下(100 kPa、200 kPa)的雙軸壓縮模擬。從試驗(yàn)樣本的摩爾強(qiáng)度線可以看出,不同圍壓條件下的應(yīng)力摩爾圓基本相切于一直線,故本次試驗(yàn)樣本的內(nèi)摩擦角可取24.3°。

3 土工袋的力學(xué)強(qiáng)度

為了研究土工袋的力學(xué)性能,只需在原試驗(yàn)樣本上加一圈袋子小顆粒,如圖6所示,假設(shè)顆粒級配及樣本內(nèi)摩擦角不變。袋子顆粒間保留一定距離,允許袋子在外荷載作用下發(fā)生伸縮變形;且假設(shè)間距相等,為袋子顆粒粒徑的2.5倍。該土工袋試樣模型中,土顆粒數(shù)為1 199,袋子顆粒數(shù)為396,袋子顆粒粒徑為0.1 cm,其密度設(shè)為0.8 g/m3。土顆粒的計算參數(shù)如表1所示,袋子顆粒的計算參數(shù)見表2,模擬中不考慮袋子與剛性邊間的摩擦作用。

加載方式與無袋子情況相同。需要注意的是,在計算主應(yīng)變時,將土工袋簡化為理想的矩形結(jié)構(gòu),即利用上下剛性邊的相對位移計算土工袋的軸向應(yīng)變,利用左右剛性邊的相對位移計算土工袋的側(cè)向應(yīng)變。

表2 土工袋試樣DEM計算參數(shù)

圖6 土工袋試樣雙軸壓縮模擬

圖7為土工袋在圍壓50 kPa條件下的偏應(yīng)力-軸向應(yīng)變和體應(yīng)變-軸向應(yīng)變關(guān)系曲線。從圖7可看出,土工袋的偏應(yīng)力隨著軸向應(yīng)變的增長逐漸增大,直到土工袋發(fā)生破壞后,迅速下降;而土工袋的體應(yīng)變在加載初期也發(fā)生微小的剪縮現(xiàn)象,隨后快速膨脹,到后期增長速率變小,這與無袋子情況下的土體變形規(guī)律一致。不同的是,土工袋破壞時的偏應(yīng)力(680 kPa)遠(yuǎn)大于無袋子的情況(70 kPa),并且土工袋破壞時的總體變形明顯大于無袋子的情況。袋子具有可延性,到后期袋子的延伸速率下降,逐漸達(dá)到袋子的抗拉強(qiáng)度。從圖中還可以看出,隨著軸向應(yīng)變的增長,偏應(yīng)力有好幾處(紅圈標(biāo)記)發(fā)生一定的下降后重新上升,這是由于雖然土工袋仍未破壞,但此時袋內(nèi)的土體早已達(dá)到極限狀態(tài),袋內(nèi)土體產(chǎn)生局部應(yīng)力集中,部分土顆粒發(fā)生旋轉(zhuǎn),在袋子的側(cè)向約束下重新達(dá)到新的平衡狀態(tài),而在此過程中體應(yīng)變并無明顯變化。

圖7 土工袋的應(yīng)力-應(yīng)變關(guān)系曲線

圖8(a)為土工袋破壞前某一狀態(tài)下(軸向應(yīng)變ε1=3%時)相鄰袋子顆粒間的張力大小沿土工袋周邊的分布情況,可見,袋子張力分布并不是均勻的,在袋子不同的位置張力大小不同,但基本上在一個平均值附近上下波動。對所有袋子顆粒間的張力取平均值可以得到袋子平均張力的變化情況,圖8(b)為模擬過程中袋子的平均張力隨著軸向應(yīng)變的變化。袋子的平均張力并不呈線性增長趨勢,這是由于在變形過程中袋子周長的伸長與土工袋軸向應(yīng)變并不呈線性關(guān)系。袋子破壞時袋子平均張力接近12.5 kN/m。

圖8 土工袋袋子張力的變化

根據(jù)劉斯宏等[2]對土工袋加筋原理的闡述,將袋子張力等效為作用于土工袋軸向及側(cè)向的附加應(yīng)力(2T/B,2T/H),強(qiáng)度公式可表示為

(3)

圖9 袋內(nèi)土體主應(yīng)力之間的關(guān)系

式中:σ1為土工袋大主應(yīng)力;σ3為土工袋小主應(yīng)力;T為袋子的張力;B為土工袋的寬度;H為土工袋的高度;KP為被動土壓力系數(shù),KP=(1+sinφ)/(1-sinφ),φ=24.3°。

取圖7中紅圈內(nèi)的大小主應(yīng)力,再結(jié)合圖8中對應(yīng)軸向應(yīng)變下的平均張力,可以得到袋內(nèi)土體主應(yīng)力之間的關(guān)系如圖9所示,大小主應(yīng)力之間基本呈線性關(guān)系,且斜率近似與KP一致,驗(yàn)證了理論強(qiáng)度式(3)的合理性。

假設(shè)土工袋在壓縮過程中體積不變,可以推得土工袋的軸向應(yīng)變ε1與側(cè)向應(yīng)變ε3之間的關(guān)系為ε3=ε1/(ε1-1)。圖10顯示DEM數(shù)值模擬得到的土工袋軸向應(yīng)變與側(cè)向應(yīng)變之間的關(guān)系(黑色實(shí)線),發(fā)現(xiàn)在加載初期土工袋基本上滿足體積不變的假設(shè),但當(dāng)軸向應(yīng)變超過2%以后,真實(shí)的側(cè)向應(yīng)變絕對值大于假設(shè)體積不變得到的值,即土工袋的體積發(fā)生了變化。DEM數(shù)值模擬中的ε1與ε3的關(guān)系更接近于多項(xiàng)式表達(dá)式(如圖10紅色虛線所示):

ε3= 0.007 85ε13-0.138 34ε12-

0.800 65ε1-0.049 87

(4)

圖12 土工袋內(nèi)部顆粒的細(xì)觀分析

圖10 土工袋軸向應(yīng)變與側(cè)向應(yīng)變的關(guān)系

通過幾何換算,袋子的應(yīng)變εbag可以表示為

(5)

式中:m為初始狀態(tài)下土工袋的寬度與高度的比值。

將式(5)代入T=kεbag(其中k為袋子接觸模型的抗拉系數(shù)),并與式(4)一起代入式(3)可以推導(dǎo)出適合此雙軸壓縮試驗(yàn)的應(yīng)力-應(yīng)變關(guān)系:

(6)

式中:H0、B0分別為土工袋的初始高度和寬度。圖11為DEM數(shù)值模擬結(jié)果與理論計算值的對比,可見推導(dǎo)公式與模擬結(jié)果吻合較好。

圖11 DEM數(shù)值模擬結(jié)果與推導(dǎo)公式的對比

4 土工袋細(xì)觀結(jié)構(gòu)變化

圖12為不同加載時刻的土工袋細(xì)觀分析圖。ε1=0時為圍壓50 kPa的情況,因?yàn)轭w粒處于等向壓縮狀態(tài),土工袋的內(nèi)部應(yīng)力比較均勻,力鏈群形成環(huán)狀分布,表明此時外荷載由土顆粒均勻承擔(dān)。當(dāng)上下剛邊以一定速度壓縮土工袋時,土工袋被壓扁,內(nèi)部應(yīng)力重新分布,ε1=2.8%時可以觀察到內(nèi)部應(yīng)力出現(xiàn)局部集中化,力鏈方向大致為大主應(yīng)力的方向;而土顆粒的位移較小且較亂,由于此時偏應(yīng)力不大,在雙軸壓縮的情況下,土顆粒大致都往內(nèi)部移動,導(dǎo)致土工袋在加載初期出現(xiàn)微小剪縮現(xiàn)象,同時袋子4個邊角上的土顆粒有向外擠出的現(xiàn)象。當(dāng)ε1=5.1%時,內(nèi)部應(yīng)力向中部集中,顆粒位移加大,顆粒基本上向土工袋兩側(cè)跑動,中部的顆粒位移較小;當(dāng)ε1=11.1%時,土體內(nèi)部出現(xiàn)滑移帶,顆粒有明顯向兩側(cè)發(fā)生大位移的趨勢;隨著軸向應(yīng)力的繼續(xù)增加,內(nèi)部顆粒應(yīng)力重新調(diào)整,直至土工袋內(nèi)部基本上都處于強(qiáng)力鏈的狀態(tài),在軸向應(yīng)變達(dá)到12.6%時發(fā)生破壞。在加載后期,顆粒的運(yùn)動也呈對稱性地向兩端擴(kuò)散,出現(xiàn)明顯的滑移帶,直至土工袋發(fā)生破壞時,土顆粒有明顯向右下側(cè)滑移的趨勢(圖12(i)中圓圈處),這是由于袋子就在此處發(fā)生了破壞(如圖12(j)所示),對袋內(nèi)土體失去約束作用,導(dǎo)致顆粒向外跑動,同時發(fā)現(xiàn)此時強(qiáng)力鏈開始從右下側(cè)消失,表明袋子破壞是土工袋達(dá)到極限強(qiáng)度的直接原因。

5 結(jié) 論

a. 土工袋破壞時的偏應(yīng)力遠(yuǎn)大于無袋子土體破壞時的偏應(yīng)力,且土工袋破壞時的總體變形明顯大于無袋子土體破壞時的總體變形。

b. 袋子張力分布并不是均勻的,在袋子不同的位置張力大小不同,但基本上在一個平均值附近上下波動。

c. 土工袋雙軸壓縮DEM數(shù)值模擬得到的應(yīng)力-應(yīng)變關(guān)系與理論計算值吻合較好,軸向應(yīng)變與側(cè)向應(yīng)變的關(guān)系更加接近于某個多項(xiàng)式函數(shù)。

d. 在豎向荷載作用下,土工袋內(nèi)部顆粒間力鏈由環(huán)狀均勻分布向大主應(yīng)力方向發(fā)展,顆粒位移呈對稱性向兩側(cè)擴(kuò)散趨勢,且袋子破壞是土工袋達(dá)到極限強(qiáng)度的直接原因。