朱 磊 沈才華 李東彪 李雪松 陳 偉 郭金勇

(1.中交南京交通工程管理有限公司, 南京 211800;2.河海大學(xué) 土木與交通學(xué)院, 南京 210098)

分子模擬方法是以計(jì)算機(jī)程序和模擬算法為工具在分子層面來(lái)研究物質(zhì)的各種理化性質(zhì),再利用統(tǒng)計(jì)學(xué)總結(jié)出模擬體系的宏觀物理量.分子動(dòng)力學(xué)方法(molecular dynamic,MD)是其中一種非常重要、應(yīng)用廣泛的分子模擬方法[1].早在1957 年,Alder和Wainwright采用硬球模型,利用分子動(dòng)力學(xué)研究氣體和液體的狀態(tài)方程,之后,學(xué)者們開(kāi)始了大量的分子動(dòng)力學(xué)研究[2].近年來(lái)分子動(dòng)力學(xué)廣泛應(yīng)用于材料科學(xué)、無(wú)線(xiàn)電電子學(xué)、有機(jī)化學(xué)等領(lǐng)域[3-5].分子動(dòng)力學(xué)模擬方法可以通過(guò)構(gòu)造比較理想的分子構(gòu)型模型,定量地再現(xiàn)真實(shí)固體中所發(fā)生的動(dòng)態(tài)過(guò)程,但鮮有人應(yīng)用分子動(dòng)力學(xué)模擬SiO2材料脆性開(kāi)裂過(guò)程.劉青松[6]基于Reax FF 反應(yīng)力場(chǎng),在微納尺度下,應(yīng)用大規(guī)模分子動(dòng)力學(xué)方法,模擬了非晶態(tài)SiO2塊體結(jié)構(gòu)和非晶態(tài)SiO2薄片孔洞結(jié)構(gòu)的單軸拉伸變形-斷裂過(guò)程;研究了微觀結(jié)構(gòu)及尺寸對(duì)非晶態(tài)SiO2力學(xué)特性的影響,揭示了結(jié)構(gòu)的微觀變形機(jī)理.韓強(qiáng)等[7]運(yùn)用分子動(dòng)力學(xué)方法模擬了含有中心裂紋的扶手椅型石墨烯的拉伸破壞行為.譚旺[8]運(yùn)用分子動(dòng)力學(xué)方法對(duì)SiC/Al復(fù)合材料的界面力學(xué)性能和界面結(jié)構(gòu)進(jìn)行模擬研究分析,模擬結(jié)果顯示塑性變形過(guò)程和斷裂機(jī)制為:首先位錯(cuò)從界面發(fā)射,并向Al內(nèi)部運(yùn)動(dòng),遇到界面發(fā)生反射,隨后發(fā)生滑移,并且有發(fā)生旋轉(zhuǎn)和晶體重構(gòu)的現(xiàn)象,滑移后有空洞生成,最后空洞增大形成裂紋而最終斷裂.楊利等[9]用分子動(dòng)力學(xué)方法,研究了預(yù)制邊界裂紋并加入Nb的單晶γ-TiAL合金的裂紋擴(kuò)展過(guò)程,分析了裂紋在1K 溫度下的擴(kuò)展行為.馬磊等[10]采用分子動(dòng)力學(xué)方法,結(jié)合Tersoff勢(shì)函數(shù),模擬了α-SiO2晶體的力學(xué)性能,并研究了溫度對(duì)α-SiO2力學(xué)性能的影響.可見(jiàn),隨著分子模擬方法理論的完善,越來(lái)越多的學(xué)者開(kāi)始采用分子動(dòng)力學(xué)對(duì)不同材料的力學(xué)行為進(jìn)行模擬探究,特別為揭示宏觀現(xiàn)象下的微觀機(jī)理提供了一種很好的途徑.考慮不同工程荷載特點(diǎn)下裂紋擴(kuò)展過(guò)程的研究較少,因此本文采用分子動(dòng)力學(xué)模擬技術(shù)探究巖石材料中常見(jiàn)的SiO2材料脆性開(kāi)裂的過(guò)程,建立斷裂因子、微觀參數(shù)、能量等裂紋擴(kuò)展過(guò)程中微宏觀的關(guān)系,為揭示固體材料裂紋擴(kuò)展的微觀機(jī)理研究提供參考,這對(duì)促進(jìn)建筑材料的脆性開(kāi)裂過(guò)程微觀機(jī)理研究具有重要意義.

