呂世明,沈才華,郭 丹,陳 偉,李雪松,唐 凱

(1.保利長大海外工程有限公司,廣東 廣州 510623；2.河海大學(xué) 巖土力學(xué)與堤壩工程教育部重點實驗室,江蘇 南京 210098)

1 前 言

分子動力學(xué)(molecular dynamic,MD)是一種非常重要、應(yīng)用廣泛的分子模擬方法[1]。該方法通過分子間作用力改變分子的坐標(biāo)和動量,得到整個體系的動力學(xué)信息[2]。早在1957年,Alder和Wainwright[3]采用分子動力學(xué)方法來研究氣體和液體的狀態(tài)方程,此后,學(xué)者們開展了大量的分子動力學(xué)研究工作。1995年,Plimpton[4]證明了分子動力學(xué)并行運算的高效性,開發(fā)了功能強大的分子動力學(xué)軟件LAMMPS(large-scale atomic molecular massively parallel simulator),該軟件目前廣泛用于各領(lǐng)域。Zhou等[5]采用分子動力學(xué)模擬方法對裂紋擴展過程進行了研究,結(jié)果表明裂紋在擴展初期擴展速度快速增大,最后達到穩(wěn)定。Holland等[6]利用分子動力學(xué)軟件對Si的裂紋尖端擴展速度和能通量進行了比較。Huang等[7]用Pearson,MEAM 和Tersoff勢函數(shù)分別對SiC中的缺陷點或簇的形成能以及遷移能進行研究。國內(nèi)采用該研究方法逐漸增多,但對巖石的研究較少。何滿潮等[8]采用分子動力學(xué)法研究高嶺石在不同溫度下拉伸力學(xué)性能的溫度效應(yīng)。孫曉慶[9]對在不同的系統(tǒng)溫度、加載速率的條件下Si中裂紋擴展進行分子動力學(xué)模擬計算。孟筠青等[10]運用分子動力學(xué)模擬法,構(gòu)建了包含7 000 個趙莊煤大分子的基體仿真模型,模擬結(jié)果與試驗結(jié)果較為吻合。分子動力學(xué)中勢函數(shù)是核心,因此陳強等[11]對常用勢函數(shù)進行了歸類總結(jié)評價。LJ勢函數(shù)應(yīng)用最廣[12],但主要用于惰性氣體或液相轉(zhuǎn)變的分子模擬研究[13]。郭宗標(biāo)等[14]基于Tersoff勢函數(shù)研究了SiC 沿著固定三個方向的初始裂紋的擴展力學(xué)行為。巖石作為一種天然材料,內(nèi)部存在大量微裂紋損傷,采用分子動力學(xué)方法研究原生裂紋的擴展以及在不同狀態(tài)下的變化過程,對于研究巖體的失穩(wěn)破壞具有重要意義。

2 基于微觀裂紋模型的建立

LAMMPS是一款開源軟件,免費而且功能強大,可以很好的模擬材料分子尺度的力學(xué)行為,但目前針對巖石材料的研究很少,主要原因是目前沒有與之一一對應(yīng)的、應(yīng)用成熟的勢函數(shù)可以使用,所以本研究針對巖石的主要組成成分SiO2,建立了分析巖石基材的裂紋擴展微觀機理研究的方法。

2.1 勢函數(shù)與模型

勢函數(shù)根據(jù)多體作用的復(fù)雜程度可分為對勢和多體勢,SiO2其分子結(jié)構(gòu)的研究相對較多,普遍認為采用多體勢模擬較合適,常用的多體勢有Feuston-Garofalini(FG)勢、Vashishta、Vessal 等。本 實 驗SiO2的模擬采用由兩體勢和三體勢組成的Vashishta勢函數(shù)[15],表達式為:

Vashishta勢函數(shù)的參數(shù)表見表1。

表1 Si-O-O 參數(shù)表[15]Table 1 Si-O-O parameter table[15]

使用VMD 軟件的建立SiO2模型(其中SiO2分子構(gòu)型見圖1),建模方法見圖2,試樣尺寸為x×y×z=199.16?×99.56?×6.948?。

圖1 二氧化硅空間構(gòu)型 (a)分子構(gòu)型；(b)空間構(gòu)型Fig.1 Spatial configuration of silica(a)molecular configuration；(b)spatial configuration

圖2 計算模型示意圖 (a)SiO2 模型示意圖；(b)模型實例Fig.2 Schematic diagram of calculation model (a)schematic diagram of silica model；(b)model examples

