高英俊， 羅志榮， 鄧芊芊， 黃禮琳， 林 葵

（1.廣西大學(xué)廣西有色金屬及特色材料加工重點(diǎn)實(shí)驗室，南寧 530004；2.廣西大學(xué)物理科學(xué)與工程技術(shù)學(xué)院，南寧 530004）

韌性材料的微裂紋擴(kuò)展與分叉的晶體相場模擬

高英俊1，2， 羅志榮1，2， 鄧芊芊1， 黃禮琳1， 林 葵1

（1.廣西大學(xué)廣西有色金屬及特色材料加工重點(diǎn)實(shí)驗室，南寧 530004；2.廣西大學(xué)物理科學(xué)與工程技術(shù)學(xué)院，南寧 530004）

采用晶體相場方法研究韌性單晶材料在雙軸拉伸條件下微裂紋擴(kuò)展與分叉的演化過程，分析應(yīng)變、溫度、初始裂口形狀等因素對裂紋擴(kuò)展和分叉的影響.結(jié)果表明：對于簡單的單向拉伸，應(yīng)變需要達(dá)到一定的臨界值，裂紋擴(kuò)展才會啟動.對于二組相互垂直的雙軸拉伸作用，當(dāng)應(yīng)變達(dá)到臨界值后，裂紋擴(kuò)展過程中會發(fā)生分叉現(xiàn)象.溫度越高，裂紋擴(kuò)展越快且分叉越多.裂紋在擴(kuò)展過程中，體系能量不斷降低，當(dāng)裂紋出現(xiàn)分叉時，體系能量降低得更快.在裂紋擴(kuò)展過程中，有時是會在裂尖處前方出現(xiàn)微小的空洞，類似在裂紋尖端前方出現(xiàn)位錯發(fā)射情況，這些微小的空洞逐漸擴(kuò)大連成裂紋.本文所得結(jié)果與相關(guān)模擬結(jié)果和實(shí)驗結(jié)果符合.

晶體相場方法；微裂紋擴(kuò)展；應(yīng)變；韌性材料

0 引言

材料在介觀和宏觀尺度上的性能在很大程度上由其內(nèi)部的微觀缺陷，例如，空位、空洞、位錯、晶界和微裂紋所決定.這些缺陷起因于在原子尺度發(fā)生的復(fù)雜非平衡動力學(xué)過程.就目前所采用的多尺度模擬方法來說，對這些微觀特征結(jié)構(gòu)的實(shí)時演化過程的建模和模擬是一個重大挑戰(zhàn).相場方法是當(dāng)今研究微觀結(jié)構(gòu)演化的強(qiáng)有力的數(shù)值計算方法［1－2］.傳統(tǒng)的相場方法［3－7］是建立在平衡態(tài)均勻場的基礎(chǔ)上的，忽略了許多由原子的周期排列結(jié)構(gòu)產(chǎn)生的物理特性，難以反映晶體學(xué)結(jié)構(gòu)特性以及原子尺度的行為信息，因而無法從根本上闡明微觀組織演化過程中原子尺度上的動力學(xué)機(jī)理.雖然，目前可用分子動力學(xué)方法［8－9］模擬預(yù)測這些微結(jié)構(gòu)演化動力學(xué)過程，且能準(zhǔn)確地描述這些現(xiàn)象，但該方法適用的特征時間尺度主要在原子振動時間尺度（10－12到10－14秒），不易拓展到擴(kuò)散時間尺度（10－6秒），并且，分子動力學(xué)要求加載的應(yīng)變率很大，通常達(dá)到107～1010s－1，與實(shí)際工業(yè)的應(yīng)變速度102s－1相差甚遠(yuǎn)［10，11］.最近，Elder等［11－13］基于密度泛函理論提出了晶體相場（phase-field-crystal，PFC）模型.該模型給出了新的自由能函數(shù)形式，引入的序參量為局域密度場，它將無序相（例如液相等）的密度場定義為一常數(shù)，將固相的密度場表示成周期性函數(shù)（波）的形式，進(jìn)而通過周期的原子密度函數(shù)表現(xiàn)晶體的晶格結(jié)構(gòu).這樣的周期結(jié)構(gòu)的密度場就很自然地與彈性效應(yīng)、晶粒取向和位錯的運(yùn)動等由周期結(jié)構(gòu)產(chǎn)生的物理特性緊密地關(guān)聯(lián)起來［13］.PFC模型既可以描述晶體學(xué)結(jié)構(gòu)特性以及原子尺度的行為，又可以揭示時間尺度為10－6秒量級的原子、缺陷運(yùn)動的行為特征.由此可以預(yù)見，PFC方法在今后的微結(jié)構(gòu)演化的研究中將展現(xiàn)其強(qiáng)大的優(yōu)勢，能夠廣泛地應(yīng)用到材料微結(jié)構(gòu)研究的各個領(lǐng)域，為研究原子尺度組織結(jié)構(gòu)的10－6秒量級的行為特征，例如研究晶體內(nèi)部的空位、位錯的運(yùn)動、晶界的遷移和微裂紋擴(kuò)展及其對材料宏觀性能的影響，提供先進(jìn)和強(qiáng)有力的計算方法.