1 基于分子動(dòng)力學(xué)的微觀裂紋擴(kuò)展模擬方法研究

分子動(dòng)力學(xué)方法可以利用計(jì)算機(jī)模擬進(jìn)行較大體系的長(zhǎng)時(shí)間尺度的動(dòng)態(tài)問(wèn)題的研究,該方法特別適用于微納米裂紋及擴(kuò)展的研究.利用該方法不但可以觀察到裂紋擴(kuò)展中原子鍵斷裂的細(xì)節(jié)、裂紋擴(kuò)展的形式,還可計(jì)算出原子的位置和速度進(jìn)而確定模擬中原子的運(yùn)動(dòng)軌跡、裂紋的擴(kuò)展速度、體系的總勢(shì)能和動(dòng)能等物理量[11].

LAMMPS(Large-scale Atomic/Molecular Massively Parallel Simulator)是一款開(kāi)源軟件,由美國(guó)能源部的桑迪亞國(guó)家實(shí)驗(yàn)室(Sandia National Laboratories)研發(fā),該軟件免費(fèi)且功能強(qiáng)大,不足之處是沒(méi)有像其他商業(yè)軟件的前處理與后處理用戶(hù)界面,使用起來(lái)可能帶來(lái)不便,因此LAMMPS的建模方法有:使用LAMMPS 的代碼建模;借助其他軟件進(jìn)行建模.

LAMMPS對(duì)所含成分比較少的物質(zhì)的納觀行為有較好的模擬效果,但是巖石是復(fù)雜混合體多孔結(jié)構(gòu),目前沒(méi)有與之一一對(duì)應(yīng)的勢(shì)函數(shù)可以使用,所以本文采用巖石材料中最常見(jiàn)且分子構(gòu)型較簡(jiǎn)單的二氧化硅(SiO2)來(lái)分析巖石基材的裂紋擴(kuò)展過(guò)程,探索建立工程材料的納觀機(jī)理分析方法.

分子動(dòng)力學(xué)模擬方法的原理是通過(guò)原子間的相互作用勢(shì),按照經(jīng)典牛頓運(yùn)動(dòng)定律求出原子軌跡及其演化過(guò)程.因此,建立合適的勢(shì)能函數(shù)是進(jìn)行分子動(dòng)力學(xué)模擬的第一步,也是最關(guān)鍵的一步.勢(shì)能函數(shù)的正確與否,直接關(guān)系到模擬結(jié)果的精確性和可靠性[12].

本文SiO2的模擬采用的是vashishta勢(shì)函數(shù)[13],表達(dá)式為:

表1 Si-O-O 參數(shù)表[12]

二氧化硅(SiO2)的模擬采用的建模方法是第二種方法,使用VMD 軟件建立SiO2模型,該軟件為一種分子可視化程序,是Visual Molecular Dynamics的縮寫(xiě),由University of Illinois at Champaign的Theoretical Biophysics Group研發(fā)(具體分子構(gòu)型如圖1所示),建模過(guò)程如圖2所示,尺寸為:x·y·z=199.16×99.56×6.948?.

圖1 二氧化硅空間構(gòu)型

圖2 計(jì)算模型示意圖

2 NVE系綜裂紋擴(kuò)展行為分析

