溫振川，寧歲婷，張　睿，高英俊

(廣西大學物理科學與工程技術(shù)學院，廣西南寧　530004)

?

納米級缺口起裂臨界條件的晶體相場模擬*

溫振川，寧歲婷，張睿，高英俊**

(廣西大學物理科學與工程技術(shù)學院，廣西南寧530004)

(College of Physical Science and Technology，Guangxi University，Nanning，Guangxi，530004，China)

摘要：【目的】研究初始缺口大小對裂紋起裂的影響。【方法】采用晶體相場(phase-field-crystal，PFC)模型研究單軸拉伸應(yīng)力作用下，材料裂紋起裂的動態(tài)演化過程，并結(jié)合自由能分析裂紋起裂與擴展過程中，不同缺口尺寸對裂紋起裂的影響?！窘Y(jié)果】沿x方向施加恒應(yīng)變，初始缺口y方向尺寸較小時，改變?nèi)笨趚方向尺寸，缺口并不起裂；當y方向尺寸達到一定值時，初始缺口會起裂。當y方向尺寸接近起裂臨界值時，x方向尺寸的增大會延緩缺口的起裂。當x方向尺寸改變量較小時，裂紋并不起裂，當增大到x方向尺寸臨界值時，缺口開始起裂并擴展。【結(jié)論】在應(yīng)力大小一定的情況下，缺口起裂與應(yīng)力作用方向垂直的缺口尺寸密切相關(guān)，只有與應(yīng)變作用方向垂直的尺寸達到臨界值，缺口才會發(fā)生起裂。

關(guān)鍵詞：晶體相場模型臨界缺口尺寸缺口起裂

0　引言

【研究意義】材料在生產(chǎn)制造和實際應(yīng)用的過程中，受加載作用下容易產(chǎn)生微裂紋，這些微裂紋在一定的條件下會發(fā)生擴展，從而影響材料的壽命[1]。早期的研究是從材料斷裂的宏觀表象出發(fā)，建立了斷裂力學理論和模型，來解釋裂紋萌生和擴展的相關(guān)行為[2]。但實際上，裂紋的萌生與擴展[3]在納米尺度就開始了。因此，深入研究材料的斷裂特征，就不能只停留在給出表象的宏觀尺度上，而應(yīng)當從納觀尺度分析裂紋萌生與擴展機理?！厩叭搜芯窟M展】材料在介觀和宏觀尺度上的性能，在很大程度上由其內(nèi)部的微觀缺陷所決定[4]，例如：空洞、空位、位錯、晶界以及微裂紋等。這些缺陷起因于在原子尺度發(fā)生的復雜非平衡動力學過程[5]。就目前所采用的多尺度模擬研究方法來說，對這些微觀特征結(jié)構(gòu)的實時演化過程的建模和模擬是一個重大的挑戰(zhàn)。最近，Elder等[5-7]基于密度泛函理論提出了晶體相場(phase-field-crystal，PFC)模型。該模型給出了新的自由能函數(shù)形式，引入周期序參量作為局域密度場，它將無序相(液相等)的密度場定義為一常數(shù)，將固相的密度場表示成周期性函數(shù)(波)的形式，進而通過周期的原子密度函數(shù)表現(xiàn)晶體的周期結(jié)構(gòu)。這樣的密度場就很自然地與彈性效應(yīng)、晶粒取向和位錯的運動等由周期性結(jié)構(gòu)產(chǎn)生的物理特性緊密的關(guān)聯(lián)起來[8]。目前，PFC方法已有許多具體的應(yīng)用，例如：模擬位錯滑移和攀移[9]、相結(jié)構(gòu)轉(zhuǎn)變[10]、異質(zhì)外延生長[11]、納米晶粒長大[12]等?！颈狙芯壳腥朦c】材料微裂紋擴展的研究目前主要集中在材料受到外部載荷累計改變作用下的疲勞裂紋尖端擴展行為以及裂紋分叉特性[13-14]，而應(yīng)用先進的PFC模型，研究初始缺口尺寸大小對裂紋起裂影響的工作還鮮有報道?！緮M解決的關(guān)鍵問題】應(yīng)用PFC方法模擬拉伸應(yīng)變作用下[15]，納米級裂紋起裂的動態(tài)演化過程，詳細分析預制中心缺口尺寸大小對裂紋起裂的影響。

