江舒棋，趙紅華，張 超，張大帥，王小紅

(1. 工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗室，大連 116024；2. 大連理工大學(xué)工程力學(xué)系，大連 116024；3. 湖南大學(xué)土木工程學(xué)院，長沙 410082)

蒙脫石是膨脹土的重要組分之一，是膨脹土吸水膨脹性能的決定因素[1]。非飽和膨脹土的強(qiáng)度變化會導(dǎo)致滑坡、泥石流等地質(zhì)災(zāi)害的發(fā)生。膨脹土地基中的復(fù)雜樁-土相互作用，在地基吸水膨脹后對單樁扭矩承載特性的影響機(jī)理值得研究[2]。當(dāng)土體含水率較低時，毛細(xì)作用和短程吸附作用占據(jù)主導(dǎo)地位，即土體處于高吸力狀態(tài)。研究高吸力狀態(tài)下的持水特性(土水特征曲線)具有重要意義。非飽和土力學(xué)中最重要的基礎(chǔ)本構(gòu)關(guān)系之一就是土水特征曲線，它描述了土體基質(zhì)吸力與土體含水量之間的關(guān)系[3]。土水特征曲線目前基本上都是通過物理模型或室內(nèi)試驗得出[4]。

其中，基質(zhì)吸力的測量一直是被研究關(guān)注的實(shí)驗難題，如何準(zhǔn)確測量出土體的吸力值是研究的重點(diǎn)，目前的測量方法主要有壓力板法、濾紙法和張力計法等，也有學(xué)者對膨脹土采用電阻率結(jié)構(gòu)參數(shù)預(yù)測基質(zhì)吸力[5]。白福清等[6]介紹了濾紙法量測基質(zhì)吸力的原理，利用NaCl溶液建立“雙圈”No.203型濾紙的總吸力率定曲線，測試并給出了含水率介于16%～30%的膨脹壓實(shí)土總吸力和基質(zhì)吸力的土水特征曲線。劉星志[7]采用多種實(shí)驗方法，以非飽和紅土和花崗巖殘積土為研究對象，找出最適合的數(shù)學(xué)擬合模型建立全吸力范圍的土水特征曲線。國外的不少學(xué)者也發(fā)展了測量較高基質(zhì)吸力的不同方法，如采用吸附臺[8]、新型陶瓷板等[9]。

Morancho等[10]研究了重度擊實(shí)的膨潤土在高吸力變化條件下的力學(xué)行為。Lu等[11]研究了土體在高吸力范圍內(nèi)的持水機(jī)理和滯回現(xiàn)象，提出高吸力范圍內(nèi)吸附作用對吸力貢獻(xiàn)起著主導(dǎo)的作用。

Khorshidi等[12]考慮土壤顆粒與水分子相互作用，陽離子與水分子相互作用，提出土水特征曲線理論干燥模型來量化基質(zhì)吸力，預(yù)測的土體在零含水量時的最低基質(zhì)勢能(即最大基質(zhì)吸力455 MPa～1100 MPa)較為合理，得出結(jié)論：描述非飽和土力學(xué)和水力特性的本構(gòu)理論不僅要考慮表面張力，還要考慮范德華力和陽離子水化能。

文獻(xiàn)表明國內(nèi)外對于高吸力土體的持水特性和持水機(jī)理的研究目前尚處于起步階段，在高吸力測量上實(shí)現(xiàn)比較困難，這是由于極低含水量狀態(tài)下對基質(zhì)吸力的測量方法和試驗精度都提出了較高的挑戰(zhàn)。雖然有許多研究土水特征曲線的室內(nèi)試驗和預(yù)測方法，但在實(shí)際工程中由于非飽和土的高離散性，基質(zhì)吸力的結(jié)果與室內(nèi)試驗有較大差異。針對上述問題，許多研究者想要探索一條區(qū)別于物理模型和室內(nèi)實(shí)驗不同的計算模擬方法。近年來計算機(jī)的高速發(fā)展，使計算機(jī)數(shù)值模擬在各行各業(yè)的應(yīng)用中大放異彩，細(xì)觀和分子尺度模擬得到快速發(fā)展[13 ? 15]，有學(xué)者采用分子動力學(xué)方法對結(jié)合材料界面起裂的過程進(jìn)行模擬，提出相應(yīng)的界面破壞準(zhǔn)則[16]。分子動力學(xué)方法也被應(yīng)用于研究黏土的水-土相互作用和納米力學(xué)性質(zhì)[17 ? 19]。Zhang等[20]針對非飽和土的最低勢能(或最高吸力)，從熱動力學(xué)原理出發(fā)，提出了基于系統(tǒng)自由能的勢能計算理論框架，模擬了二氧化硅礦物晶體表面單層水分子膜吸附過程，得出二氧化硅礦物表面的最低勢能可以達(dá)到?2.0 GPa。這為計算機(jī)水平高速發(fā)展的現(xiàn)在提供了一種新的預(yù)測土水特征曲線的方法。

