高溫下鈣蒙脫石膨脹特性的分子動力學模擬*

楊亞帆 王建州 商翔宇 王濤 孫樹瑜

1) (中國礦業大學 深部巖土力學與地下工程國家重點實驗室，徐州 221116)

2) (阿卜杜拉國王科技大學 物理科學與工程學部，圖瓦 23955-6900)

3) (中國地質大學(武漢) 地球物理與空間信息學院，武漢 430074)

高溫下蒙脫石的膨脹特性在核廢料深部封存、二氧化碳封存及頁巖氣開發等應用中有著重要影響，但相關機理尚不明確.本工作使用分子動力學模擬為技術手段計算5 MPa 和298—500 K 等條件下，1.40—4.00 nm晶面間距(d)的一系列飽和鈣蒙脫石的膨脹壓力.以模擬所得的數值結果為依據，基于水化效應、雙電層效應和離子關聯效應等模型推演膨脹壓力隨溫度與d 的變化規律，并與相應的實驗數據進行對比.模擬結果表明，當d 較小時，因為高溫會弱化水化力的強度，鈣蒙脫石膨脹壓力震蕩的幅度降低，同時水化力作用的d 的范圍減小.當d 較大時，因為高溫強化離子關聯效應，膨脹壓力降低，同時雙電層力的作用的d 的范圍增加.在較高溫度和較大d 時，膨脹壓力為收縮力，阻礙膨脹.這些膨脹壓力的變化規律與前期鈉蒙脫石體系的研究類似.然而，通過對比兩種蒙脫石體系的模擬結果，發現兩種體系存在顯著的差異—鈣蒙脫石比鈉蒙脫石更難膨脹到較大的d.此模擬結果與前人實驗觀測的結果相符.我們進一步將此差異歸于鈣蒙脫石的離子關聯效應要遠大于鈉蒙脫石.有別于分子模擬中對于離子關聯效應的精確描述，連續化的Poisson-Boltzmann方程因為忽略了離子關聯效應，從而無法表達出與兩種體系模擬結果都相吻合的膨脹壓力變化規律.

1 引言

蒙脫石在核廢料深部封存、二氧化碳封存以及頁巖氣開發等應用中發揮著重要作用.例如，蒙脫石具有很強的陽離子交換能力，對核廢料等污染物質的吸附能力突出，因此被用作核廢料深部封存的阻隔材料[1-3].當地下水與蒙脫石發生接觸，水分子會進入蒙脫石微觀空隙中，受溫度和蒙脫石陽離子種類等條件影響，蒙脫石會發生不同程度的膨脹[4，5].在深部環境中，蒙脫石的膨脹往往會受到巖石壁面的約束，從而對周圍環境產生壓力[6，7].由于核廢料在衰變過程會產生大量的熱量，產生的高溫條件會影響蒙脫石膨脹壓力特性，從而影響封存的穩定性.因此，研究高溫下蒙脫石膨脹壓力的特性對優化核廢料封存等過程至關重要.

根據層間陽離子種類的不同，常見的蒙脫石有鈉蒙脫石和鈣蒙脫石.文獻中存在較多關于這兩種蒙脫石的膨脹特性以及溫度對膨脹壓力影響的實驗研究[8-13].鈉蒙脫石在與液態水接觸后可膨脹至原體積的20 倍[8，9]，而鈣蒙脫石則膨脹至1.9 nm的晶面間距后難以繼續膨脹[10，11].升高溫度，鈉蒙脫石的膨脹壓力增加，而鈣蒙脫石的膨脹壓力降低[12，13].這些結果表明，鈉蒙脫石和鈣蒙脫石的膨脹壓力特性之間存在較大差異，然而，導致這些差異的微觀機理尚不明確.

