基于分子動力學模擬的鈣基蒙脫石表面潤濕性研究

項林語， 李長冬,2*， 李浩林， 孟杰， 黃德崴

(1.中國地質大學(武漢) 工程學院， 武漢 430074； 2.湖北巴東地質災害國家野外科學觀測研究站， 武漢 430074)

黏土礦物作為地表土體的重要組成部分，對其結構性質具有十分關鍵的影響。在降雨入滲和地表徑流等作用下，親水性的黏土礦物吸水膨脹，容易導致土體結構破壞。當土體中的水含量進一步增加時，由于強親水性的黏土礦物表面的水化層增厚，顆粒間的排斥力和空間位阻效應會導致顆粒難以碰撞聚結[1]，從而在水中呈現為穩定的分散懸浮狀態，不易沉降分離[2]。而河流和雨水的水動力較強，易攜帶著大量的泥沙流失，進而導致水土流失、土地荒漠化、水庫淤積和旱澇等災害頻發。因此，研究水在黏土礦物表面的潤濕性對進一步了解土體結構破壞和水土流失等現象至關重要。

而蒙脫石作為土體中主要的黏土礦物之一，具有較強的陽離子交換能力和水化膨脹特性，因此其與水的作用機制廣受關注。崔家慶等[3]開展了不同摻雜比例的蒙脫石的等溫吸水試驗，通過對相關吸附參數的分析，揭示了土體遇水弱化的機理。劉國強等[4]基于含蒙脫石組分的煤泥水的過濾試驗，分析了不同抑制劑對蒙脫石的水化膨脹特性的影響和作用機理，為煤泥水的處理提供了方法。鐘宜航等[5]對提純后的蒙脫石開展了吸附試驗，研究了接觸時間、溫度、固液比、初始濃度、pH和離子強度等因素對蒙脫石吸附能力的影響，并討論了不同的等溫吸附模型對于該吸附過程的適用性。為了進一步探究蒙脫石水化膨脹機制，黃緩緩[6]通過蒙脫石自然沉降試驗和過濾試驗，分析了蒙脫石在不同種類離子溶液中的體積膨脹程度和過濾速度，以分析金屬陽離子種類對蒙脫石水化膨脹的抑制性能。陶源清等[7]對浸泡在不同濃度的NaCl溶液中的蒙脫石樣品進行了熱重分析，通過對導數熱重分析(derivative thermogravimetry, DTG)曲線的對比，研究了鹽濃度對蒙脫石的去水化作用的影響，并分析了其作用機制。盡管已有大量的研究對蒙脫石的親水性進行了分析，但針對蒙脫石礦物表面與水分子之間原子尺度的吸附性機制仍有待進一步探究。

由于蒙脫石的粒徑極其微小，采用傳統的方法單獨開展試驗研究比較困難，且受多種因素影響導致結果存在較大誤差。目前，分子動力學模擬已廣泛用于研究分子或原子水平的試驗過程[8]，為探究分子間的相互作用機制提供了方法。采用分子動力學模擬研究了水分子在鈣基蒙脫石表面的潤濕性，提供了一種能夠精確計算蒙脫石表面的接觸角的方法，通過對不同時刻液滴形態快照、接觸角、徑向分布函數(radial distribution function, RDF)和原子密度分布的分析，從原子尺度解釋了蒙脫石表面親水的原因，以期為土體結構破壞和水土流失等現象的研究提供理論依據。

1 模型構建與模擬方法

1.1 模型構建

蒙脫石是一種典型的2∶1型層狀鋁硅酸鹽礦物。晶體結構由兩層硅氧四面體和一層夾在其間的鋁氧八面體組成，硅氧四面體和鋁氧八面體通過氧原子連接。鈣基蒙脫石晶體屬于單斜晶系，C2/m空間群，對稱型為L2PC，表面電荷密度為-0.124 5 C/m2。3個晶軸方向的軸長分別為：a=0.523 nm，b= 0.906 nm，c= 1.530 nm；3個軸間(a和b、a和c、b和c)的夾角分別為：α=90°、β=99°、γ=90°[9]。模型分子式為

