毛細水作用下非飽和土壓縮過程的微觀非連續(xù)變形數值分析

李 強 ，李同錄 ，李 華 ，沈 偉 ，李 萍 ，張常亮

（1.長安大學地質工程與測繪學院, 陜西 西安 710054；2.黃土高原水循環(huán)與地質環(huán)境教育部野外科學觀測研究站, 甘肅 正寧 745399；3.中國人民武裝警察部隊工程大學裝備管理與保障學院, 陜西 西安 710086）

微結構演化是土體變形破壞的本質，掌握其動態(tài)演化規(guī)律是深入認識土體宏觀變形破壞機制的關鍵。壓汞技術（MIP）、掃描電鏡（SEM）、CT掃描系列是當前研究土體微觀變形的主要實驗技術。一些學者利用MIP或SEM對荷載作用前后的黏土[1]、混合黏土[2]和黃土[3-4]觀測發(fā)現，荷載主要改變較大的粒間孔隙。但這2種觀測技術對土樣都有擾動，不能追蹤觀測荷載作用前后某一位置微結構的具體變化。Yu等[5]結合μ-CT，嘗試設計出了微型土樣加載裝置，觀察土樣變形過程中顆粒的位置變化。但觀察時需暫停加載，將土樣取出，因此對土樣造成較大擾動。同樣的問題也存在于MIP和SEM中。普通CT觀測室較大，可以在加載過程中同步觀察微結構的變化。利用普通CT對黃土[6]、膨脹土[7]等研究發(fā)現，隨著荷載的增加，土體孔隙逐漸分布均勻，不規(guī)則大孔隙演變成圓形小孔隙，小孔隙不壓縮。但普通CT分辨率相對較低，很難捕獲到清晰的微結構特征。綜合來看，前人對土體變形前后的微結構變化研究較多，認識較深，但對變形過程中微結構演化尚缺乏直接的觀測手段，現有的微觀試驗技術很難克服這一問題。因此，將數值方法引入土的微觀結構分析，是目前較為現實的一種途徑。

土在微觀上是由顆粒構成的離散體系，可以用離散介質力學方法研究其微觀力學性質。非連續(xù)變形分析方法[8]（Discontinuous Deformation Analysis，DDA）是常用的離散介質力學分析方法之一。該方法起初是用于宏觀塊體的分析，例如塊狀結構巖體的滑坡穩(wěn)定性[9]、圍巖變形與破壞[10]和壩基滑移穩(wěn)定性等[11]。近幾年，一些學者初步嘗試將DDA擴展到分析土體微觀力學性質中。例如張國新等[12]利用規(guī)則的多邊形塊體構建砂土DDA 模型，研究了砂土應力-應變關系；郭培璽等[13]將多邊形顆粒隨機填充，構建了粗粒料的 DDA 模型，研究粗粒料的力學特性；郭龍驍等[14]通過人為構建干燥黃土DDA模型，模擬了黃土的單向固結試驗。然而這些研究都基于原有DDA模型，只考慮了粒間摩擦、咬合等作用，并沒有考慮毛細作用的影響。對于非飽和土，毛細作用是基質吸力的主要來源，也是工程問題關注的重點。巖土工程領域也提出一些考慮毛細作用的離散介質分析方法[15-16]，但這些方法一般在相同大小球形或圓形顆粒之間添加一個等大的或者與距離有關的應力表征毛細作用，這很難有效地考慮毛細水分布狀態(tài)對非飽和土力學性質的影響。

為此，李強等[17]提出了考慮粒間毛細水作用的DDA方法，能夠計算不同含水率下的粒間毛細水分布和相應毛細力，可用于模擬靜態(tài)條件下非飽和土的持水特性。在此基礎上，本文考慮進一步擴展該分析方法，將其用于研究荷載作用下非飽和土微結構的動態(tài)演化機制。

本文參考黃土骨架粒的大小及形態(tài)建立了理想的非飽和土骨架結構模型，并設置自動加載系統，模擬了非飽和土側限壓縮試驗。根據模擬結果，分析了壓縮過程中孔隙比、吸力和飽和度的變化規(guī)律；通過追蹤孔隙的演變得到了不同含水率下孔隙的變形特征。

1 土粒間毛細水算法

1.1 算法思路

在非飽和土中，當相互距離較近的顆粒被毛細水黏連起來，將會形成一個相互作用的系統，此系統中包含毛細吸力和表面張力2 個作用力[18]。所以在顆粒間施加毛細力時，關鍵在于如何將這2個力施加于土粒之間。