目前，PFC方法已經(jīng)有一些具體的應(yīng)用，例如，模擬位錯滑移和攀移［10］、相結(jié)構(gòu)轉(zhuǎn)變［14］、異質(zhì)外延生長［15］、納米晶粒長大［16］等.PFC方法最新的研究前沿和發(fā)展方向［17－21］之一主要集中在：將加工變形的外界作用因素引入PFC模型中，建立形變和位錯動力學(xué)的PFC模型，用于超塑性、疲勞、斷裂、蠕變等行為和性能的機(jī)理研究.盡管傳統(tǒng)相場法和其它方法對材料的裂紋擴(kuò)展研究已有較多的報道［22－25］，但是應(yīng)用先進(jìn)的PFC方法開展微裂紋擴(kuò)展和分叉的研究工作還少見有報道.本文應(yīng)用PFC方法模擬在雙軸拉伸應(yīng)變作用下納米級微裂紋擴(kuò)展與分叉的演化過程，詳細(xì)分析應(yīng)變、溫度、初始缺口形狀等因素對微裂紋擴(kuò)展和分叉的影響.

1 PFC模型與方法

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

與傳統(tǒng)的相場模型不同，PFC模型的自由能密度函數(shù)是通過對原子密度場變量取極小值得到的［11］.因此，該模型能夠揭示晶體學(xué)結(jié)構(gòu)特性以及原子尺度的行為.對于固態(tài)金屬材料，其原子的位置不依賴于時間，呈規(guī)則排列.因此，要求相場變量必須能夠反映原子周期性排列的特征.引入周期性相場變量，其局域位置的最大值對應(yīng)于原子的位置；另一方面，液相中的原子位置隨時間隨機(jī)變化，取時間平均可看成為常量.符合這兩方面要求的相場變量定義，可用保守的原子密度場函數(shù)作為相場變量，其表達(dá)式［11］可寫成

式中，等號右邊第一項反映晶格原子的周期排列結(jié)構(gòu)特征，第二項反映液相的原子無規(guī)則運(yùn)動，平均值為一個常量.此時體系無量綱的自由能函數(shù)［11］可以寫成

式中，γ是反映體系溫度的參數(shù)，?2為Laplace算子.該能量模型自洽地包含了晶體結(jié)構(gòu)的物理特征，例如，多晶取向、彈性效應(yīng)、塑性變形等.

1.2 二維體系的極小自由能密度函數(shù)

其中，ρ0和γ為體系自由能函數(shù)的兩個重要參數(shù).由體系的極小自由能密度函數(shù)，可以計算并畫出體系不同相結(jié)構(gòu)區(qū)的相圖.對于二維體系，相區(qū)有均勻無序相和固相，且固相又有三角格子相和條狀相兩種，文獻(xiàn)［11］已給出三角格子相和條狀相的自由能函數(shù)形式.利用這些相的自由能函數(shù)，再按照平衡相圖的計算方法得到的二維相圖如圖1（a）所示.其中，L、T和S分別代表均勻無序相、三角格子固相和條狀固相.