本文基于文獻(xiàn)[20]提出的基質(zhì)勢能(或基質(zhì)吸力)計算理論框架，運(yùn)用分子動力學(xué)方法和元動力學(xué)方法，模擬鈉基蒙脫石和鈣基蒙脫石吸附不同水分子數(shù)的過程，從分子尺度采用高斯函數(shù)疊加統(tǒng)計偏置勢能，得到不同吸附狀態(tài)下的體系自由能波動曲線，基于系統(tǒng)自由能的變化將自由能勢能轉(zhuǎn)化為基質(zhì)吸力，建立高吸力范圍內(nèi)的土水特征曲線，為鈉基蒙脫石和鈣基蒙脫石最高吸力值提供分子尺度上的參考值，比較分析了兩種蒙脫石基質(zhì)吸力在不同含水率范圍內(nèi)的細(xì)微差別；對比實(shí)驗方法和理論模型預(yù)測方法，為高吸力范圍內(nèi)膨潤土的土水特征曲線提供了分子尺度下研究的案例。

1 基質(zhì)吸力與自由能勢能的轉(zhuǎn)換

Zhang等[20]對二氧化硅表面最大基質(zhì)吸力的計算框架提供了吸力與自由能勢能轉(zhuǎn)換的理論依據(jù)，土的基質(zhì)勢能可以采用孔隙水與純水在單位體積下的自由能差值來表示，因此根據(jù)Noy-Meir等[21]的研究，基質(zhì)勢能Ψ(kPa或kJ/m3)可以表示為：

此外，土中水的化學(xué)勢μw可以表示為：

式中：F/(kcal/mol)為自由能(一般為Gibbs自由能)；nw為水分子的摩爾數(shù)；Xi代表所有其他的狀態(tài)變量(如溫度、壓力)，這些狀態(tài)變量隨系統(tǒng)狀態(tài)變化而變化。

根據(jù)式(1)和式(2)，可以得到基質(zhì)勢能新的表達(dá)式[12]：

式中： ΔF1為緊密吸附狀態(tài)與解除吸附狀態(tài)水的自由能差值； ΔF2為純水與解除吸附狀態(tài)水的自由能差值，也就是純水的水化自由能，根據(jù)文獻(xiàn)[22]，已知 ΔF2為6.3 kcal/mol。

因此，最重要的是計算緊密吸附狀態(tài)與解除吸附狀態(tài)水的自由能差值 ΔF1。自由能計算有很多種方法，包括傘狀采樣[23]、自由能擾動[24]和熱力學(xué)積分[25]以及多元動力學(xué)(Metadynamics)[26]等。

2 Metadynamics原理

多元動力學(xué)(Metadynamics)是Parrinello等[26]提出的一個計算多維自由能勢能面的方法，通過對研究體系不斷地施加額外的高斯型的偏置勢能，使目標(biāo)分子脫離自由能勢阱，推動研究體系的自由能勢能面從低自由能區(qū)向高自由能區(qū)域移動，探索整個反應(yīng)空間，達(dá)到避免重復(fù)采樣和增強(qiáng)采樣的目的，最后通過施加的偏置勢能函數(shù)來求解自由能勢能面數(shù)據(jù)[27 ? 28]。在n個自由度空間構(gòu)造的偏置勢能函數(shù)V(s,t)可表示為[25]：

