李 強(qiáng) 李同錄 喬志甜 張常亮 李 萍 沈 偉 郭龍驍

(①長安大學(xué)地質(zhì)工程與測繪學(xué)院， 西安 710054， 中國)(②黃土高原水循環(huán)與地質(zhì)環(huán)境教育部野外科學(xué)觀測研究站， 正寧 745399， 中國)(③博洛尼亞大學(xué)地質(zhì)和環(huán)境科學(xué)學(xué)院， 博洛尼亞 40121， 意大利)(④九州大學(xué)土木與結(jié)構(gòu)工程學(xué)院， 福岡 819-0395， 日本)

0 引 言

非飽和土是一個固-液-氣三相體系，土-水相互作用是決定其力學(xué)性質(zhì)的內(nèi)在因素(Fredlund et al.，1993； 陳鴻賓等， 2019; 高英等， 2019)。而土-水相互作用是一個復(fù)雜的物理化學(xué)過程，包括結(jié)合水的吸附作用以及毛細(xì)水的表面張力作用。其中：毛細(xì)作用是基質(zhì)吸力的重要來源之一，相對飽和土，基質(zhì)吸力增加了非飽和土顆粒間的有效應(yīng)力，從而提高了土的抗剪強(qiáng)度，對非飽和土性質(zhì)有重要影響(邢鮮麗等， 2015； 蔣秀姿等， 2018； 李威等， 2018； 李同錄等， 2019)。

由于毛細(xì)作用是一種微觀現(xiàn)象，其本質(zhì)需要從微觀角度去探索。早期對非飽和土中毛細(xì)作用的微觀研究基于理想的接觸球模型(Fisher， 1926)。這種方法將土顆粒簡化為等大或不等大的光滑球體，將其抽象為立方體或四面體的結(jié)構(gòu)形式的土粒骨架，在土粒骨架間施加毛細(xì)水作用，分析土中的毛細(xì)水分布、表面張力、球間引力等。這種簡化模型有助于定性地理解毛細(xì)水在非飽和土中的作用。隨著數(shù)值計算能力的提高，為利用數(shù)值模擬方法來分析土粒間毛細(xì)作用提供了可能。土體在微觀上是不連續(xù)的，因此目前主要采用離散元法(DEM)和非連續(xù)變形分析法(DDA)等離散介質(zhì)模型來模擬其微觀力學(xué)行為。根據(jù)顆粒形狀，離散元(DEM)又可分為塊體離散元和球體離散元2類。目前基于塊體離散元對非飽和土的微觀力學(xué)性質(zhì)研究較少，大多基于球體離散元。對于球體離散元，顆粒的建模方式一般分為兩種，第1種是使用單一的圓盤或球形，另一種是多個圓盤通過黏結(jié)力進(jìn)行黏結(jié)形成類多邊形(周偉等， 2009； 李楊等， 2017)。這種方法能夠有效地模擬顆粒的破裂過程，但形成的顆粒與實際顆粒的體積有一定差距，會導(dǎo)致計算時產(chǎn)生誤差。在利用球體離散元模擬非飽和土中的毛細(xì)作用時，通常采用圓盤或球形顆粒構(gòu)建土骨架，在相同大小顆粒接觸處施加相同大小的應(yīng)力以考慮毛細(xì)力的作用(Jiang et al.，2004； Shen et al.，2016)。這種算法的主要問題在于未考慮不同含水率下毛細(xì)水的分布狀態(tài)，而毛細(xì)水的分布狀態(tài)與土體力學(xué)性質(zhì)有直接關(guān)系。DDA最初用于模擬塊狀結(jié)構(gòu)巖體的變形破壞過程(鄔愛清等， 1997； 周少懷等， 2000； 焦玉勇等， 2007)，由于土體中顆粒多呈多棱角的不規(guī)則形狀，所以近年來開始用于分析土體的微觀力學(xué)作用。其中張國新等(2000)利用規(guī)則的多邊形構(gòu)建微觀砂土DDA模型，研究了砂土微觀變形機(jī)理。郭培璽等(2008)通過將假想多邊形隨機(jī)填充的方式構(gòu)建了粗粒料DDA模型，基于此模型研究了粗粒料的力學(xué)特性及顆粒運(yùn)動規(guī)律。郭龍驍?shù)?2017)將實際黃土顆粒的形狀通過人工擺放的方式構(gòu)建黃土微觀模型，模擬了微觀黃土單向固結(jié)實驗。這些研究雖然對微觀模型做了一些改進(jìn)，并表明了DDA用于土粒微觀力學(xué)作用研究的可行性。但都基于傳統(tǒng)的二維DDA算法，尚未考慮毛細(xì)水作用。實際問題大多是三維，所以近幾年開始發(fā)展三維DDA算法(Zhang et al.，2016； Meng et al.，2019)，但由于其復(fù)雜的接觸關(guān)系造成自由度的大幅增加，使得三維DDA算法進(jìn)展緩慢，目前還很不成熟。因此，本課題組近年來開始研究考慮非飽和土毛細(xì)水的二維DDA算法(Guo et al.，2019)。