分子模擬被較為廣泛地應用于蒙脫石的膨脹特性及相關機理研究[1，14-22].Akinwunmi 等[14]使用分子動力學模擬方法，通過測量附加在蒙脫石上彈簧的形變計算了300 K 和1×105Pa 下鈉蒙脫石和鈣蒙脫石的膨脹壓力，模擬結果表明:當晶面間距較小時，鈣蒙脫石的膨脹壓力要大于鈉蒙脫石的膨脹壓力;當晶面間距較大時，鈣蒙脫石的膨脹壓力要小于鈉蒙脫石的膨脹壓力.Akinwunmi 等[15]使用了同樣的模擬設置研究了等體積條件下，溫度對鈉蒙脫石的膨脹壓力的影響.他們發現高溫會增加膨脹壓力.但是，Akinwunmi 等[14，15]使用的分子模擬設置無法準確描述較小晶面間距下水化力的震蕩分布[1].然而，只有準確地描述水化力的震蕩分布才能準確計算蒙脫石的平衡晶面間距[1，17，18]，因此Akinwunmi 等[14，15]采用的模擬設置無法預測蒙脫石的穩定膨脹狀態.Honorio 等[16]使用巨正則蒙特卡洛(GCMC)分子模擬方法研究了溫度對鈉蒙脫石膨脹壓力的影響.GCMC 方法可以相對準確地計算較小晶面間距下水化作用.他們發現當晶面間距較小時，高溫會降低鈉蒙脫石膨脹壓力的振幅.然而，GCMC 方法的誤差較大，無法準確地描述較大晶面間距下的膨脹壓力.近期，Yang等[1]改進了分子動力學模擬中測量膨脹壓力的體系設置，新的設置可以準確預測不同晶面間距下的膨脹壓力.此外，他們研究了高溫下鈉蒙脫石的膨脹壓力特性:當晶面間距較小時，高溫會弱化水化力，從而降低膨脹壓力的振幅;當晶面間距較大時，高溫會增強離子關聯效應，使膨脹壓力降低.盡管如此，文獻[1]未對鈉蒙脫石的膨脹自由能進行分析，也未對比鈉蒙脫石與鈣蒙脫石膨脹特性的差異.此外，文獻[1]中缺少高溫下鈣蒙脫石的膨脹特性及機理的分子模擬研究.

針對上述研究的不足之處，本文使用分子動力學模擬方法研究鈣蒙脫石的膨脹特性及機理，采用文獻[1]中的模擬設置計算5 MPa 和298—500 K等條件下，1.40—4.00 nm 晶面間距的一系列飽和鈣蒙脫石的膨脹壓力，對比鈣蒙脫石的膨脹特性與文獻[1]報道的鈉蒙脫石膨脹特性.我們計算的蒙脫石的膨脹壓力可以準確地預測蒙脫石的穩定晶面間距，與實驗測量的數據吻合[10，11]，而Akinwunmi等[14，15]報道的膨脹壓力曲線因存在較大誤差無法預測蒙脫石的穩定晶面間距.同時，基于水化效應、雙電層效應和離子關聯效應等模型推演膨脹壓力隨溫度與晶面間距的變化規律.此外，還對比了經典Poisson-Boltzmann(PB)方程預測的雙電層力與分子模擬計算的膨脹壓力.

2 計算方法

2.1 測量膨脹壓力的體系設置

本工作采用了文獻[1]提出的改進的測量膨脹壓力的體系設置，這一方法可以準確地測量不同晶面間距的蒙脫石的膨脹壓力，現做簡要描述.如圖1所示，模擬體系的中部為飽和的兩層蒙脫石層(saturated clay).蒙脫石層左右兩側被水環境(water environment)包裹.為了防護鈣離子擴散到水環境中，蒙脫石層邊緣處設置了兩個僅與鈣離子作用的隱式的壁面(implicit wall).水環境外側被兩個原子活塞(piston)包裹.其中，左側活塞固定不動，右側活塞被施加的外力作用以控制體系的環境壓力.原子活塞由一層鋪滿模擬盒子yz平面的原子組成，每個原子與水分子發生短程作用.此外，活塞的每個原子受到額外的施加在z方向的外力fext-PA/N，其中P為環境壓力、A為模擬盒子在xy平面的面積以及N為右側活塞包含的原子數目.模擬盒子采用周期性邊界條件，活塞外側為真空區域(vacuum)用于減弱體系與z方向周期鏡像的作用.

圖1 研究飽和鈣蒙脫石膨脹壓力的分子體系;蒙脫石原子類型與顏色關系:橙:鈣離子，紅:氧，白:氫，黃:硅，藍:鋁，粉:鎂;體系中的水分子不作具體展示Fig.1.Molecular system for studying the swelling pressure of Ca-montmorillonite in water.The color code for clay atoms:orange:Ca ion，red:O，white:H，yellow:Si，cyan:Al，pink:Mg.Water molecules are not explicitly shown.

