非飽和土粒間毛細(xì)作用的微觀不連續(xù)變形分析*

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

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

0 引 言

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

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

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

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

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

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

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

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

(1)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖5 樣品粒度分布曲線

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

圖6 鏡下粗顆粒形狀及多邊形簡(jiǎn)化后形狀

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

圖7 微觀黃土骨架DDA模型

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

表1 模型參數(shù)表

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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