1.3 演化動力學(xué)方程與數(shù)值化處理

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

式中，ζ為高斯隨機(jī)噪聲項，具有零平均值．本文不需考慮ζ的作用.

圖1 （a）PFC模型的二維相圖［11］；（b）（a）中方框的放大圖Fig.1 （a）2D phase diagram of PFCmodel（b）Enlarged map of the rectangle block in（a）

為求解復(fù)雜的動力學(xué)方程式（7），還必須將動力學(xué)方程在時間和空間上進(jìn)行離散化處理，即采用數(shù)值求解的辦法.在本文的數(shù)值求解中，采用顯型Euler迭代公式［26］

式中，Δx為空間步長，j和n分別代表i的最近鄰格點(diǎn)與次近鄰格點(diǎn).

1.4 施加應(yīng)變的方法

在受到拉伸作用后，樣品中所有格點(diǎn)的坐標(biāo)位置都會發(fā)生相應(yīng)的變化.因此，可通過格點(diǎn)的坐標(biāo)位置變換來反映應(yīng)變的作用所引起的樣品變形.設(shè)變形前的樣品中格點(diǎn)位置的坐標(biāo)用（x，y）表示，受拉伸作用引起變形后，格點(diǎn)位置的坐標(biāo)用（x′，y′）表示．樣品在x方向和y方向分別受到應(yīng)變εx和εy的拉伸作用下，變形前后坐標(biāo)變換關(guān)系為x′＝x／（1＋εx）和y′＝y(tǒng)／（1＋εy）．由于應(yīng)變作用引起樣品變形，其內(nèi)部原子密度函數(shù)變?yōu)棣眩▁′，y′） ＝ρ（x／（1＋εx），y／（1＋εy）），因此，在考慮應(yīng)變作用引起變形的情況下，原子密度函數(shù)表達(dá)式（4）將變換為ρ（x′，y′），體系的能量也隨之發(fā)生相應(yīng)的變化，但控制演化的動力學(xué)方程的形式保持不變．將樣品變形后的原子密度函數(shù)ρ（x′，y′）代入動力學(xué)方程（7）中，按照數(shù)值化計算步驟進(jìn)行數(shù)值求解，然后進(jìn)行可視化處理［5］，即可得到裂紋擴(kuò)展的演化細(xì)節(jié)圖.本文模擬不涉及具體材料的物性參數(shù)，所用參數(shù)均已無量綱化處理［11］，并將連續(xù)空間離散為四方格子，采用周期性邊界條件.計算區(qū)域網(wǎng)格為1 024×512 gp（gp表示格子點(diǎn)數(shù)），空間步長設(shè)為Δx′＝Δy′＝π／3，時間步長Δt＝0.05，其它參數(shù)在下節(jié)中介紹.

1.5 初始裂口的設(shè)置

根據(jù)圖1給出的相圖，選取體系自由能密度函數(shù)的溫度參數(shù)γ＝－1.0，原子密度參數(shù)ρ0＝0.49（圖1（b）中的A點(diǎn)），位于液相與三角格子固相共存區(qū).本文模擬單晶樣品，其原子的排列方向與x軸夾角θ為15°.在樣品的中心位置，挖出一個缺口作為初始裂紋的位置.為了便于比較，設(shè)置了兩種初始裂口，分別為20×10 gp的長方形缺口和直徑為16 gp的圓形缺口，缺口面積近似相等，其局部放大圖如圖2所示.缺口處周圍的原子環(huán)境如圖2（a）和2（c）所示，圖2（b）和2（d）為缺口周圍原子密度分布圖.缺口處的參數(shù)設(shè)為γ＝－1.0，ρ0＝0.79（圖1（b）中的B點(diǎn)），其正好處在液相與三角格子固相共存區(qū)附近的液相區(qū).選取該點(diǎn)的參數(shù)的好處在于有利于三角晶相結(jié)構(gòu)向裂紋結(jié)構(gòu)轉(zhuǎn)變.