模擬盒子在x，y和z三個方向的尺寸分別為11.15，3.66 和28.00 nm.蒙脫石模型為懷俄明鈣蒙脫石，晶格模型分子式為 Ca0.375[Si7.75Al0.25][Al3.5Mg0.5]O20[OH]4.其表面電荷密度為—0.1245 C/m2.一個蒙脫石層在x，y和z三個方向的尺寸分別為0.66，3.66 和12.66 nm.在模擬過程中，通過在模擬過程中不更新蒙脫石層原子的位移來維持兩個蒙脫石層的晶面間距(d-spacing)d固定在1.40—4.00 nm 范圍內.研究的溫度條件為298，400 和500 K，所在溫度范圍在用于穩定性分析的黏土水化相圖的溫度范圍內[16].取決于溫度的不同，體系中水分子的個數在17700—24200 范圍內以確保水環境區域在z方向的長度大于3.00 nm.此外，為了減弱蒙脫石層受其在x方向的周期鏡像的影響，蒙脫石層與其x方向的鏡像的距離均大于7.0 nm (見圖1).

當體系達到平衡狀態時，通過監測每一個蒙脫石層在垂直方向的力計算其絕對值的平均值得出平均受力.膨脹壓力可以由平均受力除以蒙脫石層在yz平面的面積計算求得.

2.2 分子動力學模擬細節

模擬工作基于LAMMPS[23]開源代碼.蒙脫石的力場參數為ClayFF[24].水分子采用SPC 模型[25].水分子的鍵和角由SHAKE 算法[26]固定.體系采用的勢能函數為

其中rij為i原子與j原子之間的距離、εij和σij分別為i原子與j原子之間能量參數和尺度參數、qi和qj分別為i原子與j原子的電荷以及ε0為絕對介電常數，詳細參數見表1.不同類型原子間的Lennard-Jones 參數由Lorentz-Berthelot 混合規則計算:σij(σi+σj)/2，.這些力場參數已經在課題組的前期工作中進行了驗證[1，17，18].短程范德華作用的截斷半徑為1.2 nm.庫倫作用的短程部分的截斷半徑為1.2 nm，在截斷半徑外的長程部分由PPPM 方法[27]計算，其受力精度設置為1×10—5.模擬系綜采用的是NPT 系綜.溫度控制方法為Nosé-Hoover 法[28]，阻尼系數為100 個時間步長.除了右側活塞加了外力來控制水環境的壓力，沒有加任何控壓條件.采用步長為2 fs 的velocity Verlet 算法[29]來計算原子運動軌跡.體系的平衡模擬時間為至少20 ns.當體系達到平衡后，繼續模擬的50—88 ns 被用來計算相關性質.

表1 水、蒙脫石和離子的Lennard-Jones 參數和電荷Table 1.Lennard-Jones parameters and partial charges of water，montmorillonite，and ion.

3 結果與討論

3.1 鈣蒙脫石的膨脹特性

圖2 展示的是在298 K 和5 MPa 下，飽和蒙脫石(MMT)在1.40—4.00 nm 晶面間距(d)范圍內的膨脹壓力的變化曲線(數據見表2).當d小于3.00 nm 時，鈣蒙脫石的膨脹壓力隨著d的增加震蕩衰減.膨脹壓力的峰值位于d=1.40，1.80，2.10和2.60 nm 處，對應的壓力值分別為577.94，22.75，3.31 和0.84 MPa.膨脹壓力的谷值位于d=1.70，1.95 和2.30 處，對應的壓力值分別為—6.80，—4.02和—0.53 MPa.當d大于3.00 nm 時，鈣蒙脫石的膨脹壓力降低至0 左右.在納米尺度下，親水帶電荷的固體表面的膨脹壓力主要由水化力和雙電層力決定[30，31].水化力在較小d時對膨脹壓力起主導作用，而雙電層力在較大d時對膨脹壓力起主導作用[30，31].膨脹壓力的震蕩分布是典型的水化作用的結果[30]，因此，當d小于3.00 nm 左右時，膨脹壓力主要取決于水化力.當d大于3.00 nm 左右，膨脹壓力曲線無震蕩表現，此時雙電層力決定了此處的膨脹壓力.

圖2 飽和蒙脫石在298 K 和5 MPa 下膨脹壓力曲線;鈉蒙脫石的數據源于文獻[1]Fig.2.Swelling pressure curves of saturated montmorillonite at 298 K and 5 MPa.Data for Na-montmorillonite is taken from Ref.[1].

