二維固結過程中土體各向異性影響的研究★

裴美娟，楊 鵬，李樂晗

(蘭州石化職業技術大學土木工程學院，甘肅 蘭州 730060)

1 概述

土體的固結過程是巖土工程中重要的研究課題，對于實際工程而言，土體固結過程伴隨著土體體積的變化，進而導致地基的沉降，會對工程建設造成較大的影響。再者，巖土體具有多相性質，土顆粒之間的孔隙往往被水所填充，而在外力作用下，土體會產生超孔隙水壓力，超孔隙水壓力的產生會對土體的強度造成強烈的影響，會導致土體發生剪切破壞，而太沙基的滲流固結原理指出了固結排水過程中孔壓的變化規律，可以說固結過程不僅僅體現了土體體積的變化，而且更為重要的是能夠描述排水過程中土顆粒之間有效應力的變化規律。

土體的固結過程受到了眾多學者的關注與研究，不同土體性質有較大差異，其固結排水特性不同，由于土體滲透性具有較多的影響因素，因此，不同因素下，其滲透、固結性質會表現出差異。鄧岳保等[1]研究了溫度效應下土體的固結蠕動耦合行為，張玉國等[2]給出了外荷載變化條件下組合樁復合地基的固結計算方法，指出不同邊界條件會導致孔壓消散的速率及分布出現差異。王潔等[3]建立了豎井的非線性固結模型，并引入了井阻變化的影響。當土體處于欠固結狀態時，其固結特性與正常固結土體不同[4]，張志龍等[5-6]基于分段線性方法分析了欠固結土體的大變形非線性固結過程，發現該過程中土體的孔壓消散具有明顯的非線性性質。賀建清等[7]則研究了有機質土體的次固結行為，有機質含量會對固結特性產生影響，有機質的增加會造成次固結系數的上升。循環荷載作用下，汪磊等[8]研究了循環荷載下非飽和土的固結過程，并采用變換方法給出了解析解。土體性質是具有復雜性的，不同因素影響下都具有不同的特性，因此，一些學者也開展了復雜條件下的固結行為研究[9-11]，這些研究為工程應用提供了技術支持。

目前，對于固結過程中土體各向異性的研究較少，黃土受到成因、氣候等因素的影響，具有強烈的豎向節理，豎向滲透性往往比水平滲透性大得多，一維滲流固結理論無法描述出二維空間中的各向異性表現，多數的研究基于各向同性假設開展，無法準確刻畫出黃土的固結過程，因此，開展各向異性下二維固結理論的研究，對于評價黃土等特殊土的固結特性具有重要的意義。本文通過Galerkin法進行控制方程的離散，獲得其有限元形式，并基于Python語言編程開發專有的有限元求解器，完成各向異性下的二維固結計算，能夠有效反映出不同各向異性程度下超孔壓的變化，并對比分析了不同異性程度下及排水邊界條件二維固結模型的孔壓消散規律以及分布特性。

2 控制方程的離散

2.1 控制方程及邊界條件

在二維空間中，考慮豎向滲透系數kv與水平滲透系數kh不同，同時理論模型依然滿足水流連續性原理，孔隙水的排出依然滿足達西定律，這樣豎向滲流速度與水平滲流速度可以表示為：

(1)

(2)

其中，u(x,y,t)為孔隙水壓力函數，是關于空間坐標和時間的函數；γw為水的重度；vh，vv分別為水平向滲流速度和豎向滲流速度。

而土體體積的變化完全由孔隙水的排出量決定，即有式(3)：

(3)

其中，εv為土體體積應變。

引入土體體積壓縮系數,同時由太沙基有效應力原理得到二維滲流固結控制方程,如式(4)所示。

(4)

其中，Ch為水平固結系數,Ch=kh/mvγw;Cv為豎向固結系數,Cv=kv/mvγw;mv為土體體積壓縮系數。

方程(4)為各向異性條件下的二維滲流固結控制方程,該方程為標準的二維拋物型偏微分方程,屬于非穩態物理場問題,需要在邊界條件以及初始條件給出的基礎上,才能夠完成求解。

2.2 控制方程的離散

為了對控制方程進行求解,考慮到空間維度的復雜性,采用有限元方法進行偏微分方程的求解。基于Galerkin法進行控制方程的離散,獲得其有限元形式。首先,單元選擇三節點三角形單元,此時,單元內任意點處的場函數值可以通過節點值進行插值表達,見式(5)。

u(x,y,t)=[N]·[u(t)]

(5)

其中，[N]為形函數行向量;[u(t)]為t時刻單元3個節點處的函數值,為列向量。

針對方程(4),構建其加權殘差,得到對應的等效積分強形式,見式(6)。

(6)

利用分部積分思路或者格林公式進行式(6)的處理,可以進一步得到等效積分弱形式,見式(7)。

(7)

由于形函數是空間坐標的一次函數,因此求導后為常數,僅與單元節點的坐標值有關,同時,將式(5)代入到式(7)后,可以得到空間離散后的形式,見式(8)。

(8)

由于關于空間的積分是可以直接進行計算,這樣將空間項積分記作[B]矩陣,而含有時間項的積分記作[S]矩陣,則得到空間離散后的形式，見式(9)。

(9)

由于該問題屬于非穩態問題,利用后向差分進行時間項的離散,得到式(10)。

(10)