本文在二維DDA算法的基礎(chǔ)上，提出了一種用于非飽和微觀結(jié)構(gòu)力學(xué)分析的拓展DDA算法。該算法可算出不同飽和度時水分在土顆粒之間的分布形態(tài)及彎液面半徑，再用Young-Laplace方程計算出粒間毛細(xì)負(fù)壓，用毛細(xì)負(fù)壓和表面張力表征毛細(xì)水對土顆粒的作用，這為進(jìn)一步開展非飽和土各種力學(xué)行為的分析奠定了基礎(chǔ)?；谠撍惴?，構(gòu)建了一個中粗顆粒微觀模型，對不同含水率下土體中的毛細(xì)水分布及基質(zhì)吸力進(jìn)行計算，獲得模型土的土水特征曲線(SWCC)。

1 非飽和土粒間毛細(xì)力的計算

當(dāng)土處于非飽和狀態(tài)時，土粒之間主要由水聯(lián)結(jié)，在土中形成基質(zhì)吸力。當(dāng)土特別干時，如黃土含水率一般在7%以下，粒間聯(lián)結(jié)以結(jié)合水為主(王鐵行等， 2014)； 當(dāng)含水率增高時，結(jié)合水膜增厚，結(jié)合水聯(lián)結(jié)減弱，土粒間以毛細(xì)水作用為主。天然狀態(tài)下伏于地下的土都偏濕，以毛細(xì)水作用為主，因此暫先不考慮結(jié)合水的作用，只考慮土粒間的毛細(xì)水聯(lián)結(jié)。

1.1 基本假設(shè)與算法流程

通常所說的毛細(xì)現(xiàn)象就是指液體在毛細(xì)管中的移動，表面張力和彎液面內(nèi)外壓差是造成這一現(xiàn)象的本質(zhì)因素(高世橋等， 2010)。其中表面張力是水體在毛細(xì)管中移動的起因，內(nèi)外壓差是其動力來源。

在巖土工程領(lǐng)域，通常將土壤孔隙理想化為毛細(xì)管，并將彎液面的內(nèi)外壓差稱為基質(zhì)吸力，其與交界面幾何形狀的關(guān)系滿足Young-Laplace方程，在二維情況下表現(xiàn)形式為：

(1)

式中：ua為孔隙氣壓；uw為孔隙水壓；ua-uw為基質(zhì)吸力；Ts為表面張力；θ為土水間的接觸角；r液面曲率半徑。

式(1)中，表面張力Ts與接觸角θ在溫度和材料等條件確定的情況下均為常量，所以基質(zhì)吸力的大小主要取決于彎液面的曲率半徑r。由于實際土體中毛細(xì)作用復(fù)雜，為了建立模型，對復(fù)雜條件需進(jìn)行適當(dāng)簡化，本文的模型基于如下假設(shè)：

(1)土顆粒表面為理想表面，忽略增濕和減濕接觸角的不同； 不考慮礦物表面特征對接觸角的影響，假定所有土顆粒的接觸角都相等。

(2)不考慮滲流時間和滲流路徑，認(rèn)為土體系統(tǒng)內(nèi)任意時刻毛細(xì)水分布均處于平衡狀態(tài)。平衡狀態(tài)的條件是：水壓、氣壓處處相等； 水和氣在液面處凝結(jié)和蒸發(fā)平衡； 所有彎液面曲率半徑相等。