表2 在不同溫度和晶面間距下，分子動力學模擬計算的鈣蒙脫石的膨脹壓力及預測誤差Table 2.Swelling pressures of Ca-montmorillonite at different temperatures and d-spacings and the corresponding standard deviations obtained from molecular dynamics simulations.

3.1.1 水化作用

水化作用可以通過水在蒙脫石層所形成的孔中的密度分布研究.圖3 展示的是不同d的孔中水分子在垂直于蒙脫石層方向的一維密度分布曲線.在所選取的d范圍內，孔中形成不同數量的水層.當d從1.40 nm 增加到3.00 nm，對應的孔中水的層數從二層增加到八層.當d=3.00 nm，孔中部的水的密度分布與自由水的密度分布幾乎相同，說明水化作用已經較小.對于孤立的蒙脫石表面(見圖3(a))，水化作用在其表面產生了四層水層，這也表明當孔中水層數超過八層時，位于孔中心部分的水的性質和自由水相似，水化作用幾乎無影響.此外，當d較大時，靠近蒙脫石表面的第一層和第二層水層的密度分布的形狀隨著d的變化幾乎不變，這說明蒙脫石對前兩層水有較強的水化作用.然而，靠近蒙脫石表面的第三層和第四層水層的密度分布的形狀隨著d的變化發生較大的變化.這些水層的密度分布隨d的變化主要受減弱的水化作用、孔中水層間的相互作用以及離子密度分布的變化等因素的影響.

圖3 在298 K 和5 MPa 下，水在孤立鈣蒙脫石表面(d →∞)的平衡密度分布(a)和鈣蒙脫石層間水在不同晶面間距下的平衡密度分布((b)—(l));在((b)—(l))圖中，虛線(點線)為平移后的(a)圖中的密度分布;x=0 對應孔的中心Fig.3.Equilibrium density profiles of water near one quasi-isolated clay surface (a) and inside Ca-montmorillonite pore with various d-spacings ((b)—(l)) at 298 K and 5 MPa.In ((b)—(l))，the shifted density profile from (a) are plotted as dashed(dotted) lines.x=0 corresponds to the center of the pore.

水化力對膨脹壓力的影響可以通過水的密度分布來解釋.當d從1.95 nm(圖3(g))降低到1.70 nm(圖3(d))，孔中的水的層數從五層降為四層.在孔中部水層的擠出過程中，首先由于空間約束作用孔中的水層受到擠壓，對應的膨脹壓力隨d的降低而增加.當d等于1.80 nm(圖3(e))時，擠壓作用最強，對應的的膨脹壓力達到峰值.繼續降低d，中部的水層被逐步擠出，在強水化力的作用下，靠近兩個壁面的四層水層會維持水層結構，因此中部區域會形成負壓并對壁面產生一個收縮力，對應的膨脹壓力隨d的降低而降低.當d等于1.70 nm(圖3(d))時，收縮力最強，對應的的膨脹壓力達到谷值.類似的，當孔中的水層數從六層降低到五層，對應的膨脹壓力峰值和谷值分別位于d=2.10 和1.95 nm.此時的峰谷值的絕對值要遠小于對應于水層數從五層降為四層的峰谷值.這是因為隨著孔的d的增加插入孔中部的水受到的水化作用越弱，中心部分的水分子的性質更像自由水.對于d大于2.30 nm的情況，水化力進一步減弱.值得注意的是，當d從1.60(圖3(c))降低到1.40 nm(圖3(b))時，水層數從四層降低到二層，對應的膨脹壓力單調增加，未發現谷值.這是因為靠近壁面的水層受到強烈的水化作用，具體表現在對應的水層的峰值要明顯高于較大d時的情況.

3.1.2 雙電層作用

通過分析離子密度分布和對比PB 方程預測的結果來研究雙電層作用.圖4 展示的是不同d的孔中鈣離子在垂直于蒙脫石層方向的一維密度分布曲線.鈣離子吸附到蒙脫石表面附近形成內層水合物和外層水合物，這與文獻報道的結果相吻合[1，17，18，32].內層水合物的陽離子與固體表面間不存在水分子，而外層水合物的陽離子完全被水分子包圍[33].靠近蒙脫石表面的第一個鈣離子峰的形狀隨著d的改變變化不大，這表明內層水合物與表面緊密結合.遠離蒙脫石表面的鈣離子會形成外層水合物，并且其密度分布隨著d的變化發生較大的改變.總體來說，雖然部分鈣離子形成了內層水合物，孔的中部仍然存在較多鈣離子，這導致了孔中心區域有較高的離子濃度.因為水環境中不含離子，滲透效應會產生排斥的雙電層力[30].但是滲透效應會被源于帶負電的固體表面與帶正電的孔中鹽溶液之間的靜電效應減弱[30].使用分子模擬對每個效應進行單獨研究較為困難，因此，我們使用經典的PB 方程研究雙電層力.

