含微生物氣泡高飽和砂土循環三軸剪切CFD-DEM模擬

蔣明鏡,沈振義,劉 俊,楊 濤,譚亞飛鷗

(1. 天津大學 建筑工程學院土木系, 天津 300072; 2. 天津大學 水利工程仿真與安全國家重點實驗室, 天津 30007;3. 同濟大學 地下建筑與工程系, 上海 2000921; 4. 同濟大學 土木工程防災重點實驗室, 上海 200092;5. 上海理工大學 環境與建筑學院, 上海 200093)

0 引言

砂土液化是指土體受到振動作用,孔隙水在壓力作用下力圖向上排出,土顆粒受重力影響向下沉落,使土顆粒處于懸浮狀態,導致土體抗剪強度完全消失的現象。砂土液化現象是地震災害的主要形式之一,通常會引起建筑物基礎的不均勻沉降及結構的破壞,甚至造成嚴重災害和人員傷亡。如1975年的海城地震、1976年的唐山地震和2008年的汶川地震都引起了大區域砂土的液化,許多建筑物由于基礎的不均勻沉降遭到不同程度的破壞。目前,減少地基液化主要有換填、強夯地基土的物理處理方法和使用人工化學材料的化學處理方法,但物理方法能量效率低下、化學材料污染環境。因此,探尋新的提高砂土地基抗液化能力的方法具有重要意義。微生物灌漿技術作為一種防止砂土液化的方法,為地基處理提供了新的選擇。

國內外學者對砂土液化的研究,大多集中在室內試驗和數值模擬。在室內試驗方面,Xia和Hu[1]對四組不同飽和度和五組施加不同反壓的試樣進行了循環三軸試驗。結果表明:飽和度對抗液化能力有明顯的影響,反壓也對抗液化能力有影響。但此試驗是通過反壓技術來提高試樣的飽和度,未能準確探究飽和度對液化的影響。Okamura和Soga[2]通過控制飽和度、初始圍壓和初始孔隙壓力壓三個參數,進行了一系列三軸試驗。結果表明:飽和度對抗液化強度有顯著影響。Tsukamoto[3]對粉砂進行了不排水抗剪強度的研究。結果表明:部分飽和和非飽和粉砂的抗液化能力較飽和粉砂明顯提高。He等[4]通過反硝化反應生成的氮氣制取含微生物氣泡的高飽和砂土試樣,進行了振動模擬臺試驗。結果表明:微生物產生的氣泡可以提高土體的抗液化程度。然而,對于微生物氣泡改善土體工程力學性質的微觀機理認識尚淺,有必要對其微觀加固機理進行詳細研究。

針對室內試驗在試樣制取耗時長、成本高和重復性差上的不足,許多學者開始采用數值模擬方法對砂土液化微觀機理進行研究。離散單元法最早由Cundall和Strack[5]提出,其用微觀信息解釋了離散顆粒材料的宏觀力學特性。隨后發展出考慮顆粒粗糙接觸的二維抗轉動模型[6]。目前,顆粒模型進入三維時代,考慮抗轉動和扭轉的三維模型[7]也被提出。計算流體動力學是流體力學的一個分支,運用特定的數值方法和算法求解與流體流動相關的問題[8]。Tsuji等[9]1993年首次將CFD嵌入DEM源程序,模擬固體顆粒中氣體流動行為。結果表明:氣流速度在氣泡出現和壓力的波動頻率改變時與理論解較為吻合。Hoomans等[10]用該方法模擬氣泡的形成,并采用硬球模型代替軟球模型,對氣體和顆粒之間的相互作用采用N-S方程描述,模擬結果與方程求得的結果能很好吻合。Wei等[11]采用離散單元法對不排水循環荷載作用下的砂土液化過程進行了數值模擬,通過對配位數等的變化的分析,探究了砂土液化機理。蔣明鏡和張望城[12]在傳統的CFD流體運動控制方程和大量室內試驗結果的基礎上,引入了考慮流體弱可壓縮性的流體狀態方程,建立了流體密度與壓力時間的關系,利用FISH語言將CFD與PFC2D離散元計算耦合起來,用以對飽和砂土的固結不排水剪切試驗進行模擬。對于微生物巖土技術抗液化能力的試研究,許多學者的研究重點依舊集中于生物膠結,對微生物氣泡的重視程度不足,鮮有微生物氣泡數值模擬的成果。因此,對含有生物氣泡非飽和砂土開展CFD與DEM耦合的液化分析很有必要。

