閆盛熠 ，賀 鑫 ，洪望兵 ，李長冬* ，孟 杰 ，付國斌

(1.中國地質(zhì)大學(武漢)工程學院，武漢 430074；2.浙江華東建設工程有限公司，杭州 310000)

黏土礦物在中國廣泛分布于各類土壤、巖石之中，因其粒度小且具有高比表面特性而在工程地質(zhì)中發(fā)揮著重要的作用。其中蒙脫石(MMT)是地表、近地表最常見的黏土礦物之一，它是一種層狀的硅鋁酸鹽黏土礦物，其每個晶片層由兩層硅氧四面體中間夾一層鋁氧八面體組成，而因其極易產(chǎn)生晶格取代效應，導致結(jié)構(gòu)層產(chǎn)生負電荷，故層間常吸附有K+、Ca2+、Na+、Mg2+等正電荷以保持結(jié)構(gòu)為中性[1-3]。

分子動力學模擬是一門結(jié)合了多學科的研究材料化學、物理性質(zhì)的有力方法。該方法基于傳統(tǒng)牛頓力學來模擬分子體系在不同構(gòu)成狀態(tài)下的運動，從而進一步計算體系的相關參數(shù)以及宏觀性質(zhì)。與傳統(tǒng)的研究方法相比[4]，分子動力學模擬技術通常用于更好地研究原子之間的物理和化學相互作用，并獲得原子尺度上的主要性質(zhì)，為更深層次的研究提供可靠的依據(jù)[5-7]。而工程巖土體的宏觀變化往往是由其微觀結(jié)構(gòu)的改變引發(fā)的[8]，因此使用分子動力學模擬技術來對實際工程中的黏土礦物進行微觀尺度的研究就具有非常重要的研究意義。

近些年來，分子動力學模擬技術的發(fā)展為學者們從不同方面對黏土礦物進行研究提供了新的途徑[9-10]。對黏土礦物水化特征的分子動力學模擬研究多集中于不同條件下蒙脫石的結(jié)構(gòu)特征、力學特征的變化規(guī)律，結(jié)果表明含水量及溫壓對蒙脫石晶層間離子分布及力學特征有著顯著影響[11-13]。而在改性蒙脫石領域，分子動力學模擬多用來探究蒙脫石在不同替代陽離子下的物理力學性質(zhì)和形態(tài)結(jié)構(gòu)特征[14-16]。除此之外，還有部分學者基于分子動力學模擬探究了黏土礦物顆粒在不同環(huán)境條件下的水化機理，并構(gòu)建了相關模型[17-19]。已有的蒙脫石分子動力學研究所采用的模型尺寸各不相同，而模型尺寸的差異對計算產(chǎn)生的影響尚未明確，最終可能會導致蒙脫石力學性質(zhì)的模擬結(jié)果出現(xiàn)偏差。但對于蒙脫石晶胞尺寸效應的分子動力學研究還很缺乏，因此需要進一步探究尺寸效應對蒙脫石模擬計算的影響。

尺寸效應的影響因素包含有兩個方面——尺寸及形態(tài)，又由于室內(nèi)試驗手段在研究蒙脫石微觀性質(zhì)方面較為困難，因此，采用LAMMPS軟件以及Paracloud超算云平臺，并利用分子動力學模擬技術研究鈉基蒙脫石晶胞在不同尺寸和形態(tài)下的性質(zhì)變化。最終對模擬結(jié)果進行分析，利用灰色關聯(lián)度理論評價尺寸與形態(tài)對蒙脫石力學性質(zhì)的影響程度，為進一步研究蒙脫石微觀結(jié)構(gòu)對其性質(zhì)的影響提供一定的參考依據(jù)。

1 蒙脫石力學特征分子動力學模擬

1.1 蒙脫石分子動力學模型構(gòu)建