圖2 初始裂口的局部放大圖（εx和εy分別為x方向和y方向的應(yīng)變.）Fig.2 Local enlarged zone of initial crack rim（εxandεyare strains in x and y direction.）

2 結(jié)果與分析

2.1 應(yīng)變對裂紋的影響

設(shè)置初始裂口為長方形，如圖2（a）所示.圖3為樣品在y方向單軸拉伸作用下經(jīng)過5.0×105ts（ts表示時間步數(shù)）演化的裂紋擴(kuò)展情況，其中y方向的應(yīng)變分別為εy＝7%、7.5%和8%.由圖3（a）可見，對于單軸拉伸情況，當(dāng)εy＜7.5%時，裂紋不擴(kuò)展；當(dāng)εy≥7.5%時，裂紋才開始擴(kuò)展．由此可見，裂紋沿x方向擴(kuò)展的臨界應(yīng)變約為7.5%.

圖3 5.0×105ts時裂紋擴(kuò)展情況（εx＝0%）Fig.3 Crack propagation at t＝5.0×105ts andεx＝0%

設(shè)置初始裂口為長方形，如圖2（a）所示.樣品同時受到x和y方向的雙軸拉伸作用，設(shè)沿y方向的應(yīng)變保持εy＝10%，沿x方向的應(yīng)變分別為εx＝0%、3%、5%、7%、7.5%和8%，經(jīng)過5.0×105ts演化的裂紋擴(kuò)展分布如圖4所示.由圖4（a）可見，當(dāng)εx＝0%時，裂紋擴(kuò)展方向沿著原子間距最寬的方向進(jìn)行，與x方向的夾角為15°，個別地方出現(xiàn)轉(zhuǎn)折（圖4（a1）），然后又恢復(fù)到原來的前進(jìn)方向；圖4（a1）給出了圖4（a）的裂紋波折區(qū)域的放大圖.由圖4（a1）可見，裂紋的邊緣與原子的排列方向平行，在轉(zhuǎn)折處轉(zhuǎn)折方向與原方向夾角約為60°.當(dāng)εx＜5%時，裂紋不出現(xiàn)分叉，存在鋸齒邊緣，裂紋擴(kuò)展的形狀像一把“劍”.當(dāng)εx＝5%時（圖4（c）），在裂紋前端轉(zhuǎn)折處，開始出現(xiàn)微小空洞（見圖4（c）中圓形區(qū)域）.當(dāng)εx＞7%時，裂紋在擴(kuò)展中開始出現(xiàn)分叉，εx越大，分叉越多．例如，當(dāng)εx＝8%時（圖4（f）），由主裂紋分叉為子裂紋，子裂紋又分叉為孫子裂紋，形成了多枝狀結(jié)構(gòu)，且裂紋越來越寬.由圖4中還可以看到在裂紋周圍出現(xiàn)灰白波紋分布，這是在裂紋周圍形成的原子位移的細(xì)微變化，反映了裂紋周圍彈性應(yīng)變場分布發(fā)生的畸變.圖5是體系能量隨時間的變化曲線.由此可見，應(yīng)變εx越大，能量降低就越快，表明裂紋擴(kuò)展也越快，分叉也越多．當(dāng)εx較小時，能量變化曲線近似為直線，表明這一過程中，裂紋擴(kuò)展過程不分叉．裂紋要出現(xiàn)分叉，施加的εx存在一個臨界值，在本文設(shè)定的條件下，由圖5可見，該εx的臨界值接近7%．當(dāng)εx＞7% 時，曲線下降速率明顯加快，表明開始出現(xiàn)裂紋顯著分叉．由圖4的裂紋變化情況，可以看到裂紋演化的形狀變化特征：開始時，施加的εx較小，裂紋形狀像一把“劍”；隨著應(yīng)變εx增加，裂紋出現(xiàn)分叉擴(kuò)展，逐漸長成“蝦”狀；最后，當(dāng)εx≥8%時（圖4（f）和（g）），出現(xiàn)多個分支結(jié)構(gòu)，裂紋形狀變成了“蟹”狀.文獻(xiàn)［28］也給出了類似的現(xiàn)象.