1　PFC模型與方法

1.1體系的自由能密度函數(shù)

與傳統(tǒng)的相場模型不同，PFC模型采用具有周期性的局域原子密度場作為相場變量[6]。因此，PFC模型能夠反映晶體周期結(jié)構(gòu)特性以及原子尺度的行為。對于固態(tài)金屬材料，其原子的位置不依賴于時間，呈現(xiàn)規(guī)則排列。文獻[6]給出了既能夠描述晶體周期結(jié)構(gòu)，又能描述體系均勻結(jié)構(gòu)的密度函數(shù)，其表達式可寫成

(1)

式中，G為倒格子，r為空間位置矢量，an,m為Fourier系數(shù)。式(1)中等號右邊第一項反映晶格原子的周期性排列特征；第二項ρ0反映均勻無序相的原子密度分布平均值，是一個常量。此時體系無量綱的自由能函數(shù)F可以寫成

(2)

(3)

式中：A0是一個特定的常數(shù)，反映固相原子密度周期結(jié)構(gòu)的振幅。由能量函數(shù)取極小值得到：波矢q=2π/a，A0表達式為

(4)

式中，ρ0和γ為體系能量函數(shù)的兩個重要參數(shù)。由體系的能量密度函數(shù)可以計算并畫出體系的相圖。對于二維體系，相區(qū)有均勻無序相和晶體相，且晶體相又有六角格子相和條狀相兩種。文獻[6]已給出均勻相和條狀相的能量密度函數(shù)表達式，按照平衡相圖的計算方法可得到二維PFC相圖如圖1所示。

S、T、L分別表示條狀相、三角相、液相；圖(b)為圖(a)中虛線框的放大圖

S,T and L represent stripe phase,triangular phase and liquid phase,respectively；(b)is enlarged image of the dashed box of (a)

圖1PFC模型的二維相圖[6]

Fig.1Two-dimensional phase diagrams of PFC model[6]

1.2動力學方程與數(shù)值化處理

保守的原子密度場變量的演化可用與時間相關(guān)的Cahn-Hilliard動力學方程[6]描述

(5)

該動力學方程是六階非線性高階偏微分方程，在坐標空間求解十分復雜和繁瑣。因此，為了提高求解動力學方程(5)的效率，需要轉(zhuǎn)換到波矢空間求解。這里采用半隱式傅里葉法[6,17]求解，對式(5)變換到傅里葉空間[18]，整理可得

(6)

1.3拉伸應(yīng)變的施加以及參數(shù)設(shè)置

對樣品采用等應(yīng)變(等位移)單軸拉伸加載。在受到拉伸作用后，樣品中所有格點的坐標位置都會發(fā)生相應(yīng)的變化。因此，通過格點坐標位置的變化來反映應(yīng)變作用所引起的樣品變形。設(shè)變形前的樣品中格點位置的坐標用(x,y)表示，施加等應(yīng)變作用引起變形后，格點位置的坐標用(x′,y′)表示。樣品在y方向受到應(yīng)變εy的拉伸作用。變形前后坐標變換關(guān)系為x′=x和y′=y/(1+εy)。由于等應(yīng)變作用引起樣品變形，其內(nèi)部原子密度函數(shù)變?yōu)棣?x′,y′)=ρ(x,y/(1+εy))，因此，在考慮等應(yīng)變作用引起變形的情況下，原子密度函數(shù)表達式(4)將變換為ρ(x′,y′)，體系的能量也隨之發(fā)生相應(yīng)的變化，但控制演化的動力學方程的形式保持不變。將樣品變形后的原子密度ρ(x′,y′)代入動力學方程(5)中，采用半隱式傅里葉變換方法整理得到(6)，對其求解并進行可視化處理，即可得到裂紋擴展的演化圖。