3 不同影響因素作用下對裂紋擴展行為的影響

3.1 加載速率對裂紋擴展行為的影響

微正則系綜(NVE系綜),在10 K 溫度下,研究了加載速率分別為10?/ps和1?/ps時,對Ⅰ型裂紋擴展行為的影響。圖3(a)顯示:當(dāng)加載速率為10?/ps時,作用一段時間后,空位形成。圖3(b)顯示:繼續(xù)作用一段時間后預(yù)制裂紋尖端和空位兩個不同位置分別形成一個新的裂紋尖端,稱為“裂紋1(即預(yù)制裂紋擴展的裂紋尖端crack tip)”和“裂紋2(即空位vacancy)”,隨著載荷能量的進一步傳播,在兩個裂紋尖端形成能量重分布的“競爭”態(tài)勢,最終裂紋1先開始起裂,裂尖1被拉開,從而產(chǎn)生能量釋放并影響裂紋2周圍能量場,這可能會抑制裂紋2的擴展使得裂紋最終沿著裂尖1向前擴展。

圖3 加載速率為10?/ps裂紋擴展行為 (a)空位形成；(b)裂紋沿裂紋尖端擴展Fig.3 Crack growth behavior at a loading rate of 10?/ps (a)vacancy formation；(b)crack propagates along the crack tip

當(dāng)加載速率為1?/ps時,從圖4(a)可以看出,初始裂紋尖端出現(xiàn)空位,但在低加載速率下裂紋并未沿著空位擴展,而是從裂紋尖端向前擴展,原因是雖然空位處原子間結(jié)合能較低,但慢速加載導(dǎo)致原子有充足的時間重新調(diào)整位置,降低勢能,提高整體穩(wěn)定性,空位得到填充,導(dǎo)致裂紋沿尖端擴展。

圖4 加載速率為1?/ps裂紋擴展行為 (a)空位形成；(b)裂紋沿裂紋尖端擴展Fig.4 Crack growth behavior at loading rate of 1?/ps (a)vacancy formation；(b)crack propagates along the crack tip

如圖5所示,圖2(a)中圓2區(qū)域在裂紋擴展后能量變化稍低于圓1、3區(qū)域,裂紋穿過圓2區(qū)域用時也較長,該區(qū)域原子有更多的時間調(diào)整位置,降低勢能,提高穩(wěn)定性,所以圓2區(qū)域最后總能量偏低?？傮w勢能、動能,A 區(qū)域平均應(yīng)力以及總體溫度的變化與加載速率為10?/ps時的計算結(jié)果規(guī)律基本相同,但模擬發(fā)現(xiàn)加載速率會引發(fā)裂紋擴展路徑的變化,這與試驗結(jié)果基本一致。

圖5 加載速率對裂紋擴展的影響 (a)不同加載速率圓1能量變化；(b)不同加載速率圓2能量變化；(c)不同加載速率圓3能量變化；(d)不同加載速率系統(tǒng)總能量變化；(e)A 區(qū)域平均應(yīng)力變化；(f)系統(tǒng)溫度變化Fig.5 Effect of loading rate on crack growth (a)energy change of circle 1 under different loading rates；(b)energy change of circle 1 under different loading rates；(c)energy change of circle 3 under different loading rates,(d)total energy change of the system with different loading rates；(e)change of average stress in area A；(f)system temperature change

3.2 力加載對裂紋擴展行為的影響

NVE系綜,10 K 溫度下,將圖2(a)中速度加載換成力加載,并將圓2、3 位置進行調(diào)整,調(diào)整后如圖6(a)所示。應(yīng)力采用維里應(yīng)力計算方法,如式(4)所示:

式中:V為需要計算應(yīng)力部分的體積；i,j為原子i周圍的j原 子；r iα,r jα為i原 子 和j原 子 在α方 向 的 位置；f ijβ為β方向j原子對i原子的力；m i為i原子的質(zhì)量；v iα和v iβ為i原子在α和β方向的速度。

調(diào)整加載力的大小,當(dāng)加載應(yīng)力＜1.793 GPa時,裂紋不擴展；當(dāng)加載應(yīng)力≥1.793 GPa時,裂紋擴展。初始裂紋長度為9.5?,根據(jù)式(5)和式(6)可計算出起裂應(yīng)力強度因子與裂紋擴展驅(qū)動力:

力加載計算的裂紋擴展應(yīng)力強度因子比速度加載情況下大4.9倍,裂紋擴展驅(qū)動力增大了24倍。裂紋一旦開裂,迅速擴展,由于裂紋擴展速度過快,而且是端部加載,導(dǎo)致模型縱向斷裂,如圖6(d)所示,究其原因,是圓2、3部分的水平應(yīng)力大于豎直應(yīng)力(圖6(f),(g))；斷裂后,斷裂的部分受力快速運動,因此系統(tǒng)溫度急劇升高,系統(tǒng)動能也急劇升高,如圖6(h)、(i)所示。

圖6 力加載對裂紋擴展行為的影響 (a)力加載模型示意圖；(b)裂紋沿裂紋尖端擴展(5 ps)；(c)縱向起裂(25 ps)；(d)縱向裂紋擴展(26 ps)；(e)圓1水平應(yīng)力與豎直應(yīng)力；(f)圓2水平應(yīng)力與豎直應(yīng)力；(g)圓3水平應(yīng)力與豎直應(yīng)力；(h)溫度變化；(i)動能、勢能、總能量變化；(j)三個圓區(qū)域能量變化Fig.6 Effect of force loading on crack growth behavior (a)schematic diagram of force loading model；(b)crack propagation along crack tip(5 ps)；(c)longitudinal crack initiation(25 ps)；(d)longitudinal crack growth(26 ps)；(e)horizontal stress and vertical stress of circle 1；(f)horizontal stress and vertical stress of circle 2；(g)horizontal stress and vertical stress of circle 3；(h)temperature change；(i)change of kinetic energy,potential energy and total energy；(j)energy change of three circular regions

力加載方式顯示,當(dāng)力小于某閾值時,微裂紋不會擴展,所以理論上脆裂蠕變變形不會發(fā)生,這與巖石力學(xué)的長期強度理論基本一致,但數(shù)值模擬的計算裂紋起裂應(yīng)力遠大于長期強度值,主要是因為數(shù)值模擬的是微觀的、完整的分子結(jié)構(gòu)應(yīng)力狀態(tài)。本研究模擬的微觀裂紋未考慮實際巖石微環(huán)境存在缺陷的影響,而現(xiàn)實中,巖石暴露在空氣中,巖石裂紋尖端含有水或者電離的羥基OHˉ,SiO2的Si—O鍵與之發(fā)生反應(yīng),見式(7):

圖7 應(yīng)力腐蝕下的G-R 曲線Fig.7 G-R curves under stress corrosion

3.3 不同溫度下力加載對裂紋擴展行為的影響

NVE系綜下,對平均加載應(yīng)力1.793 GPa的系統(tǒng)溫度設(shè)置分別為1、2、3 K,研究力加載情況下,溫度對裂紋擴展行為的影響。如圖8(a)～(c)所示,對比這三個裂紋擴展的行為,發(fā)現(xiàn)1、3 K 溫度下,沿著+y方向分叉的裂紋位置在沿著ˉy方向分叉裂紋的左邊,2 K溫度下,二者的位置相反。三個裂紋擴展均表現(xiàn)出初始慢,之后越來越快的特點,水平向裂紋最終分叉導(dǎo)致裂紋不再沿著+x方向擴展。不同溫度下,水平裂紋的最終擴展長度和平均擴展速率如圖8(d)、(e)所示。從圖可見,不論裂紋最終擴展長度還是平均擴展速率,從1 K 到3 K 均是增加。

圖8 不同溫度下力加載裂紋擴展行為 (a)1 K 溫度下裂紋擴展行為；(b)2 K 溫度下裂紋擴展行為；(c)3K 溫度下裂紋擴展行為；(d)溫度與裂紋總長的關(guān)系；(e)溫度與裂紋平均擴展速率的關(guān)系Fig.8 Crack growth behavior under force loading at different temperatures (a)crack growth behavior at 1 K；(b)crack growth behavior at 2 K；(c)crack growth behavior at 3 K；(d)relationship between temperature and crack length；(e)relationship between temperature and average crack growth rate

由此可見,力加載作用下,裂紋擴展行為不僅與力的大小有關(guān),還與溫度有關(guān),溫度影響原子運動,從而影響勢函數(shù),原子之間作用力改變,裂紋擴展行為隨之改變。

3.4 模型厚度對裂紋擴展行為的影響