采用Wyoming型鈉基蒙脫石晶體[20]，所構(gòu)建的模型基本條件如下：單斜晶系，空間群為C2/m，對稱型為L2PC結(jié)構(gòu)。其晶層間依靠范德華力和靜電力聯(lián)結(jié)，在a軸和b軸方向延伸，在c軸方向堆疊。由于蒙脫石存在晶格取代效應，所以其八面體中的Al3+常易被Mg2+、Fe2+、Fe3+等異價陽離子置換，這也是蒙脫石晶層電荷的主要來源之一。研究中首先構(gòu)建了8個單位晶胞(4a×2b×1c)且含有2層水分子的超晶胞片層，晶格參數(shù)a=0.523 nm，b=0.906 nm，c=1.53 nm，α=γ=90°，β=99°。模型中每32個Si4+中有1個被Al3+取代，每8個Al3+中有1個被Mg2+取代，因此晶格取代效應產(chǎn)生的層間電荷為-0.75e/元胞(e為電子)，由層間Na+離子平衡。最終建立的水化蒙脫石晶格模型分子式為Na0.75(Si7.75Al0.25)(Al3.5Mg0.5)O20(OH)4·nH2O。所建模型如圖1所示。

圖1 鈉基蒙脫石分子模型Fig.1 Molecular model of Na-montmorillonite

1.2 模擬內(nèi)容

1.2.1 力場與水分子模型

力場本質(zhì)上是一種描述系統(tǒng)中各原子勢能與其位置關系的函數(shù)，它使用不同參數(shù)的組合來計算粒子系統(tǒng)的勢能。分子動力學模擬可利用不同的力場來模擬對應的微觀情況，并通過計算得到研究對象在不同條件下的宏觀物理力學性質(zhì)。因此，所選力場的精度、適用性以及測定某一方面性質(zhì)的能力就決定了分子動力學模擬結(jié)果的科學性、有效性。

選取Clayff力場[21]來模擬鈉基蒙脫石在尺寸效應影響下的力學特征。Clayff力場多用于黏土礦物的模擬計算，它認為分子總勢能是鍵結(jié)與非鍵結(jié)相互作用的結(jié)果。Clayff力場使用的水分子模型為柔性SPC模型[22]，其相關參數(shù)如表1所示。

表1 柔性SPC水分子模型參數(shù)Table 1 Parameters of flexible SPC water molecular model

1.2.2 結(jié)構(gòu)弛豫及力學計算

基于LAMMPS軟件并結(jié)合Paracloud超算云平臺(并行云計算)對水化鈉基蒙脫石進行了大規(guī)模的分子動力學模擬計算。先對蒙脫石晶胞進行擴胞處理，然后在LAMMPS軟件中對其進行弛豫處理，最后計算蒙脫石晶胞的相關力學參數(shù)。模擬中設定邊界條件為三維周期性邊界條件，時間步長設置為0.1 fs。

利用LAMMPS軟件可以計算得到所建晶胞的彈性剛度常數(shù)Cij與柔度系數(shù)Sij，然后依據(jù)Voigt-Reuss-Hill近似理論[23]可計算得到研究材料的體積模量K、剪切模量G、彈性模量E以及泊松比ν。

Voigt-Reuss-Hill近似理論中的Voigt理論假設材料各組分具有相同的應變，且與多晶體的平均應變值相等，因此可用材料的彈性剛度常數(shù)Cij來計算多晶體的體積模量KV和剪切模量GV[24]：

(1)

(2)

而Reuss理論則認為材料各組分承受了相同的應力，且這個應力與多晶體所受應力相等，因此可用材料的柔度系數(shù)Sij來計算多晶體的體積模量KR和剪切模量GR[25]：

(3)

(4)

Hill理論則根據(jù)彈性極限原理證明出了Voigt理論和Reuss理論得到的結(jié)果分別為多晶體力學參數(shù)的上限和下限，因此二者的算術平均值可近似當作多晶體的平均模量，這即是Voigt-Reuss-Hill近似[26]：

(5)

(6)

根據(jù)式(5)和式(6)可計算得到多晶體的彈性模量E和泊松比ν：

(7)

(8)

式中：KV為多晶體在Voigt理論下的體積模量；GV為多晶體在Voigt理論下的剪切模量；KR為多晶體在Reuss理論下的體積模量；GR為多晶體在Reuss理論下的剪切模量；KH為多晶體在Hill理論下的體積模量；GH為多晶體在Hill理論下的剪切模量；E為多晶體的彈性模量；ν為多晶體的泊松比。

2 模擬結(jié)果與分析

2.1 分子建模