式中：ω為高斯峰的高度；si,t為t時刻反應(yīng)坐標(biāo)第i維的值，時間t正比于加峰頻率； σsi為控制高斯峰寬度的變量。與反應(yīng)坐標(biāo)相關(guān)的偏置勢能V與自由能F轉(zhuǎn)換關(guān)系如下：

其中，C為只與時間有關(guān)的常數(shù)。Metadynamics構(gòu)建偏置勢函數(shù)需要選取合適的集成變量(原子間距離或配位數(shù))作為其反應(yīng)坐標(biāo)，反應(yīng)坐標(biāo)的合適選取對于結(jié)果的影響非常巨大。此外，從上述偏置勢能和自由能計算公式可以看出，模擬時間t趨于∞時，Metadynamics中添加的偏置勢能并不收斂于真實(shí)的自由能數(shù)據(jù)，而是在其附近震蕩，因此不能直接使用得到的自由能數(shù)據(jù)，而需要選擇使用自由能變化值作為后續(xù)計算數(shù)據(jù)。

3 模擬細(xì)節(jié)與方法

本研究采用國內(nèi)外廣泛應(yīng)用的Wyoming蒙脫石作為研究對象，其化學(xué)結(jié)構(gòu)式為M0.75(Si7.75Al0.25)(Al3.5Mg0.5)O20(OH)4，其中M指代陽離子。Wyoming蒙脫石結(jié)構(gòu)為單斜晶系，C2/m空間群，對稱型為L2PC，晶格常數(shù)α=90°，β=99°，γ=90°，a≈0.523 nm，b≈0.906 nm，c值隨層間含水量的變化而變化，在未吸附水分子時c≈0.96 nm，含1層水時c≈1.25 nm，含2層水時c≈1.55 nm，含3層水時c≈1.85 nm。本次研究的對象是由8個單晶胞按4a×2b×1c組成的超晶胞，a、b方向尺寸約為2.1×1.8 nm，空間群P1。取代規(guī)則是在硅氧四面體結(jié)構(gòu)中，每32個Si中有一個被Al替代；鋁氧八面體中，每8個Al被一個Mg替代，且相鄰的原子不能同時被取代。類質(zhì)同象產(chǎn)生的電荷由層間陽離子來平衡，故層間應(yīng)有6個一價陽離子或3個二價陽離子。本次模擬涉及到的陽離子有一價的Na離子和二價的Ca離子。圖1為模擬體系的原子結(jié)構(gòu)圖，紫色粒子為一價的Na離子，粉色粒子為二價的Ca離子。

圖1 吸附16個水分子的鈉基蒙脫石和鈣基蒙脫石晶胞結(jié)構(gòu)Fig.1 The structure of sodium montmorillonite and calcium montmorillonite adsorbing 16 water molecules

首先在Material studio軟件中按照上述結(jié)構(gòu)參數(shù)建立蒙脫石結(jié)構(gòu)模型，接下來分別建立吸附16個、32個、48個、64個、96個水分子的蒙脫石結(jié)構(gòu)模型，對應(yīng)的含水量在模擬體系中分別為4.9%、9.8%、14.7%、19.6%、29.5%。將上述各結(jié)構(gòu)模型數(shù)據(jù)導(dǎo)入到大型并行分子動力學(xué)軟件LAMMPS中，利用軟件對原子間相互作用勢等參數(shù)進(jìn)行設(shè)定。Cygan等[29]提出的黏土力場ClayFF對蒙脫石和高嶺石等黏土礦物結(jié)構(gòu)體系具有很好的可移植性，忽略了黏土礦物晶體中鋁氧八面體和硅氧四面體的共價鍵作用，這些粒子間只有庫侖力和范德華力的相互作用，共價鍵只存在于水分子和八面體羥基中。在ClayFF中，總勢能Etotal包括原子間庫侖能Ecoul、原子間范德華能Evdw以及成鍵原子間(水分子和羥基)的鍵伸縮能Ebond和鍵角彎曲能Eangle：

