楊亞帆 王建州 商翔宇 王濤 孫樹(shù)瑜

1) (中國(guó)礦業(yè)大學(xué) 深部巖土力學(xué)與地下工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，徐州 221116)

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

3) (中國(guó)地質(zhì)大學(xué)(武漢) 地球物理與空間信息學(xué)院，武漢 430074)

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

1 引言

蒙脫石在核廢料深部封存、二氧化碳封存以及頁(yè)巖氣開(kāi)發(fā)等應(yīng)用中發(fā)揮著重要作用.例如，蒙脫石具有很強(qiáng)的陽(yáng)離子交換能力，對(duì)核廢料等污染物質(zhì)的吸附能力突出，因此被用作核廢料深部封存的阻隔材料[1-3].當(dāng)?shù)叵滤c蒙脫石發(fā)生接觸，水分子會(huì)進(jìn)入蒙脫石微觀空隙中，受溫度和蒙脫石陽(yáng)離子種類等條件影響，蒙脫石會(huì)發(fā)生不同程度的膨脹[4，5].在深部環(huán)境中，蒙脫石的膨脹往往會(huì)受到巖石壁面的約束，從而對(duì)周圍環(huán)境產(chǎn)生壓力[6，7].由于核廢料在衰變過(guò)程會(huì)產(chǎn)生大量的熱量，產(chǎn)生的高溫條件會(huì)影響蒙脫石膨脹壓力特性，從而影響封存的穩(wěn)定性.因此，研究高溫下蒙脫石膨脹壓力的特性對(duì)優(yōu)化核廢料封存等過(guò)程至關(guān)重要.

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

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

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

2 計(jì)算方法

2.1 測(cè)量膨脹壓力的體系設(shè)置

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

圖1 研究飽和鈣蒙脫石膨脹壓力的分子體系;蒙脫石原子類型與顏色關(guān)系:橙:鈣離子，紅:氧，白:氫，黃:硅，藍(lán):鋁，粉:鎂;體系中的水分子不作具體展示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三個(gè)方向的尺寸分別為11.15，3.66 和28.00 nm.蒙脫石模型為懷俄明鈣蒙脫石，晶格模型分子式為 Ca0.375[Si7.75Al0.25][Al3.5Mg0.5]O20[OH]4.其表面電荷密度為—0.1245 C/m2.一個(gè)蒙脫石層在x，y和z三個(gè)方向的尺寸分別為0.66，3.66 和12.66 nm.在模擬過(guò)程中，通過(guò)在模擬過(guò)程中不更新蒙脫石層原子的位移來(lái)維持兩個(gè)蒙脫石層的晶面間距(d-spacing)d固定在1.40—4.00 nm 范圍內(nèi).研究的溫度條件為298，400 和500 K，所在溫度范圍在用于穩(wěn)定性分析的黏土水化相圖的溫度范圍內(nèi)[16].取決于溫度的不同，體系中水分子的個(gè)數(shù)在17700—24200 范圍內(nèi)以確保水環(huán)境區(qū)域在z方向的長(zhǎng)度大于3.00 nm.此外，為了減弱蒙脫石層受其在x方向的周期鏡像的影響，蒙脫石層與其x方向的鏡像的距離均大于7.0 nm (見(jiàn)圖1).

當(dāng)體系達(dá)到平衡狀態(tài)時(shí)，通過(guò)監(jiān)測(cè)每一個(gè)蒙脫石層在垂直方向的力計(jì)算其絕對(duì)值的平均值得出平均受力.膨脹壓力可以由平均受力除以蒙脫石層在yz平面的面積計(jì)算求得.

2.2 分子動(dòng)力學(xué)模擬細(xì)節(jié)

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

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

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

3 結(jié)果與討論

3.1 鈣蒙脫石的膨脹特性

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

