基于DEM-CFD耦合方法的砂土滲透破壞強度劣化模型研究

邱松楠,黎曉冬

(1.國家開放大學,北京 100039;2.北京工業大學 城市與工程安全減災省部共建教育部重點實驗室,北京 100124;3.南京理工大學 機械工程學院,江蘇 南京 210094)

滲透破壞是水利工程防災防護領域首要考慮的問題。在我國重大涉水項目工程中,時常出現因滲透破壞導致的管涌、地基沉降、邊坡失穩等諸多問題,以雅魯藏布江引水入疆工程為例,由于地理地質環境因素的限制,在砂土地基下進行隧道開挖、水渠建設需考慮滲透沖刷的影響。目前針對滲透破壞引發的土石壩潰壩的研究方法主要包括試驗研究和數值研究。其中試驗研究可根據場地條件劃分為現場試驗和室內試驗,苗發盛等[1]基于滲透-環剪試驗研究了滑坡土在不同滲透條件下強度弱化的性質。李軍紅[2]研究了非均質多孔介質在滲流耦合因素下細觀結構變異特征及宏觀強度劣化規律,并分析了細觀結構變異對土滲透性和強度特性的影響。由于現場試驗在場地選取、費用、時間成本及風險控制等方面存在較大困難,而基于相似理論的室內試驗缺乏可靠的土工試驗機,整體試驗的操作難度高且供排水系統尚不完善。因此,數值模擬基于其成本低、靈活性大、可重復性強的特點,已被廣泛應用于管涌破壞、潰壩致災等水利工程領域。

數值模擬方法中,離散元顆粒流-網格流體耦合模型方法(DEM-CFD)可以較好的描述特定工況下研究區域內的顆粒分布情況,通過迭代計算速度場與壓力場的變化,獲得下一階段顆粒的變化趨勢。顆粒流方法很好的解決了固-固、流-固的多相耦合作用過程,相比于有限元程序顆粒流模擬滲透破壞工程具有很大的優勢。諸多研究結果表明:采用離散元顆粒流-網格流體耦合模型方法能夠模擬從低雷諾流到高雷諾流范圍很大的流體運動,DEM-CFD模型方法模擬多孔介質在不同雷諾數下的流體滲流和顆粒移動過程的結果真實度較高。倪小東等[3]利用顆粒流數值方法分析了管涌發展過程中模型內部細觀變化規律,并于室內砂槽模型試驗結果進行對比,證明了顆粒流程序研究管涌問題的合理性。周健等[4]對比了顆粒流模型試驗與小尺寸細觀模型試驗結果,進一步驗證了顆粒流方法的可行性和合理性。Sterpi[5]通過比對重塑粉砂試樣在不同的時間間隔內涌砂量的變化,建立了涌砂量隨時間和水力梯度的變化關系,并將滲透破壞細顆粒遷移模型與用于無側限滲流分析的有限元程序和應力分析程序相耦合,分析了滲透破壞現象對土體和周圍結構的影響。Maeda等[6]建立了二維離散元模型對高滲透過程的簡化模擬,分析了滲透破壞過程中土體變形和顆粒重分布的相互作用,指出滲透破壞降低了土體強度,導致土體出現畸形應變。Hicher[7]基于微觀結構法分析了顆粒的應力-應變關系,發現在高應力條件下,土體產生較大變形,滲透破壞后的土體內摩擦角所獲得的數值結果與現場堤壩上觀察的破壞模式相一致。此外,離散元計算方法也已經成功應用于砂土的三軸剪切試驗等分析中[8-12],可以有效的解釋滲透破壞發生前后砂土力學性能和細觀結構的變化。

基于以上離散元顆粒流-網格流體耦合模型方法的研究成果,擬建立DEM-CFD滲透破壞數值模型,并通過分析累計涌砂量的變化,推導滲透破壞強度劣化模型,并結合試驗結果與前人研究結果對本文劣化模型的適用性進行驗證。

1 基于DEM-CFD耦合的數值方法

1.1 求解方程