陳寶等[27]、朱贊成等[28]、談云志等[29]研究表明蒙脫石的比表面積為550～650。因此擬研究各晶胞方向上等比例擴胞的試驗共7組，其擴大倍數(shù)分別為1～7倍；同時設計5 000倍數(shù)下，比表面積在0～700范圍內(nèi)共5組模擬試驗，最終選取尺寸為A組：20a×250b×1c；B組：100a×2b×25c；C組：100a×50b×1c；D組：500a×2b×5c；E組：500a×10b×1c。其中，a=0.523 nm，b=0.906 nm，c=1.53 nm。不同形態(tài)的各組晶胞如圖2所示。

圖2 不同擴展形態(tài)的蒙脫石分子模型Fig.2 Models of montmorillonite cell with different expansion shapes

2.2 模擬結(jié)果

根據(jù)分子動力學模擬計算得到了不同尺寸效應參數(shù)下蒙脫石晶胞的彈性剛度常數(shù)，將所得結(jié)果代入式(1)～式(8)，所得出蒙脫石晶胞力學參數(shù)變化特征如圖3所示。

由圖3(a)、圖3(b)可知，隨著擴大倍數(shù)的增大，蒙脫石晶胞的力學參數(shù)產(chǎn)生了明顯的變化。其中，體積模量、剪切模量及彈性模量皆呈先減后增的趨勢，而泊松比則與之相反，且分別在擴大倍數(shù)為2時達到最小、最大值。這一現(xiàn)象表明蒙脫石體積的增大顯著提升了蒙脫石晶胞的剛度，同時降低了其在受力之后的橫向變形量，而剪切模量、彈性模量出現(xiàn)負值則是因為晶體體積的增大及模擬過程中進行的能量最小化處理導致所建分子模型結(jié)構(gòu)的力學穩(wěn)定性降低[30]，同時減弱了蒙脫石晶體的延展性，增強了材料的脆性，使材料易發(fā)生脆性斷裂[31]。

圖3 蒙脫石力學參數(shù)模擬結(jié)果Fig.3 The simulation results of montmorillonite mechanical Parameters

由圖3(c)、圖3(d)可知，隨著不同形態(tài)的變化，除體積模量外，其余力學參數(shù)整體皆呈減小的趨勢。其中，C組的各項力學參數(shù)達到最小值。這表明形態(tài)的改變使蒙脫石晶體原子間的結(jié)合形式發(fā)生了變化。

2.3 關聯(lián)度分析

為了探討晶胞尺寸與形態(tài)這兩個影響因子對蒙脫石力學性質(zhì)的影響程度，利用灰色關聯(lián)度理論[32]來對模擬結(jié)果進行分析，分析結(jié)果見表2。

由表2可知，晶胞尺寸對蒙脫石的剪切模量等力學參數(shù)影響更大，而晶胞形態(tài)則對體積模量的影響更大。這表明蒙脫石晶胞體積的增大對各原子之間的鍵合強度影響較深，而其在三軸應力下抵抗體積變形的能力則與各軸方向上的尺寸密切相關。

表2 蒙脫石力學參數(shù)與各因素間的關聯(lián)度Table 2 Correlation degree between mechanical parameters of montmorillonite and various factors

3 討 論

研究了分子動力學模擬中不同尺寸參數(shù)的選取對鈉基蒙脫石晶體模擬結(jié)果的影響。在分子動力學模擬過程中，對蒙脫石晶體模型進行能量最小化處理，模擬過程中原子所受總作用力與總勢能如圖4所示。

由圖4(a)、圖4(b)可知，不同擴大倍數(shù)的蒙脫石晶體原子在模擬過程中的總勢能呈指數(shù)增長[33]，而其所受總作用力則具有較好的線性變化特征。當蒙脫石晶體擴大2倍時，其內(nèi)部原子所受總作用力的線性擬合結(jié)果較好。這表明由于能量最小化對不同尺寸蒙脫石晶體內(nèi)部結(jié)構(gòu)影響程度不同，而不同晶格尺寸對原子勢能、原子應力、原子位移的影響也不相同[34]，在特定溫度下，不同尺寸晶體的變形機制對礦物的強度有決定性影響[35]。因此，能量最小化調(diào)整不同尺寸晶體原子結(jié)構(gòu)，同時也改變了蒙脫石晶體原子間能量分布的方式，最終導致蒙脫石的變形破壞形式和儲能-耗能機制產(chǎn)生變化，使得自身的壓縮性和橫向變形能力以及抵抗剪切應變的能力顯著下降。

圖4 不同尺寸參數(shù)蒙脫石原子總作用力與總勢能Fig.4 The total force and potential energy of montmorillonite atoms with different size parameters