為提出該算法，先作以下假設：①只考慮粒間毛細力；②忽略不同顆粒土-水的接觸角差異，即所有接觸角都相同；③假定體系處于平衡狀態(tài)，即模型內部沒有水力梯度；④土體中的顆粒均為凸多邊形。

算法的思路是：先給定含水率ww，假定該含水率下土顆粒之間的毛細水彎液面半徑為rw，基于假設①～③，所有彎液面半徑都相同，若求得該半徑，根據Young-Laplace方程，可計算毛細吸力為：

式中：ua——孔隙氣壓力/kPa；

uw——孔隙水壓力/kPa；

Ts——表面張力/（mN·m-1）；

R——彎液面半徑/μm；

α——接觸角/（°）。

根據毛細水在顆粒上的浸潤面積，計算作用在顆粒表面毛細吸力的合力。將該合力與表面張力疊加，嵌入DDA的運動方程組中。

計算某一確定含水率的彎液面半徑rw時可用黃金分割法逼近。首先假定彎液面初始半徑區(qū)間為[a,b]，利用圓心軌跡交匯法分別確定a和b對應的所有土粒之間的彎液面。彎液面和所連接的土粒邊界構成封閉的毛細水分布區(qū)。統計毛細水分布區(qū)的面積，得到與a和b分別對應的毛細水含水率wa和wb。若wb＞ww＞wa，可采用黃金分割法逼近ww。若wa＞ww或wb＜ww，調整半徑區(qū)間[a，b]的范圍，重新求解。最終得出與給定含水率ww對應的毛細水彎液面半徑rw。

該方法構建的毛細水分布在顆粒接觸或接近的區(qū)域，毛細水呈凹透鏡狀，見圖1。在這個透鏡形區(qū)域內毛細水產生毛細吸力(ua-uw)和表面張力Ts2 個作用力。在二維平面上，將這2個力對土顆粒的作用進行簡化。其中，作用在土粒邊界上的毛細吸力合力按其作用段分段計算，第i段毛細吸力的作用長度為Li，其合力fsi為：

圖1 單顆粒上毛細作用力計算Fig.1 Calculation of the capillary force on the single particle

該力的作用點位于浸濕段中點，方向指向該段的外法線。表面張力Ts僅作用于固、液、氣三相交界點處，在溫度等條件確定的情況下，其值為常量，方向沿彎液面的切線由兩端指向中間。而后分別將2 個力在笛卡爾坐標系下分解到x、y兩個方向：(Tsx，Tsy), (fsix,fsiy)，將分解后的力分別以顆粒為單元求和，再嵌入到DDA的運動方程中。

先利用理想模型驗證毛細水算法的有效性。構建2 個半徑同為10 μm的圓盤，上面圓盤固定，下面圓盤可移動，初始狀態(tài)下2 圓盤不接觸，見圖2。分別用理論解和以上算法計算不同間距下的粒間毛細力。

圖2 不同水量下兩圓盤模型中的毛細水分布Fig.2 Capillary water in the double-disk model at different water content

不考慮圓盤的重力。在二維平面上，當粒間間距為d時，計算不同填充角δ對應的毛細水彎液面半徑r[15]：

單位厚度水量Vw為：

粒間毛細力(ua-uw)為：

利用式（3）—（5）計算得，間距d為0.50，0.25，0.05 μm的粒間毛細力，計算結果見圖3。

圖3 不同間距下的毛細力理論和模擬計算結果Fig.3 Capillary pressure calculated with theory and simulation

毛細力作用的數值模型中，顆粒的重度、接觸角與理論計算一致，均設為0，模型所需的其余參數（密度、彈性模量、泊松比、內摩擦角、黏聚力）與表1中的土顆粒一致。該方法計算結果見圖3。理論解與模型結果吻合很好，說明毛細作用模型是可靠的，還可將其進一步擴展到土體的微結構演化分析中。

2 非飽和土壓縮過程模擬

2.1 骨架結構生成

參考黃土骨架顆粒的形態(tài)和大小，建立理想的具有大孔隙結構的非飽和土體微觀骨架模型。圖4為用于構建模型的實測黃土粒度分布曲線。由圖4可知，該土樣的骨架大多以粉粒為主，這部分顆粒粒徑集中分布在20～30 μm，因此本文只選取該范圍粗顆粒建立理想土骨架模型。

圖4 樣品粒度分布曲線Fig.4 Particle size distribution curve

