多晶金屬材料斷裂過程數值模擬方法

孫兆瓏， 任克亮， 紀華， 丁麗宏， 張亞婷， 王婷， 楊佳， 李佶敏

(寧夏大學物理與電子電氣工程學院， 銀川 750021)

生產生活中廣泛使用的金屬材料，大多是由不規則形狀的眾多晶粒組成的多晶體，晶粒之間的邊界稱為晶粒邊界或晶粒間界，簡稱晶界[1]。在多晶體的形成過程中，晶界區域總會包含一些雜質原子，這些雜質原子會削弱或強化晶界處的結合力,進而影響材料的整體性能[2]。

對于多晶材料來說，晶界結構的強弱與分布情況是影響其材料性能的重要因素之一，針對這一問題，學者們分別從實驗、數值模擬等手段對其進行研究分析。路君等[3]通過總結前人實驗結果，討論分析了晶粒大小與各類多晶材料力學性能之間的關系，給出了Hall-Petch公式的適用范圍。辛存[4]通過晶體塑性有限元法(crystal plasticity finite element modeling，CPFEM)對Voronoi多晶體模型收拉伸載荷的情況進行了模擬，發現晶界處會產生應力集中、接觸壓力劇變等現象。戴澍等[5]將晶粒簡化為正六邊形，利用有限元法對多晶材料晶界處的應力情況進行了模擬分析，給出了晶界上正應力分布情況與晶界和外載荷方向夾角的關系。楊靜等[6]通過有限單元法中的內聚力單元模擬多晶材料的晶界屬性，對多晶材料的疲勞損傷破壞進行了模擬，預測了裂紋的萌生擴展情況。楊丹[7]用分子動力學方法探究了不同晶粒尺寸、孿晶含量等因素對納米多晶Al的力學性能與變形機制的影響。盧敏[8]分別通過分子動力學方法、有限元方法對NiTi合金的沿晶和穿晶斷裂行為進行了模擬分析，結果表明，晶粒尺寸越小，材料的斷裂韌性越好。

在目前較為成熟的研究手段中，通過具體實驗能夠最客觀地獲取晶粒結構與各項材料參數之間的關系，但成本較高。用有限元方法模擬時，能夠較好地得到材料內部的應力分布情況，但在模擬裂紋擴展時，要在裂紋尖端設置奇異單元、加密網格，求解應力強度因子，確定裂紋擴展方向，并同時要求裂紋頂點和裂紋邊與有限元的節點單元邊重合，相應的網格必須重新劃分，使計算的效率和精度同時降低。同時，多晶材料隨機多變的晶粒組織構型也對網格劃分造成了一定的難度。利用分子動力學方法對多晶材料進行分析時，可以得出微觀尺度下原子的行為和受力、能量物理參數，材料的晶體結構，并能夠給出外界載荷作用下微裂紋的演化過程，但分析過程需要消耗大量的計算資源。

Silling[9]于2000年提出了近場動力學(peridynamic)方法。該方法是一種基于非局部思想的無網格方法，通過積分方程求解代替微分方程求解，對解決損傷、斷裂等不連續、多尺度問題時具有明顯的優勢，目前已廣泛運用于各類材料的破壞分析中[10-11]。Askari等[12]通過近場動力學微彈性本構模型對多晶硅破壞過程進行了模擬，其結果相對傳統數值方法更接近實驗現象。伊航[13]基于近場動力學理論對功能梯度材料的破壞過程進行了模擬，分析了材料梯度對試件起裂與完全破壞之間的關系。姜曉偉等[14]基于常規態型近場動力學，引入基于能量的失效判斷準則，對復合材料的分層擴展過程進行模擬，其結果能夠較好地捕捉到“指甲蓋形”的分層前沿，與試驗結果達到了較好的一致性。孫杰等[15]通過在傳統鍵型近場動力學本構模型中引入非局部長程作用力強度尺寸效應的核函數修正項，部分消除了傳統鍵型近場動力學模型對泊松比的限制，并對多裂紋脆性板的破壞過程進行模擬分析，給出了初始裂紋數量、位置和方向對結構破壞形式、擴展路徑和承載能力的影響規律。

因此，現針對多晶金屬材料內部微結構的多樣性和損傷破壞過程的復雜性，基于Voronoi圖基本思想，在Fortran環境下對多晶材料的晶粒組織結構進行模擬，建立晶粒大小、均勻度可控的晶粒組織模型。結合近場動力學方法在處理不連續、多尺度等問題的優勢，對外界載荷作用下多晶材料內部裂紋擴展及結構破壞過程進行模擬，對晶界強度、晶粒大小等因素影響下多晶材料的斷裂性能進行研究。

1 多晶材料破壞分析的近場動力學模型

1.1 近場動力學理論