圖4 5.0×105ts時裂紋擴(kuò)展情況（εy＝10%）Fig.4 Crack propagation asεy＝10%at t＝5.0×105ts

圖6給出了圖4（b）中虛線方形區(qū)域內(nèi)裂紋尖端B的擴(kuò)展速率曲線.由該圖可見，裂紋沿直線擴(kuò)展過程中，速率較快，在轉(zhuǎn)折處，裂紋擴(kuò)展速率減小；從轉(zhuǎn)折到恢復(fù)直線擴(kuò)展，速率有一個加速突變的過程，然后又恢復(fù)到近似勻速運(yùn)動過程.在轉(zhuǎn)折處，出現(xiàn)速率振蕩現(xiàn)象［29］，見圖6中速率的峰和谷.由圖5的變化趨勢和圖6速率的變化特點(diǎn)可以預(yù)測，當(dāng)應(yīng)變εx大于臨界值后繼續(xù)增加，則裂紋擴(kuò)展速度加快，導(dǎo)致分叉更快、更多，將使得裂紋擴(kuò)展速率出現(xiàn)更復(fù)雜的振蕩情況.類似文獻(xiàn)［29］指出的裂紋擴(kuò)展速率超過臨界值后，發(fā)生裂紋分叉，裂紋擴(kuò)展速率出現(xiàn)振蕩現(xiàn)象.

圖5 不同εx對應(yīng)的能量曲線（εy＝10%）Fig.5 Energy curves atεy＝10%and differentεx

圖6 在圖4（b）中虛線方形區(qū)域內(nèi)裂紋尖端B的裂紋擴(kuò)展速率振蕩曲線Fig.6 Oscillations of crack propagation in Fig.4（b）rectangle block at crack tip B

2.2 溫度參數(shù)對裂紋擴(kuò)展的影響

設(shè)置初始裂口為長方形，如圖2（a）所示，在裂口處取ρ0＝0.79，溫度參數(shù)分別取γ＝－1.0、－0.95和－0.90（分別對應(yīng)圖1（b）中的B點(diǎn)、B1點(diǎn)和B2點(diǎn)）得到3類樣品.對樣品施加的應(yīng)變?yōu)棣舩＝8% 和εy＝10%.圖7給出了演化4.0×10 ts時不同溫度參數(shù)取值情況下的裂紋擴(kuò)展形貌圖.由圖7（c）可見，γ＝－0.9對應(yīng)溫度較高的情況，裂紋擴(kuò)展較快，分叉也較多；γ＝－1.0對應(yīng)溫度較低的情況，裂紋擴(kuò)展較慢，分叉也較少．圖8給出了三種溫度參數(shù)值情況下的體系能量隨時間的變化曲線．由圖可見，體系能量隨著演化時間增加而減?。茫剑?.9時，能量曲線減小最快；γ＝－1.0時，能量減少最慢.在最初的演化時間里，能量曲線隨時間線性地減小，經(jīng)過一段時間后，曲線開始向下彎曲，能量降低有加速趨勢，表明從此時開始，裂紋開始分叉，隨后的裂紋擴(kuò)展過程中，出現(xiàn)更多的裂紋分叉；裂紋擴(kuò)展分叉越多，體系能量降低得越快.能量曲線降低的本質(zhì)原因是由于彈性應(yīng)變能的釋放要比裂紋界面能（表面能）引起的能量增加要大得多，使得體系總能量隨著裂紋長大不斷降低，裂紋得以擴(kuò)展.

圖7 4.0×105ts時取不同溫度參數(shù)值γ對應(yīng)的裂紋擴(kuò)展形貌圖Fig.7 Crack propagation with different temperature parameterγat t＝4.0×105ts

2.3 初始裂口形狀及周圍原子密度的影響