Ca0.375[Si7.75Al0.25][Al3.5Mg0.5]O20[OH]4·nH2O[10]。

為了模擬液滴在鈣基蒙脫石表面的潤濕過程，選取了干燥的Mt (001)表面。將模型上層真空層厚度設置為100 ?，建立了Mt (001)表面液滴的初始結構，直徑約為82 ?，并構建了209.2 ? × 217.44 ?× 201 ?的仿真模擬盒。在模擬的初始狀態，將液滴放置于基底表面的正上方。如圖1所示。

在蒙脫石表面潤濕性的分子動力學模擬中，分子數量的選擇會影響模擬的效果：分子數量過多會降低計算速度，而分子數目過少則會影響液滴的形狀，甚至難以形成穩定的輪廓[11]。因此，本次模擬選擇了合適的分子數量建模，在完成分子動力學馳豫過程使得系統穩定之后，體系中包含了9 171個水分子，液滴輪廓較為清晰穩定。

紅色球體、藍色球體分別代表氫原子和氧原子;基底上綠色的球體 代表鈣離子圖1 鈣基蒙脫石分子模型Fig.1 Molecular model of Ca-montmorillonite

1.2 模擬計算方法

使用 LAMMPS[12]軟件開展了分子動力學模擬，采用周期性邊界條件(periodic boundary condition, PBC)，原子之間的力可以跨越邊界。這種邊界的優點是可以用少量的分子模擬較大的液滴，從而大大節省了計算資源，并且可以保證模型中有明顯的固-液-氣三相線和接觸角。此外，PPP周期性邊界條件還可以避免線張力對接觸角的影響[13]。

采用正則系綜(canonical ensemble, NVT)[14]對體系進行分子動力學模擬，即粒子數N恒定，體積V恒定，溫度T恒定。采用擴展的簡單點電荷(extended simple point charge, SPC/E)模型模擬該系統中的水分子，并計算了水分子能量和溫度等物理參數的變化。在SPC/E模型中，假定水分子為剛性分子，且O—H鍵長固定為1 ?，H—O—H鍵角固定為109.47°[15]，此模型可以再現水分子在界面處的物理動態行為，而且水的密度、擴散系數等物理參數與實驗數據的一致性較好。使用Nosé-Hoover恒溫器[16]模擬300 K的溫度，采用SHAKE算法計算模型中的鍵和角度。所有相互作用的截斷半徑設置為9 ?，截斷半徑之外的長程靜電相互作用采用基于傅里葉的Ewald求和方法：PPPM方法(particle particle particle mesh)來計算，大大降低了計算的誤差，精度可達10-4。根據Lorentz Berthelot組合規則確定不同原子和離子相互作用的L-J(Lennard-Jones)勢參數。選取Velocity-Verlet算法[17]求解原子運動方程，模擬總步數為5×106步，設置步長為1.0 fs，足以使系統達到相對穩定狀態。在模擬過程中，基底始終被固定在初始位置。此外，在納米尺度下，重力對液滴的影響較小，因此研究中忽略了重力的影響。

采用CLAYFF力場來描述原子間的相互作用，其被廣泛用于描述黏土礦物-水界面的物理過程。選擇L-J勢來描述范德華相互作用[18]，然后用靜電力對其進行擴充，其勢能函數表達式為

(1)

式(1)中：i和j為任意兩個原子；qi和qj分別為i和j的電荷；rij為i和j之間的距離；C為庫倫相互作用常數；σij為位勢為零時(即平衡位置)i和j之間的有限距離；εij為勢阱深度；ε0為相對介電常數。

力場參數如表1[19]所示。采用OVITO軟件包實現三維模型可視化，可以觀察液滴形狀的變化，獲得分子的坐標參數。