初始模型溫度為10 K,初始裂紋長(zhǎng)度10?,采用消除A 和B交界面的作用模擬初始裂縫,即不考慮初始裂縫寬度和形狀.位移加載速度為10?/ps,模擬分析裂紋擴(kuò)展過(guò)程.模擬結(jié)果顯示:在拉伸過(guò)程中,能量逐漸向裂紋尖端聚集,裂紋尖端出現(xiàn)應(yīng)力集中,分子之間吸引力較小的部位慢慢被拉開(kāi),直至勢(shì)函數(shù)作用距離之外,尖端開(kāi)裂.如圖3(a)所示,在加載10 ps左右,裂紋尖端出現(xiàn)空位,顯然,該空位沒(méi)有出現(xiàn)在裂紋正前方,是二氧化硅的構(gòu)型決定的,即原子的空間排布導(dǎo)致空位處的原子之間勢(shì)能較低.加載到20 ps,如圖3(b)所示,預(yù)制裂紋尖端及空位各自形成一個(gè)裂紋尖端,即“裂紋1”和“裂紋2”(空位處),隨后,在兩個(gè)裂紋尖端的“競(jìng)爭(zhēng)”中,裂尖之間的原子被拉開(kāi)直至拉斷,“裂紋1”和“裂紋2”合并后繼續(xù)向前擴(kuò)展,如圖3(c)、(d)所示.圖3(d)、(e)、(f)分別表示裂紋擴(kuò)展至圖2(c)所示的“圓1”、“圓2”和“圓3”的位置.

圖3 裂紋擴(kuò)展過(guò)程

在裂紋擴(kuò)展路徑,即拉伸方向上取11個(gè)貼近斷裂面的點(diǎn)(圖4(a)),間距為20?×10=200?.觀察這些點(diǎn)y方向位移變化(圖4(b)),判斷各點(diǎn)起裂時(shí)間,并計(jì)算裂紋擴(kuò)展的大致速度以及裂紋長(zhǎng)度變化(圖4(c)、(d)).由圖4(c)可看出,裂紋在擴(kuò)展過(guò)程中并不是勻速前進(jìn)的,而是上下震蕩.

圖4 裂紋擴(kuò)展速率與長(zhǎng)度

結(jié)合圖4(d)可看出,裂紋長(zhǎng)度曲線(xiàn)斜率在260 ps之前逐漸變小,之后曲率增大.所以速率變化的總體趨勢(shì)是先減小(260 ps之前)后增大.不計(jì)裂紋停止的時(shí)間,裂紋擴(kuò)展速率在2?/ps(200 m/s)左右,最高速率超過(guò)4?/ps(400 m/s),該計(jì)算結(jié)果與S.Chandra等人[14]計(jì)算的Al在1K 溫度下擴(kuò)展速率320 m/s比較接近,所以本模型計(jì)算結(jié)果基本合理可靠.

在加載過(guò)程中,圖2(c)的A 區(qū)域的平均應(yīng)力如圖5(a)所示,整個(gè)過(guò)程應(yīng)力變化為先減小后略微增大,從260 ps左右開(kāi)始再次減小.應(yīng)力的計(jì)算采用的是維里應(yīng)力[14],維里應(yīng)力表征原子間相互作用力,可正可負(fù),表達(dá)式為:

式中:V是需要計(jì)算應(yīng)力部分的體積;i,j表示原子i周?chē)膉原子;r iα,r jα表示i原子和j原子在α方向的位置;f ijβ表示β方向j原子對(duì)i原子的力;m i表示i原子的質(zhì)量;v iα和v iβ分別表示i原子在α和β方向的速度.

初始應(yīng)力較高的原因是,在裂紋擴(kuò)展前,A 區(qū)域的右下角處于預(yù)制裂紋尖端的附近,該處應(yīng)力集中,從而導(dǎo)致A 區(qū)域總體平均應(yīng)力較大.裂紋擴(kuò)展后,裂紋尖端逐漸遠(yuǎn)離A 區(qū)域,因此該區(qū)域應(yīng)力被釋放后下降,且應(yīng)力分布逐漸均勻.應(yīng)力峰值時(shí)刻對(duì)應(yīng)裂紋加速擴(kuò)展時(shí)刻.圖5(b)裂紋應(yīng)力強(qiáng)度因子按照斷裂力學(xué)應(yīng)力強(qiáng)度因子公式來(lái)計(jì)算:

式中;KⅠ是裂紋應(yīng)力強(qiáng)度因子;Y是與裂紋形狀和加載方式有關(guān)的系數(shù)[15],這里取1.12;σ是A 區(qū)域平均應(yīng)力,取圖5(a)中的數(shù)據(jù);c是裂紋總長(zhǎng)度,取4(d)中的數(shù)據(jù).應(yīng)力強(qiáng)度因子計(jì)算結(jié)果如圖5(b)所示,應(yīng)力強(qiáng)度因子逐漸升高,20 ps左右裂紋起裂,此時(shí)裂紋應(yīng)力強(qiáng)度因子可能達(dá)到了斷裂韌性(根據(jù)圖5(a)計(jì)算的A 區(qū)域平均最大應(yīng)力為380 MPa,初始裂紋長(zhǎng)度為10,代入公式(5)可獲得其斷裂韌度KⅠC=0.02 MPa·m1/2),導(dǎo)致裂紋開(kāi)裂.裂紋驅(qū)動(dòng)力的表達(dá)式為:

其中:GⅠ是裂紋驅(qū)動(dòng)力;B是模型寬度;-?U/?c是指系綜總勢(shì)能隨著裂紋擴(kuò)展降低的量;E是彈性模量.20 ps左右裂紋起裂.據(jù)資料可得[16]:二氧化硅彈性模量約為70 GPa,計(jì)算得到起裂裂紋驅(qū)動(dòng)力為0.005 7 J/m2.

如圖6所示,裂紋失穩(wěn)擴(kuò)展對(duì)應(yīng)圖中的切點(diǎn),c＞c0時(shí)有?C/?c＞?R/?c,即裂紋擴(kuò)展動(dòng)力增長(zhǎng)率高于阻力增長(zhǎng)率,裂紋失穩(wěn)擴(kuò)展.本文理論曲線(xiàn)取自《應(yīng)力強(qiáng)度因子手冊(cè)》[17].由圖5(c)可看出,理論曲線(xiàn)與計(jì)算結(jié)果在虛線(xiàn)之前吻合較好,虛線(xiàn)之后開(kāi)始偏離,可能是加載方式以及模型尺寸效應(yīng)導(dǎo)致.當(dāng)裂紋驅(qū)動(dòng)力GⅠ達(dá)到最大值,對(duì)應(yīng)圖5(b)中應(yīng)力強(qiáng)度因子的KⅠmax處,從該處開(kāi)始,應(yīng)力強(qiáng)度因子下降,亦即裂紋驅(qū)動(dòng)力下降,有?G/?c＜?R/?c,裂紋停止擴(kuò)展,此時(shí)對(duì)應(yīng)圖6中G曲線(xiàn)的頂點(diǎn).在裂紋驅(qū)動(dòng)力(阻力)曲線(xiàn)圖中,可大致分為Ⅰ、Ⅱ、Ⅲ3個(gè)階段,Ⅰ階段?G/?c＜?R/?c,裂紋阻力大于驅(qū)動(dòng)力,裂紋未擴(kuò)展,Ⅱ階段?G/?c＞?R/?c,裂紋驅(qū)動(dòng)力大于阻力,裂紋擴(kuò)展,Ⅲ階段?G/?c＜?R/?c,裂紋擴(kuò)展動(dòng)力不足,停止擴(kuò)展.

圖5 A 區(qū)域應(yīng)力與裂紋應(yīng)力強(qiáng)度因子

圖6 裂紋擴(kuò)展動(dòng)力-阻力曲線(xiàn)

在裂紋擴(kuò)展過(guò)程中,圖2(c)中3個(gè)圓形裂紋尖端區(qū)域能量變化如圖7(a)所示,3個(gè)區(qū)域變化趨勢(shì)大致相同,在裂紋擴(kuò)展的瞬間能量突然升高,隨后維持穩(wěn)定,而且NVE 系綜無(wú)其他能量耗散項(xiàng),因此可認(rèn)為能量升高是因?yàn)?個(gè)區(qū)域裂紋斷裂后的表面儲(chǔ)存有表面能.從圖7(b)中可看出,C 區(qū)域系統(tǒng)動(dòng)能變化不明顯,可以認(rèn)為3個(gè)圓形裂紋尖端區(qū)域能量升高的部分均是表面能,每個(gè)圓形區(qū)域裂紋開(kāi)裂截面面積為41.688?2,開(kāi)裂后能量增長(zhǎng)約為0.6 e V,所以表面能為0.23 J/m2.計(jì)算過(guò)程是在極低溫度(10 K)下進(jìn)行,而且斷裂過(guò)程裂紋尖端始終是尖銳的,所以模型高度脆性.Lawn[18]給出的高度脆性材料的表面能在10-1～10之間,本模型計(jì)算的表面能結(jié)果在該范圍內(nèi),說(shuō)明本方法模擬分析結(jié)果基本合理.