圖4 在298 K 和5 MPa 下，鈣離子在孤立鈣蒙脫石表面(d →∞)的平衡密度分布(a)和鈣蒙脫石層間鈣離子在不同晶面間距下的平衡密度分布((b)—(l));在((b)—(l))圖中，虛線(點線)為平移后的(a)圖中的密度分布;x=0 對應孔的中心Fig.4.Equilibrium density profiles of Ca ion near one quasi-isolated clay surface (a) and inside Ca-montmorillonite pore with various d-spacings ((b)—(l)) at 298 K and 5 MPa.In ((b)—(l))，the shifted density profile from (a) are plotted as dashed(dotted) lines.x=0 corresponds to the center of the pore.

對于由帶負電荷的固體表面與只含陽離子的水溶液所組成的夾縫孔，PB 方程預測滲透效應比靜電效應更強，因而產生一個正的膨脹壓力[30]

其中w為孔的寬度，kB為玻爾茲曼常數，T為溫度以及ρion，w(x0)是陽離子在孔中心平面的密度.一維PB 方程有如下形式[30]:

其中ψ為靜電勢;q為陽離子化合價以及εb為水的介電常數.對于固體表面帶負電，且孔中僅有陽離子的情況，一維PB 方程的解為[30]

其中K為常數，其表達式為[30]

根據電荷平衡的邊界條件，K可從下式求得[30]

其中σ為固體的表面電荷.

對于現有的體系，PB 方程需要四個參數:表面電荷(σ—0.1245C/m2)、鈣離子化合價(q2)、孔寬度w以及孔中溶劑的介電常數εb.孔的寬度的取值為水密度分布圖中最靠近固體表面的峰值的距離加上一個水分子的直徑(近似為0.25 nm).介電常數由自由水的介電常數近似[1，31，34]，我們在前期工作[1]中已對不同條件下的自由水的介電常數進行了計算，具體見表3.

表3 水在5 MPa 和不同溫度下的介電常數和耦合參數Table 3.Water dielectric constants and coupling parameters at 5 MPa and various temperatures.

圖2 中的黑虛線展示的是從PB 方程計算得到的雙電層力.預測的雙電層力在所選取的d的范圍內均為正值.當d從3.00 nm 增加到4.00 nm時，雙電層力從0.36 MPa 降低到0.17 MPa.而在相同的d范圍下，分子動力學模擬預測的膨脹壓力接近0.因為懷俄明蒙脫石的表面電荷較高以及鈣離子化合價較高，PB 方程所高估的雙電層力可能是源于對離子關聯效應的忽略[31，35].根據強關聯理論，離子關聯效應會降低雙電層力，且離子關聯越強對雙電層力的減弱效應越明顯[36，37].離子關聯效應的大小可以用耦合參數(coupling parameter)來描述:

其中lBe2/(4πε0εbkBT)為Bjerrum 長度.當耦合參數Ξ遠小于1 時，PB 方程可以準確地描述雙電層力.隨著Ξ增加，PB 方程對雙電層力的高估更加明顯.例如，當Ξ3.06，雙電層力被高估約30%[31];當Ξ20，分子模擬預測的雙電層力可為負值[37]，而PB 方程預測的雙電層力不可能成為負值.表3列出了計算的結果，在298 K 和5 MPa 下，鈣蒙脫石體系的耦合參數為30.38，遠大于1.這說明了對于鈣蒙脫石體系，離子關聯效應對雙電層力起著決定作用.

3.1.3 膨脹自由能

單位面積的蒙脫石膨脹自由能可通過膨脹壓力Ps對d的積分計算[38，39]:

其中A為蒙脫石面積;d0=1.60 nm 為選取的參考晶面間距.圖5 中黑實線展示的是鈣蒙脫石在298 K 和5 MPa 下的膨脹自由能曲線.在d=1.68，1.92 和2.27 nm 處存在局部最低值，對應的孔中的水的層數為四、五和六層.其中當d=1.92 nm 時，膨脹的能壘最大，說明孔中水層數為五層時，蒙脫石較為穩定.這與實驗觀測到的鈣蒙脫石在膨脹到1.9 nm 晶面間距時較穩定的現象相吻合[10，11].