本文基于CFD-DEM耦合程序,運用PFC3D5.0商業軟件,模擬了不同生物氣泡含量的高飽和砂土的固結不排水循環剪切。在模擬過程中,通過跟蹤試樣的力鏈分布、液化振次、孔壓比、軸向應變和配位數在加載過程中的變化情況,從宏微觀角度分析了微生物處理砂土固結不排水循環剪切作用下的試樣力學特性。將模擬結果與室內試驗相對比,探究微生物氣泡對砂土抗液化性能的影響及其作用機理。

1 CFD-DEM流-固耦合計算基本原理

目前,流固耦合的方法主要有3種。筆者采用的CFD-DEM耦合方法考慮了不排水試驗中的邊界移動、含氣泡流體的可壓縮性、以及土體試樣在循環荷載作用下加卸載剛度的不同。此方法適用于微生物氣泡引起流體壓縮系數變化的不排水循環三軸模擬。

1.1 顆粒運動方程

顆粒運動方程與經典DEM方法一致,均遵循牛頓第二定律:

(1)

(2)

1.2 流體運動控制方程

(1) 連續性方程

根據質量守恒定律推導得出的流動控制方程稱為連續方程。單位時間內,通過體積單元的質量凈通量等于體積元內質量的變化量,即

(3)

式中:εf為一個流體計算單元中的孔隙率;ρf為流體密度,▽為拉普拉斯算子;uf和ug分別為流體運動速度和顆粒運動速度。

(2) 動量方程

動量方程遵循牛頓第二定律,其代表的物理意義是作用于控制體上力的總和等于控制體的質量乘以控制體運動時的加速度。其可由公式(4)表達:

-εf▽P+▽·(εfμf▽uf)

(4)

式中:P為流體壓力;μf為流體動力粘滯系數。

(3) 能量方程

將熱力學第一原理應用于流體模型的控制方程稱為能量方程。其代表的物理意義是控制體內的能量變化率等于流入控制體內的凈熱流量加上體積力和表面力對控制體做功的功率?？捎晒?5)表達:

▽·(Puf)+εfρfFuf

(5)

(4) 液體狀態方程

在不考慮溫度和壓力耦合的情況下,能量方程是滿足的。而只有動量方程和連續性方程是無法閉合求解的。為了形成一個封閉的方程組,本文引入了流體的狀態方程。一般情況下,水的可壓縮性很低,可以忽略不計。而對于含氣泡的非飽和液體,其可壓縮性就會隨著氣泡含量的增加而提高。根據液體增壓前后質量保持不變的性質,得到如下液體壓縮系數的表達式:

(6)

本文通過液體的壓縮系數建立起CFD的狀態方程,如下:

ρf-ρf0=ρfC(p-p0)

(7)

式中:ρf0為流體初始密度;p0為流體初始壓力。

微視頻也是隨著網絡發展而產生的。我們都知道現在社會的生活節奏是相當的快，人們一天幾乎沒有很好的時間去享受自己的空閑時間，也沒有時間去利用自己的娛樂時間。長的視頻需要花大量的時間去看，而且看過之后心里也總是在想著劇情，容易擾亂自己的思想，不能專心工作。也正是因為這樣的情況，我們的段視頻應運而生，它不需要花費很多的時間，簡單地幾分鐘就可以瀏覽完[4]。那我們的商務英語翻譯課程就可以制作成簡短的視頻，如果內容長，我們就可以做成相關系列的東西，便于人們利用閑碎的時間去學習，只需簡短的幾分鐘，便可以輕松搞定。再有閑暇時間還可以拿出來重新，溫習一下，讓知識得到鞏固的效果。這就是利用做簡單的微視頻的好處。

1.3 CFD網格移動與網格內孔隙比計算

在模擬過程中,試樣顆粒被六道剛性墻包圍,上下墻的速度保持不變,側墻的速度根據伺服系統得到。因此,CFD邊界條件采用速度邊界,即靠近邊界的液體速度與墻的速度一致。加載過程中,由于試樣形態發生變化,CFD的邊界網格需要與墻同時移動,內部網格則采用簡單的均勻變形方式移動,以x方向為例,內部網格的運動控制方程為:

(8)

1.4 耦合軟件及流程

本文DEM計算通過PFC5.0進行,安裝在Windows系統下,采用FISH語言和Python語言混合編程。而CFD計算通過OpenFOAM進行,安裝在Linux系統下,采用Python語言編程。兩者通過基于TCP的sockets編程實現CFD-DEM的數據交換。耦合計算流程示意圖如圖1所示。

圖1 CFD-DEM耦合流程示意圖Fig.1 Process schematic of CFD-DEM coupling