力場中水分子為SPC結(jié)構(gòu)模型，通過SHAKE算法控制水分子的鍵長與鍵角[30]。設(shè)置截斷距離1.25 nm，庫侖能采用Ewald加和長程力算法，范德華能采用LJ12-6勢。體系設(shè)置周期性邊界條件。采用NVT系綜控制體積和溫度，設(shè)置溫度300 K，采用Nose-Hoover算法控制溫度，Verlet速度算法進(jìn)行運(yùn)動方程積分，模擬的時間步長設(shè)置為1 fs。

在自由能勢能面計算部分，選擇水分子質(zhì)量中心到蒙脫石表面(001面)的距離D為集成變量，Metadynamics方法的計算精度和收斂速度很大程度上取決于偏置勢函數(shù)，也即決定偏置勢函數(shù)的參數(shù)高斯峰高度和寬度以及加峰頻率。因此，選擇合適的偏置勢參數(shù)具有重要的影響。本次模擬設(shè)置高斯峰的寬度為0.01 nm, 高斯峰的高度為1 kcal/mol，偏置因子為6，溫度300 K，加峰頻率間隔為1000個時間步，即每皮秒加峰一次。最后，每100次高斯峰堆積計算一次自由能。由于Metadynamics需要在非常大的時間尺度內(nèi)對自由能勢能面采樣才能獲得比較真實(shí)可靠的數(shù)據(jù)。根據(jù)多次擬合結(jié)果比較，設(shè)置模擬總時長為100 ns，并采集數(shù)據(jù)。

4 結(jié)果與討論

4.1 自由能勢能變化與基質(zhì)勢能結(jié)果

圖2和圖3分別為鈉基蒙脫石和鈣基蒙脫石自由能勢能隨水分子質(zhì)量中心距蒙脫石表面距離的變化曲線，不同顏色的系列線代表吸附不同水分子數(shù)的鈉基蒙脫石和鈣基蒙脫石。從圖中可以看出，在元動力學(xué)模擬初始階段，體系自由能勢能先逐漸減小，直至最低值；接著，隨水分子質(zhì)量中心距蒙脫石表面的距離增加，即水分子遠(yuǎn)離蒙脫石表面，水分子與蒙脫石表面相互作用勢增加，勢能逐漸增加。當(dāng)水分子質(zhì)量中心距蒙脫石表面距離大于15 ?時，4種含水量狀態(tài)下勢能均趨于定值，并在此定值附近振蕩。此時認(rèn)為蒙脫石表面對水分子的吸附作用可以忽略，即解吸狀態(tài)(desorbed state)[20]。

對比分別吸附16個水分子、32個水分子、48個水分子、64個水分子和96個水分子數(shù)的蒙脫石體系勢能曲線，可以發(fā)現(xiàn)：隨含水量增加，觀察到體系自由能值在減小。且隨著蒙脫石含水量的增加，達(dá)到的自由能勢能面波谷越平緩，說明含水量越高的蒙脫石體系越容易克服自由能勢壘，表明其狀態(tài)越穩(wěn)定，這與一直以來的實(shí)驗值和理論結(jié)果吻合。此外，從圖中可以看出不同含水量的蒙脫石的自由能變化曲線的波谷所在位置大致相同，該波谷對應(yīng)五種含水量狀態(tài)的蒙脫石吸附水分子的勢能最低點(diǎn)，即水分子緊密吸附狀態(tài)。將圖2和圖3中的相關(guān)數(shù)據(jù)統(tǒng)計整理并加以計算得到表1，可以看出，水分子距離蒙脫石表面的距離在9 ?附近時，體系達(dá)到勢能最低點(diǎn)的狀態(tài)。其中，對于鈉基蒙脫石，吸附16個水分子的蒙脫石結(jié)構(gòu)達(dá)到勢能最低點(diǎn)的橫坐標(biāo)的距離較其他含水量，距離蒙脫石表面最近。并且，任意含水量狀態(tài)下的鈣基蒙脫石達(dá)到勢能最低點(diǎn)的橫坐標(biāo)的距離均較鈉基蒙脫石稍遠(yuǎn)。