圖5 飽和蒙脫石在298 K 和5 MPa 下的膨脹自由能曲線Fig.5.Swelling free energy curves of saturated montmorillonite at 298 K and 5 MPa.

3.1.4 對比鈉蒙脫石膨脹特性

課題組的前期工作[1]在相同條件下研究了鈉蒙脫石體系的膨脹壓力.鈣蒙脫石的膨脹壓力與鈉蒙脫石的膨脹壓力特性有一定相似性，但也存在較大差異.下面將做具體討論.

圖2 對比了298 K 和5 MPa 下，鈉蒙脫石和鈣蒙脫石的膨脹壓力.當d較小時，水化力對膨脹壓力影響較大.可以發現鈣蒙脫石膨脹壓力的振幅要大于鈉蒙脫石膨脹壓力的振幅，并且鈣蒙脫石水化力對膨脹壓力起主導作用的d的范圍(d小于3.00 nm)比鈉蒙脫石的(d小于2.30 nm)更大.這些現象的發生主要是因為鈣離子的水化能更低[40]，這導致了更強的水化作用.當d較大時，雙電層力對膨脹壓力影響較大.此時，分子模擬預測的鈣蒙脫石的膨脹壓力要小于鈉蒙脫石，這與Akinwunmi等[14]的模擬結果相符.因為忽略了離子關聯效應，PB 方程預測的膨脹壓力均大于分子模擬的結果.鈉蒙脫石在298 K 和5 MPa 下的耦合參數值為3.81[1]，遠小于鈣蒙脫石的耦合參數，這說明鈣蒙脫石的離子關聯效應更強，對雙電層力的減弱效果也更強.這解釋了在較大d時，鈣蒙脫石的膨脹壓力要低于鈉蒙脫石膨脹壓力這一現象.

圖5 對比了鈉蒙脫石和鈣蒙脫石的膨脹自由能曲線，兩者有較大差異.鈉蒙脫石在d=1.60 nm附近存在局部最低值.當d大于1.70 nm，鈉蒙脫石的膨脹自由能隨d的增加單調降低.這表明如果鈉蒙脫石膨脹到d=1.70 nm 后，它會自發的繼續膨脹.而鈣蒙脫石的自由能曲線存在多個局部最低值(見上節)，且當d大于3.00 nm 時，自由能曲線幾乎不隨d的改變發生變化.因此，當d較大時，鈣蒙脫石相對于鈉蒙脫石更難膨脹，這與其他學者通過實驗方法觀測到的現象相吻合[10，11].

3.2 高溫下鈣蒙脫石的膨脹特性

圖6 展示的是在5 MPa 環境壓力和不同溫度(298，400 和500 K)下不同晶面間距的飽和鈣蒙脫石的膨脹壓力(數據見表2).根據d的不同，溫度對膨脹壓力的影響較為復雜.下面將根據d的范圍不同做具體討論.

圖6 飽和鈣蒙脫石在不同溫度和5 MPa 下膨脹壓力曲線Fig.6.Swelling pressure curves of saturated Ca-montmorillonite at various temperatures and 5 MPa.

3.2.1 高溫下的水化作用

當d較小時，水化力決定了膨脹壓力，高溫會降低膨脹壓力的振幅.具體而言，若低溫時膨脹壓力為正，高溫會降低膨脹壓力;若低溫時膨脹壓力為負，高溫會增加膨脹壓力.這一現象在d位于峰值和谷值位置時較為明顯.這和課題組的前期工作[1]研究的高溫下鈉蒙脫石體系的結果類似，同時和其他學者在不同體系的模擬結果相吻合[16，41].高溫降低膨脹壓力的振幅是因為高溫會弱化水化效應.高溫對的水化力的弱化可以通過水的密度分布來理解.圖7 展示的是不同溫度下水分子在不同d的孔中的密度分布.溫度對水分子的密度分布有較大的影響.總體而言，高溫會使固體表面的水分子脫附，降低孔中水密度分布的峰值.孔中水的水化結構受到破壞，對應的水化力降低.然而，當d=1.70 nm 時(298 K 的谷值位置)，隨著溫度從298 K升高到400 K，膨脹壓力從—6.80 MPa 升高到10.74 MPa.繼續升溫到500 K，膨脹壓力增加至23.48 MPa(見圖6 中插圖).此現象不符合高溫降低膨脹壓力振幅這一規律.這個特例可以從圖7(c)所展示的密度分布來理解.觀察圖7(c)可以明顯地發現，隨著溫度從298 K 升高到400 K，孔中水的層數從四層增加到了五層.繼續升溫到500 K，孔中心部分的水層更明顯.因為高溫會弱化水化效應，靠近固體壁面的水層的峰值隨溫度升高而降低.水層峰值的降低意味著其結構更不穩定，因此，高溫下孔的中心區域更容易形成新的水層，因而增加了水化力.

