非飽和土中彎液面形態(tài)與液橋力的分子動(dòng)力學(xué)模擬

駱莉莎,孫天佑,申志福,周 峰

(1.江蘇開放大學(xué) 建筑工程學(xué)院,江蘇 南京 210036;2.南京工業(yè)大學(xué) 交通運(yùn)輸工程學(xué)院,江蘇 南京 211800)

自Roscoe等[1]提出飽和黏土的劍橋模型以來(lái),土力學(xué)逐漸形成了以宏觀試驗(yàn)為依據(jù),連續(xù)介質(zhì)力學(xué)與宏觀本構(gòu)理論為基礎(chǔ),數(shù)值模擬為工具的實(shí)用體系,在巖土工程實(shí)踐中發(fā)揮了重要作用[2];而工程實(shí)踐中的非飽和土的理論研究則相對(duì)滯后得多。Bishop[3]提出的非飽和土有效應(yīng)力公式認(rèn)為非飽和土的強(qiáng)度取決于基質(zhì)吸力(即單變量理論)。Fredlund等[4]對(duì)非飽和土的應(yīng)力表征進(jìn)行了系統(tǒng)研究,提出了雙變量理論。沈珠江[5]提出非飽和土有效應(yīng)力原理需要做根本的修正,之后學(xué)者們紛紛提出多種具有一定微觀機(jī)制的非飽和土應(yīng)力理論[6-8]。以上研究都是從宏觀角度試圖描述非飽和土的力學(xué)特性。非飽和土的力學(xué)性質(zhì)遠(yuǎn)比飽和土復(fù)雜,歸根到底是微觀上氣-液-固三相的復(fù)雜相互作用。非飽和土液橋動(dòng)態(tài)演化及顆粒間液橋力鏈分布規(guī)律是揭示非飽和土宏、微觀力學(xué)特性機(jī)制的關(guān)鍵[9]。為此,需首先系統(tǒng)認(rèn)識(shí)非飽和土微觀尺度(納米尺寸)的力學(xué)規(guī)律;然而,天然非飽和黏土中片間毛細(xì)水接觸角、液橋力的定量測(cè)試存在極大困難,目前多以定性分析為主,這在一定程度上阻礙了非飽和土力學(xué)的進(jìn)一步發(fā)展。

得益于計(jì)算機(jī)軟硬件的飛速發(fā)展,從分子、原子層面對(duì)各類巖土體礦物與其他物質(zhì)之間的相互作用進(jìn)行模擬,成了巖土力學(xué)研究的重要手段。在巖土力學(xué)研究中,分子動(dòng)力學(xué)模擬較其他研究手段有其獨(dú)特的價(jià)值。分子動(dòng)力學(xué)模擬可以反映試驗(yàn)中難以直接觀測(cè)的納米尺度力學(xué)規(guī)律;分子動(dòng)力學(xué)模擬具有廣泛適用性,可同時(shí)反映巖土體微觀物理過(guò)程、力學(xué)過(guò)程、化學(xué)過(guò)程等,而無(wú)需對(duì)各個(gè)過(guò)程分別研究再試圖整合形成完整的規(guī)律表達(dá)。正是利用了這些優(yōu)勢(shì),何滿潮等[10]采用分子動(dòng)力學(xué)對(duì)高嶺石進(jìn)行單軸拉伸模擬,從微觀角度還原高嶺石變形及破壞的過(guò)程。Abduljauwad等[11]通過(guò)宏觀膨脹試驗(yàn)、微觀樣品研究以及分子水平建模3個(gè)階段，開發(fā)了土體納米級(jí)本構(gòu)模型,用于研究非飽和狀態(tài)下蒙脫土的膨脹率和膨脹壓力。Amarasinghe等[12]通過(guò)分子動(dòng)力學(xué)模擬獲得了非飽和土層間彎液面形態(tài),但其僅采用幾個(gè)彎液面瞬時(shí)形態(tài)確定接觸角,結(jié)果不夠精確。Zhang等[13]采用分子動(dòng)力學(xué)模擬獲得了非黏土礦物間液橋力。目前尚無(wú)從分子動(dòng)力學(xué)模擬角度對(duì)各種固相礦物(包括黏土礦物和非黏土礦物)間液橋形態(tài)、液橋力進(jìn)行系統(tǒng)研究,將所模擬的微觀力學(xué)規(guī)律與宏觀試驗(yàn)現(xiàn)象相結(jié)合的分析更是少有。