(3)土體中所有土顆粒為凸形顆粒。

基于上述假設(shè)，用Young-Laplace方程可計算出粒間毛細(xì)水的基質(zhì)吸力。本算法是在原有DDA算法(Shi， 1992)的基礎(chǔ)上進(jìn)行拓展，拓展部分的算法如圖1所示。

圖1 拓展DDA算法中粒間毛細(xì)力計算流程

1.2 顆粒間毛細(xì)水量的計算

土粒間的毛細(xì)水分布在顆粒接觸處和距離最近處。傳統(tǒng)的DDA算法可以識別所有土粒間的接觸點和接近點，對于任意接觸點，其鄰近有兩個或多個土粒，首先判定兩個土粒間毛細(xì)水的分布和面積，在此基礎(chǔ)上可判定多個土粒間毛細(xì)水的分布和面積。

兩個土粒間毛細(xì)水的分布是由曲率半徑相等的兩個彎液面和兩個土粒的邊界所限定。由于土粒位置固定，若給定彎液面的曲率半徑r、確定出兩個彎液面的圓心點，就可算出毛細(xì)水的分布區(qū)域。確定彎液面圓心可采用圓心軌跡交匯法(圖2)，具體計算過程如下：

圖2 雙顆粒系統(tǒng)中毛細(xì)水分布示意圖

(1)以選定的某一顆粒A為對象，假定顆粒A表面任意一點為液面的起始點SA1，過SA1作顆粒邊的外法線，將外法線逆時針旋轉(zhuǎn)接觸角θ，在旋轉(zhuǎn)后的線上確定與SA1的距離為r的一點OA1(圖2a)。

(2)使線段SA1OA1長度與外法線夾角大小不變，沿顆粒A表面逆時針移動一周，端點OA1在顆粒A外圍形成一條平行于顆粒表面直邊的軌跡線。

(3)該軌跡線在角點處存在缺口，故顆粒角點處的軌跡需要單獨處理。以角點為圓心，以r為半徑，作一段圓弧，將缺口處的軌跡線連接起來(圖2c)。按此方法處理后，顆粒A外圍將形成一條完整的圓心軌跡線。

(4)當(dāng)相鄰兩個顆粒之間的最近距離小于2r時，才可能形成毛細(xì)水彎液面。故在距顆粒A表面距離為2r的范圍內(nèi)尋找相鄰顆粒B。在顆粒B上隨機(jī)選取一點SB1，過SB1做顆粒表面線的外法線，將外法線順時針旋轉(zhuǎn)接觸角θ，并確定與SB1距離為r的一點OB1，按如上相同的做法使SB1OB1順時針繞著顆粒B表面一周，做出相鄰顆粒B的圓心軌跡線。

(5)顆粒A和顆粒B的圓心軌跡線存在兩個交點，分別以這兩個交點為圓心，以A顆粒的邊為起點，B顆粒的邊為終點，逆時針畫出兩個圓弧。這兩個當(dāng)中有一個是凹液面，另一個是凸液面。凸液面不合理，應(yīng)去除。

(6)做兩個液面與A顆?；駼顆粒交點處的液面切線，設(shè)該切線在液體一方與固液交界線之間的夾角為β。當(dāng)βi>90°時，液面為凸液面，將其刪除； 當(dāng)βi<90°時，為凹液面，將其保留。

(7)將顆粒A和顆粒B的角色置換，重復(fù)步驟(1)～(6)則可確定另一側(cè)的凹液面(圖2b)。將兩側(cè)凹液面組合，則確定了毛細(xì)水分布區(qū)(圖2d)。

圖3 多顆粒中毛細(xì)水彎液面形式

在圖3的3個土顆粒中，共有3個毛細(xì)區(qū)， 6個液面。當(dāng)土中含水量較低時，各液面相互之間不會產(chǎn)生重疊(圖3a)。當(dāng)土中含水量逐漸增多時，彎液面C2、C3、C5會出現(xiàn)重疊區(qū)域，在未重疊區(qū)域內(nèi)會出現(xiàn)一個封閉的氣泡(圖3b)。在實際土體中，這個氣泡在二維平面上應(yīng)當(dāng)呈圓形，并非曲邊三角形，但這并不影響其中含水量的計算以及對土性質(zhì)的分析。