本研究模擬不涉及具體材料的物性參數(shù)，所用參數(shù)均已無量綱化處理，并將連續(xù)空間離散為四方格子網(wǎng)格，采用周期性邊界條件。計算網(wǎng)格為512Δx×512Δy，空間步長Δx=Δy=π/3，時間步長Δt=0.05。根據(jù)單模近似下的二維相圖(圖1)，本研究選用無量綱密度為ρ0=0.49的三角相表征晶體相，對應(yīng)的溫度參量取為γ=-1.0，如圖1b中A點所示。為模擬中心裂紋，在樣品中心設(shè)置一定大小的缺口作為初始裂口，其參數(shù)設(shè)置為ρ0=0.79，γ=-1.0，如圖1b中B點所示。之所以選取該點是因為其對應(yīng)液相與三角相共存區(qū)域附近的原子密度，有利于三角相向裂紋結(jié)構(gòu)轉(zhuǎn)變。原子排列取向角θ=1°。

2　結(jié)果與分析

2.1缺口x方向起裂的臨界尺寸

在施加相同應(yīng)變情況下，不同的初始缺口尺寸會影響裂紋的萌生與擴展。如圖2所示，本研究施加拉應(yīng)變方向為單x方向，應(yīng)變量為ε=10%。固定缺口y方向長度為10個格子，改變?nèi)笨趚方向的尺寸，分別取10，12，14，16，18，20，22個格子(如表1樣品Ai所示)。圖3a～d為實驗條件下，樣品A6缺口形貌圖，從圖3可以看出，整個過程中初始缺口沒有起裂。圖4給出了不同尺寸初始缺口樣品的自由能曲線。由于采取恒應(yīng)變施加的方法，從該圖中可觀察到，在整個過程，自由能曲線呈水平直線，體系能量保持常量，這表明裂紋沒有起裂。由此可以推斷，缺口y方向長度固定且較小時，改變?nèi)笨趚方向的尺寸，不產(chǎn)生缺口起裂；但是隨著x方向尺寸的增大，同等應(yīng)變作用下，整個體系的自由能升高。

圖2施加單軸拉伸應(yīng)變的示意圖

Fig.2Effect diagram of single biaxial stress

表1樣品參數(shù)及起裂情況

Table 1Sample parameters and crack initiation

樣品Sample缺口尺寸NotchsizeΔx×Δy起裂CrackinitiationA110×10NA212×10NA314×10NA416×10NA518×10NA620×10NA722×10NB110×10NB210×12NB310×14NB410×16YB510×18YB610×20YB710×22YC110×14NC212×14YC314×14YC416×14YC518×14YC620×14YC722×14Y

注：N表示沒有起裂，Y表示發(fā)生起裂

Note:N means no crack initiation，Y means crack initiation

t表示時間步長t is time step size圖3　樣品A6缺口形貌Fig 3　The morphology of crack of sample A6

圖4樣品Ai體系的自由能變化曲線

Fig.4Free energy change curve of the evolution of sample Ai

2.2缺口y方向起裂的臨界尺寸

固定缺口x方向尺寸為10個格子，改變初始缺口y方向的尺寸，分別取10，12，14，16，18，20，22個格子(如表1樣品Bi所示)。圖5給出了不同尺寸缺口樣品的自由能曲線，從圖中可以觀察到，當y方向缺口尺寸小于等于14個格子時，自由能幾乎是一條水平的直線，反映了整個體系能量基本沒有變化，也表明缺口沒有起裂。但隨著y方向缺口尺寸增大，不同樣品能量升高；當y方向尺寸大于14個格子時，自由能曲線在時間步長大約為25萬的時候開始降低，這意味著缺口在該時刻開始起裂，體系自由能得到釋放。對于y方向缺口尺寸不同的起裂樣品，起裂時間與曲線降低速率變化相一致。圖6為樣品B6在x方向的應(yīng)變10%作用下的裂口起裂擴展演化圖。從圖6a～b中能看出，缺口沒有發(fā)生明顯的變化，圖6c中缺口出現(xiàn)明顯的細小裂縫，對應(yīng)圖4自由能曲線能量開始降低的時刻，說明缺口已經(jīng)開始起裂，是由于起裂擴展過程中釋放應(yīng)變能，導致自由能曲線降低。圖7顯示了樣品B1～B7在時間步長40萬步時的裂紋形貌，圖7a～c缺口隨著選取缺口參數(shù)的不同而有細微的變化，但無明顯起裂情況，符合圖4中自由能曲線1，2和3對應(yīng)的能量變化規(guī)律；圖7d～g這4個形貌圖中，裂紋的形貌基本一致，符合圖4中自由能曲線4，5，6和7對應(yīng)的能量變化規(guī)律。該結(jié)果與格里菲斯強度理論[19]基本一致，即裂紋擴展的應(yīng)變大小與垂直應(yīng)變作用方向的缺口尺寸密切相關(guān)[20]。