圖2 飽和蒙脫石在298 K 和5 MPa 下膨脹壓力曲線;鈉蒙脫石的數(shù)據(jù)源于文獻(xiàn)[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 在不同溫度和晶面間距下，分子動(dòng)力學(xué)模擬計(jì)算的鈣蒙脫石的膨脹壓力及預(yù)測(cè)誤差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 水化作用

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

圖3 在298 K 和5 MPa 下，水在孤立鈣蒙脫石表面(d →∞)的平衡密度分布(a)和鈣蒙脫石層間水在不同晶面間距下的平衡密度分布((b)—(l));在((b)—(l))圖中，虛線(點(diǎn)線)為平移后的(a)圖中的密度分布;x=0 對(duì)應(yīng)孔的中心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.

水化力對(duì)膨脹壓力的影響可以通過(guò)水的密度分布來(lái)解釋.當(dāng)d從1.95 nm(圖3(g))降低到1.70 nm(圖3(d))，孔中的水的層數(shù)從五層降為四層.在孔中部水層的擠出過(guò)程中，首先由于空間約束作用孔中的水層受到擠壓，對(duì)應(yīng)的膨脹壓力隨d的降低而增加.當(dāng)d等于1.80 nm(圖3(e))時(shí)，擠壓作用最強(qiáng)，對(duì)應(yīng)的的膨脹壓力達(dá)到峰值.繼續(xù)降低d，中部的水層被逐步擠出，在強(qiáng)水化力的作用下，靠近兩個(gè)壁面的四層水層會(huì)維持水層結(jié)構(gòu)，因此中部區(qū)域會(huì)形成負(fù)壓并對(duì)壁面產(chǎn)生一個(gè)收縮力，對(duì)應(yīng)的膨脹壓力隨d的降低而降低.當(dāng)d等于1.70 nm(圖3(d))時(shí)，收縮力最強(qiáng)，對(duì)應(yīng)的的膨脹壓力達(dá)到谷值.類似的，當(dāng)孔中的水層數(shù)從六層降低到五層，對(duì)應(yīng)的膨脹壓力峰值和谷值分別位于d=2.10 和1.95 nm.此時(shí)的峰谷值的絕對(duì)值要遠(yuǎn)小于對(duì)應(yīng)于水層數(shù)從五層降為四層的峰谷值.這是因?yàn)殡S著孔的d的增加插入孔中部的水受到的水化作用越弱，中心部分的水分子的性質(zhì)更像自由水.對(duì)于d大于2.30 nm的情況，水化力進(jìn)一步減弱.值得注意的是，當(dāng)d從1.60(圖3(c))降低到1.40 nm(圖3(b))時(shí)，水層數(shù)從四層降低到二層，對(duì)應(yīng)的膨脹壓力單調(diào)增加，未發(fā)現(xiàn)谷值.這是因?yàn)榭拷诿娴乃畬邮艿綇?qiáng)烈的水化作用，具體表現(xiàn)在對(duì)應(yīng)的水層的峰值要明顯高于較大d時(shí)的情況.

3.1.2 雙電層作用

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

圖4 在298 K 和5 MPa 下，鈣離子在孤立鈣蒙脫石表面(d →∞)的平衡密度分布(a)和鈣蒙脫石層間鈣離子在不同晶面間距下的平衡密度分布((b)—(l));在((b)—(l))圖中，虛線(點(diǎn)線)為平移后的(a)圖中的密度分布;x=0 對(duì)應(yīng)孔的中心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.

對(duì)于由帶負(fù)電荷的固體表面與只含陽(yáng)離子的水溶液所組成的夾縫孔，PB 方程預(yù)測(cè)滲透效應(yīng)比靜電效應(yīng)更強(qiáng)，因而產(chǎn)生一個(gè)正的膨脹壓力[30]

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

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

其中K為常數(shù)，其表達(dá)式為[30]

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

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

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

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