表1 力場參數[19]Table 1 Force field parameters[19]

2 結果與分析

2.1 體系平衡狀態

當體系中的原子之間相互作用時，體系的溫度和能量等參數逐漸趨于穩定，當這些參數的數值變化范圍在波動容許的誤差范圍內時，可以將此作為評價熱力學平衡狀態的依據[20]。如圖2所示，在 5 ns 的范圍內，溫度始終在300 K的設定值上下波動，最大誤差不超過1.3%。勢能在1 ns前急劇下降，幅度可達到總下降幅度的69.4%，然后勢能降低的速率逐漸減慢，在4.5 ns后，體系勢能平均變化量不超過1.0%。

分析模擬過程中，液滴的形態變化也是確定平衡狀態的一種方法。圖3展示了不同時刻液滴在鈣基蒙脫石表面的潤濕情況的正視圖。隨著模擬時間的推移，液滴的結構逐漸被破壞，使得水分子在蒙脫石基底表面擴散。水分子的動能降低，減少了鈣離子和水中氧的結合。而蒙脫石表面的氧捕獲了水分子中的氫，從而限制了納米液滴在蒙脫石表面的擴散。因此，在1 ns之前，液滴形態變化較大，接觸角變化也較為迅速，而在4 ns之后，液滴在蒙脫石表面的形態變化較不明顯，接觸角值變化不大。綜上所述，這些參數在4.5 ns之后的變化幅度都較小。因此，在模擬的終止時間5 ns時，該體系已經達到了熱力學平衡狀態。

圖2 物理參數隨時間變化曲線Fig.2 Curve of physical parameters with time

圖3 不同時刻下的液滴形態Fig.3 Droplet morphology at different times

2.2 接觸角計算

接觸角是表征黏土礦物表面的潤濕性的最直接的參數，其被定義為氣-液界面在氣-液-固三相交點處所作切線與固-液交界線的夾角。當水在固體表面的接觸角小于90°時，通常認為該固體表面是親水的；反之，當接觸角大于90°時，則認為固體表面是疏水的[21]。

當系統達到平衡后，將液滴的質心坐標表示為(x,y,z)。如圖4(a)所示，基于原點坐標(x,y, 0)，將空間直角坐標系中的點用柱狀坐標(r,z)來表示。z= 0 ?對應于Mt (001)基底表面，z軸指向基底外表面的法線方向，r為點和z軸(即徑向)之間相應的距離。在網格劃分中，Δr和Δz分別為r方向和z方向上的網格微元，仿真箱被劃分為Δr= 1 ?，Δz= 1 ?的圓柱箱，如圖4(b)所示。

計算仿真箱中點的密度，得到液滴的密度分布等值線圖[圖5(a)]。可以明顯地看出，液滴密度等高線的曲線形態近似為一個圓，且氣-液界面存在明顯的過渡層。根據等值線圖所顯示的液滴密度輪廓線，選擇ρ(r,z) = 0.50 g/cm3的輪廓線作為液滴與空氣的交界面，并提取等高線的點集數據。

z和r分別為基底外表面的法線方向和徑向的數值；Δr和Δz 分別為r方向和z方向上的網格微元圖4 柱坐標系示意圖Fig.4 Schematic diagram of cylindrical coordinate system

ρ為液滴密度圖5 接觸角計算Fig.5 Contact angle calculation

如圖5(b)所示，采用最小二乘法對氣-液界面的點集進行擬合。擬合的公式為圓方程，可表示為

(2)

式(2)中：a′和b′分別為圓心坐標在r和z方向上的值；Rcir為擬合圓的半徑。