近場動力學理論的基本思想是在某一空間域R內，將實體模型離散成一系列包含著周圍實體模型的體積、密度等信息的物質點。這些物質點分別與近場域Hx(近場域半徑δ)內的其他物質點之間存在相互影響，這種影響被稱為“鍵”[11]，鍵的強度通過力密度f進行描述，如圖1所示。

u′、u為相鄰物質點位移矢量；x′、x為相鄰物質點坐標位置矢量；t為當前時刻圖1 物質點間矢量關系Fig.1 Vector relations between points of matter

f=f(u′-u,x′-x,t)

(1)

根據牛頓第二定律可得其運動方程

(2)

在近場動力學方法中，為了描述材料的破壞，需要定義“鍵”的伸長率s和極限伸長率s0。在式(2)中令x′-x=ξ，u-u′=η，其中ξ表示相鄰物質點的相對位置矢量，η表示相鄰物質點的相對位移矢量，則伸長率s可表示為

(3)

“鍵”的極限伸長率s0，其大小與該材料楊氏模量E、斷裂時的臨界能量釋放率Gc以及近場域半徑δ的選取有關，具體表達式[16]為

(4)

式(4)中：Gc與材料的斷裂韌性KIC存在以下關系：

(5)

在近場動力學分析中，通過引入一個與時間t相關的標量函數μ(ξ,t)來描述鍵的破壞，其形式為

(6)

當物質點之間鍵的伸長率s超過極限伸長率s0時，則定義“鍵”斷裂，此后該物質點對之間不再產生相互作用，此時μ(ξ,t)的值為0，反之為1。

在材料破壞分析中，通過求解某一物質點與其相鄰物質點間的破壞鍵的數量與初始鍵總數的比值來定義物質點x在t時刻的損傷度，其函數形式為

(7)

1.2 多晶材料晶界影響

多晶材料由很多個單晶顆粒與包圍單晶顆粒的晶界組成，由于晶界偏聚等現象致使材料在晶粒內部和晶界附近有不同的強度系數，材料在外界載荷作用下將發生晶間或穿晶破壞。為了考慮晶界對多晶材料斷裂性能的影響，通過引入晶界強度系數β[12]，并建立關系式為

(8)

對于同一晶粒內的物質點之間的“鍵”的破壞情況仍用式(6)判別；若兩物質點屬于不同晶粒，即結合式(8)，采用以下布爾函數[式(9)]判別。

(9)

相應的物質點的損傷判定則定義為

(10)

1.3 基于Voronoi圖的多晶材料幾何模型

在基于近場動力學方法研究多晶材料的斷裂性能時，模型結構中必須包含材料的微觀結構特點。Voronoi圖能夠基于一組可控的特定點將對平面進行分隔，能很好地反映晶體的微觀結構，目前被廣泛用于復雜問題的處理中[4,18]。用Voronoi方法來建立多晶體微觀結構圖，在保證晶粒隨機分布的同時，使得晶格內的每一點到晶粒中心點位置距離最短，在近場動力學分析時，把同一晶格內的所有物質點的屬性設定與中心點相同。

盡管Voronoi圖可以在Arcgis、Sketchup等軟件內自動生成等，但為了建立基于Voronoi圖的多晶材料破壞分析的近場動力學模型，需要在Fortran環境下編寫生成Voronoi圖的程序結構，具體過程如下。

(1)通過隨機抽樣并結合幾何模型尺寸，確定的一系列隨機點坐標位置(xi,yi)作為構建泰森多邊形的特征點，同時賦予這些特征點一個不同的顏色作為區分。根據Voronoi圖定義，這些點代表著晶粒的中心位置。

(2)在幾何模型尺寸范圍內以均勻分布的形式生成近場動力學物質點。

(3)計算每個物質點與所有特征點之間的距離，篩選出與當前物質點距離最短的特征點，并賦予相同的顏色。

(4)對求解域內的所有物質點重復步驟(3)，可將多晶材料的實體模型轉化為可用近場動力學方法分析的Voronoi多晶材料模型。

通過控制特征點的個數及最小間隔，就可以控制晶體材料中晶粒大小及均勻度，如圖2為設置不同特征點個數時所對應晶粒大小變化情況。其中不同顏色用于區別每個獨立的晶粒組織，其目的主要表征在顏色不同的色塊間存在晶界。

圖2 不同晶粒個數的晶體結構圖Fig.2 Crystal structure diagram of different number of grains

2 數值計算

為了研究多晶金屬材料的斷裂性能，以7xxx系高強度鋁合金為研究對象，材料參數如表1所示。

表1 材料參數Table 1 Material parameters

其中計算模型結構為長×寬=1 mm×1 mm，厚度0.002 mm的正方形板；在近場動力學分析中,取物質點間距Δx=0.002 mm, 近場域半徑δ=3.015Δx, 物質點的數量是250 000，時間步長Δt=0.25 ns，在模型中心預制長度為0.2 mm中心穿透裂紋并在豎直方向施加40 m/s速度載荷后，分別對不同晶界強度系數和不同晶粒大小所對應的模型的斷裂過程進行分析。