圖2 中的黑虛線展示的是從PB 方程計(jì)算得到的雙電層力.預(yù)測(cè)的雙電層力在所選取的d的范圍內(nèi)均為正值.當(dāng)d從3.00 nm 增加到4.00 nm時(shí)，雙電層力從0.36 MPa 降低到0.17 MPa.而在相同的d范圍下，分子動(dòng)力學(xué)模擬預(yù)測(cè)的膨脹壓力接近0.因?yàn)閼讯砻髅擅撌谋砻骐姾奢^高以及鈣離子化合價(jià)較高，PB 方程所高估的雙電層力可能是源于對(duì)離子關(guān)聯(lián)效應(yīng)的忽略[31，35].根據(jù)強(qiáng)關(guān)聯(lián)理論，離子關(guān)聯(lián)效應(yīng)會(huì)降低雙電層力，且離子關(guān)聯(lián)越強(qiáng)對(duì)雙電層力的減弱效應(yīng)越明顯[36，37].離子關(guān)聯(lián)效應(yīng)的大小可以用耦合參數(shù)(coupling parameter)來(lái)描述:

其中l(wèi)Be2/(4πε0εbkBT)為Bjerrum 長(zhǎng)度.當(dāng)耦合參數(shù)Ξ遠(yuǎn)小于1 時(shí)，PB 方程可以準(zhǔn)確地描述雙電層力.隨著Ξ增加，PB 方程對(duì)雙電層力的高估更加明顯.例如，當(dāng)Ξ3.06，雙電層力被高估約30%[31];當(dāng)Ξ20，分子模擬預(yù)測(cè)的雙電層力可為負(fù)值[37]，而PB 方程預(yù)測(cè)的雙電層力不可能成為負(fù)值.表3列出了計(jì)算的結(jié)果，在298 K 和5 MPa 下，鈣蒙脫石體系的耦合參數(shù)為30.38，遠(yuǎn)大于1.這說(shuō)明了對(duì)于鈣蒙脫石體系，離子關(guān)聯(lián)效應(yīng)對(duì)雙電層力起著決定作用.

3.1.3 膨脹自由能

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

其中A為蒙脫石面積;d0=1.60 nm 為選取的參考晶面間距.圖5 中黑實(shí)線展示的是鈣蒙脫石在298 K 和5 MPa 下的膨脹自由能曲線.在d=1.68，1.92 和2.27 nm 處存在局部最低值，對(duì)應(yīng)的孔中的水的層數(shù)為四、五和六層.其中當(dāng)d=1.92 nm 時(shí)，膨脹的能壘最大，說(shuō)明孔中水層數(shù)為五層時(shí)，蒙脫石較為穩(wěn)定.這與實(shí)驗(yàn)觀測(cè)到的鈣蒙脫石在膨脹到1.9 nm 晶面間距時(shí)較穩(wěn)定的現(xiàn)象相吻合[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 對(duì)比鈉蒙脫石膨脹特性

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

圖2 對(duì)比了298 K 和5 MPa 下，鈉蒙脫石和鈣蒙脫石的膨脹壓力.當(dāng)d較小時(shí)，水化力對(duì)膨脹壓力影響較大.可以發(fā)現(xiàn)鈣蒙脫石膨脹壓力的振幅要大于鈉蒙脫石膨脹壓力的振幅，并且鈣蒙脫石水化力對(duì)膨脹壓力起主導(dǎo)作用的d的范圍(d小于3.00 nm)比鈉蒙脫石的(d小于2.30 nm)更大.這些現(xiàn)象的發(fā)生主要是因?yàn)殁}離子的水化能更低[40]，這導(dǎo)致了更強(qiáng)的水化作用.當(dāng)d較大時(shí)，雙電層力對(duì)膨脹壓力影響較大.此時(shí)，分子模擬預(yù)測(cè)的鈣蒙脫石的膨脹壓力要小于鈉蒙脫石，這與Akinwunmi等[14]的模擬結(jié)果相符.因?yàn)楹雎粤穗x子關(guān)聯(lián)效應(yīng)，PB 方程預(yù)測(cè)的膨脹壓力均大于分子模擬的結(jié)果.鈉蒙脫石在298 K 和5 MPa 下的耦合參數(shù)值為3.81[1]，遠(yuǎn)小于鈣蒙脫石的耦合參數(shù)，這說(shuō)明鈣蒙脫石的離子關(guān)聯(lián)效應(yīng)更強(qiáng)，對(duì)雙電層力的減弱效果也更強(qiáng).這解釋了在較大d時(shí)，鈣蒙脫石的膨脹壓力要低于鈉蒙脫石膨脹壓力這一現(xiàn)象.

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

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

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

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