可以看出，彎液面曲率半徑r的大小與土中的含水率Vw之間存在相關(guān)關(guān)系。若人為給定含水率Vw，需通過反復(fù)迭代計算出r的值，具體算法如下：

(1)先假定一個彎液面初始半徑區(qū)間[ra，rb]，按上面所述的算法分別確定ra和rb所對應(yīng)的所有土粒之間的彎液面，彎液面和所連接的土粒邊界構(gòu)成封閉的毛細(xì)水分布區(qū)。隨即將毛細(xì)水分布區(qū)的像素指定為特定且與模型中其他任何要素不同的顏色屬性。最后通過統(tǒng)計整個模型中毛細(xì)水像素所代表的面積，就可以計算毛細(xì)水的含量Va和Vb，這種方法優(yōu)點在于計算簡單，精度相對較高，并且不會重復(fù)計算毛細(xì)水的重疊部分。

(2)若計算出的Vb>Vw>Va，采用黃金分割算法逼近Vw。若Va>Vw或Vb

1.3 毛細(xì)作用力計算

在實際土體中，毛細(xì)水作用在土顆粒間的力包括凹透鏡形水體內(nèi)的毛細(xì)負(fù)壓(ua-uw)作用在土顆粒上的力和兩側(cè)的表面張力Ts。其中毛細(xì)負(fù)壓(ua-uw)是作用在土顆粒與水的接觸面上的面力，表面張力Ts是作用在三相界線上的線力。據(jù)此，該算法在二維平面上進(jìn)行了相應(yīng)簡化，將毛細(xì)負(fù)壓作用在被毛細(xì)水浸濕的顆粒邊界上，方向指向浸濕邊界的外法線方向(圖3)。計算時，將毛細(xì)負(fù)壓(ua-uw)作用在土顆粒上的力簡化為(ua-uw)×L的力，記為fc，其中，L為毛細(xì)水浸濕段長度，毛細(xì)負(fù)壓(ua-uw)在r已確定的情況下，可由式(1)獲得，作用點位于浸濕段中點，方向指向顆粒表面的外法線(圖4)。

將表面張力Ts僅作用于固、液、氣三相交界點處，在溫度等條件確定的情況下，其值大小為一常量，方向沿彎液面的切線由兩端指向中間(圖4)。

圖4 單顆粒上毛細(xì)作用力計算

在直角坐標(biāo)系中，將所有力都分解到x軸和y軸方向上，再求和。同時還必須計算所有毛細(xì)力產(chǎn)生的力矩，以顆粒的形心為轉(zhuǎn)動中心，計算其合力矩(圖4)。將合力及合力矩添加到DDA的荷載向量中，對每個顆粒都做如此計算，整個土粒間的毛細(xì)作用就施加在計算模型中了。值得注意的是，封閉氣泡中的吸力和表面張力未直接作用于土顆粒上，故不考慮氣泡周圍的力。

2 黃土微觀結(jié)構(gòu)模型的構(gòu)建

現(xiàn)以黑方臺L1黃土樣為參考，建立其理想土骨架結(jié)構(gòu)DDA模型，來研究黃土中的土-水相互作用。該土樣的干密度為 1.32 g·cm-3，比重為2.71，孔隙比為1.06，液、塑限分別為27.5%和18.3%，圖5為其粒度分布曲線。

圖5 樣品粒度分布曲線

考慮到計算效率，模型的顆粒數(shù)量不宜太多。根據(jù)圖5， 20～30μm是顆粒組分相對較多的部分，本文只選取該范圍內(nèi)的粉粒建立理想黃土微觀骨架模型。在放大100倍下顯微鏡下觀察了大量的黃土中該范圍內(nèi)粉粒散粒的形狀，圖6a為其中部分顆粒，發(fā)現(xiàn)大部分粗顆粒多呈不規(guī)則多邊形或扁平狀。將鏡下觀察的顆粒進(jìn)行矢量化，由于矢量化后的形狀包含的角點數(shù)量較多，為保證計算效率，需簡化。簡化原則為當(dāng)多個角點位于同一邊或近似在同一邊時，僅保留兩個角點(喬志甜等， 2020)。簡化后的多邊形的邊數(shù)盡量不大于8(圖6b)。