通過對氣-液界面的最小二乘擬合，可以確定a′、b′和Rcir的值。此外，在擬合曲線時剔除了距蒙脫石表面小于11.45 ?的數據點，以消除界面上因吸附層等導致密度顯著波動的因素造成的影響。通過以上方法計算可得，在300 K溫度下，水在鈣基蒙脫石表面的接觸角為18.03°，這與Ethington[22]和Ballah等[23]試驗測得的數據一致性較好，也說明了鈣基蒙脫石表面的親水性。

2.3 RDF分析

徑向分布函數可以從原子尺度來表征分子團聚的結構和動態特性[24]。為了研究純水在鈣基蒙脫石表面的潤濕性，計算并繪制了Ca-Owater(水分子中的氧原子)和Ca-Hwater(水分子中的氫原子)的徑向分布函數，即以Ca為中心原子的一定范圍內出現氧原子和氫原子的概率分布函數。

如圖6所示，300 K溫度下，Ca-Owater的徑向分布函數呈現出兩個峰：在r= 2.445 ?處出現一個第一個明顯的峰，在該處Owater原子之間的距離較小，在范德華力、靜電力等作用下，水分子排列得比較緊密；隨著距離的增大，在r= 4.545 ?處出現了一個不顯著的峰，但第二個峰的峰值遠小于第一個峰，水分子的密度和液滴結構的有序度降低。Ca-Hwater的徑向分布函數也呈現了類似的雙峰特性：在r= 3.135 ?處出現了第一個十分顯著的峰，而在r= 4.995 ?處出現了一個較低的峰。這種雙峰特征表明了在Ca2+陽離子周圍形成了兩個水化層。

圖6 鈣基蒙脫石表面的徑向分布函數Fig.6 Radial distribution function of Ca-montmorillonite surface

通過對比Ca-Owater和Ca-Hwater的徑向分布函數曲線，可以看出Ca-Owater的第一峰位比Ca-Hwater的第一峰位出現得更早，這說明了鈣離子周圍出現氧原子的概率大于氫原子，水中的氧原子更易與鈣離子結合形成水合鈣離子，這是導致蒙脫石表面親水的原因之一。Ca2+與水中的氧原子之間的距離更短，表明水中的氧原子指向Ca2+，而氫原子背離Ca2+方向。此外，隨著距離的增大，導致分子間作用力的減弱，水分子密度和水分子排列的有序度減小，進而導致兩條曲線的RDF值的差值呈現出逐漸減小的趨勢。

2.4 水分子密度分析

水化層中的水分子與基底間的相互作用會影響接觸角的大小，進而改變表面的親水性，且水化層的結構是影響潤濕性的關鍵因素，因此水化層的性質在固體表面潤濕性的研究中起到了十分重要的作用。

為了研究鈣基蒙脫石表面的水化層性質，繪制了水中氧原子和氫原子沿著z軸方向的密度分布曲線，如圖7所示。在z= 8.889 ?處，氧原子密度曲線呈現一個明顯的主峰，代表著明確的第一水化層，這種顯著的水化分層現象也說明了蒙脫石表面的親水性。第一主峰的較高的峰值證實了鈣基蒙脫石表面和水分子之間存在著較強的相互作用，也表明了水分子在界面層的密集程度。與氧原子密度分布曲線的特征不同，在第一水化層內，氫原子密度曲線在z= 8.222 ?和z= 9.333 ?處存在兩個相距較近的峰。如圖7所示，在第一水化層內的氫原子和氧原子的密度遠遠高于上層水分子中的原子密度，且隨著與基底距離的增大，水化分層現象消失。這是由于水分子間距離增大，分子間相互作用力減小，液滴的結構變得松散，從而導致水分子密度降低。這一現象與距一定距離不再對水層起作用的理論相符合，同時也表明了第一水化層的水化性質的研究尤為重要。

圖7 沿z軸方向不同原子的密度分布Fig.7 Density profile of different atoms along the z-axis

3 討論

3.1 蒙脫石表面親水性的原因分析