為此,本文采用分子動(dòng)力學(xué)模擬,對(duì)非飽和土中普遍存在的黏土礦物、非黏土礦物間液橋形態(tài)與液橋力進(jìn)行研究,探討各主要因素的影響規(guī)律,修正了利用Young-Laplace方程計(jì)算液橋力的理論方法,利用納米尺度模擬結(jié)果初步揭示了一些非飽和土宏觀試驗(yàn)現(xiàn)象的微觀機(jī)制。

1 固相礦物間液橋的分子動(dòng)力學(xué)模擬

一般而言,黏土中可同時(shí)包含黏土片和粒狀非黏土礦物顆粒,物理原型示意圖如圖1所示。由圖1可知:毛細(xì)水可存在于黏土片之間、非黏土礦物顆粒之間以及黏土片-非黏土礦物顆粒之間。本文選擇正長(zhǎng)石(KSi3Al3O12H12)及α-石英(SiO2)兩種原生礦物作為非黏土礦物的代表,選擇高嶺石(Al4[Si4O10](OH)8,不考慮同晶置換)作為黏土礦物的代表(蒙脫石需考慮晶層間水分子的嵌入及毛細(xì)水中溶解的離子,建模復(fù)雜;為簡(jiǎn)化起見,本文未選擇蒙脫石)。本文僅模擬圖1中黏土片、非黏土礦物顆粒表面平行的情形,為降低計(jì)算量,固相部分僅模擬圖1中陰影區(qū)域。試算結(jié)果表明:進(jìn)一步增大模擬的固相礦物厚度和寬度不影響本文結(jié)論。首先,單獨(dú)構(gòu)建各固相礦物和液相水的分子動(dòng)力學(xué)模型;隨后,將二者組合構(gòu)建同種或異種礦物間毛細(xì)水的分子動(dòng)力學(xué)模型,進(jìn)而模擬液橋彎液面形態(tài)。本文采用開源程序LAMMPS計(jì)算,得到南京工業(yè)大學(xué)高性能計(jì)算中心的計(jì)算支持。

圖1 分子動(dòng)力學(xué)模型對(duì)應(yīng)的物理原型(未按比例繪制)

1.1 單晶胞與固相礦物模型

本文所用固相礦物的晶體模型已多次被國(guó)內(nèi)外學(xué)者引用并證實(shí)[14-16],晶胞參數(shù)如表1所示,單晶胞模型如圖2(a)—2(c)所示。將單晶胞按照3a×31b×c的方式擴(kuò)展構(gòu)建正長(zhǎng)石薄片,按照4a×95b×c的方式擴(kuò)展構(gòu)建α-石英薄片,按照4a×45b×c的方式擴(kuò)展構(gòu)建高嶺石薄片,以形成尺寸相近的片狀礦物顆粒,其中,a、b、c分別為單晶胞寬度(x)、長(zhǎng)度(y)、厚度(z) 3個(gè)維度的擴(kuò)展,前綴數(shù)值代表單晶胞擴(kuò)展份數(shù)。以高嶺石為例,擴(kuò)展形成的薄片尺寸為2.06 nm×40.25 nm×0.74 nm。

表1 單晶胞參數(shù)

1.2 水分子與液相模型

本文采用在黏土模擬中廣泛應(yīng)用的簡(jiǎn)單點(diǎn)電荷(SPC)水分子模型,其原子核與電子軌道可簡(jiǎn)化為點(diǎn)電荷,通過(guò)調(diào)整H—O—H鍵角與O—H鍵長(zhǎng)以平衡水分子偶極矩,具體參數(shù)如圖2(d)所示。建立3 nm×27 nm×27 nm的周期性空間,空間中隨機(jī)生成若干上述柔性SPC水分子，并弛豫平衡以形成“水盒子”,即形成液相模型。

圖2 單晶胞及SPC水分子模型

1.3 固相礦物-毛細(xì)水組合模型