式(10)即為離散后的有限元形式，為了實現數值求解，基于Python語言開發數值模塊完成計算，采用面向對象的編程范式，構建控制方程實例的類，而單元剛度矩陣的計算以及整體剛度矩陣的集成作為類的方法進行設計，邊界條件、初始條件作為屬性進行輸入，邊界條件的引入在每一次迭代計算時進行更新。

3 數值結果分析

3.1 算例設置及參數

為了對比各向異性對于土體固結的影響，設置不同組別下的參數組，具體信息見表1。控制豎向固結系數保持不變，而水平向固結系數為豎向固結系數乘以比例因子，進行各項異性的表達，其中參數組別A中，豎向固結系數等于水平向固結系數，為各向同性，而其他組別，隨著比值r的降低，水平向固結系數是隨之降低的，典型的黃土含有大量的豎向裂縫，豎向滲透性較強，而水平滲透性低。

表1 固結系數取值

模型長高為5 m×3 m，如圖1所示。模型具有4個邊界，其中邊界1和邊界2為不透水邊界，為第二類齊次邊界條件，頂部地面與大氣連通，孔壓始終為0，邊界4為第一類齊次邊界條件。為了考慮邊界條件對于各向異性的影響，邊界3分別考慮為不透水和透水兩種情況，即模型1中邊界3為不透水邊界，而模型2中邊界3與外界連通，孔壓為0。采用自編程開發的有限元程序進行計算，網格劃分如圖2所示，模型節點個數為1 834，模型單元數量為3 506，孔壓場為標量場，節點處自由度為1。計算總時間為1e4 s，時間步數為1 000步，模型初始孔隙水壓力為20 kPa。

3.2 結果分析與對比

圖3為模型1在參數A下末時刻的孔隙水壓力分布云圖，在該邊界條件下，當水平向固結系數降低到參數E時，其孔壓分布云圖見圖4，對比發現，兩者區別較小，且分布規律較為一致。在模型1中，由于左側邊界為不透水邊界，而滲透系數在水平向的降低，主要影響到水平向滲透性的降低，但是由于邊界的限制，兩側邊界不透水，因此，在該種條件下土體滲透性的各向異性對于固結過程的影響是較小的。圖5為參數A結果與參數E結果的偏差孔壓分布圖，可以看到，各向異性條件下偏差孔壓的幅值較小，在0.02 kPa之內。圖6為模型在5組參數下的平均孔隙水壓力變化曲線，可以看到各向異性在此種邊界條件下，對于模型孔壓消散的影響較小，與各向同性下的計算結果偏差不大，平均孔壓曲線基本重合。

當邊界3為透水邊界，其邊界條件為第一類齊次邊界，圖7為模型2在參數A條件下第500時間步的孔隙水壓力云圖，圖8為在參數E下的孔壓分布，可以看出，隨著水平向滲透性的降低，孔壓的消散速度也變慢，圖8孔壓幅值要高于圖7的，反映出邊界條件會對各向異性的描述產生一定程度的影響。

為了對比該邊界條件下，水平向固結系數降低所帶來的影響，分別給出了模型2的參數A與參數E的第100步、第500步及結束時的偏差孔壓分布云圖，分別見圖9～圖11。在第100步時，固結系數的降低使得模型整體孔隙水壓力偏大，并且在遠離排水邊界的位置幅值較大，最大達到了6 kPa。隨著固結過程的繼續發生，孔隙水的持續排出，孔壓偏差較大的位置逐步向中間位置移動，而當達到末時刻時，偏差孔壓主要集中在右側邊界部分。造成這一種變化規律，主要是受到了右側邊界性質的影響，并且土體水平向固結系數的降低，表示土體水平方向排水能力的降低，在邊界條件一定的情況下，土體的各向異性性質會造成固結過程產生差異，尤其是在邊界透水情況下，這一差異更為明顯。圖12給出了在不同參數組下模型固結過程中平均孔隙水壓力隨時間變化曲線，可以看出，隨著水平向固結系數的降低，模型的排水速度逐漸降低，表明在此時邊界條件下，各向異性的性質對于模型的整體的孔壓分布是有明顯影響的。

綜上，對于二維空間中的土體固結過程而言，土體的各向異性是會對固結過程中的孔壓變化產生一定程度影響的，但是還需要考慮邊界條件的影響，如果該方向邊界不排水，那么該方向的滲透性對于孔壓消散的影響較小。當兩個方向均含有排水邊界時，各向異性表現得較為明顯，與各向同性的結果相比，孔壓偏差主要出現在遠離排水邊界位置處，并且伴隨著固結的進行逐步向排水邊界發展。

4 結論

為了探究各向異性條件下的土體固結特性研究，針對二維固結模型開展了研究，基于Galerkin法進行了控制方程的離散，并利用Python語言進行了有限元程序的開發，通過自編程的程序計算分析了不同異性程度下二維模型的固結過程，并與各向同性的結果進行了對比，主要研究結論如下：

1)各向異性條件下，土體固結過程中超孔隙水壓力的變化規律與邊界條件有關，當該方向邊界為不排水時，該方向滲透性(固結系數)的變化對于孔壓的分布影響較小，且模型整體平均孔壓變化較小。

2)當兩個方向均有排水邊界時，土體各向異性對于固結過程影響明顯，與各向同性結果有明顯區別，偏差孔壓主要出現在不排水邊界附近，并隨著固結的進行逐漸向排水邊界側移動。