3.2.1 高溫下的水化作用

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

圖7 在5 MPa 下，水在孤立鈣蒙脫石表面(d →∞)的平衡密度分布(a)和鈣蒙脫石層間水在不同晶面間距下的平衡密度分布((b)—(l));在((b)—(l))圖中，虛線(點(diǎn)線)為平移后的(a)圖中的密度分布.x=0 對(duì)應(yīng)孔的中心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.

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

3.2.2 高溫下的雙電層作用

當(dāng)d較大時(shí)，雙電層力決定了膨脹壓力，高溫降低膨脹壓力，使膨脹壓力變?yōu)樨?fù)值.這一規(guī)律與實(shí)驗(yàn)研究的結(jié)果相符合[12，13].此外，隨著d的增加，膨脹壓力趨近于0.因?yàn)楹雎粤穗x子關(guān)聯(lián)效應(yīng)，PB 方程預(yù)測(cè)的雙電層力與分子模擬計(jì)算的結(jié)果有較大差距:雙電層力為正值，表現(xiàn)為膨脹力，且當(dāng)d較大時(shí)，受溫度的影響較小(見(jiàn)圖6).為了估計(jì)溫度對(duì)離子關(guān)聯(lián)效應(yīng)的影響，計(jì)算了不同溫度下的耦合參數(shù).表3 列出的計(jì)算結(jié)果表明，耦合參數(shù)隨著溫度的升高而增加，當(dāng)溫度為500 K 時(shí)，耦合參數(shù)升高至80.84.這說(shuō)明離子關(guān)聯(lián)效應(yīng)隨溫度的升高而增加，因此，其對(duì)雙電層力的減弱效應(yīng)隨溫度升高而增強(qiáng).這解釋了溫度對(duì)分子模擬計(jì)算的膨脹壓力的影響.類似地，鈉蒙脫石的雙電層力也受溫度升高而降低[1].例如，d=2.60 nm 及5 MPa 下，溫度從298 K 升高到500 K，鈉蒙脫石的耦合參數(shù)從3.81 升高到10.13，對(duì)應(yīng)的膨脹壓力從1.09 MPa降低到0.11 MPa[1].因?yàn)殁}蒙脫石的耦合參數(shù)遠(yuǎn)大于鈉蒙脫石的耦合參數(shù)，這導(dǎo)致了鈣蒙脫石在高溫下的膨脹壓力為負(fù)值，表現(xiàn)為收縮力.

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

3.2.3 高溫下的膨脹自由能