NVE 系綜下,溫度為10 K,初始裂紋長度為10?,加載速度為10?/ps,模型厚度由6.948? 增大到34.74?,對模型厚度變化影響進行研究,結(jié)果發(fā)現(xiàn),模型厚度會改變裂紋擴展行為。從圖9(a)、(b)可以看出在擴展的最后階段,裂紋尖端偏離了原來的方向,有一層原子“粘”在下表面,裂紋沿著新的表面擴展,分析認為產(chǎn)生這一現(xiàn)象的原因主要是裂紋尖端始終自洽尋找斷裂韌度或者擴展阻力小的地方擴展,厚度方向的增加使裂紋擴展過程中裂紋尖端的“弱鍵”增加,使得擴展行為發(fā)生了細微的變化。從圖9(d)可以看出,厚度增加后,裂紋擴展過程A 區(qū)域(見圖6(a))平均拉應(yīng)力增加。厚度增加后,系統(tǒng)內(nèi)原子數(shù)量增加,所以系統(tǒng)動能勢能和總能量變化均增大,裂紋擴展形成的表面增大(圖9(c)),擴展后三個圓區(qū)域內(nèi)能量增幅變大,如圖9(e～g)所示。假設(shè)裂紋形狀系數(shù)不變時,模型厚度增大的同時也導(dǎo)致開裂過程中所需應(yīng)力強度因子計算值增大(圖9(h)),分析認為這是由裂紋厚度與分子結(jié)構(gòu)的相對尺度大小決定的,即厚度與分子結(jié)構(gòu)尺度相當(dāng)時,存在宏觀力學(xué)上的尺度效應(yīng),此時的裂紋形狀系數(shù)理論上應(yīng)該不同。

圖9 不同厚度裂紋擴展行為 (a)模型正視圖；(b)模型左視圖；(c)系統(tǒng)能量變化；(d)A 區(qū)域應(yīng)力變化；(e)不同厚度圓1能量變化；(f)不同厚度圓2能量變化；(g)不同厚度圓3能量變化；(h)不同厚度應(yīng)力強度因子計算值Fig.9 Crack growth behavior of different thickness (a)model front view；(b)left view of the model；(c)system energy change；(d)stress change in area A；(e)energy change of circle 1 with different thickness；(f)energy change of circle 2 with different thickness；(g)energy change of circle 3 with different thickness；(h)calculated values of stress intensity factors for different thickness

4 結(jié) 論

采用分子動力學(xué)模擬技術(shù)可以很好的模擬分析加載速率、力加載方式、不同溫度下力加載方式及模型厚度對裂紋擴展行為的微觀影響規(guī)律,對揭示固體材料的裂紋擴展微觀機理和探究裂紋控制技術(shù)具有重要理論指導(dǎo)意義。通過對超低溫SiO2脆性材料的裂紋擴展規(guī)律研究得出如下主要結(jié)論:

1.加載速率會引發(fā)裂紋擴展路徑的變化。力加載方式顯示,力小于某閾值時,微裂紋不會擴展,數(shù)值模擬的計算裂紋起裂應(yīng)力遠大于實際巖石的長期強度值,主要是由于實際巖石存在缺陷并受到復(fù)雜環(huán)境因素的影響。

2.由于溫度影響原子運動,從而影響勢函數(shù),原子之間作用力改變,裂紋擴展行為隨之改變,表現(xiàn)出裂紋總長和裂紋平均擴展速率隨溫度的增加而增加。

3.裂紋模型厚度增加后,系統(tǒng)內(nèi)原子數(shù)量增加,所以系統(tǒng)動能、勢能和總能量變化均增大。裂紋擴展形成表面增大,擴展后三個圓區(qū)域內(nèi)能量增幅變大,模型厚度增大的同時導(dǎo)致開裂過程中所需平均應(yīng)力強度因子計算值增大,這主要是由假設(shè)的裂紋形狀系數(shù)不變不符合實際情況導(dǎo)致,總體上模擬的計算結(jié)果和裂紋擴展規(guī)律基本合理。

分子動力學(xué)模擬技術(shù)為探究巖石類工程材料的脆裂蠕變微裂紋擴展機理提供了一種有效途徑,當(dāng)然目前勢函數(shù)的選擇和參數(shù)的確定具有一定經(jīng)驗性,對不同分子結(jié)構(gòu)和不同環(huán)境下勢函數(shù)選擇的研究還有待進一步深入。