非飽和土中包括固-液-氣三相,理想情況下分子模型體系中應(yīng)包括氣體分子。但實(shí)際模擬中發(fā)現(xiàn),空氣分子在模型體系中的數(shù)量占比非常少,本文按照文獻(xiàn)中的慣常處理方法,即在模擬中不考慮空氣分子,這種處理不影響模擬結(jié)論[12]?，F(xiàn)以高嶺石-毛細(xì)水-高嶺石組合模型為例,表述建模過(guò)程。構(gòu)建2.06 nm(x方向)×56.00 nm(y方向)×24.00 nm(z方向)的三向周期性邊界空間,將上、下兩片片狀高嶺石礦物平行放置于模型空間中間位置,兩片之間凈距(d)為4.3 nm(足夠大以避免顆粒表面水化層的重疊[16])。兩片之間嵌入切割后的“水盒子”(尺寸為lx×ly×lz),lx=2.06 nm,為“水盒子”在x方向的尺寸;ly=16.00 nm,為“水盒子”在y方向的尺寸;lz=4.00 nm,為“水盒子”在z方向的尺寸。x方向上周期性模型空間尺寸、固相礦物尺寸與“水盒子”尺寸一致；y方向上固相礦物顆粒與周期性邊界的距離約為8 nm，以消除周期性邊界對(duì)計(jì)算結(jié)果的影響[13];固相礦物與“水盒子”之間的距離為0.15 nm，以防止水分子與固相礦物原子重疊[13]。高嶺石-毛細(xì)水模型的初始構(gòu)型如圖3所示。

圖3 高嶺石-毛細(xì)水模型初始構(gòu)型

1.4 力場(chǎng)

本文采用Cygan等[17]提出的CLAYFF力場(chǎng)進(jìn)行計(jì)算。在CLAYFF力場(chǎng)中,相互作用能由4個(gè)部分組成,分別為庫(kù)倫能、范德華能、鍵長(zhǎng)伸縮勢(shì)能、鍵角扭曲勢(shì)能,其中,鍵長(zhǎng)伸縮勢(shì)能與鍵角扭曲勢(shì)能構(gòu)成力場(chǎng)中鍵的相互作用能(Ub),僅存在于水分子與黏土羥基中,Ub計(jì)算式為

(1)

式中:k1、k2分別為鍵長(zhǎng)和鍵角計(jì)算常數(shù);r0、θ0分別為平衡時(shí)的鍵長(zhǎng)與鍵角;rij、θijk分別為鍵長(zhǎng)與鍵角,下標(biāo)ij、ijk表示原子i、j與k的不同組合;“bond”表示對(duì)所有鍵的鍵長(zhǎng)伸縮勢(shì)能求和；“angle”表示對(duì)所有鍵的鍵角扭曲勢(shì)能求和。

力場(chǎng)內(nèi)其他相互作用能由庫(kù)倫能與范德華能描述,庫(kù)倫能表示點(diǎn)電荷之間的靜電相互作用,范德華能常用Lennard-Jones(L-J)函數(shù)[18]描述,二者構(gòu)成力場(chǎng)中非鍵的相互作用能(Unb),Unb計(jì)算式為

(2)

式中:e為電子電荷,ε0為真空介電常數(shù),qi、qj分別為原子i、j的電荷,D0,ij、R0,ij為原子對(duì)i與j的等效范德華能參數(shù)。

D0,ij與R0,ij計(jì)算式分別為

(3)

(4)

式中:D0,i、R0,i為原子i的范德華能參數(shù),D0,j、R0,j為原子j的范德華能參數(shù)。

本文所有力場(chǎng)模型參數(shù)見文獻(xiàn)[17]。

1.5 液橋彎液面形態(tài)模擬

對(duì)初始構(gòu)型進(jìn)行能量最小化處理后,在恒原子數(shù)、恒溫、恒體積系統(tǒng)下模擬毛細(xì)水形態(tài),采用Nose-Hoover法控制體系的溫度為298 K,體積不變。依據(jù)文獻(xiàn)[9-12]結(jié)果和本研究模擬嘗試,L-J相互作用的截?cái)喟霃皆O(shè)置為10 ?,靜電相互作用采用顆粒-顆粒與顆粒-網(wǎng)格(PPPM)算法處理,精度為0.000 1,每2步檢查一次近鄰表,時(shí)間步為1.0×10-15s。