通過繪制放大100倍的黃土粗顆粒形狀，并用多邊形將其擬合，建立黃土粗顆粒形態(tài)庫。利用落體投放法建立模型，具體步驟為：①建立內徑為540 μm×400 μm，盒壁為50 μm的樣盒模型；②從顆粒庫中隨機抽取350個在20～30 μm的顆粒，并將其隨機掉落于540 μm×400 μm區(qū)域內，形成的微觀結構模型見圖5（a）。

計算三維模型的孔隙率n3d和孔隙比e：

該模型的初始二維孔隙率n2d為28.0%，孔隙比e為0.390。利用式（6）（7）可以估算出三維孔隙率n2d和孔隙比e分別為48.0%和0.910。

2.2 參數確定

土顆粒和樣盒參數見表1。土顆粒一般是物源區(qū)巖石風化后的產物，所以密度、彈性模量、泊松比等基礎物理力學參數可以參考巖石中（花崗巖）的相應參數。本文旨在探討毛細作用對土體的影響，沒有考慮細粒的作用，黏聚力設為0 kPa。摩擦角根據直剪實驗確定。

表1 模型參數表Table 1 Parameters in the model

樣盒對土樣起到約束作用，因此為剛性體，物理參數參考銅的相應參數，如表1。此外，為避免盒壁與土粒的摩擦等非土顆粒間作用力對土體變形的影響，將樣盒表面設置為光滑不浸濕，即分別設置內摩擦角、黏聚力為0°，0 kPa，接觸角為179°，表面張力為0 mN/m。

2.3 加載系統的構建

骨架模型建立以后，進一步施加毛細作用。構建的模型土樣初始含水率為5%、10%、18%。利用毛細水算法確定模型中相應含水率下的毛細水面積及分布范圍，并計算出相應毛細力，施加于土顆粒之間。各含水率下毛細水在土中的分布狀態(tài)如圖5（b）—（d），對應的飽和度分別為18.0%、36.0%、68.0%。

加載控制系統的實現方法為：在上述模型上部設置一加載板，加載板向土樣逐級加載。加載板與樣盒的特性參數、厚度、親水性均相同，見表1。加載板與樣盒之間為光滑接觸，只能在豎直方向運動。在加載板上均勻布置4個荷載點，荷載點的間隔為加載板寬度的1/5，用于給土樣加載，如圖5所示。在加載板上布置3個監(jiān)測點記錄監(jiān)測點垂直位移，監(jiān)測點間隔為1/4加載板寬度。

加載控制系統的關鍵在于如何確定施加下一級荷載的時間節(jié)點。實現方法為：首先賦予系統預計施加的各級荷載值，荷載值可任意給定，但需按加載順序逐級遞增。開始加載后，監(jiān)測點反饋位移信息。在當級荷載下，若監(jiān)測點反饋的位移量穩(wěn)定（時步增加100步，位移量小于0.01 μm）且不是最終荷載，則自動施加下一級荷載；若位移未穩(wěn)定，維持當級荷載，直至穩(wěn)定；當位移穩(wěn)定且為最后一級荷載，結束壓縮過程。每級荷載作用穩(wěn)定后，自動輸出當級荷載下的土樣變形量。本次試驗的荷載序列為1 kPa→6.25 kPa→12.5 kPa→25 kPa→50 kPa→100 kPa→200 kPa→400 kPa→800 kPa。

在土顆粒之間增加毛細作用力，結合加載系統，模擬非飽和土側限壓縮試驗。在壓縮過程中，土顆粒將重新排列，孔隙結構變形，毛細水重新分布。在每一時步，模型將重新計算毛細水彎液面半徑及毛細力，并更新荷載矩陣，而后進入下一增量步。同時，模型將實時輸出圖像，通過圖像可直觀地觀察孔隙結構形態(tài)及毛細水分布狀態(tài)。圖6為整個數值試驗的計算流程。

圖6 非飽和土壓縮試驗模擬的計算流程Fig.6 Flow chart of simulation for the unsaturated soil compression test

3 模擬結果分析

3.1 e-lgp壓縮曲線

實驗過程中測量點的豎向位移即模型土的豎向變形S轉換為二維孔隙比e：

式中：e0——初始孔隙比；

H——初始模型高度/μm。

根據式（6）（7）估算出三維孔隙比，三維孔隙比e與豎向荷載p之間的關系曲線（e-lgp）如圖7。