圖6 鏡下粗顆粒形狀及多邊形簡化后形狀

從中隨機(jī)選取350個大小為20～30μm不同形狀的顆粒用于生成土的模型。然后將選取的顆粒隨機(jī)丟落于內(nèi)徑為0.54mm×0.40mm的樣盒模型中，土中顆粒會在重力的作用下達(dá)到相對穩(wěn)定狀態(tài)，形成的模型作為初始模型(圖7)，相應(yīng)的初始孔隙比e為0.39，小于實際黃土的孔隙比。在模型生成過程中未考慮空氣阻力和靜電引力的影響，還難以形成與實際土一致的孔隙結(jié)構(gòu)，對此還有待進(jìn)一步研究。

圖7 微觀黃土骨架DDA模型

由于黃土中20～30μm的粉粒大多為石英，密度可參考石英礦物的密度。此外黃土顆粒最初來自物源區(qū)的巖石，模型使用的彈性模量E可以參考巖石中(花崗巖)的相應(yīng)參數(shù)。此外，本文只考慮毛細(xì)作用，故將黃土顆粒之間的黏聚力設(shè)為0。土樣盒也為一組DDA單元，相對于黃土顆粒，樣盒為剛體，故將樣品盒的彈性模量設(shè)為黃土顆粒彈模的20倍。樣盒與土粒之間的內(nèi)摩擦角、黏聚力均設(shè)置為0。具體參數(shù)如表1。

表1 模型參數(shù)表

3 毛細(xì)作用下土水特征曲線模擬結(jié)果

利用本文拓展的DDA算法，給定一系列的含水率，可計算出該土樣的飽和度及相應(yīng)的毛細(xì)水彎液面半徑，用Young-Laplace方程(1)可計算出不同半徑相應(yīng)的基質(zhì)吸力，計算結(jié)果如表2所示，由此得到飽和度與基質(zhì)吸力之間的關(guān)系，即SWCC曲線(圖8)。圖9給出了其中6種飽和度工況下的毛細(xì)水分布。

表2 不同飽和度下彎液面半徑及吸力值

圖8 SWCC曲線模擬結(jié)果與實測結(jié)果

圖8中的實測數(shù)據(jù)來自黑方臺L1原狀黃土壓力板儀實驗結(jié)果。由圖8中的DDA模擬結(jié)果可以看出，SWCC曲線存在著明顯的邊界效應(yīng)區(qū)(0～8kPa)、過渡區(qū)(8～110kPa)、殘余區(qū)(>110kPa)。在邊界效應(yīng)區(qū)，飽和度變化基本穩(wěn)定，基質(zhì)吸力在該范圍內(nèi)隨飽和度的下降緩慢增加。該階段土中水相分布狀態(tài)如圖9a，土體內(nèi)部大部分孔隙均被水體充滿，只有少量氣泡存在于水中。當(dāng)飽和度進(jìn)一步下降，SWCC將處于過渡區(qū)，邊界效應(yīng)區(qū)和過渡區(qū)之間的界限值即進(jìn)氣值，模擬試驗獲得進(jìn)氣值約為8kPa，對應(yīng)的飽和度為90.3%。在進(jìn)氣值附近，土體中的水相分布狀態(tài)如圖9b，此時土體中的氣泡開始擴(kuò)大并出現(xiàn)少量獨立的彎液面，土體開始進(jìn)入非飽和狀態(tài)。圖9c、圖9d是完全處于過渡階段時土體中水相分布狀態(tài)，土體中氣體開始占據(jù)大部分孔隙，獨立的彎液面數(shù)量明顯增多。此時，飽和度下降速度較快，基質(zhì)吸力的增幅相對較大。隨飽和度進(jìn)一步降低，SWCC將處于殘余區(qū)。過渡區(qū)與殘余區(qū)同樣存在一個界限吸力值，即殘余值，該值大小約為110kPa，對應(yīng)的飽和度6.31%。在該飽和度附近土體中水的分布狀態(tài)如圖9e，此時土體中的水僅分布于顆粒的接觸位置及較小的孔隙中。圖9f是完全處于殘余階段時土體中水的分布狀態(tài)，土體中大孔隙均被氣體占據(jù)，極少量的水分布在顆粒的接觸位置，土體中水量較少，基質(zhì)吸力較高。并且在該階段，基質(zhì)吸力隨飽和度緩慢下降而迅速增加。