圖8 展示的是不同溫度下鈣蒙脫石在不同晶面間距下的膨脹自由能曲線.從圖8 可以明顯地看出，隨著溫度的升高膨脹自由能曲線向下遷移.值得注意的是當(dāng)溫度提高到400 K 以上時(shí)，自由能曲線的局部最低點(diǎn)的個(gè)數(shù)從三個(gè)降低到兩個(gè).溫度為400 K 時(shí)，自由能局部最低值位于d=1.93 和2.18 nm 處.溫度為500 K 時(shí)，自由能局部最低值位于d=1.95 和2.11 nm 處.298 K 時(shí)，位于d=1.68 nm 附近的局部最低點(diǎn)在高溫下消失，這源于3.2.1 中介紹的當(dāng)d=1.70 nm 時(shí)，溫度對(duì)膨脹壓力的反常影響(即高溫下膨脹壓力的谷值為正).向下遷移的自由能曲線及消失的局部最低點(diǎn)共同表明，當(dāng)溫度升高時(shí)，鈣蒙脫石更容易膨脹到d=1.9 nm 處.此外，隨溫度從298 K 上升到500 K，位于2.2 nm 附近的局部最低點(diǎn)的位置的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āng)d較大時(shí)，溫度對(duì)自由能曲線有較大的影響.當(dāng)溫度從298 K 升高到400 K 時(shí)，自由能隨d增加的變化規(guī)律從基本上不隨d的改變而變化轉(zhuǎn)變成隨d的增加而增加.繼續(xù)升溫到500 K，自由能隨d增加的增幅最大.這是源于3.2.2 中介紹的溫度升高所導(dǎo)致的離子關(guān)聯(lián)效應(yīng)的強(qiáng)化，因此，負(fù)的膨脹壓力降低(收縮力增強(qiáng)).這表明當(dāng)鈣蒙脫石膨脹到2.2 nm 附近，升高的溫度會(huì)限制蒙脫石繼續(xù)膨脹.

4 結(jié)論

本文使用分子動(dòng)力學(xué)模擬研究了5 MPa 和298—500 K，1.40—4.00 nm 晶面間距的一系列飽和鈣蒙脫石的膨脹壓力特性，基于水化效應(yīng)、雙電層效應(yīng)和離子關(guān)聯(lián)效應(yīng)等模型推演膨脹壓力隨溫度與d的變化規(guī)律，并與相應(yīng)的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比.主要結(jié)論有以下幾點(diǎn):

1) 在298 K 和5 MPa 下，分子模擬預(yù)測(cè)飽和鈣蒙脫石在晶面間距為1.9 nm 處較穩(wěn)定，且鈣蒙脫石比鈉蒙脫石更難膨脹到較大的晶面間距，這些模擬結(jié)果與文獻(xiàn)中實(shí)驗(yàn)觀測(cè)的結(jié)果相符.將兩種蒙脫石膨脹程度的差異歸因于鈣蒙脫石的離子關(guān)聯(lián)效應(yīng)遠(yuǎn)大于鈉蒙脫石所導(dǎo)致的更低的膨脹壓力.

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

3) 當(dāng)晶面間距大于約3.0 nm 時(shí)，飽和鈣蒙脫石的膨脹壓力主要由雙電層力決定.根據(jù)強(qiáng)關(guān)聯(lián)理論，溫度從298 K 增至500 K，離子關(guān)聯(lián)效應(yīng)增強(qiáng)，對(duì)應(yīng)的鈣蒙脫石膨脹壓力降低，此外，雙電層力對(duì)膨脹壓力起主導(dǎo)作用的晶面間距范圍增加至少0.7 nm.較高溫時(shí)，膨脹壓力表現(xiàn)為收縮力并阻礙膨脹.

4) 有別于分子模擬對(duì)于離子關(guān)聯(lián)效應(yīng)的精確描述，連續(xù)化的Poisson-Boltzmann 方程因忽略了離子關(guān)聯(lián)效應(yīng)，在溫度于298—500 K 范圍內(nèi)，預(yù)測(cè)的膨脹壓力均為膨脹力，從而無(wú)法準(zhǔn)確表達(dá)鈣蒙脫石膨脹壓力的變化規(guī)律.

以上結(jié)論揭示了高溫下鈣蒙脫石的膨脹特性及相關(guān)機(jī)理，有助于優(yōu)化相關(guān)黏土材料在核廢料處理等涉及到高溫條件的過(guò)程中的應(yīng)用.此外，本項(xiàng)研究?jī)H涉及蒙脫石與純水體系，未探討鹽溶液對(duì)蒙脫石膨脹特性的影響.實(shí)驗(yàn)研究表明鹽溶液對(duì)蒙脫石的膨脹特性狀態(tài)有較大影響[10，11]，為探究鹽溶液的影響機(jī)制，計(jì)劃開(kāi)展受鹽溶液影響下蒙脫石膨脹特性的分子模擬研究.