采用顆粒流程序進行滲透過程的研究,首先要確保模擬滲流過程與多孔介質中的流動相適應。在多孔介質中, 流體在雷諾數較低時應符合達西定律, 當雷諾數超過一定值時, 流體流動偏離達西定律, 流速和水力梯度則成非線性關系。利用PFC3D內置FISH語言可定義流固兩相的壓力梯度方程和作用力方程,求解不可壓縮流體中兩相介質的Navier-Stokes連續方程和運動方程。流體的作用力施加在砂土顆粒上,砂土顆粒的運動和孔隙率等也間接影響流體運動的參數及其狀態,二者充分考慮流固耦合作用,具體的求解方法分為液體相方程和固體相方程[13-15]。

1.2 接觸模型

顆粒間的接觸剛度模型為線性接觸模型(The Linear Contact Model),服從廣義胡克定律,為方便計算和模擬理想狀態下純砂土的滲流過程,將模擬過程中固體顆粒間粘結力設置為0,而顆粒間的滑動模型則采用摩擦滑動模型(The Slip Model)。兩種接觸模型概念如圖1所示。

圖1 顆粒切向、法向接觸模型

1.3 迭代程序

具體迭代求解過程如圖2所示。

圖2 顆粒流數值模型迭代程序

飽和砂土體固液耦合反應方程采用PFC3D中CFD(Fixed Coarse-Grid Fluid Model)模型,通過顯式差分方式求解。在i-1時步中,顆粒位置決定顆粒間疊合量,從而確定粒間接觸力的大小,將接觸力與拖曳力相加,計算顆粒加速度,確定i時步顆粒的位移。流體計算基于交錯網格法進行離散計算,流域內的顆粒移動引起土體孔隙率的變化,導致速度場、應力場分布發生改變,流體對顆粒的拖曳力發生改變, 影響顆粒在i時步的位移, 進而引起孔隙率的改變。

2 滲透破壞的DEM-CFD數值模型

2.1 砂體模型

此次模擬選定的顆粒粒徑范圍在0.075 mm～4 mm,其中0.075 mm～0.5 mm的細顆粒含量(FC)為30%,0.5 mm～1 mm的細顆粒含量為0,1 mm～2 mm的細顆粒含量為20%,2 mm～3 mm的細顆粒含量為30%,3 mm～4 mm的細顆粒含量為30%。試樣的尺寸高度80 mm,底面直徑39.1 mm。因為細顆粒數目較多,計算效率普遍較低,因此對流體網格劃分單元。

為了能真實模擬砂樣的實際狀態,砂樣的生成按實際試驗步驟進行。通過FISH函數在給定的模型空間內隨機生成一定密實狀態的基料,循環消除內力;經過足夠長的時間后,砂樣內部的不平衡力基本消散,砂樣處于相對自然沉積狀態。砂樣內部共有20 000個顆粒,如圖3所示。為了增加計算時步,提高計算精度,可放大顆粒半徑,應用Jiang等[10]提出的相似理論。在數值模擬分析中采用相似系數為n=1000,即顆粒放大1 000倍,各直徑的顆粒數目由Jiang等[10]的方法計算。在砂樣飽和后,進行初始耦合計算,使顆粒體達到靜水壓力狀態。施加邊界條件:在砂樣頂部施加100 kPa的豎向壓力,保持底部壓力為0。根據彈性梁在純軸向荷載和純切向荷載作用下的分析,得到顆粒法向和切向剛度數值,具體參數見表1。

表1 數值模型參數

圖3 模擬砂樣

未發生滲透破壞的砂樣(設為對照組)在50 kPa的圍壓下固結,固結后設置顆粒間的摩擦因數為0.5,然后調整圍壓大小為試驗設定數值。在剪切過程中,未發生滲透破壞的砂樣和滲透破壞后的砂樣分別在各自恒定圍壓下,按照恒定1%/s的應變速率進行豎向壓縮。墻和顆粒間的摩擦因數取值用于消除顆粒與墻的邊界效應。

2.2 試驗方案