由2.1，2.2的結(jié)果可推斷，在單x方向施加固定的應(yīng)變量作用下，裂紋擴展主要是受初始缺口y方向尺寸變化的影響。只有y方向缺口尺寸達到臨界值，缺口才會起裂。

圖5　樣品Bi裂紋演化過程的體系自由能曲線Fig.5　Free energy curve of sample Bi

t表示時間步長

t is time step size

圖6x方向應(yīng)變10%作用下樣品B6的裂紋擴展演化

Fig.6The morphology evolution of crack propagation of sample B6under 10% of strain at x direction

圖7時間步長40萬步時樣品B1～B7的裂紋形貌

Fig.7The morphology of the sample B1～B7at the time step of 400 thousands

2.3缺口y方向臨界尺寸與起裂的關(guān)系

設(shè)置初始缺口y方向尺寸為14個格子，x方向尺寸分別取10，12，14，16，18，20，22個格子(如表1樣品Ci所示)，研究在單x方向施加固定的應(yīng)變量作用下，缺口x方向尺寸變化是否對樣品裂紋擴展情況產(chǎn)生影響。圖8給出了不同x方向尺寸樣品裂紋演化的自由能曲線。從圖中可以觀察到，當x方向尺寸取10個格子時，曲線是一條比較平滑的直線，這意味著整個樣品體系的能量幾乎沒有變化，裂紋沒有起裂；當x方向尺寸取大于10個格子時，自由能曲線均有明顯的下滑過程，能量顯著降低，但曲線走向并不是在同一時刻開始降低，這表明不同的x尺寸對應(yīng)的能量開始降低的起始步長不同。隨著缺口x方向尺寸增大，能量降低時間不斷延緩，缺口起裂時間不斷滯后。由此可知，當缺口y方向尺寸達到裂紋擴展臨界尺寸值時，x方向尺寸增加，會延緩缺口的起裂。

圖8樣品Ci裂紋演化過程的體系自由能曲線

Fig.8Free energy curve of evolution of crack propagation of sample Ci

從圖9中可以看出，當時間步長t達到15萬步(圖9b)時，初始缺口還沒有出現(xiàn)明顯的擴展，缺口邊緣已經(jīng)出現(xiàn)了原子層局部滑移。當時間步長t達到25萬步(圖9c)時，缺口已經(jīng)出現(xiàn)了明顯的裂紋擴展，并且隨著時間步長的增加，裂紋擴展越明顯。圖10顯示了樣品C1～C7在時間步長40萬步時的形貌，圖10a對應(yīng)圖8中的自由能曲線1，整個過程初始裂口沒有起裂；圖10b和圖10c中樣品的裂紋形態(tài)基本一致，分別對應(yīng)圖8中的自由能曲線2和3；圖10d～f中3個樣品的裂紋形態(tài)一致，但是不同于圖10b～c樣品裂紋擴展情況，圖10d，10e和10f中的樣品分別對應(yīng)圖8中的自由能曲線4，5和6；圖10g裂紋擴展情況落后于圖10其他已經(jīng)起裂樣品，與圖8中其對應(yīng)的自由能曲線7降低的拐點落后于其他自由能曲線的情況相符。由此可推斷，在施加單方向恒應(yīng)變及缺口y方向尺寸達到裂紋擴展臨界尺寸值的條件下，x方向尺寸改變量較小時，裂紋并不起裂；當增大到臨界值時，缺口開始起裂并擴展。

圖9單x方向10%應(yīng)變作用下樣品C6的裂紋擴展演化

Fig.9The morphology evolution of crack propagation of sample C6under 10% of strain in uniaxial x direction

圖10時間步長為40萬步時樣品C1～C7的形貌

Fig.10The morphology of the sample C1～C7at the time step of 400 thousands

3　結(jié)論