在模擬過(guò)程中，固定固相原子位置,水分子在力場(chǎng)作用下運(yùn)動(dòng)至平衡狀態(tài),形成彎液面。以正長(zhǎng)石-毛細(xì)水模型為例,總模擬時(shí)間超過(guò)2.0×10-9s后達(dá)到平衡,即溫度和勢(shì)能達(dá)到穩(wěn)定(圖4)。

圖4 正長(zhǎng)石-毛細(xì)水模型中溫度及勢(shì)能隨時(shí)間變化曲線

1.6 模擬算例

本文通過(guò)11個(gè)算例分析了礦物種類、固相礦物間距、含水量、固相礦物排列取向4個(gè)因素對(duì)彎液面形態(tài)、接觸角、液橋力的影響,分子動(dòng)力學(xué)模擬結(jié)果如表2所示。由于高嶺石上、下表面基團(tuán)不同,為做區(qū)分,表2中高嶺石指羥基面在上、硅氧面在下的高嶺石模型,反轉(zhuǎn)高嶺石指硅氧面在上、羥基面在下的高嶺石模型。表2中不區(qū)分上下片時(shí),即指上、下兩個(gè)片狀固相礦物取向相同;因正長(zhǎng)石和α-石英礦物上、下面基團(tuán)相同,故不做區(qū)分。算例1—7模擬黏土礦物之間液橋作用,算例8、9模擬非黏土礦物之間液橋作用,算例10、11模擬黏土-非黏土礦物之間液橋作用。

表2 模擬算例與結(jié)果

2 液橋形態(tài)與基質(zhì)吸力

以算例1為例，分析分子動(dòng)力學(xué)模擬結(jié)果。圖5(a)—5(c)為算例1中高嶺石-毛細(xì)水體系在達(dá)到平衡后的3個(gè)瞬時(shí)構(gòu)型。根據(jù)各態(tài)歷經(jīng)假說(shuō)理論,體系達(dá)到平衡狀態(tài)后的瞬時(shí)構(gòu)型仍隨時(shí)間動(dòng)態(tài)變化,足夠長(zhǎng)時(shí)間的模擬后體系將經(jīng)歷所有可能的構(gòu)型,因此根據(jù)一段時(shí)間內(nèi)的平均原子位置確定的彎液面形態(tài)具有宏觀物理意義。本文獲取彎液面形態(tài)的具體方法:在yz平面均勻劃分出4 800個(gè)等面積的正方形小區(qū)域,統(tǒng)計(jì)平衡狀態(tài)后5.0×10-10s內(nèi)各小區(qū)域中水分子的平均密度。圖5(d)為高嶺石-毛細(xì)水體系中的水分子平均密度分布，r為彎液面曲率半徑，D為固相礦物(包括強(qiáng)吸附水)間距，θ1—θ4為4個(gè)位置的接觸角。由圖5(d)可知:彎液面邊緣有一層密度較低(<0.9 g/cm3)的區(qū)域,這是由于水分子在彎液面邊緣處有蒸發(fā)形成水汽的趨勢(shì)。而高嶺石表面存在一層密度波動(dòng)很大的水分子層,這是由于高嶺石表面基團(tuán)親水性強(qiáng),水分子受到吸引而形成強(qiáng)吸附水層。實(shí)際上,由表2可知:本文模擬的固相礦物表面都存在厚度不一的強(qiáng)吸附水層,厚度約為0.4～0.6 nm。在宏觀土力學(xué)中,對(duì)于粉粒、砂粒而言,該厚度非常小而不需考慮;對(duì)于黏土而言,該厚度與黏土片間距在一個(gè)數(shù)量級(jí),其作用不可忽略。

圖5 高嶺石-毛細(xì)水體系模擬結(jié)果(算例1)