表1 不同含水量鈉基、鈣基蒙脫石基質(zhì)勢能與勢能最低狀態(tài)Table 1 The matric potential and the lowest potential energy state of sodium and calciummontmorillonites with different water content

圖2 鈉基蒙脫石吸附不同水分子數(shù)自由能勢能變化曲線Fig.2 Free-energy landscapes of sodium montmorillonite adsorption with different water molecules

圖3 鈣基蒙脫石吸附不同水分子數(shù)自由能勢能變化曲線Fig.3 Free-energy landscapes of calcium montmorillonite adsorption with different water molecules

圖4是表1各項基質(zhì)勢能隨含水量變化的總結(jié)曲線，兩種顏色的系列線分別代表鈉基蒙脫石、鈣基蒙脫石基質(zhì)勢能隨含水量的變化曲線?？梢钥闯觯汕笆龌|(zhì)吸力與自由能勢能轉(zhuǎn)換公式計算得到各含水量狀態(tài)下蒙脫石的基質(zhì)吸力，能夠明顯發(fā)現(xiàn)，隨著含水量升高基質(zhì)吸力變小，且呈現(xiàn)在低含水量范圍基質(zhì)吸力急劇減小的特點(diǎn)。在含水率小于20%時，鈉基蒙脫石的基質(zhì)吸力略小于鈣基蒙脫石，含水率超過此臨界點(diǎn)后鈉基蒙脫石基質(zhì)吸力略大于鈣基蒙脫石。

圖4 高吸力下鈉基、鈣基蒙脫石的土水特征曲線Fig.4 Soil-water characteristic curve of sodium and calcium montmorillonites at high suction range

Khorshidi等[12]主要研究了陽離子對土體基質(zhì)勢能的影響，并基于離子價態(tài)及陽離子與水分子的距離建立了基質(zhì)吸力的計算公式，與實(shí)驗值對比驗證了最大基質(zhì)吸力與離子類型的密切關(guān)系。據(jù)其理論公式，鈣離子相關(guān)礦物的最大基質(zhì)吸力約為1.11 GPa，而鈉離子相關(guān)礦物的最大基質(zhì)吸力約為0.555 GPa。此結(jié)果與本文得到的在低含水率狀態(tài)下，鈉基蒙脫石的基質(zhì)吸力小于鈣基蒙脫石的基質(zhì)吸力結(jié)果一致。

水分子在礦物表面的吸附不僅有范德華力相互作用，還有表面極化分子之間存在的靜電相互作用。因此，分子模擬的結(jié)果往往比較理想化，同時實(shí)驗過程中可能出現(xiàn)的人為的或環(huán)境造成的不可避免的誤差，往往會使結(jié)果偏小。

為保證成果的說服力，嘗試了多種理論數(shù)學(xué)模型擬合進(jìn)行驗證。由于Fredlund & Xing模型嚴(yán)格適用于吸力值在0 GPa～1 GPa的土水特征曲線[31]，本文采取分子模擬方法的特殊導(dǎo)致所得到的吸力值要大于1 GPa，因此不適宜采用Fredlund & Xing模型進(jìn)行擬合?，F(xiàn)選取土壤持水曲線van Genuchten模型[32]進(jìn)行曲線擬合驗證，表達(dá)式如下：

式中：θ為含水率；θr和θs分別為殘余含水量和飽和含水量；ψ為吸力值；α、n、m為經(jīng)驗擬合參數(shù)。

圖5和圖6分別為van Genuchten模型擬合曲線與模擬數(shù)據(jù)的對比，可以看出，擬合曲線較全面地覆蓋上文模擬得到的數(shù)據(jù)點(diǎn)，效果較好，增加了本文所采用方法計算土水特征曲線的合理性。

圖5 VG模型擬合曲線與鈉基蒙脫石的模擬數(shù)據(jù)Fig.5 VG model fitting curve and Na-MMT simulation data

圖6 VG模型擬合曲線與鈣基蒙脫石的模擬數(shù)據(jù)Fig.6 VG model fitting curve and Ca-MMT simulation data

4.2 模擬盒子對計算結(jié)果的影響