在沿x方向施加恒定10%應(yīng)變量條件下，缺口y方向尺寸遠小于臨界值時，改變?nèi)笨趚方向尺寸并使其增大，缺口不起裂；固定缺口x方向尺寸，改變?nèi)笨趛方向尺寸并使其增大，當缺口y方向尺寸增大到臨界值時，開始起裂并擴展；當缺口y方向尺寸接近臨界值時，則改變?nèi)笨趚方向尺寸時，當x方向尺寸較小時，裂紋并不起裂，隨著x增大，當達到x尺寸臨界值時，缺口開始起裂并擴展。表明在施加恒應(yīng)變作用下，缺口起裂與應(yīng)力作用方向垂直的缺口尺寸密切相關(guān)，只有與應(yīng)變作用方向垂直的尺寸達到臨界值，缺口才會發(fā)生起裂。

參考文獻：

[1]衣林,陳躍良,徐麗,等.金屬材料疲勞微裂紋的萌生與擴展研究[J].飛機設(shè)計,2012,32(2):63-67. YI L,CHEN Y L,XU L,et al.Research on initiation and propagation of fatigue micro-crack for metal materials[J].Aircraft Design,2012,32(2):63-67.

[2]鐘群鵬,周煜,張崢.裂紋學[M].北京:高等教育出版社,2015. ZHONG Q P,ZHOU Y,ZHANG Z.Crack Science[M].Beijing:Higher Education Press,2015.

[3]高英俊,羅志榮,黃禮琳,等.韌性材料的微裂紋擴展和連通的晶體相場模擬[J].中國有色金屬學報,2013,23(07):1892-1899. GAO Y J,LUO Z R,HUANG L L,et al.Phase field crastal modeling for microcrack propagation and connecting of ductile materials[J].The Chinese Journal of Nonferrous Metals,2013,23(07):1892-1899.

[4]劉曉波,徐慶軍,劉劍.鋁裂紋擴展行為的分子動力學模擬[J].中國有色金屬學報,2014,24(6):08. LIU X B,XU Q B,LIU J.Molecular dynamics simulation of crack propagation behavior of aluminum[J].The Chinese Journal of Nonferrous Metals,2014,24(6):08.

[5]ELDER K R,KATAKOWSKI M,HAATAJA M,et al.Modeling elasticity in crystal growth[J].Physical Review Letters,2002,88(24):245701.

[6]ELDER K R,GRANT M.Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals[J].Physical Review E,2004,70(5Pt1):051605.

[7]ELDER K R,PROVATAS N,BERRY J,et al.Phase-field crystal modeling and classical density functional theory of freezing[J].Physical Review B,2007,75(6):064107.

[8]BERRY J,GRANT M,ELDER K R.Diffusive atomistic dynamics of edge dislocations in two dimensions[J].Physical Review E,2006,73(3Pt1):031609.

[9]楊濤,陳錚,董衛(wèi)平.應(yīng)力誘發(fā)雙位錯組亞晶界湮沒的晶體相場模擬[J].金屬學報,2011,47(10):1301-1306. YANG T,CHEN Z,DONG W P,et al.Phase field crystal simulation of stress-induced annihilation of sub-grain boundary with double-array dislocation[J].Acta Metallurgica Sinica,2011,47(10):1301-1306.

[10]羅志榮,黃世葉,茹謝辛,等.晶體相場法模擬大角度晶界的變形過程[J].廣西科學,2013,20(4):311-315. LUO Z R,HUANG S Y,RU X X,et al.Phase-field-crystal modeling for deformation process of high-angle grain boundary[J].Guangxi Sciences,2013,20(4):311-315.

[11]高英俊,羅志榮,黃創(chuàng)高,等.晶體相場方法研究二維六角相向正方相結(jié)構(gòu)轉(zhuǎn)變[J].物理學報,2013,62(5):050507-050510. GAO Y J,LUO Z R,HUANG C G,et al.Phase-field-crystal modeling for two-dimensional transformation from hexagonal to square structure[J].Acta Physica Sinica,2013,62(5):050507.

[12]YU Y M,BACKOFEN R,VOIGT A.Morphological instability of heteroepitaxial growth on vicinal substrates:A phase-field crystal study[J].Journal of Crystal Growth,2011,318(1):18-22.