為消除水汽界面及強(qiáng)吸附水層對(duì)彎液面接觸角獲取的影響,本文將礦物表面的強(qiáng)吸附水視為固相的一部分,而將上、下兩部分強(qiáng)吸附水之間，平均密度為0.9～1.1 g/cm3的區(qū)域稱為液橋。以液橋與強(qiáng)吸附水接觸位置的彎液面切線傾角作為接觸角(θ),取上片的左、右兩接觸角平均值作為上部接觸角,采用相同方法得到下部接觸角。各算例得到的平衡后某瞬時(shí)構(gòu)型、彎液面形態(tài)及接觸角均列于表2中。

本文模擬結(jié)果與采用靜態(tài)滴落法獲得的黏土和非黏土礦物毛細(xì)水接觸角試驗(yàn)值如表3所示。由表3可知:本文模擬的接觸角與試驗(yàn)值有一定的差異,這可能是由于試驗(yàn)測(cè)試的是單個(gè)平面上的毛細(xì)水形態(tài),而本文模擬的是兩礦物平板間的毛細(xì)水形態(tài)。此外,正長(zhǎng)石、α-石英與毛細(xì)水的接觸角分別與Zhang等[13]模擬得到的34°、29°接近;而高嶺石與毛細(xì)水的接觸角與Song等[22]模擬的另一種黏土礦物葉蠟石與毛細(xì)水的接觸角范圍25°～34°基本吻合(葉蠟石與高嶺石表面基團(tuán)相似),這說(shuō)明了本文模擬結(jié)果的可靠性。

表3 毛細(xì)水接觸角模擬值與試驗(yàn)值

Young-Laplace方程是描述液橋彎液面兩側(cè)空氣與水的壓力差(即狹義的基質(zhì)吸力)與彎液面曲率半徑之間關(guān)系的控制方程,如式(5)所示。

(5)

式中:ua為氣壓;uw為水壓;σ為水的表面張力,298 K時(shí)σ=71.97×10-3N/m;R1、R2分別為彎液面沿切面的兩個(gè)垂直方向的曲率半徑。

由于本文模擬未考慮空氣分子,故式(5)中ua=0。在本文模擬的情形中,彎液面沿x方向(周期性擴(kuò)展方向)無(wú)限擴(kuò)展,對(duì)應(yīng)的曲率半徑無(wú)窮大;將oyz平面上的液橋彎液面假設(shè)為半徑為r的圓弧,則式(5)可轉(zhuǎn)變?yōu)?/p>

(6)

r與D的幾何關(guān)系為

(7)

根據(jù)式(6)和(7)計(jì)算所得的平均基質(zhì)吸力(表2)與分子動(dòng)力學(xué)模擬的液橋中水分子集合的壓力基本一致,而與基于所有水分子集合(液橋+強(qiáng)吸附水)的壓力相差很大，這說(shuō)明式(6)和(7)適合描述受固相礦物作用很小、相對(duì)自由的液橋特性。

算例1—3保持毛細(xì)水體積不變而改變lz,可近似反映非飽和土中含水量不變而孔隙比增大(或干密度降低)、飽和度降低的工況,結(jié)果表明:在相同含水量下,彎液面接觸角隨lz(或孔隙比)的增大而增大,基質(zhì)吸力隨lz(或孔隙比)的增大而減小。對(duì)多種非飽和土進(jìn)行試驗(yàn)，結(jié)果表明:在“含水量(質(zhì)量分?jǐn)?shù)或體積分?jǐn)?shù))-吸力”空間中,不同孔隙比試樣的土水特征曲線可存在交叉;當(dāng)吸力<吸力閾值時(shí),相同含水量下的吸力隨孔隙比的增大而增大;而當(dāng)吸力≥吸力閾值時(shí),相同含水量下的吸力隨孔隙比的增大而減小;吸力閾值因土的種類、礦物成分、孔隙比等因素不同而在很大范圍內(nèi)發(fā)生變化[23-30]。算例1—3的結(jié)果從微觀角度驗(yàn)證了王鐵行等[29]提出的“密度變化引起土體孔隙尺寸和薄膜水曲率半徑發(fā)生變化”機(jī)制猜想。因此,研究非飽和土的土水特征曲線、確定基質(zhì)吸力時(shí)，應(yīng)充分考慮土體密實(shí)狀態(tài)的影響。

