分子動(dòng)力學(xué)模擬中的鈉基蒙脫石尺寸參數(shù)選取及其影響

閆盛熠 ，賀 鑫 ，洪望兵 ，李長(zhǎng)冬* ，孟 杰 ，付國(guó)斌

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

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

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

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

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

1 蒙脫石力學(xué)特征分子動(dòng)力學(xué)模擬

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

采用Wyoming型鈉基蒙脫石晶體[20]，所構(gòu)建的模型基本條件如下：?jiǎn)涡本?，空間群為C2/m，對(duì)稱型為L(zhǎng)2PC結(jié)構(gòu)。其晶層間依靠范德華力和靜電力聯(lián)結(jié)，在a軸和b軸方向延伸，在c軸方向堆疊。由于蒙脫石存在晶格取代效應(yīng)，所以其八面體中的Al3+常易被Mg2+、Fe2+、Fe3+等異價(jià)陽(yáng)離子置換，這也是蒙脫石晶層電荷的主要來(lái)源之一。研究中首先構(gòu)建了8個(gè)單位晶胞(4a×2b×1c)且含有2層水分子的超晶胞片層，晶格參數(shù)a=0.523 nm，b=0.906 nm，c=1.53 nm，α=γ=90°，β=99°。模型中每32個(gè)Si4+中有1個(gè)被Al3+取代，每8個(gè)Al3+中有1個(gè)被Mg2+取代，因此晶格取代效應(yīng)產(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 力場(chǎng)與水分子模型

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

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

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

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

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

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

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

(1)

(2)

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

(3)

(4)

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

(5)

(6)

根據(jù)式(5)和式(6)可計(jì)算得到多晶體的彈性模量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。因此擬研究各晶胞方向上等比例擴(kuò)胞的試驗(yàn)共7組，其擴(kuò)大倍數(shù)分別為1～7倍；同時(shí)設(shè)計(jì)5 000倍數(shù)下，比表面積在0～700范圍內(nèi)共5組模擬試驗(yàn)，最終選取尺寸為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 不同擴(kuò)展形態(tài)的蒙脫石分子模型Fig.2 Models of montmorillonite cell with different expansion shapes

2.2 模擬結(jié)果

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

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

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

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

2.3 關(guān)聯(lián)度分析

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

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

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

3 討 論

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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