圖9 不同飽和度下毛細(xì)水的分布狀態(tài)

綜上所述，模擬結(jié)果符合典型土水特征曲線(McQueen et al.，1974)的形態(tài)。飽和度相對較高的情況下，模擬結(jié)果和實測曲線較為吻合。原因在于：此時，毛細(xì)作用是基質(zhì)吸力的主要來源。毛細(xì)水DDA算法以及模型雖對復(fù)雜條件進(jìn)行了一些簡化，但并未脫離實際情況。模型中使用的顆粒粒徑分布以及其形狀，都是參照實際黃土構(gòu)建的。并且，構(gòu)建的理想模型基本包含了實際土體中主要的孔隙結(jié)構(gòu)。所以，在毛細(xì)作用起主要作用的區(qū)間內(nèi)，能夠較為準(zhǔn)確地計算出基質(zhì)吸力。但隨著飽和度下降，與實測曲線有所偏差。這一差距主要是因為：實測的土水特征曲線中的吸力應(yīng)包含毛細(xì)水的毛細(xì)作用和結(jié)合水的吸附作用兩部分(Tuller et al.，2005)，一般在殘余含水率之前，毛細(xì)作用起控制作用，之后吸附作用起控制作用。此外，真實黃土中存在一定量的黏粒，黏粒往往會吸附較多的結(jié)合水，這類吸附作用主要來源于范德華力、電場力和水合作用，而非毛細(xì)作用。并且在非飽和土中，接觸角并非是個常數(shù)，會隨外界條件的變化而變化。一般在飽和度降低的時候，接觸角會變小(Yang et al.，2012； 楊松等， 2015)，進(jìn)而造成基質(zhì)吸力變大。這些影響因素中，未考慮結(jié)合水的吸附作用是導(dǎo)致低含水率模擬結(jié)果與實際結(jié)果誤差的主要原因，所以建立合理的結(jié)合水非連續(xù)變形算法是作者們將要進(jìn)一步研究的問題。

總的來說，本文使用的模型僅考慮了20～30μm的粗粒，且未考慮吸附作用和接觸角的變化，因此，隨著飽和度下降，差距也隨之明顯。但從整體上來看，該算法能夠較為合理地考慮毛細(xì)作用，并能計算不同飽和度條件下的毛細(xì)水分布狀態(tài)。

4 結(jié) 論

本文在只考慮毛細(xì)水的作用下，基于Young-Laplace方程，將毛細(xì)水作用力嵌入已有DDA算法中，利用擴(kuò)展后的DDA算法對不同飽和度工況下的土體中毛細(xì)水分布及基質(zhì)吸力進(jìn)行了模擬，結(jié)果表明：

(1)擴(kuò)展的DDA算法能夠通過反復(fù)迭代的方式確定不同飽和度下彎液面的半徑，并能夠利用圓心軌跡交匯法確定毛細(xì)水彎液面與顆粒的搭接位置，進(jìn)而計算出顆粒間毛細(xì)水分布范圍。通過給毛細(xì)水區(qū)域上色的方式，仿真模擬不同飽和度毛細(xì)水在土體中的分布狀態(tài)。

(2)擴(kuò)展的DDA能夠利用彎液面半徑準(zhǔn)確計算出土中的基質(zhì)吸力，并根據(jù)浸濕區(qū)域大小，將基質(zhì)吸力與表面張力一起作用于土顆粒之上，以此考慮毛細(xì)作用對土體力學(xué)性質(zhì)的影響，這為非飽和土宏觀力學(xué)行為的微觀分析提供了新途徑。

(3)選取黃土中含量最多粗顆粒建立黃土微觀模型，計算出的土水特征曲線與典型的土水特征曲線形態(tài)一致。并且在高飽和度條件下，與實測曲線吻合性較好。這表明擴(kuò)展的DDA算法能夠合理地考慮非飽和土粒間毛細(xì)作用。