圖7 模擬壓縮曲線Fig.7 Simulated compression curves

從圖7的結果來看，不同含水率下的e-lgp曲線可分為平緩和陡降2個階段，且含水率高的試樣，同一荷載下的孔隙比小。平緩段和陡降段的轉折點對應的荷載被稱為側限壓縮屈服應力（Psc）[19]，模擬得到的Psc隨含水率的增大而減小。這些規(guī)律說明基質吸力有利于土體結構穩(wěn)定，不易發(fā)生變形?，F有同工況下的試驗曲線[20-21]也表現出與模擬結果同樣的規(guī)律。這說明，數值實驗能夠反映壓縮過程中土體宏觀變形特征。

3.2 壓縮過程中飽和度和基質吸力的變化

分析在壓縮過程中，因顆?；?，孔隙變形，引起飽和度和吸力的變化。飽和度Sr和基質吸力ψ隨荷載p增加的變化情況見圖8。

根據圖8，隨荷載的增加，Sr-p呈上升趨勢，ψ-p曲線呈下降趨勢。飽和度Sr的上升幅度隨初始含水率的增加而增加。說明土體結構的變形與含水率有關。當含水率相對較低時，土體抵抗變形的能力較強，壓縮終止后土體結構變形很小。隨含水率的逐漸增加，土樣的變形量逐漸增大，這與圖7反映的結果一致。

圖8 壓縮過程中飽和度和基質吸力的變化Fig.8 Variation of saturation and matric suction during compression

對于ψ-p曲線，壓縮結束后，各含水率條件下基質吸力ψ降幅均較小，其中10%含水率下的降幅最小。原因在于：基質吸力在壓縮過程中受到毛細作用與飽和度2個因素的控制，且這2個因素作用效果相反。隨著土樣壓縮，孔隙半徑減小，毛細作用增強，基質吸力增大；另一方面由于孔隙體積變小，飽和度增大，導致基質吸力減小。由壓縮曲線（圖7）可知，當含水率較低為5%時，壓縮終止后，土體的孔隙比，或孔隙體積變化不大，所以毛細作用引起的基質吸力變化不大；而含水率較高為18%時，孔隙體積變化大，但由于土樣在18%的含水率下已接近飽和，毛細作用對基質吸力的影響很微弱。因此這2種狀態(tài)下，毛細作用影響微弱，基質吸力主要隨壓縮過程中飽和度的增加而減小。當含水率為10%時，毛細作用對基質吸力的增強作用與飽和度對基質吸力的減弱作用都有較大的影響，兩者相互抵消，因此基質吸力變化量小于5%與18%的土樣。但3種含水率下基質吸力總體呈下降趨勢，說明飽和度相較于毛細作用而言，是影響基質吸力的主要因素。

總的來說，隨荷載的增加，土樣飽和度上升，基質吸力下降，飽和度漲幅隨土樣含水率增加而增加，基質吸力在土樣干燥及近飽和時降幅較大。

3.3 微觀結構變形分析

數值模擬可以直觀地觀察土壤微觀結構的變化過程，并可以追蹤特定孔隙的在不同含水率狀態(tài)時的變形過程，這是實驗技術難以達到的。圖9展示了初始土樣（含水率0%）的孔隙分布與5%、10%、18%含水率下壓縮終止時的孔隙分布情況。并追蹤觀察了頂（孔隙a）、中（孔隙b）、底（孔隙c）部3個典型孔隙在不同含水率下的變化。

根據圖9，在荷載作用下，土樣的變形隨含水率的上升有顯著增加，且含水率較低時，顆粒接觸方式包括大量的點接觸和面接觸，見圖9（a），隨著含水率增高，土體在壓縮后點接觸數量明顯減少，以面接觸為主，見圖9（b）。這說明，含水率越高，土壤越易被壓縮。

圖9 初始孔隙分布與不同含水率條件下壓縮終止時孔隙分布Fig.9 Model pore distribution at the end of compression under different water content

對于上部孔隙a，在土樣含水率為5%時，孔隙a在壓縮終止時產生了變形，孔隙面積被壓縮，但孔壁顆粒位移較小，孔隙仍保持初始形態(tài)；含水率增大到10%和18%時，孔隙a均被壓扁，孔壁顆粒發(fā)生了明顯的位移，孔壁坍塌，孔壁顆粒進入孔隙中，孔隙由1個大孔分解成幾個小孔隙。