設(shè)置的初始裂口為圓形，如圖2（c）所示.同時施加的應(yīng)變?yōu)棣舩＝8%和εy＝10%，溫度參數(shù)γ＝－1.0．為研究初始裂口處的原子密度的影響，裂口處的初始原子密度ρ0設(shè)置了兩種值，分別為0.79（圖1（b）中的B點(diǎn)）和0.90（圖1（b）中的C點(diǎn)）.圖10（a）和（b）分別給出了這兩種ρ0參數(shù)下演化5.0×105ts時的裂紋分叉形貌圖.從圖10可見，ρ0＝0.90時，裂紋分叉較多，分枝較豐富，沿x方向和沿y方向都發(fā)展較快，且發(fā)展較均衡；而ρ0＝0.79時，相對而言裂紋沿x方向的發(fā)展較快，而沿y方向擴(kuò)展較慢．由圖9能量曲線變化可見，體系能量減小較快的是ρ0＝0.90情況．因為初始裂口處原子密度ρ0＝0.90比ρ0＝0.79偏離相圖的晶相與均勻無序相共存區(qū)更遠(yuǎn)，不利于三角格子固相的穩(wěn)定，而有利于裂紋的擴(kuò)展.由圖10還可見，在裂紋的擴(kuò)展過程中，多個裂紋的子分枝出現(xiàn)鋸齒折線結(jié)構(gòu).在分叉的地方，有的是先在分叉處前方出現(xiàn)微小的空洞，類似在裂紋尖端前出現(xiàn)位錯發(fā)射，然后這些微小的空洞逐漸連成裂紋.這是由于裂紋尖端處存在非線性彈性與塑性變形作用而引起的現(xiàn)象，文獻(xiàn)［30］也觀察到類似的結(jié)果.

圖8 取三種不同溫度參數(shù)值（γ＝－1.0，－0.95和－0.90）對應(yīng)的能量曲線Fig.8 Energy curves with different parameter（γ＝ －1.0， －0.95 and －0.90）

圖9 初始圓形裂口處ρ0分別為0.79和0.90時的能量曲線Fig.9 Energy curves asρ0＝0.79 and 0.90 at circular crack rim

圖10 圓形初始裂口在應(yīng)變εx＝8%和εy＝10%情況下的裂紋形貌圖（5.0×105ts）Fig.10 Crack propagation asεx＝8%andεy＝10%at circular crack rim

圖10（a）與圖4（f）所對應(yīng)的設(shè)置條件，除了初始裂口形狀不同，其它條件相同.對比兩圖可以發(fā)現(xiàn)，圓形裂口的裂紋沿x方向的擴(kuò)展占有較明顯的優(yōu)勢，而長方形缺口的裂紋較均勻的向四周擴(kuò)展.主要原因是PFC方程的高度非線性和初始裂口形狀不同將會導(dǎo)致裂紋周圍的原子作用環(huán)境不同，可見，初始裂口形狀，對裂紋擴(kuò)展也有明顯的影響［24］.

3 結(jié)論

采用PFC方法能夠在原子空間尺度和擴(kuò)散時間尺度較好地研究分析模擬韌性材料的裂紋擴(kuò)展和分叉的詳細(xì)過程和細(xì)節(jié).主要結(jié)論如下：

1）裂紋的擴(kuò)展方向沿著原子間距最寬的方向進(jìn)行，施加的應(yīng)變存在一個臨界值，超過該臨界值，裂紋擴(kuò)展才得以啟動；對于裂紋分叉，需要x和y方向的應(yīng)變達(dá)到臨界應(yīng)變值才能出現(xiàn)裂紋分叉.

2）體系的溫度參數(shù)取值較小時，則裂紋發(fā)展較快，分叉也較多.裂紋在擴(kuò)展過程中，體系能量在不斷下降；當(dāng)裂紋出現(xiàn)分叉時，體系能量減少得更快.裂紋的子分枝出現(xiàn)鋸齒折線結(jié)構(gòu)，可看成裂紋擴(kuò)展過程的振蕩現(xiàn)象.在裂紋擴(kuò)展過程中分叉的地方，有時是先在裂尖處前方出現(xiàn)微小的空洞，類似在裂紋尖端前方出現(xiàn)位錯發(fā)射，然后逐漸連成裂紋.裂口處原子密度對裂紋擴(kuò)展有明顯影響.