2.1 不同晶界強度系數時的斷裂特性

為了研究晶界強度對多晶材料斷裂性能的影響，在計算模型內隨機生成大小形狀不同的600個晶粒后(其中平均晶粒直徑d=46 μm)，對晶界強度系數β分別取1、2、4和0.5時多晶材料的斷裂過程進行模擬研究，得到了如圖3～圖6所示的損傷云圖。

從圖3～圖6所示的損傷云圖中，通過觀察物質點的損傷度變化情況，可較為清晰的反映裂紋擴展的動態過程，當β=1時，晶界強度與晶內強度一致，可認為模型是均質材料，晶界的存在不產生影響。此時隨著加載，模型中裂紋呈對稱形式沿水平方向兩端擴展，直到2 950時間步完全貫通模型；當β=2、β=4時，隨著加載，裂紋的擴展會受到晶界的阻礙，在晶界處會產生小幅偏轉及分叉，由于晶界分布的隨機性導致裂紋不再呈現對稱性，表現為穿晶破壞。與β=1相比，如果取相同加載時間，裂紋在β=2、β=4時擴展尺寸較短，結構完全破壞所需時間較長；當β=0.5時，材料在晶界處強度較低，隨著加載，裂紋會沿晶界擴展，表現為晶間破壞，與β≥1情形相比，裂紋擴展速度最快，材料完全破壞用時短。

圖3 β=1時裂紋擴展過程Fig.3 Damage cloud of the model when β=1

圖4 β=2時裂紋擴展過程Fig.4 Damage cloud of the model when β=2

圖5 β =4時裂紋擴展過程Fig.5 Damage cloud of the model when β=4

圖6 β =0.5時裂紋擴展過程Fig.6 Damage cloud of the model when β=0.5

在近場動力學方法中，通過判定某物質點與其相鄰物質點間“鍵”的破壞來模擬裂紋擴展。同時，也可以用每個物質點的損傷度進行累計來反映模型的破壞情況。如圖7所示，在晶界強度系數取不同值時，加載初期，材料損傷度均為水平直線，表示材料還未出現新的損傷。隨著加載的繼續，損傷度開始增加，對應裂紋開始擴展，擴展速率較慢，當繼續加載，裂紋進入快速擴展期。晶界強度系數不同，裂紋起裂時間不同。當β<1，裂紋開始擴展時間較早，β≥1時裂紋開始擴展時間較晚。

圖7 不同晶界強度系數累計損傷度隨時間變化Fig.7 Cumulative damage curves of different grain boundary strength coefficients over time

2.2 不同晶粒大小時的斷裂特性

為了研究晶粒大小對多晶材料斷裂性能的影響，在計算模型中，分別對晶粒平均直徑取65、46、38、33 μm 4種情況的斷裂性能進行了研究。由于實驗中7xxx系鋁合金多以準解理斷裂為主[19]，即脆性穿晶斷裂，本文選取其晶界強度系數β=4，對比分析了4種晶粒大小對應的損傷特性。

如圖8所示，隨著晶粒平均直徑的減小，晶界成分在模型中占比提高，對裂紋擴展的阻礙作用也在增強,在相同速度載荷下，裂紋完全貫穿模型所需時間增加,即材料的抗斷裂性能越好。

圖8 不同晶粒尺寸累計損傷度隨時間變化Fig.8 Cumulative damage curves of different grain sizes over time

在相同速度載荷作用下，通過求解不同晶粒大小所對應的材料完全破壞時間步，建立完全破壞時間步與晶粒平均直徑d倒數的平方的關系曲線如圖9所示，結果發現晶體材料的抗斷裂性能與晶粒平均直徑倒數的平方呈線性關系，且晶粒細化有助于提高多晶材料的力學性能。

圖9 晶粒尺寸對模型完全破壞時間的影響Fig.9 Effect of grain size on complete failure time of the model

3 結論

以近場動力學理論為基礎，引入Voronoi圖方法建立多晶材料的晶體結構圖，以Voronoi圖的特征點的個數來控制晶體的大小，通過晶界強度系數來區別晶粒內部和晶界處的強度關系，建立了描述晶體結構破壞的近場動力學模型。對7xxx系鋁合金的拉伸破壞過程進行模擬，得到以下結論。

(1)通過控制晶界強度系數β，該模型能夠模擬多晶材料穿晶斷裂或晶間斷裂兩種不同的斷裂模式，同時模擬結果表明材料在穿晶斷裂時的抗斷裂性能要優于晶間斷裂的抗斷裂性能。

(2)通過對不同晶粒大小模型的裂紋擴展模擬發現，在穿晶斷裂模式下，完全破壞時間步與晶粒平均直徑d倒數的平方呈正相關線性關系，即晶粒細化有助于提高材料的抗斷裂性能。