鈉基、鈣基蒙脫石低含水率狀態下基質勢能的多元動力學研究

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

(1. 工業裝備結構分析國家重點實驗室，大連 116024；2. 大連理工大學工程力學系，大連 116024；3. 湖南大學土木工程學院，長沙 410082)

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

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

Morancho等[10]研究了重度擊實的膨潤土在高吸力變化條件下的力學行為。Lu等[11]研究了土體在高吸力范圍內的持水機理和滯回現象，提出高吸力范圍內吸附作用對吸力貢獻起著主導的作用。

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

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

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

1 基質吸力與自由能勢能的轉換

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

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

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

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

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

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

2 Metadynamics原理

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

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

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

3 模擬細節與方法

本研究采用國內外廣泛應用的Wyoming蒙脫石作為研究對象，其化學結構式為M0.75(Si7.75Al0.25)(Al3.5Mg0.5)O20(OH)4，其中M指代陽離子。Wyoming蒙脫石結構為單斜晶系，C2/m空間群，對稱型為L2PC，晶格常數α=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。取代規則是在硅氧四面體結構中，每32個Si中有一個被Al替代；鋁氧八面體中，每8個Al被一個Mg替代，且相鄰的原子不能同時被取代。類質同象產生的電荷由層間陽離子來平衡，故層間應有6個一價陽離子或3個二價陽離子。本次模擬涉及到的陽離子有一價的Na離子和二價的Ca離子。圖1為模擬體系的原子結構圖，紫色粒子為一價的Na離子，粉色粒子為二價的Ca離子。

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

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

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

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

4 結果與討論

4.1 自由能勢能變化與基質勢能結果

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

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

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

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

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

圖4是表1各項基質勢能隨含水量變化的總結曲線，兩種顏色的系列線分別代表鈉基蒙脫石、鈣基蒙脫石基質勢能隨含水量的變化曲線。可以看出，由前述基質吸力與自由能勢能轉換公式計算得到各含水量狀態下蒙脫石的基質吸力，能夠明顯發現，隨著含水量升高基質吸力變小，且呈現在低含水量范圍基質吸力急劇減小的特點。在含水率小于20%時，鈉基蒙脫石的基質吸力略小于鈣基蒙脫石，含水率超過此臨界點后鈉基蒙脫石基質吸力略大于鈣基蒙脫石。

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

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

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

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

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

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

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

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

4.2 模擬盒子對計算結果的影響

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

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

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

5 結論與展望

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

關于模擬體系設置的盒子大小：通過案例驗證發現c方向的大小會影響層間水分子躍過勢能壁壘達到穩定值時，與體系勢能最低點處的差值。進而影響采用Metadynamics方法計算基質勢能結果的準確性，因此，采用此方法時，需要考慮設置合適的體系模擬盒子大小。

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