算例1、4、5通過(guò)改變ly而增大毛細(xì)水體積,其余參數(shù)不變,可近似反映非飽和土中含水量增大而孔隙比不變、飽和度增大的工況,結(jié)果表明:當(dāng)固相礦物間距(或孔隙比)不變時(shí),毛細(xì)水體積(或含水量)的增大不會(huì)引起彎液面接觸角和基質(zhì)吸力的變化,此結(jié)果對(duì)應(yīng)土水特征曲線的陡降段。本文模擬的毛細(xì)水孤立存在于平行固相礦物薄片之間,相當(dāng)于理想化的單孔徑非飽和黏土微觀結(jié)構(gòu),Marinho[31]針對(duì)單孔徑假設(shè),從理論上推導(dǎo)出的土水特征曲線的確存在陡降段;而實(shí)際非飽和土為多孔徑材料,將各單孔徑對(duì)應(yīng)的土水特征曲線“疊加組合”,即可得與實(shí)際一致的緩變型土水特征曲線。本文模擬結(jié)果為基于孔徑分布預(yù)測(cè)非飽和土土水特征曲線提供了理論依據(jù)，而實(shí)際狀態(tài)下的礦物顆粒間取向隨機(jī)排列,如何反映這種排列對(duì)非飽和土土水特征曲線的影響還需進(jìn)一步研究。

算例1—5反映的高嶺石間彎液面接觸角變化規(guī)律與Zhang等[13]所得α-石英間彎液面接觸角的變化規(guī)律一致,說(shuō)明這些規(guī)律對(duì)非飽和土中的黏土礦物和非黏土礦物具有一定普適性。

算例1、6、7反映了高嶺石礦物的相對(duì)面原子類型(即片狀顆粒排列取向)對(duì)毛細(xì)水彎液面形態(tài)的影響,其中，算例6中羥基面相對(duì),而算例7中硅氧面相對(duì)(圖2(c))。由于Si—O鍵中的O原子與水分子中H原子形成的氫鍵平均鍵能稍小于羥基H—O鍵中的O原子與水分子中H原子形成的氫鍵平均鍵能,故算例6中羥基面相對(duì)時(shí)的毛細(xì)水彎液面接觸角大于算例7中硅氧面相對(duì)時(shí)的毛細(xì)水彎液面接觸角,而前者基質(zhì)吸力小于后者基質(zhì)吸力。黏土片排列取向的差異對(duì)非飽和土宏觀物理、力學(xué)特性的影響程度取決于實(shí)際土體中這些排列取向的相對(duì)比例,影響程度與相對(duì)比例的定量關(guān)系還需進(jìn)一步探究。

算例1、8、9反映了礦物類型對(duì)彎液面接觸角和基質(zhì)吸力的影響,結(jié)果表明:正長(zhǎng)石接觸角大于石英接觸角;而高嶺石接觸角居于正長(zhǎng)石接觸角與石英接觸角之間或者低于石英接觸角,這取決于高嶺石礦物間的排列取向;同時(shí)，相應(yīng)基質(zhì)吸力也有一定差異。

算例10、11給出了黏土-非黏土礦物之間液橋彎液面的形態(tài)與基質(zhì)吸力,結(jié)果表明:異種固相礦物間接觸角明顯大于任一單種固相礦物間接觸角。雖然這一結(jié)論尚無(wú)試驗(yàn)結(jié)果支撐(也很難進(jìn)行試驗(yàn)驗(yàn)證),但鑒于本文模擬方法與整體結(jié)果的可靠性,該結(jié)論可作為研究非黏土礦物含量影響非飽和黏土力學(xué)特性的微觀力學(xué)機(jī)制的參考[9]。

3 修正液橋力計(jì)算方法

分子動(dòng)力學(xué)模擬可直接計(jì)算固相礦物與毛細(xì)水間作用力,即液橋力(Fz)。算例1中高嶺石-毛細(xì)水體系中的上片高嶺石受到的動(dòng)態(tài)液橋力與時(shí)間平均液橋力隨時(shí)間變化曲線如圖6所示。

圖6 高嶺石-毛細(xì)水體系上片高嶺石受到的動(dòng)態(tài)液橋力與時(shí)間平均液橋力隨時(shí)間變化曲線