三維視圖下試樣的示意圖如圖4所示,按照選定粒徑范圍的顆粒組成,四周設為剛性不透水邊界,下砂面采用多孔底板模擬可蝕顆粒自由流出邊界,上砂面模擬試驗室的透水石。在上砂面施加大小為100 kPa的壓力邊界條件,在10 m/s的水力梯度條件下進行DEM-CFD耦合運算,使試樣最終達到穩定滲流狀態。

圖4 模擬試樣三維示意圖

滲流過程中在試樣下砂面向下20 mm位置處設置承載板,用于計算流出顆粒數目,即累計涌砂量n。計算完成后進行三軸剪切試驗過程中,將下砂面的多孔底板替換為墻構件模擬的透水石,滲流過程邊界和三軸剪切邊界三維視圖如圖5所示。詳細試驗方案及參數設置見表2。

表2 試驗方案

圖5 滲流過程邊界和三軸剪切邊界三維視圖

3 結果分析

3.1 不同時刻細顆粒的運移情況

滲流過程邊界和三軸剪切邊界三維圖如圖6所示,圖6展示出土體內部細顆粒的運移情況,在試樣內部多個位置同時出現細顆粒運移現象,滲流后期試樣下部出現裸露的粗顆粒骨架。累計涌砂量隨時間呈拋物線型變化,在滲流時間為0.5 h達到最大。該現象產生的原因是由于試樣細顆粒含量較多且施加較低滲透坡降造成的。由圖6可以看出,試樣上部和底部顆粒的運移情況較中部顆粒更加明顯,試樣中部孔隙基本被完全充滿,只有下游出口附近的細顆粒運移流失后,才能給上游側的細顆粒運移提供空間,但由于大量細顆粒同時出現運移,可能造成局部淤堵,進而導致分層運移的情況。

3.2 滲透破壞前后的力鏈變化情況

滲透破壞后力鏈的變化情況主要取決于細顆粒在土體應力傳遞結構中所起的作用。Thevanayagam等[16]根據土體細顆粒含量大小,將土體應力傳遞方式分為3種情況。第1種情況,當細顆粒含量較少時,細顆粒松散分布在粗顆粒形成的孔隙中,幾乎不參與土體內部應力傳遞;第2種情況,當細顆粒含量較多時,部分細顆粒開始承擔應力,顆粒間形成弱力鏈,但仍以骨架顆粒為主傳遞應力;第3種情況,當細顆粒含量達臨界值,顆粒足夠充滿孔隙時,粗細顆粒同時參與土體內部應力傳遞。上述模擬中CFD試樣屬于情況2,其中部分細顆粒參與應力傳遞,形成了數量眾多、錯綜復雜的弱力鏈,隨著細顆粒不斷流失,該情況下結構破壞方式為弱力鏈逐漸減少,強力鏈失去支撐,進而造成骨架應力傳遞結構失穩。

滲透破壞前后的力鏈變化情況如圖7所示,力鏈中線條的粗細代表接觸力的大小,粗線條表示強力鏈,細線條表示弱力鏈。由圖7中CFD 0可看出,試樣滲透破壞前的強弱力鏈差別不明顯,而滲透破壞后,試樣上部弱力鏈數量減少,同時強力鏈明顯變粗,說明土體內部應力傳遞結構發生了明顯改變。

圖7 滲透破壞前后的力鏈變化情況

3.3 滲透破壞前后的應力-應變情況

滲透破壞后模型應力應變特性與力鏈變化之間存在內在聯系,模型內部的結構性損傷必然引起峰值強度的降低,三軸剪切模擬試驗結果如圖8所示,由圖8可以看出,DEM-CFD耦合運算方法下滲透破壞過程沒有改變砂樣的應力應變線的形狀,不同滲流時間下的砂樣應變軟化性依然明顯。與滲透破壞前模擬砂樣的試驗結果相比,軟化后臨界強度值下降明顯,原因是由于模擬球狀顆粒在三軸應力條件下快速達到臨界狀態時,大量力鏈斷裂,承載能力快速下降,從而導致顆粒滑移,強度降低。與此同時,峰值強度也出現較為明顯的變化,在CFD 0組試驗中峰值強度482 kPa,而在CFD 9組試驗中峰值強度為427 kPa,隨著滲透破壞時間的增加,峰值強度劣化度逐漸增加。