2 微生物處理過程的時效性

筆者采用文獻[16]中的氣泡生成方法,制備得到微生物處理2天、3天和5天后的含不同氣泡含量的高飽和砂土試樣。在將試樣中的氣泡體積與飽和度轉換為數值模擬中的流體壓縮系數時,不考慮水表面張力的影響,采用劉侃等[17]提出的含氣孔隙流體壓縮系數的表達式,如下所示:

(9)

式中:下標0代表初始狀態;Va為氣泡體積;Vw為水的體積;uw為孔隙水壓;Pa為標準大氣壓;ua為氣泡壓力;Sr為飽和度;C0為飽和情況下的壓縮系數。最終的參數計算結果如表1所示。

3 循環三軸CFD-DEM耦合模擬

為研究砂土抗液化能力,進行砂土固結不排水循環三軸室內試驗時,制取試樣后,施加圍壓使試樣處于等應力狀態,通過激振系統在豎向施加循環荷載,直至試樣達到破壞狀態。但室內試驗無法獲取試樣的力鏈分布、配位數等微觀數據,數值模擬可以直觀試樣破壞的整個過程,便于探究砂土液化的宏微觀機理。

表1 數值模擬試樣微生物處理過程時效性的參數

3.1 模擬步驟

含微生物氣泡的高飽和砂土循環三軸模擬分成樣、預壓、固結、循環剪切四個過程。前三個過程采用DEM完成,最后的循環剪切過程采用CFD-DEM耦合。

(1) 成樣。為與室內試驗[18]對比,本文離散元模擬對象為渥太華50-70砂,級配采用Ottawa 50-70砂級配,初始孔隙比為0.73。模型微觀參數如表2所列,離散元試樣的制備采用分層欠壓法[19],分10層制備,每層2 000顆粒。

表2 模型微觀參數

(2) 預壓。制得均勻試樣后,固定四道側墻,通過伺服系統對上下墻施加大小為12.5 kPa的應力進行預壓。

(3) 固結。預壓穩定后,通過伺服系統對試樣進行等向固結。固結試樣的平均有效主應力p=100 kPa,偏應力q=0。

(4) 循環剪切。此階段中加入了CFD模塊進行耦合模擬,CFD參數借鑒文獻[20]。剪切時DEM時步被固定為1×10-6s/cycle,DEM與CFD的交換頻率為25(即DEM計算25次,CFD計算1次),CFD時步為2.5×10-5s/cycle。液體的初始密度為1 000 kg/m3,壓縮系數為4.5×10-7。網格數目為6×6×12 (長×寬×高),網格邊界與DEM中六道墻體的坐標一致。循環加載前施加100 kPa反壓。通過上下墻加載進行剪切,當偏應力達到應力幅值qcyc時,上下墻反向繼續運動,當偏應力在相反方向達到應力幅值qcyc則再次反向。側墻仍然采用伺服系統控制,伺服應力為總圍壓減去孔壓后的有效應力。

3.2 模擬方案

為研究不同氣泡含量的高飽和砂土的抗液化強度,選取了微生物處理2天、3天及5天的高飽和試樣進行對比,詳細工況見表3。其中,循環應力比為循環荷載作用下的最小主應力與最大主應力的比值,其表征平均應力與偏應力的聯合作用。為方便描述,對每個試樣進行了編號,例如,US2代表微生物處理2天的砂土試樣,US2-1/2/3/4代表對微生物處理2天砂土試樣的不同循環應力比循環三軸試驗。

表3 微生物處理高飽和砂土循環三軸數值模擬工況列表

4 DEM模擬結果

4.1 力鏈分布

圖2為US2-3試樣在循環剪切過程中特征時刻的力鏈演化。圖中力鏈越粗代表顆粒接觸力越大,力鏈顏色同時也代表接觸力的大小,紅色代表的接觸力最大,接下來依次是黃色、綠色,藍色代表接觸力最小。如圖所示,隨著加載的進行,試樣顆粒接觸力越來越小。試樣經過等向固結后,內部力鏈分布較為均勻。試樣在受壓時,豎向受力為主導,內部力鏈主要沿豎向分布。試樣在受拉狀態時,側向受力為主導,內部力鏈主要沿水平向分布。模擬結果表明:在加載過程中,孔隙水壓力累積,承受了很大一部分力的作用,導致土骨架的有效傳力減小。當試樣中的有效應力為零時,便達到了初始液化狀態。試樣內部力鏈的分布方向隨拉壓狀態的轉換而轉換,沿主要受力方向分布。