[13]高英俊,王江帆,羅志榮,等.晶體相場方法模擬納米孿晶結(jié)構(gòu)[J].計算物理,2013,30(4):577-581. GAO Y J,WANG J F,LUO Z R,et al.Nano-twin structure simulation with phase field crystal method[J].Chinese Journal of Computational Physics,2013,30(4):577-581.

[14]HU S,CHEN Z,PENG Y Y,et al.Modeling and simulation of microcrack propagation behavior under shear stress using phase-field-crystal[J].Computational Materials Science,2016,121:143-150

[15]GOMEZ H,NOGUEIRA X.An unconditionally energy-stable method for the phase field crystal equation[J].Computer Methods in Applied Mechanics and Engineering,2012,249-252:52-61.

[16]ELDER K R,ROSSI G,KANERVA P,et al.Patterning of heteroepitaxial overlayers from nano to micron scales[J].Physical Review Letters,2012,108(22):226102.

[17]盧成健,蔣麗婷,王玉玲,等.晶體相場法模擬小角度晶界的位錯結(jié)構(gòu)及其演化[J].廣西科學,2013,20(4):316-320. LU C J,JIANG L T,WANG Y L,et al.Dislocation structure evolution in low angle grain boundary[J].Guangxi Sciences,2013,20(4):316-320.

[18]毛鴻,羅志榮,黃世葉,等.材料裂紋擴展分叉機理的晶體相場法研究[J].廣西科學,2015(5):499-505. MAO H,LUO Z R,HUANG S Y,et al.Phase-field-crystal modeling for crack propagation and branch of materials[J].Guangxi Sciences,2015(5):499-505.

[19]張俊善.材料強度學[M].黑龍江:哈爾濱工業(yè)大學出版社,2004. ZHANG S J.Strength of Materials[M].Heilongjiang:Harbin Institute of Technology Press,2004.

[20]師俊平,胡義鋒,羅峰.裂紋起裂方向及擴展路徑的數(shù)值模擬[J].力學季刊,2013,34(02):233-239. SHI J P,HU Y F,LUO F.Numerical simulation of crack initiation directions and its propagation path[J].Chinese Quarterly of Mechanics,2013,34(2):233-239.

(責任編輯：陸雁)

Phase-field-crystal Simulation of Critical Conditions for the Crack Initiation of Nanoscale Notch

WEN Zhenchuan，NING Suiting，ZHANG Rui，GAO Yingjun

Key words：phase-field-crystal model，critical notch size，notch initiation

Abstract：【Objective】To study the effect of initial crack size on crack initiation.【Methods】The dynamic of process crack propagation under the uni-axial tension deformation in single crystals is simulated by the phase-field-crystal (PFC) model.The free energy is used to analyze the influence of different initial crack size on the crack initiation in the process of crack initiation and crack propagation.【Results】When the initial crack of y direction size is smaller,there is no initiation of crack under the uniaxial tension in the x axis,although the notch size in x direction changes.when the y direction size reaches a certain value,the initial cracks will expand.The time of crack initiation is delayed with the increase of x direction size.When the initial crack of x direction size changes less,there is no initiation of crack,and the initial cracks will expand when the x direction size reaches a certain value.【Conclusion】The initiation of crack is closely related to the crack size of vertical stress direction.The crack extension occurs,if only the size of notch reaches a certain value.

收稿日期：2016-07-25

作者簡介：溫振川(1991-)，男，碩士研究生，主要從事納米材料設(shè)計與模擬實驗研究。 **通信作者：高英俊(1962-)，男，教授，博士生導師，主要從事材料納微結(jié)構(gòu)的設(shè)計與模擬實驗研究，E-mail:gaoyj@gxu.edu.cn。

中圖分類號：TG111

文獻標識碼：A

文章編號：1005-9164(2016)05-0448-06

*國家自然科學基金項目(51161003，50661001)，廣西自然科學基金重點項目(2012GXNSFDA053001)和廣西大學大創(chuàng)項目(201610593220，201610593218)資助。

廣西科學Guangxi Sciences 2016,23(5):448～453

網(wǎng)絡(luò)優(yōu)先數(shù)字出版時間：2016-11-21【DOI】10.13656/j.cnki.gxkx.20161121.014

網(wǎng)絡(luò)優(yōu)先數(shù)字出版地址：http://www.cnki.net/kcms/detail/45.1206.G3.20161121.1546.028.html