對于中部孔隙b，在土樣含水率為5%時，孔隙b未發(fā)生明顯的變形，孔隙基本維持初始大小和形態(tài)；在含水率達到10%時，孔隙b面積被壓縮減小，但孔隙形態(tài)未發(fā)生明顯的改變，孔壁未坍塌；當含水率進一步增大達到18%接近飽和時，孔隙b孔壁顆粒發(fā)生較大的位移，孔壁坍塌，孔隙分解成多個小孔。

對于底部孔隙c，在土樣含水率為5%時，孔隙基本未產生變形；含水率達到10%時，孔隙面積有微弱減小，形態(tài)不變；含水率達到18%接近飽和時，孔壁坍塌，孔隙分解成若干個小孔，原始孔隙發(fā)生徹底改變。

在壓縮時，上部孔隙（孔隙a）在含水率低時即產生變形，隨著含水率增大孔隙變形隨之增大；中部孔隙（孔隙b）在含水率較小時基本不發(fā)生變形，隨含水率增大發(fā)生變形，接近飽和時孔壁坍塌；底部孔隙（孔隙c）在含水率較低或中等時均很難發(fā)生變形，含水率接近飽和時，孔壁坍塌孔隙分解。這說明，在土體壓縮時，靠近受壓位置的孔隙優(yōu)先發(fā)生變形，遠離受壓位置的孔隙不易發(fā)生變形。曹亮等[22]利用改進的細觀結構觀測固結儀，觀察到了土樣變形過程中相似的現象；孔隙的變形形式與水分含量有關，隨含水率增大，孔隙變形逐漸由體積收縮變化演變?yōu)榭妆谔紫斗纸獬筛〉目紫丁?/p>

將土體在壓縮過程中的孔隙變形模式概化為圖10所示的概念模型，其規(guī)律為：在含水率較低的情況下，孔壁顆粒穩(wěn)定，孔隙以收縮變形為主，如圖10（b）。當含水率升高后，孔壁顆粒坍塌進入孔隙內部，以分解變形為主，如圖10（c）。Xie等[3]和Wang等[4]通過實測試驗對土體壓縮過程中孔隙、顆粒形態(tài)的統計表明，高含水率下加載過程中，土顆粒會失去穩(wěn)定，顆粒將重新排列，導致孔隙大小、形態(tài)和顆粒方向都發(fā)生了變化。這實際上也從側面證明了土體在高含水率下發(fā)生了分解變形，孔壁坍塌，較大的孔隙分解為小孔隙。

圖10 不同含水率下孔隙變形的主要形式Fig.10 Main forms of pore deformation under different water content

不同含水率下變形模式的差異本質上與基質吸力有關?；|吸力的存在，增加了土顆粒之間的有效應力，從而增強土顆粒之間的摩擦力，使得土顆粒之間更難發(fā)生滑移，利于維持土體結構，土體孔隙以收縮變形為主?；|吸力隨含水率的增加而減小，維持土顆粒穩(wěn)定的能力將隨之減弱，土顆粒之間易于滑移，孔隙的變形形式逐漸轉為分解變形為主。在基質吸力很低，土體接近飽和時，土樣極易發(fā)生壓縮變形，顆粒極易發(fā)生位移，孔隙產生較大變形，土體表現出整體的軟化現象，土體強度喪失。

4 結論

（1）根據實際黃土骨架顆粒分布、形狀建立了理想非飽和土模型，將毛細力作用施加于土骨架顆粒間，在整個結構上施加外荷載，模擬了非飽和土的常含水率壓縮過程。獲得的壓縮曲線與同工況下的試驗結果規(guī)律一致，表明本文的數值試驗能夠反映土體的宏觀變形特征，即隨荷載增加，孔隙比逐漸減小，并且同一荷載下，含水率越高，土壤的壓縮量越大。

（2）對壓縮過程中土水作用模擬結果表明，常含水率壓縮引起土壤飽和度增大，基質吸力減小。隨含水率的升高，飽和度的增幅逐漸增大；基質吸力的下降程度受毛細作用與飽和度2個因素影響，其在干燥與近飽和時降幅較大。

（3）追蹤土體特定孔隙的演變過程，結果表明在土體壓縮時，靠近受壓位置的孔隙優(yōu)先發(fā)生變形，遠離受壓位置的孔隙不易發(fā)生變形；孔隙的變化形式與含水率有關，低含水率時，孔隙主要為體積減小的收縮變形，高含水率時，孔隙主要為孔壁坍塌、孔壁顆粒進入孔隙內部的分解變形。