與表面沒有陽離子吸引而呈電中性的高嶺石不同，蒙脫石晶體結構中的類質同象取代作用(如鋁氧八面體中的Al3+被Mg2+取代，硅氧四面體中的Si4+被Al3+取代)導致其片層結構帶有負電荷，因此可以在層間吸附陽離子。當水分子接觸到蒙脫石表面時，基底上的Ca2+受到沖擊從而給水分子足夠的能量去破壞與表面氧的化學鍵。Ca-Owater和Ca-Hwater的徑向分布函數也證實了鈣離子更加親氧，游離的Ca2+更容易被水分子中的氧捕獲，形成水合鈣離子。此外，表面離子產生的靜電場作用誘導了水在蒙脫石表面形成有序的吸附。這些都導致了蒙脫石表面呈現出明顯的親水性。

3.2 潤濕性對土體結構破壞的影響

蒙脫石的潤濕性和水化膨脹特性對土體結構破壞有著十分重要的影響。由于晶胞間微弱的聯結力，蒙脫石的解理較為發育，水分子極易進入，因此其具有較強的吸水膨脹性能。當水含量增加時，具有強親水性的蒙脫石會發生水化作用，促使雙電層擴展，進而削弱了顆粒間的物理化學吸附。這種作用導致蒙脫石遇水后顆粒間及晶層間會發生膨脹，破壞了土體的結構。水分子在蒙脫石表面吸附的過程中，減少的表面自由能部分轉化為力學破壞能，也會導致土體的變形和破壞。

3.3 潤濕性對水土流失的影響

中國水土流失形勢嚴峻，尤其是在黃土高原地區，水土流失造成了極大的生態和安全威脅。因此，水土流失問題亟待解決，而探尋水土流失的成因則是水土流失治理中十分關鍵的一步。

蒙脫石的潤濕性與水土流失的過程具有緊密聯系。由于黃土的礦物組成中含有大量具有強親水性的蒙脫石，其表面的水分子在較強的吸引作用下，容易在蒙脫石表面形成一層水化膜。當蒙脫石顆粒相互靠近時，產生的水化斥力較強，顆粒易分散，阻礙了蒙脫石顆粒的絮凝沉降。由于蒙脫石顆粒表面帶有較強的電負性，顆粒間存在較強的靜電斥力，從而在水中呈現為穩定的分散懸浮狀態，導致顆粒難以聚結沉降。而黃土高原地區氣候特殊，通常在長期的干旱后出現暴雨。在暴雨和地表徑流等作用下，較強的水動力易導致水流攜帶著大量的泥沙流失，降低了土壤基質的滲透能力，削弱了土體的保水能力。這不僅造成了該地區土壤侵蝕、養分流失和土地荒漠化，淤積的泥沙還會導致河道阻塞、河床抬升，從而引發旱、澇等災害。因此，研究蒙脫石表面的潤濕性對探尋水土流失的成因及防治措施具有極其重要的意義。

4 結論

通過鈣基蒙脫石表面潤濕性的分子動力學模擬，計算了水分子在其表面的接觸角，結合了不同時刻液滴的形態、徑向分布函數和元素密度分布曲線的分析，研究了鈣基蒙脫石表面的潤濕性。得到以下結論。

(1)基于分子動力學模擬，通過對液滴密度分布的計算和液-氣界面輪廓線的擬合處理，提供了一種能夠精確計算蒙脫石表面的接觸角的方法。

(2)Ca-Owater和Ca-Hwater的徑向分布函數都呈現出雙峰特性，但Ca-Owater的第一峰位比Ca-Hwater第一峰位更早出現，表明Owater更易與鈣離子結合形成水合鈣離子。

(3)表面離子產生的靜電場作用，誘導了水在蒙脫石表面形成有序吸附，進而導致水在蒙脫石表面呈現明顯的親水性。從原子尺度揭示了蒙脫石表面親水的機制，以期為土體結構破壞和水土流失等現象的研究提供理論依據。