3）裂紋擴(kuò)展的形狀變化特征：當(dāng)施加單軸拉伸，應(yīng)變較小時，裂紋形狀像一把“劍”；同時施加x和y方向的雙軸拉伸作用，隨著應(yīng)變增加，裂紋形狀轉(zhuǎn)變成“蝦”狀；最后，裂紋出現(xiàn)多個分支結(jié)構(gòu)，裂紋形狀變成了“蟹”狀.

［1］Moelans N，Blanpain B，Wollants P.An introduction to phase-fieldmodeling ofmicrostructure evolution［J］.Calphad，2008，32：268－294.

［2］Wang Y，Li J.Phase field modeling of defects and deformation［J］.Acta Mater，2010，58：1212－1235.

［3］Chen LQ.Phase-field models formicrostructure evolution［J］.Annu Rev Mater Res，2002，32：113－140.

［4］Steinbach I.Phase-field models in materials science［J］.Modell Simul Mater Sci Eng，2009，17：73001.

［5］高英俊，羅志榮，胡項英，等.相場方法模擬AZ31鎂合金的靜態(tài)再結(jié)晶過程［J］.金屬學(xué)報，2010，46：1161－1172.

［6］羅炳池，王海鵬，魏炳波.液態(tài)三元Ni-Cu-Pb偏晶合金相分離的相場模擬 ［J］.科學(xué)通報，2009，54：7－11.

［7］王錦程，李俊杰，楊玉娟，等.定向凝固界面形態(tài)演化及其穩(wěn)定性的相場法研究 ［J］.中國科學(xué)：E輯，2009，38：16－23.

［8］Millett P C，Selvam R P，Saxena A.Molecular dynamics simulations of grain size stabilization in nanocrystallinematerials by addition of dopants［J］.Acta Meterialia，2006，54：297－303.

［9］Holm E A，F(xiàn)oiles SM.How grain growth stops：A mechanism for grain-growth stagnation in pure materials［J］.Science，2010，328：1138.

［10］楊濤，陳錚，董衛(wèi)平.應(yīng)力誘發(fā)雙位錯組亞晶界湮沒的晶體相場模擬［J］.金屬學(xué)報，2011，47：1301－1306.

［11］Elder K R，Grant M.Modeling elastic and plastic deformations in nonequilibrium processing using phase field crystals［J］. Phys Rev E，2004，70：051605－20.

［12］Elder K R，KatakowskiM，Haataja M，GrantM.Modeling elasticity in crystal growth［J］.Phys Rev Lett，2002，88：245701－5.

［13］Berry J，GrantM，Elder K R.Diffusive atomistic dynamics of edge dislocations in two dimensions［J］.Phys Rev E，2006，73：031609－16.

［14］高英俊，羅志榮，黃創(chuàng)高，等.晶體相場法研究二維六角結(jié)構(gòu)向正方結(jié)構(gòu)轉(zhuǎn)變 ［J］.物理學(xué)報，2013，62（5）：050507－010.

［15］Yu Y M，Backofen R，Voigt A.Morphological instability of heteroepitaxial growth on vicinal substrates：A phase-field crystal study［J］.JCryst Growth，2011，318：18－22.

［16］Gao Y J，Wang JF，Luo ZR，Lu Q H，Liu Y.Nano-twin structure simulation with phase field crystalmethod［J］.Chinese JComputational Physics，2013，30（4）：577－581.

［17］Stefanovic P，Haataja M，Provatas N.Phase field crystal study of deformation and plasticity in nanocrystallinematerials［J］. Phys Rev E，2009，80：046107.

［18］Chan PY，Tsekenis G，Dantzig J，et al.Plasticity and dislocation dynamics in a phase field crystalmodel［J］.Phys Rev Lett，2010，105：015502.