圖8 滲透破壞前后的土體應力-應變關系

4 DEM-CFD耦合劣化模型

4.1 基于累計涌砂量n的劣化模型

引入劣化度DR的概念,用于定量的反映土體力學性質變化規律,劣化度DR為:

(1)

式中:qmax前為滲透破壞前的峰值強度;qmax后為滲透破壞后的峰值強度。

基于上述試驗結果,模擬砂樣的峰值強度劣化度隨滲透破壞時間的增大而增大,且與累計涌砂量呈正比,即單位累計涌砂量損失過程中,強度劣化率為常數。設累計涌砂量產生n克后,砂樣峰值強度為F(n),為可微函數,令未發生滲透破壞的砂樣峰值強度為F(0),則從n克累計涌砂量到(n+Δn)克的峰值強度劣化度為:

(2)

式中:k1為單位累計涌砂量的峰值強度劣化度。

變換得:

(3)

即:

(4)

對式(4)積分可得:

(5)

聯立

(6)

并考慮實際狀況下對峰值強度劣化度的影響差別,引入修正系數α對公式進行修正,則有:

DR=1-αexp(-k1n)

(7)

由式(7)可看出,滲透破壞作用下模擬砂樣峰值強度劣化度與累計涌砂量之間呈指數函數關系。但考慮到劣化度并非完全與累計涌砂量服從正比例函數關系,采用試驗數據擬合的方式對上述DEM-CFD耦合劣化模型進行驗證。

4.2 劣化模型的驗證

Chen等[17]采用細顆粒控制試驗方法,通過使用可溶鹽顆粒代替細顆粒的方法,在控制的應力條件下將鹽溶解在水中,以實現內部滲透過程期間細顆粒的指定損失。侵蝕過程中的徑向和軸向變形均采用照相法測量,并進行了排水三軸壓縮試驗,研究結果表明損失不同數量細顆粒的土體應力應變行為中峰值摩擦角和臨界摩擦角隨著細顆粒的流失而減小。Chen的試驗結果得到了諸多試驗證明和較為廣泛的認可,代入其試驗結果驗證DEM-CFD耦合劣化模型的可靠性。

圖9與表3反映了DEM-CFD耦合劣化模型對Chen試驗結果的擬合程度,表明劣化度與累計涌砂量具有良好的相關性。結果證明,參數α的取值在1附近,參數k1在0.01附近。綜上可知, DEM-CFD耦合劣化模型公式可以在砂柱模型中進行滲透破壞過程計算,相比于雙曲線模型,該模型在峰值流量的計算結果與實測值誤差較小,在整體趨勢和計算精度方面也可得到良好的預測結果。

表3 試驗方案

圖9 計算結果與實測數據對比

5 結 論

(1) 初始細顆粒含量對模擬砂土的影響顯著,會改變土體內部結構和傳力方式。具體表現為試樣滲透破壞前的強弱力鏈差別不明顯,而滲透破壞后,試樣上部弱力鏈數量顯著減少,滲透破壞后土體的峰值強度值和臨界強度值下降明顯,且隨著滲透破壞時間的增加,峰值強度劣化度逐漸增加。

(2) 前述試驗結果和模擬結果充分證明,利用PFC3D中CFD顆粒流數值方法,是數值分析研究滲透破壞較直觀、可靠的方法,模擬結果的相對誤差也較低。

(3) 滲透破壞作用下模擬砂樣峰值強度劣化度與累計涌砂量之間呈指數函數關系,采用試驗數據擬合的方式對上述DEM-CFD耦合劣化模型進行驗證。結果表明:DEM-CFD耦合劣化模型可以在砂柱模型中進行滲透破壞過程計算,在整體趨勢和計算精度方面也可得到良好的預測結果。