而由圖4(c)、圖4(d)可知，不同形態(tài)的蒙脫石晶體所受總作用力及其消耗能量的變化趨勢正好相反，且分別在C組時達到最大、最小值。這表明不同晶體方向上的擴大尺寸對晶體性質(zhì)的影響程度并不相同，而這種變化以及能量最小化的影響讓蒙脫石的力學性質(zhì)以及儲能-耗能機制變得不太穩(wěn)定，其中C組時的蒙脫石晶體形態(tài)穩(wěn)定性最差。因此在b軸方向上的尺寸變化對蒙脫石晶體性質(zhì)的影響最為顯著，a軸方向上的尺寸變化則影響較小，這兩者皆為晶體延伸方向上的尺寸。而在c軸這一晶體堆疊方向上的尺寸變化雖然也影響蒙脫石晶體的力學性能，但是當c軸尺寸的擴大倍數(shù)較高時，可能會使所構(gòu)建的模型完全偏離實際結(jié)構(gòu)，從而導致所建蒙脫石晶體結(jié)構(gòu)的失效。因此，不對更大尺度的c軸尺寸變化作深入研究。

研究表明模型大小及形態(tài)的不同顯著改變了分子動力學模擬的結(jié)果，使研究對象的性質(zhì)產(chǎn)生了明顯的變化，這一結(jié)果與部分學者較為接近[36-38]。但是當考慮溫度變化時，尺寸對構(gòu)建模型性能的影響就變得不明確[39]，這表明溫度的改變使得模型內(nèi)部分子結(jié)構(gòu)更加復雜。因此，還需進一步研究溫度、濕度等外界因素影響下的尺寸效應對分子動力學模擬結(jié)果的影響。

闡述了尺寸參數(shù)選取對蒙脫石分子動力學模擬結(jié)果的影響，同時分析了尺寸效應對蒙脫石力學特性的影響。研究結(jié)果表明尺寸的選取不僅影響了模擬結(jié)果的準確性，也影響了模擬過程中蒙脫石模型結(jié)構(gòu)優(yōu)化的效果，而這種影響在不同溫度、濕度等外界條件下起到的作用不同。這一結(jié)果為鈉基蒙脫石的研究提供了分子尺度的信息，讓人們可以更為準確地獲取試驗研究數(shù)據(jù)，對蒙脫石的試驗研究具有較好的指導意義，從而為一些工程應用提供理論依據(jù)。

4 結(jié)論

通過對水化鈉基蒙脫石晶胞模型進行分子動力學模擬，分析了尺寸效應對蒙脫石晶胞力學性質(zhì)的影響，得到了如下結(jié)論。

(1)當晶胞形態(tài)相同時，隨著擴胞倍數(shù)的增加，蒙脫石的體積模量、剪切模量及彈性模量皆呈先減后增的趨勢，在擴大2倍的時候達到最小值；泊松比的變化則相反，在擴大2倍的時候達到最大值。這表明蒙脫石晶胞在原始形態(tài)下，當整體尺寸擴大時，其力學性質(zhì)也得到了一定的提高；而當擴大倍數(shù)為2倍時，蒙脫石晶胞抵抗剪切破壞、彈性變形的能力最差。

(2)當擴胞倍數(shù)相同時，雖然變化趨勢不同，但是蒙脫石晶胞的體積模量、剪切模量及彈性模量 和泊松比均在C組時達到最小值。這表明蒙脫石晶胞的力學性質(zhì)隨著晶胞形態(tài)的改變變得不太穩(wěn)定，而b軸這一晶體延伸方向上的尺寸變化對其性質(zhì)的影響最為顯著，a軸次之，c軸的尺寸變化則影響最小。

(3)蒙脫石晶胞的力學性質(zhì)顯然受到了尺寸效應的影響，而尺寸效應的影響因素包括2個方面——晶胞尺寸及晶胞形態(tài)?；诨疑P聯(lián)度理論計算出蒙脫石力學參數(shù)與尺寸和形態(tài)的關聯(lián)度，最終分析結(jié)果得出晶胞尺寸對蒙脫石剪切模量等參數(shù)的影響較大，而晶胞形態(tài)對蒙脫石體積模量的影響較大。

研究為蒙脫石在分子動力學模擬中分子模型尺寸的選取提供了依據(jù)，有助于進一步深化不同分子尺寸的模擬計算。