圖2 試樣US2-3加載過程中特征時刻力鏈分布演化Fig.2 The force chains distribution of sample US2-3 during loading

4.2 抗液化強度

微生物處理不同天數的高飽和砂土的抗液化強度曲線如圖3(a)所示。圖 3(b)(c)為室內試驗[1,21]中高飽和砂土循環三軸試驗抗液化強度曲線。圖2(a)中,循環應力比CSR=0.25時,微生物處理2天的試樣在2次振動后即達到液化狀態,處理3天的試樣達到液化狀態需要振動38次,而處理5天的試樣需要87次振動達到液化狀態。模擬結果表明:含微生物氣泡的高飽和砂土的抗液化強度比處理前的飽和凈砂有所提升。且處理天數越多,氣泡含量越多,試樣的抗液化強度越高。模擬結果與室內試驗相吻合。

圖3 微生物處理高飽和砂土試樣循環應力比-液化振次曲線Fig.3 Cyclic stress ratio-liquefaction vibration times curves of highly-saturated sand treated by microorganisms

4.3 孔壓比

孔壓比是指在循環應力作用下試樣的孔隙水壓力增量與側向有效固結應力的比值。圖4為循環應力比CSR=0.25時不同氣泡含量的試樣孔壓比發展曲線,同時給出飽和試樣在CSR=0.15下的孔壓比。如圖所示,微生處理后的試樣的CSR比飽和試樣大,含氣泡的試樣孔壓比曲線發展更為緩慢。飽和試樣和微生物處理2天后的試樣的孔壓比在加載后都迅速上升至1,達到初始液化狀態。試樣經過微生物處理后的第3天和第5天都沒有達到初始液化狀態,其最終的孔壓比分別在0.5和0.4附近波動,且達到累積應變破壞所需的加載時間大幅度增加。模擬結果表明:隨著微生物處理的天數增加,試樣的飽和度不斷減小,孔壓比的累積隨之變弱。

圖4 不同氣泡含量試樣孔壓比發展曲線Fig.4 Pore pressure ratio development curves of samples containing different bubbles

4.4 軸向應變

圖5所示為循環應力比為CSR=0.25時的不同微生物氣泡含量的試樣軸向應變發展曲線,圖中還給出了飽和試樣在CSR=0.15下的軸向應變曲線。從圖中可以看出,微生物處理2天后的試樣的軸向應變發展曲線與飽和試樣較類似,試樣達到初始液化狀態時的最終應變都在2%以內。微生物處理3天和5天后的試樣軸向應變不斷累積增長,直至達到5%時發生累積應變破壞。模擬結果表明:隨著微生物處理時間的增加,試樣的飽和度逐漸降低,應變累積隨之變慢。

圖5 不同氣泡含量試樣軸向應變發展曲線Fig.5 Axial strain development curves of samples with different bubble contents

4.5 力學配位數

力學配位數是指單個顆粒所含有接觸力的數目,通常要求每個接觸的法向接觸力大于零。其曲線可以表征加載過程中試樣的微觀變化。圖6為循環應力比CSR=0.25時不同微生物氣泡含量的試樣力學配位數發展曲線,同時給出了飽和試樣在CSR=0.15下的力學配位數曲線。從圖中可以看出,加載前各試樣的力學配位數都是4。微生物處理2天后的試樣與飽和試樣的力學配位數發展曲線較為相似,達到初始液化狀態時力學配位數都下降至3.5附近。而微生物處理3天和5天后的試樣不發生初始液化破壞,它們的力學配位數發展的整體規律為先上升后下降。結果表明:經微生物處理后,試樣的力學配位數發展曲線發生了變化,且當氣泡含量達到一定值后,試樣發生累積應變破壞。

圖6 不同氣泡含量試樣力學配位數發展曲線Fig.6 Mechanical coordination number development curves of samples containing different bubbles

5 結論

本文利用CFD-DEM耦合方法模擬含微生物氣泡的高飽和砂土的不排水循環三軸試驗,對比室內試驗,分析了試樣的抗液化能力,得出以下結論:

(1) 含微生物氣泡的高飽和砂土的抗液化強度較飽和凈砂高。且試樣經微生處理的天數越多,氣泡含量越多,抗液化強度越高。

(2) 高飽和砂土試樣的孔壓比和軸向應變隨著微生物氣泡含量的變化而變化,且隨著氣泡含量的增加,試樣的孔壓比的累積隨之變弱,應變累積隨之變慢,抗液化能力增強。

(3) 經微生物處理后,高飽和砂土試樣的力學配位數發生了變化,當氣泡含量達到一定值后,試樣發生的是累積應變破壞。