圖7展示了模擬時間相同、模擬盒子不同的64個水分子的鈉基蒙脫石自由能勢能變化曲線，模擬盒子在此處指代模擬體系的大小。藍(lán)色線代表模擬盒子c方向大小為58.272 ?時的自由能勢能變化曲線，紅色線代表模擬盒子c方向大小為68.272 ?時的自由能勢能變化曲線，黑色線代表模擬盒子c方向大小78.272 ?時的自由能勢能變化曲線，淺黃色線代表模擬盒子c方向大小為88.272 ?時的自由能勢能變化曲線。

從圖7中可以看出，相同體系在同樣的模擬時間內(nèi)，由于模擬盒子的大小設(shè)置不同，自由能勢能曲線并不一致。模擬盒子c方向大小在58.272 ?～78.272 ?范圍內(nèi)時，較大的體系自由能變化值更大，說明水分子達(dá)到解除吸附的狀態(tài)對模擬體系的大小有要求，這是因為在NVT系綜下，模擬盒子的大小是保持不變的，故在初始設(shè)置時要預(yù)留足夠的空間給水分子，以達(dá)到解除吸附的狀態(tài)。當(dāng)模擬盒子大小設(shè)置到一定值時，在本次模擬中指吸附64個水分子的鈉基蒙脫石c方向大小在78.272 ?及其以上時，體系自由能變化值不再隨模擬盒子的大小變化。但就達(dá)到平衡狀態(tài)時的曲線來看，c方向為88.272 ?時可以獲得更平穩(wěn)的效果。這一結(jié)果表明模擬體系的盒子大小選擇會在一定程度上影響最后的模擬效果，因此，選擇合適的模擬體系盒子大小是非常重要的。

圖7 模擬時間相同、模擬盒子不同的吸附64個水分子的鈉基蒙脫石自由能勢能變化曲線Fig.7 Free-energy landscapes of 64H2O-NaMMT in the same simulation time and different simulation box

5 結(jié)論與展望

本文基于文獻(xiàn)[12]提出的基質(zhì)吸力計算理論框架，從分子動力學(xué)模擬出發(fā)，結(jié)合應(yīng)用Metadynamics方法，模擬計算了鈉基蒙脫石和鈣基蒙脫石在低含水率狀態(tài)下的基質(zhì)吸力，建立了5%～30%含水率之間的土水特征曲線，曲線趨勢與前人結(jié)果較吻合，吸力大小值較前人實(shí)驗結(jié)果略偏大，得到的鈉基蒙脫石在含水率5%時的基質(zhì)吸力約為4.2 GPa，鈣基蒙脫石在含水率5%時的基質(zhì)吸力約為4.4 GPa。通過曲線對比，發(fā)現(xiàn)了鈉基蒙脫石與鈣基蒙脫石的基質(zhì)吸力均隨含水率的升高而指數(shù)下降，并且，在含水率小于20%時，鈉基蒙脫石的基質(zhì)吸力略小于鈣基蒙脫石，含水率超過此臨界點(diǎn)后鈉基蒙脫石基質(zhì)吸力略大于鈣基蒙脫石。

關(guān)于模擬體系設(shè)置的盒子大?。和ㄟ^案例驗證發(fā)現(xiàn)c方向的大小會影響層間水分子躍過勢能壁壘達(dá)到穩(wěn)定值時，與體系勢能最低點(diǎn)處的差值。進(jìn)而影響采用Metadynamics方法計算基質(zhì)勢能結(jié)果的準(zhǔn)確性，因此，采用此方法時，需要考慮設(shè)置合適的體系模擬盒子大小。

由于自由能的計算是非常消耗計算資源的工作，受計算條件的限制，目前所計算的案例效率較低、耗費(fèi)時間長，因此可以探索計算效率更加高效的方法來節(jié)約案例計算消耗時間，得到更加詳細(xì)的不同含水量密度體系基質(zhì)吸力模擬，建立更加完備的土水特征曲線與實(shí)驗數(shù)據(jù)對比，這對于非飽和土力學(xué)將是一個重要的研究契機(jī)。