［19］Huang Z F，Elder K R，Provatas N.Phase-field-crystal dynamics for binary systems：Derivation from dynamical density functional theory，amplitude equation formalism，and applications to alloy heterostructures［J］.Phys Rev E，2010，82：21605.

［20］Jaatinen A，Ala-Nissila T.Eighth-order phase-field-crystalmodel for two-dimensional crystallization［J］.Phys Rev E，2010，82：061602.

［21］Hakim V，Karma A.Laws of crackmotion and phase-fieldmodels of fracture［J］.JMech Phys Solid，2009，57：342－368.

［22］Spatschek R，Brener E，Karma A.Phase field modeling of crack propagation［J］.Philos Mag，2011，91：75－95.

［23］高英俊，羅志榮，黃禮琳，林葵.韌性材料的微裂紋擴(kuò)展和連通的相場方法模擬 ［J］.中國有色金屬學(xué)報，2013，23（7）：1892－1899.

［24］Huang S，Zhang S，Belytschko T，et al.Mechanics of nanocrack：Fracture，dislocation emission，and amorphization［J］.J Mech Phys Solids，2009，57：840－850.

［25］Livne A，Bouchbinder E，Svetlizky I，et al.The near-tip fields of fast cracks［J］.Science，2010，327：1359－1363.

［26］Chen L Q，YangW.Computer simulation of the domain dynamics of a quenched system with a large number of nonconserved order parameters：The grain-growth kinetics［J］.Phys Rev B，1994，50：15752－15756.

［27］Oono Y，PuriS.Computationally efficientmodeling ofordering ofquenched phases［J］.Phys Rev Lett，1987，58：836－839.

［28］Imbeni V，Kruzic J J，Marshall G W，et al.The dentin—enamel junction and the fracture of human teeth［J］.Nature Materials，2005，4：229－232.

［29］Sharon E，F(xiàn)ineberg J.Confirming the continuum theory of dynamic brittle fracture for fast cracks［J］.Nature，1999，397：333－335.

［30］Buehler MJ，Abraham F F，Gao H J.Hyperelasticity governs dynamic fracture at a critical length scale［J］.Nature，2003，426：141－146.

date：2013－08－21；Revised date：2013－11－25

Phase-Field-Crystal Modeling of Microcrack Propagation and Branching in Ductile Materials

GAO Yingjun1，2，LUO Zhirong1，2，DENG Qianqian1，HUANG Lilin1，LIN Kui1
（1．Guangxi Key Laboratory for Non-ferrousMetal and Featured Materials，Guangxi University，Nanning 530004，China；2．College of Physics Science and Engineering，Guangxi University，Nanning 530004，China）

Morphology ofmicrocrack propagation and crack branching of single crystal ductile materials under tensile are simulated with phase-field-crystalmethod.It shows that shape of notch has a strong effect on crack propagation.With uniaxial tension crack propagation appears if only strain reaches a critical value.With biaxial tension a crack branches if only strain reaches a critical value. Temperature parameter takes significant effects on crack propagation.The smaller the temperature parameter is（the temperature is higher），the faster the cracks propagate and themore the crack branches are.It is observed that free energy of system decreaseswith time and decreases faster after crack starts to branch.During crack propagation，disconnectedly isolated cavities around themain cracks can be found and these defects become new cracks.They grow up along a line and become a branching crack.Simulated results are in agreementwith other simulations and experimental results.

phase-field-crystalmethod；microcrack propagation；strain；ductilematerials

TG111.5

A

1001-246X（2014）04-0471-08

2013－08－21；

2013－11－25

國家自然科學(xué)基金（51161003）；廣西自然科學(xué)基金重點(diǎn)項目（2012GXNSFDA053001）；廣西大學(xué)廣西有色金屬及特色材料加工重點(diǎn)實(shí)驗室開放基金（GXKFJ12-01）及廣西大學(xué)科研基金（XJZ110611）資助項目

高英?。?962－），男，博士，教授，博士生導(dǎo)師，主要從事材料微結(jié)構(gòu)設(shè)計與計算模擬實(shí)驗研究，E-mail：gaoyj＠gxu.edu.cn