圖7 在5 MPa 下，水在孤立鈣蒙脫石表面(d →∞)的平衡密度分布(a)和鈣蒙脫石層間水在不同晶面間距下的平衡密度分布((b)—(l));在((b)—(l))圖中，虛線(點線)為平移后的(a)圖中的密度分布.x=0 對應孔的中心Fig.7.Equilibrium density profiles of water near one quasi-isolated clay surface (a) and inside Ca-montmorillonite pore with various d-spacings ((b)—(l)) at 5 MPa.In ((b)—(l))，the shifted density profiles from (a) are plotted as dashed(dotted) lines.x=0 corresponds to the center of the pore.

同時，隨著溫度的升高，膨脹壓力震蕩分布的d的范圍變小.如圖6 所示，當溫度為298 K 時，膨脹壓力的震蕩范圍在1.40—3.00 nm.而溫度在400和500 K 時，膨脹壓力的震蕩范圍在1.40—2.60 nm.這一現象可以通過研究水的密度分布的規律來理解.對于孤立的單個蒙脫石表面(圖7(a))，隨著溫度的升高，靠近固體表面區域形成的水層的層數從四層降低到三層.此外，遠離固體表面的水層形狀發生較大的改變.298 K 下靠近固體表面的第三層水層的厚度要明顯大于更高溫度下同樣位置的水層厚度.這些現象的發生是因為隨著溫度的升高，水分子的動能更大，因此更不容易被蒙脫石表面束縛.溫度對孔中水密度分布的影響與對孤立的單個蒙脫石表面附近水密度分布的影響相似(對比圖7(a)和圖7(i)).降低的孔中水層層數和水層厚度共同解釋了高溫下膨脹壓力震蕩范圍的降低.

3.2.2 高溫下的雙電層作用

當d較大時，雙電層力決定了膨脹壓力，高溫降低膨脹壓力，使膨脹壓力變為負值.這一規律與實驗研究的結果相符合[12，13].此外，隨著d的增加，膨脹壓力趨近于0.因為忽略了離子關聯效應，PB 方程預測的雙電層力與分子模擬計算的結果有較大差距:雙電層力為正值，表現為膨脹力，且當d較大時，受溫度的影響較小(見圖6).為了估計溫度對離子關聯效應的影響，計算了不同溫度下的耦合參數.表3 列出的計算結果表明，耦合參數隨著溫度的升高而增加，當溫度為500 K 時，耦合參數升高至80.84.這說明離子關聯效應隨溫度的升高而增加，因此，其對雙電層力的減弱效應隨溫度升高而增強.這解釋了溫度對分子模擬計算的膨脹壓力的影響.類似地，鈉蒙脫石的雙電層力也受溫度升高而降低[1].例如，d=2.60 nm 及5 MPa 下，溫度從298 K 升高到500 K，鈉蒙脫石的耦合參數從3.81 升高到10.13，對應的膨脹壓力從1.09 MPa降低到0.11 MPa[1].因為鈣蒙脫石的耦合參數遠大于鈉蒙脫石的耦合參數，這導致了鈣蒙脫石在高溫下的膨脹壓力為負值，表現為收縮力.

此外，隨著溫度的升高，雙電層力對膨脹壓力起主導作用的d的范圍增加.例如，當溫度為298 K 時，雙電層力對膨脹壓力起主導作用的d的范圍為大于3.00 nm;當溫度為500 K 時，d在2.30 nm 處為水化力的峰值.若不存在雙電層力，水化力的峰值應為正值.然而，此時的膨脹壓力值為—1.71 MPa，膨脹壓力的正負轉變表明在500 K下，雙電層力在d在2.30 nm 處已經起到主導作用，雙電層力對膨脹壓力起主導作用的d的范圍增加了至少0.7 nm.