由圖7(a)可計(jì)算出裂紋從圓1擴(kuò)展到圓2的速率v12=0.49?/ps,以及從圓2擴(kuò)展到圓3的平均速率v23=0.46?/ps.

圖7 能量變化

由裂紋尖端應(yīng)力場(chǎng)的彈性解表達(dá)式可知,距離裂紋尖端越近,應(yīng)力越大.理論上,裂紋尖端應(yīng)力趨于無(wú)窮大,這是不符合物理規(guī)律的,所以裂紋尖端應(yīng)力到底多大時(shí)裂紋開(kāi)裂是彈性力學(xué)解無(wú)法觸及的“黑匣子”.物理層面,裂紋尖端是由原子分布形成的,圖3已經(jīng)說(shuō)明了裂紋擴(kuò)展是原子之間相互作用斷開(kāi)即鍵斷裂導(dǎo)致的,所以裂紋尖端應(yīng)力達(dá)到一定程度,原子鍵被拉斷,裂紋擴(kuò)展.為了探究裂紋尖端應(yīng)力變化,得到裂紋擴(kuò)展過(guò)程,計(jì)算圖2(c)中圓1、2、3三個(gè)區(qū)域水平和豎直方向應(yīng)力變化,如圖8所示.在微觀尺度上,原子時(shí)刻都在運(yùn)動(dòng),因此應(yīng)力值是波動(dòng)的.圖8(a)是3個(gè)裂紋尖端區(qū)域水平方向應(yīng)力的變化,裂紋依次進(jìn)入3個(gè)區(qū)域,3個(gè)區(qū)域先后出現(xiàn)峰值,峰值對(duì)應(yīng)裂紋進(jìn)入該區(qū)域的時(shí)間,可見(jiàn)裂紋擴(kuò)展過(guò)程尖端水平向應(yīng)力最大可達(dá)到4 GPa,同理,8(b)表示裂紋擴(kuò)展過(guò)程中3個(gè)尖端區(qū)域豎直向應(yīng)力變化,裂紋進(jìn)入的瞬間應(yīng)力可高達(dá)6 GPa左右.對(duì)比8(a)和8(b)可發(fā)現(xiàn),水平向應(yīng)力是逐漸增加、逐漸降低,垂直向則是快速增加、裂紋擴(kuò)展后幾乎是瞬間下降至裂紋未靠近的應(yīng)力值范圍,說(shuō)明裂紋尖端在向前推移的過(guò)程中,水平向應(yīng)力受影響的原子遠(yuǎn)多于垂直向.裂紋擴(kuò)展時(shí),裂紋尖端應(yīng)力的“奇異性”是不存在的,應(yīng)力增加到6 GPa左右裂紋會(huì)裂開(kāi)并向前擴(kuò)展.

圖8 裂紋擴(kuò)展過(guò)程裂尖應(yīng)力變化

3 不同系綜對(duì)裂紋擴(kuò)展行為的影響規(guī)律

分子動(dòng)力學(xué)模擬,總是在一定的系綜下進(jìn)行的,不同系綜反映了不同的工程環(huán)境,實(shí)際工程環(huán)境非常復(fù)雜,進(jìn)行機(jī)理研究時(shí)一般假設(shè)為不同的理想環(huán)境,即用不同的系綜來(lái)描述.系統(tǒng)原子數(shù)N,體積V,能量E保持不變的環(huán)境稱(chēng)為NVE 系綜.系統(tǒng)原子數(shù)N,體積V,溫度T保持不變,且總動(dòng)量保持不變,則稱(chēng)為NVT 系綜.在NVT 系綜下進(jìn)行計(jì)算,計(jì)算結(jié)果如圖9所示.

圖9 NVT 系綜下的計(jì)算結(jié)果