由圖6可以看出:在模擬的最后0.5×10-9s時(shí)間間隔內(nèi),時(shí)間平均液橋力基本穩(wěn)定,說(shuō)明系統(tǒng)達(dá)到平衡狀態(tài);而即使在平衡狀態(tài),動(dòng)態(tài)液橋力隨時(shí)間變化仍較大,這是因?yàn)樵蛹词乖谡w平衡狀態(tài)下也會(huì)發(fā)生振動(dòng)而引起力的變化,而時(shí)間平均液橋力才是宏觀上有意義的液橋力,因此后文將對(duì)時(shí)間平均液橋力進(jìn)行分析。圖6中時(shí)間平均液橋力為負(fù)值，說(shuō)明固相礦物受到液橋的吸引作用與宏觀液橋力的作用效應(yīng)一致。

在非飽和土力學(xué)中,可從宏觀角度直接求解Young-Laplace方程,計(jì)算球體(簡(jiǎn)化的粉、砂粒)間的理論液橋力[6-7,32]；然而,本文模擬的液橋尺寸處于納米尺度,基于宏觀連續(xù)概念的Young-Laplace方程能否用于計(jì)算納米尺度的液橋力還值得探索,而Amarasinghe等[12]、Song等[22]均發(fā)現(xiàn)基于Young-Laplace方程的理論液橋力明顯大于分子動(dòng)力學(xué)模擬的液橋力,但均未給出合理解釋。本節(jié)將對(duì)此加以分析，并提出納米尺度固相礦物間液橋力理論計(jì)算的修正方法。

圖7為液橋力計(jì)算示意圖,以上片固相礦物為例,其在z方向因左、右兩處表面張力而受到的液橋力(FTz)為

注:L′y為強(qiáng)吸附水在y方向的長(zhǎng)度。

FTz=-2lxσsinθ

(8)

上片受到的由負(fù)水壓力引起的z方向液橋力(FPz)為

FPz=uwlxLy

(9)

式中:Ly為毛細(xì)水在y方向的長(zhǎng)度。

結(jié)合式(6)和(7),上片受到的總液橋力大小為

(10)

Amarasinghe等[12]、Song等[22]的傳統(tǒng)計(jì)算方法為將L′y及d作為幾何參數(shù)代入式(10),導(dǎo)致理論計(jì)算的液橋力比模擬結(jié)果大,其誤差實(shí)質(zhì)是強(qiáng)吸附水由于受到固相礦物的強(qiáng)烈作用而不再滿足Young-Laplace控制方程。本文提出的修正方法是在式(10)中用排除強(qiáng)吸附水后的液橋尺寸(即Ly和D)來(lái)計(jì)算液橋力，此處尺寸由模擬的10個(gè)瞬時(shí)構(gòu)型取均值確定,雖有一定誤差但不影響結(jié)果。取上、下片固相礦物受到的液橋力均值作為每個(gè)算例的代表液橋力,傳統(tǒng)計(jì)算方法理論值、修正計(jì)算方法理論值和分子動(dòng)力學(xué)模擬值的對(duì)比結(jié)果如圖8所示，分子動(dòng)力學(xué)模擬值是指“固相礦物+強(qiáng)吸附水”原子組與“液橋”原子組之間的作用力總和。由圖8可以看出:傳統(tǒng)計(jì)算方法得到的液橋力明顯大于模擬結(jié)果,而修正方法計(jì)算值與模擬值基本吻合。盡管排除強(qiáng)吸附水后的液橋幾何尺寸為納米尺度,但可認(rèn)為其特性仍服從宏觀的Young-Laplace控制方程,因此本文提出的液橋力修正理論計(jì)算方法合理。

圖8還給出了式(9)計(jì)算所得的基質(zhì)吸力引起的液橋力。由圖8還可以看出:液橋力主要由基質(zhì)吸力貢獻(xiàn),而表面張力的貢獻(xiàn)較小,該結(jié)果與Song等[22]的模擬結(jié)果一致。

圖8 液橋力理論計(jì)算值與模擬值對(duì)比