3.2.3 高溫下的膨脹自由能

圖8 展示的是不同溫度下鈣蒙脫石在不同晶面間距下的膨脹自由能曲線.從圖8 可以明顯地看出，隨著溫度的升高膨脹自由能曲線向下遷移.值得注意的是當溫度提高到400 K 以上時，自由能曲線的局部最低點的個數從三個降低到兩個.溫度為400 K 時，自由能局部最低值位于d=1.93 和2.18 nm 處.溫度為500 K 時，自由能局部最低值位于d=1.95 和2.11 nm 處.298 K 時，位于d=1.68 nm 附近的局部最低點在高溫下消失，這源于3.2.1 中介紹的當d=1.70 nm 時，溫度對膨脹壓力的反常影響(即高溫下膨脹壓力的谷值為正).向下遷移的自由能曲線及消失的局部最低點共同表明，當溫度升高時，鈣蒙脫石更容易膨脹到d=1.9 nm 處.此外，隨溫度從298 K 上升到500 K，位于2.2 nm 附近的局部最低點的位置的2.27 nm降低到2.11 nm.這源于3.2.1 中所介紹的高溫下孔中水層厚度的降低.

圖8 飽和鈣蒙脫石在不同溫度和5 MPa 下的膨脹自由能曲線Fig.8.Swelling free energy curves of saturated Ca-montmorillonite at various temperatures and 5 MPa.

當d較大時，溫度對自由能曲線有較大的影響.當溫度從298 K 升高到400 K 時，自由能隨d增加的變化規律從基本上不隨d的改變而變化轉變成隨d的增加而增加.繼續升溫到500 K，自由能隨d增加的增幅最大.這是源于3.2.2 中介紹的溫度升高所導致的離子關聯效應的強化，因此，負的膨脹壓力降低(收縮力增強).這表明當鈣蒙脫石膨脹到2.2 nm 附近，升高的溫度會限制蒙脫石繼續膨脹.

4 結論

本文使用分子動力學模擬研究了5 MPa 和298—500 K，1.40—4.00 nm 晶面間距的一系列飽和鈣蒙脫石的膨脹壓力特性，基于水化效應、雙電層效應和離子關聯效應等模型推演膨脹壓力隨溫度與d的變化規律，并與相應的實驗數據進行對比.主要結論有以下幾點:

1) 在298 K 和5 MPa 下，分子模擬預測飽和鈣蒙脫石在晶面間距為1.9 nm 處較穩定，且鈣蒙脫石比鈉蒙脫石更難膨脹到較大的晶面間距，這些模擬結果與文獻中實驗觀測的結果相符.將兩種蒙脫石膨脹程度的差異歸因于鈣蒙脫石的離子關聯效應遠大于鈉蒙脫石所導致的更低的膨脹壓力.

2) 當晶面間距小于約3.0 nm 時，飽和鈣蒙脫石的膨脹壓力主要由水化力決定.溫度從298 K增至500 K，水化力的強度減弱，對應的鈣蒙脫石膨脹壓力震蕩的幅度降低，此外，水化力對膨脹壓力起主導作用的晶面間距范圍減小約0.4 nm.

3) 當晶面間距大于約3.0 nm 時，飽和鈣蒙脫石的膨脹壓力主要由雙電層力決定.根據強關聯理論，溫度從298 K 增至500 K，離子關聯效應增強，對應的鈣蒙脫石膨脹壓力降低，此外，雙電層力對膨脹壓力起主導作用的晶面間距范圍增加至少0.7 nm.較高溫時，膨脹壓力表現為收縮力并阻礙膨脹.

4) 有別于分子模擬對于離子關聯效應的精確描述，連續化的Poisson-Boltzmann 方程因忽略了離子關聯效應，在溫度于298—500 K 范圍內，預測的膨脹壓力均為膨脹力，從而無法準確表達鈣蒙脫石膨脹壓力的變化規律.

以上結論揭示了高溫下鈣蒙脫石的膨脹特性及相關機理，有助于優化相關黏土材料在核廢料處理等涉及到高溫條件的過程中的應用.此外，本項研究僅涉及蒙脫石與純水體系，未探討鹽溶液對蒙脫石膨脹特性的影響.實驗研究表明鹽溶液對蒙脫石的膨脹特性狀態有較大影響[10，11]，為探究鹽溶液的影響機制，計劃開展受鹽溶液影響下蒙脫石膨脹特性的分子模擬研究.