從如9(a)可看出,C區(qū)域NVT 總能量增長(zhǎng)幅度小于NVE 系綜,分析認(rèn)為NVT 系綜可能與外界環(huán)境有能量交換.從圖9(b)可看出,兩個(gè)系綜A 區(qū)域Y方向應(yīng)力開(kāi)始時(shí)刻相似,隨后,NVT 系綜下應(yīng)力增大,NVE系綜基本不變,最后二者又都經(jīng)歷減小的過(guò)程.由圖9(c)、(d)可看出,NVE系綜下裂紋擴(kuò)展平均速率高于NVT 系綜,從圖9(e)可看出,NVT 系綜下,裂紋擴(kuò)展的應(yīng)力強(qiáng)度因子大于NVE 系綜下的計(jì)算結(jié)果,這可能是由于NVT 系綜允許與外界能量交換,能量損失較大,導(dǎo)致裂紋擴(kuò)展需要更大的裂紋應(yīng)力強(qiáng)度因子.

4 結(jié) 論

宏觀計(jì)算在一定程度上規(guī)避了裂紋尖端的力學(xué)行為,而裂紋尖端的力學(xué)行為恰恰是裂紋擴(kuò)展研究的關(guān)鍵.裂紋的起裂擴(kuò)展歸根結(jié)底是材料微觀尺度上原子之間相互作用力從有到無(wú)的過(guò)程,即鍵斷裂的過(guò)程.采用vashishta勢(shì)函數(shù),使用VMD 軟件的建立SiO2模型,使用開(kāi)源分子動(dòng)力學(xué)軟件LAMMPS 對(duì)巖石常見(jiàn)主要成分SiO2進(jìn)行微觀尺度裂紋擴(kuò)展行為研究,不僅計(jì)算結(jié)果和呈現(xiàn)規(guī)律基本合理,而且可以很好地揭示微觀系統(tǒng)對(duì)裂紋擴(kuò)展過(guò)程的影響,為探究工程材料的破壞微觀機(jī)理提供了一種新途徑.本文在建立了分子動(dòng)力學(xué)模擬技術(shù)模擬SiO2材料裂紋過(guò)程微觀機(jī)理研究方法的基礎(chǔ)上,針對(duì)NVE 和NVT 系綜下對(duì)裂紋擴(kuò)展的行為規(guī)律也做了對(duì)比,其主要結(jié)論如下:

1)在拉伸過(guò)程中,能量逐漸向裂紋尖端聚集,裂紋尖端出現(xiàn)應(yīng)力集中,勢(shì)能較小的原子之間慢慢被拉開(kāi),裂紋擴(kuò)展方向受分子構(gòu)型影響,并非直線(xiàn),尖端拉裂的過(guò)程非常短暫,計(jì)算的拉裂強(qiáng)度最高達(dá)6 GPa,而且裂紋擴(kuò)展過(guò)程中計(jì)算應(yīng)力值呈現(xiàn)波動(dòng)狀態(tài).裂紋擴(kuò)展可大致分為Ⅰ、Ⅱ、Ⅲ3個(gè)階段:Ⅰ階段?G/?c＜?R/?c,裂紋阻力大于驅(qū)動(dòng)力,裂紋未擴(kuò)展;Ⅱ階段?G/?c＞?R/?c,裂紋驅(qū)動(dòng)力大于阻力,裂紋擴(kuò)展;Ⅲ階段?G/?c＜?R/?c,裂紋擴(kuò)展動(dòng)力不足,停止擴(kuò)展.

2)不同系綜下裂紋擴(kuò)展速率是不同的,NVE 系綜下裂紋擴(kuò)展平均速率高于NVT 系綜;NVT 系綜允許與外界能量交換,裂紋擴(kuò)展的計(jì)算應(yīng)力強(qiáng)度因子大于NVE系綜的計(jì)算值.

本文方法能有效反映不同環(huán)境下裂紋微觀擴(kuò)展機(jī)理的區(qū)別,并揭示其斷裂因子和能量的變化規(guī)律,分析結(jié)論基本符合實(shí)際情況,為探究其他常見(jiàn)巖石類(lèi)材料或工程材料的微觀機(jī)理研究提供了新的途徑.