由于液橋力與非飽和土有效應(yīng)力密切相關(guān),本文算例可從微觀上解釋非飽和土應(yīng)力與強(qiáng)度的主要試驗(yàn)規(guī)律。算例1—3的結(jié)果表明:保持毛細(xì)水體積不變而增大高嶺石片間距(孔隙比)將導(dǎo)致液橋力減小,這從微觀上解釋了相同含水量非飽和黏土的抗剪強(qiáng)度隨孔隙比增大而減小的現(xiàn)象。算例1、4、5的結(jié)果表明:當(dāng)片間距不變時(shí),含水量增大可引起液橋力的小幅增大,驗(yàn)證了張鵬程等[32]猜想的黏土在一定較低土體含水量范圍內(nèi)的“濕吸力”隨含水量增大而增大。算例1、6、7的結(jié)果表明:黏土片排列取向差異對(duì)液橋力影響很小,說(shuō)明從應(yīng)力與強(qiáng)度角度不需要考慮黏土片排列取向的差異。算例1、8、9的結(jié)果表明:礦物類型對(duì)液橋力影響較大,在相同條件下,黏土礦物間液橋力大于非黏土礦物間液橋力,加之黏土比表面積明顯大于粉土、砂土比表面積,使得相同條件下的非飽和黏土有效應(yīng)力和強(qiáng)度明顯大于非飽和粉土、砂土有效應(yīng)力和強(qiáng)度。而算例10、11的結(jié)果表明:黏土礦物與非黏土礦物構(gòu)成的異種礦物間液橋力接近,甚至低于兩種相應(yīng)同種礦物組合下液橋力的較小值,該現(xiàn)象是由異類礦物間液橋形態(tài)導(dǎo)致的。

綜上所述,液橋力與礦物類型、固相礦物間距、毛細(xì)水體積有關(guān),這是宏觀上造成非飽和土有效應(yīng)力和強(qiáng)度與土的礦物成分、密實(shí)度、含水量有關(guān)的微觀機(jī)制。

4 結(jié)論與展望

本文采用分子動(dòng)力學(xué)模擬研究了正長(zhǎng)石、α-石英和高嶺石3種礦物間毛細(xì)水形態(tài)和液橋力隨含水量、固相礦物間距等因素的變化規(guī)律。

1)在分子動(dòng)力學(xué)模擬中，采用CLAYFF力場(chǎng)能很好地再現(xiàn)固相礦物間毛細(xì)水彎液面形態(tài),采用毛細(xì)水平均密度分布圖獲取彎液面接觸角合理、可靠,所模擬的接觸角與試驗(yàn)結(jié)果基本吻合。

2)在微觀上,當(dāng)毛細(xì)水體積相同時(shí),彎液面接觸角隨固相礦物間距的增大而增大,基質(zhì)吸力隨之減小;當(dāng)固相礦物間距不變時(shí),毛細(xì)水體積的增大不會(huì)引起彎液面接觸角和基質(zhì)吸力的變化;高嶺石礦物的羥基面相對(duì)時(shí)的彎液面接觸角大于硅氧面相對(duì)時(shí)的彎液面接觸角,而前者基質(zhì)吸力小于后者基質(zhì)吸力;異種固相礦物間彎液面接觸角明顯大于任一單種固相礦物間彎液面接觸角。

3)固相礦物表面存在一層強(qiáng)吸附水,排除該強(qiáng)吸附水層后剩余的毛細(xì)水(即液橋)可認(rèn)為仍服從Young-Laplace控制方程;改進(jìn)液橋力計(jì)算方法,采用液橋形態(tài)幾何參數(shù)、基于Young-Laplace方程和彎液面圓弧假定計(jì)算的液橋力理論值與分子動(dòng)力學(xué)模擬值非常接近,液橋力主要由基質(zhì)吸力貢獻(xiàn)。

因此,分子動(dòng)力學(xué)模擬可用于分析非飽和土微觀尺度的力學(xué)規(guī)律，揭示宏觀力學(xué)特性的微觀機(jī)制,是一種強(qiáng)大的研究工具。本課題后續(xù)將探索固相礦物傾斜排列、固相礦物間距動(dòng)態(tài)變化、液橋斷裂與重建等更加復(fù)雜情形下的微觀力學(xué)規(guī)律。