基于內(nèi)聚力界面的鎢合金兩相細(xì)觀結(jié)構(gòu)拉伸模擬

楊梓清

摘 要：在鎢合金粉末冶金制備工藝的基礎(chǔ)上，本研究建立了一種基于顆粒擠壓方法的鎢合金兩相結(jié)構(gòu)模型。首先在光學(xué)顯微鏡下觀察了鎢合金的微觀組織結(jié)構(gòu)，并統(tǒng)計了其顆粒分布特征，隨后利用隨機(jī)順序吸附算法生成符合實際分布的顆粒初始構(gòu)型，然后利用Python腳本程序在ABAQUS軟件中建立顆粒擠壓模型，以提高顆粒的體積分?jǐn)?shù)。當(dāng)模型模擬出的顆粒體積分?jǐn)?shù)達(dá)到鎢合金中鎢顆粒真實的體積分?jǐn)?shù)時，提取顆粒邊界輪廓，并運(yùn)用布爾運(yùn)算生成粘結(jié)相。最后通過Python腳本程序在顆粒界面生成內(nèi)聚力單元，結(jié)合有限元方法，將上述模型運(yùn)用于鎢合金的拉伸模擬，確定了內(nèi)聚力單元的參數(shù)，并與試驗結(jié)果進(jìn)行了對比驗證，討論了內(nèi)聚力參數(shù)對拉伸應(yīng)力-應(yīng)變曲線的影響。

關(guān)鍵詞：鎢合金;細(xì)觀結(jié)構(gòu);擠壓模型;內(nèi)聚力單元;有限元法;拉伸仿真

中圖分類號：TB33;TB302.3文獻(xiàn)標(biāo)識碼：A 文章編號：1003-5168（2021）15-0014-06

Abstract： Based on the powder metallurgy preparation process of tungsten alloy， this study establishes a dual-phase structure model of tungsten alloy based on particle extrusion method. First， the microstructure of the tungsten alloy is observed under an optical microscope， and its particle distribution characteristics are counted， subsequently， the random sequential adsorption algorithm is used to generate the initial configuration of the particles in line with the actual distribution， and then the particle extrusion model is established in the ABAQUS software using the Python script program to increase the volume fraction of the particles. When the particle volume fraction simulated by the model reaches the true volume fraction of tungsten particles in the tungsten alloy， the boundary contours of the particles are extracted and the Boolean operation is used to generate the binder phase. Finally， a cohesive force unit is generated on the particle interface through a Python script program， combined with the finite element method， the above model is applied to the tensile simulation of tungsten alloy， the parameters of the cohesion unit are determined， the test results are compared and verified， and the influence of cohesion parameters on the tensile stress-strain curve is discussed.

Keywords： tungsten alloy;microstructure;extrusion model;cohesion element;finite element method;tensile simulation

材料的宏觀力學(xué)特性往往會受到其微觀結(jié)構(gòu)特征的影響，特別是高比重鎢合金這種擁有復(fù)雜微觀結(jié)構(gòu)的兩相材料。已有研究表明，鎢合金的力學(xué)性能受其微觀結(jié)構(gòu)特征的影響[1-3]。比如，鎢合金的平均顆粒半徑會影響鎢合金的拉伸性能，平均顆粒尺寸越小，鎢合金的屈服強(qiáng)度和抗拉強(qiáng)度越大，表現(xiàn)出細(xì)晶強(qiáng)化效應(yīng);鎢合金內(nèi)部鎢顆粒分布不均勻，會導(dǎo)致其拉伸性能下降;鎢合金的鎢含量增加，粘結(jié)相體積相對減少，鎢-鎢連接度增加，從而導(dǎo)致其拉伸性能和沖擊性能降低。因此，在研究鎢合金復(fù)雜的力學(xué)性能時，不能忽略其微觀結(jié)構(gòu)組織特征的影響。如今，科學(xué)技術(shù)飛速發(fā)展，有限元方法已經(jīng)變得越來越成熟和完善，能夠解決許多工程問題，而對于具有兩相結(jié)構(gòu)的鎢合金，有限元仿真不應(yīng)該局限于宏觀尺度，而應(yīng)該更進(jìn)一步從細(xì)觀角度考慮其微結(jié)構(gòu)的效應(yīng)。因此，建立真實有效的鎢合金微觀結(jié)構(gòu)仿真模型具有重要意義，在細(xì)觀尺度上研究鎢合金的結(jié)構(gòu)參數(shù)可以更好地解釋鎢合金力學(xué)性能差異的原因，這對鎢合金加工制備方法的改進(jìn)有一定指導(dǎo)作用。

1 鎢合金微觀結(jié)構(gòu)

本文研究的鎢合金為95W-3.5Ni-1.5Fe合金，其化學(xué)成分組成如表1所示。根據(jù)鎢合金中每種元素的質(zhì)量分?jǐn)?shù)和密度可得，鎢合金中，鎢（W）顆粒的體積分?jǐn)?shù)為89.8%，鎳（Ni）和鐵（Fe）顆粒的體積分?jǐn)?shù)之和為10.2%。首先，將鎢合金表面拋光，并在光學(xué)顯微鏡下觀察，其微觀組織結(jié)構(gòu)如圖1（a）所示?；疑糠质墙茍A形的鎢顆粒，在粉末冶金過程中會稍微發(fā)生變形，而黑色部分為粘結(jié)相。鎢顆粒被粘結(jié)相包圍，二者之間存在明顯的界面。之后，統(tǒng)計鎢顆粒尺寸的分布特征，其直徑主要分布在10～60 [μm]，并近似服從正態(tài)分布，均值（[μ]）為28.5 [μm]，標(biāo)準(zhǔn)差（[σ]）為10.553，統(tǒng)計數(shù)目（[N]）為700，如圖1（b）所示。

2 鎢合金微觀結(jié)構(gòu)模型建立過程

鎢合金兩相微觀結(jié)構(gòu)模型的建立過程主要分三個步驟，分別是顆粒擠壓模型的建立、輪廓重建與幾何修復(fù)、內(nèi)聚力單元的插入。

2.1 顆粒擠壓模型

根據(jù)鎢合金粉末冶金的制備過程，研究人員在ABAQUS軟件中建立了二維顆粒擠壓模型，如圖2所示。首先，根據(jù)顆粒分布特征，利用RSA加密算法在MATLAB軟件中生成初始顆粒。顆粒的位置（用[x]和[y]坐標(biāo)表示）是在指定區(qū)域中隨機(jī)生成的，顆粒的尺寸（用[d]表示）也是隨機(jī)生成的，且服從上述提到的正態(tài)分布。因此，每次生成的顆粒都可以用一組隨機(jī)數(shù)（[x]、[y]和[d]）表示。如果當(dāng)前生成的顆粒不與之前已生成的顆粒重疊（兩個顆粒之間的距離大于它們的半徑和），則記下該顆粒的信息（[x]、[y]和[d]）。如果當(dāng)前生成的顆粒與已生成的顆粒重疊，則將其舍棄，并重新生成顆粒，再做判斷，直到所有顆粒的體積分?jǐn)?shù)達(dá)到要求。隨后，輸出所有顆粒的信息，編寫Python腳本程序?qū)㈩w粒信息輸入ABAQUS軟件中。接著，在ABAQUS軟件中創(chuàng)建一個可以包含所有顆粒的四邊形剛體，剛體的上邊界可以向下移動，以實現(xiàn)擠壓過程，而其他邊界則完全固定。將每個顆粒設(shè)置為彈性體，并任意給定合理的彈性模量和泊松比。網(wǎng)格劃分和接觸設(shè)置都通過Python腳本實現(xiàn)。網(wǎng)格最小尺寸為0.001 mm，網(wǎng)格劃分得越細(xì)，最終模型的輪廓越精確。顆粒之間、顆粒與剛性邊界之間都設(shè)置相互接觸。

隨著上邊界的連續(xù)向下移動，顆粒逐漸被擠壓變形，顆粒的體積分?jǐn)?shù)逐漸上升。當(dāng)所有顆粒的體積分?jǐn)?shù)達(dá)到鎢合金中鎢顆粒真實的體積分?jǐn)?shù)時，提取所有顆粒的輪廓，輸出每個顆粒邊界上所有的節(jié)點坐標(biāo)。為了保證邊界上的節(jié)點按照順序逐個輸出，在擠壓模型運(yùn)算之前，要將每個顆粒邊界上所有的節(jié)點從小到大重新編號。

2.2 輪廓重建與幾何修復(fù)

完成擠壓模型后，依次輸出變形后所有顆粒邊界上的節(jié)點坐標(biāo)，再編寫Python腳本程序?qū)⑵浒错樞蜉斎階BAQUS軟件中，以完成建模。輪廓重建的模型如圖3（a）所示。隨后，創(chuàng)建一個可以包含所有變形顆粒的矩形實體，并對所有顆粒執(zhí)行布爾運(yùn)算，以獲得粘結(jié)相的模型，如圖3（b）所示。至此，鎢合金的兩相模型已初步建立。但是，該方法中顆粒的邊界是由輸入的節(jié)點依次連接形成的，由于計算機(jī)本身的數(shù)值誤差，相鄰顆粒之間可能會存在微小的間隙，這不利于后續(xù)的網(wǎng)格劃分，并且會嚴(yán)重影響網(wǎng)格的質(zhì)量。因此，在劃分網(wǎng)格之前需要進(jìn)行幾何修復(fù)。

將上述的顆粒和粘結(jié)相模型保存為IGES文件格式，并將其導(dǎo)入Hypermesh軟件中，使用其幾何修復(fù)功能來清理重復(fù)的邊和狹窄的邊，并適當(dāng)調(diào)整幾個顆粒相交的狹小區(qū)域，以提高網(wǎng)格質(zhì)量。然后，劃分網(wǎng)格，盡管四邊形單元的求解精度高于三角形單元，但在相鄰顆粒的區(qū)域四邊形單元網(wǎng)格質(zhì)量比三角形單元差，因此使用三角形單元。接下來，為每個顆粒和粘結(jié)相創(chuàng)建單獨(dú)的單元集，以便隨后賦予材料屬性。最后，將上述模型導(dǎo)出為INP文件格式，該文件類型可以直接通過ABAQUS軟件讀取。最終在ABAQUS軟件中顯示的網(wǎng)格模型如圖4所示。

2.3 內(nèi)聚力單元的插入

在獲得了包含顆粒和粘結(jié)相的網(wǎng)格模型后，進(jìn)一步考慮顆粒邊界的存在，在顆粒-顆粒界面和顆粒-粘結(jié)相界面插入一層零厚度內(nèi)聚力單元。ABAQUS軟件可視化窗口不能直接生成零厚度的內(nèi)聚力單元，因此要編寫Python腳本來修改INP文件，實現(xiàn)內(nèi)聚力單元的插入[4]。在INP文件中，要重點修改三個部分：節(jié)點信息、單元信息以及單元集信息。為了在邊界上插入內(nèi)聚力單元，首先要將邊界上的節(jié)點分裂，邊界上的節(jié)點在顆?；蛘辰Y(jié)相集合內(nèi)出現(xiàn)幾次就需要分裂幾次，然后用新分裂出的節(jié)點替換邊界單元內(nèi)原來的節(jié)點。最后，利用分裂出來的節(jié)點創(chuàng)建內(nèi)聚力單元，如圖5所示。

3 鎢合金兩相結(jié)構(gòu)拉伸模擬

上文介紹了在ABAQUS軟件中利用Python腳本程序建立鎢合金的兩相微觀結(jié)構(gòu)仿真模型，該模型包括了鎢顆粒和粘結(jié)相，還考慮了顆粒界面的存在。接下來將利用該模型進(jìn)行拉伸仿真，與準(zhǔn)靜態(tài)下鎢合金的拉伸試驗曲線進(jìn)行對比，從而驗證該模型的可靠性。

3.1 鎢合金拉伸試驗

根據(jù)《金屬材料 拉伸試驗 第1部分：室溫試驗方法》（GB/T 228.1—2010），將95W-3.5Ni-1.5Fe合金制備成啞鈴形試樣，其尺寸如圖6所示，然后將試樣裝夾在萬能試驗機(jī)上進(jìn)行試驗，為了獲得準(zhǔn)確的應(yīng)變，在試樣中段使用引伸計測量其變形。試驗過程中采用位移驅(qū)動的方式進(jìn)行加載，加載速度為0.2 mm/min，當(dāng)觀測到試驗曲線過了屈服點之后，取下引伸計繼續(xù)加載，直到試樣被拉斷。之后，利用試驗機(jī)的力傳感器和引伸計，得到鎢合金的準(zhǔn)靜態(tài)拉伸應(yīng)力-應(yīng)變曲線。

3.2 鎢合金兩相結(jié)構(gòu)拉伸仿真模型

二維平面應(yīng)變的鎢合金兩相結(jié)構(gòu)拉伸模型如圖7所示，在模型的左右邊界施加一個沿[x]方向的位移（[Ux]），等于常數(shù)[C]，[y]方向的位移（[Uy]）為零。如上所述，該模型包括三部分：鎢顆粒、粘結(jié)相和界面。為了簡便起見，鎢顆粒的材料屬性選用純鎢的彈塑性材料參數(shù)，粘結(jié)相的材料屬性選用鐵的彈塑性材料參數(shù)，二者的材料參數(shù)均列于表2中，純鎢和鐵的塑性本構(gòu)模型如圖8所示[5]。界面給定雙線性內(nèi)聚力本構(gòu)模型，如圖9所示（[δ0]表示損傷起始點的位移，[δf]表示完全失效點的位移），該模型主要包含三個參數(shù)，即斷裂能量[Gc]、最大應(yīng)力[t0]和內(nèi)聚力剛度[k]，各參數(shù)取值如表3所示。需要說明的是，界面的內(nèi)聚力模型參數(shù)目前沒有直接的試驗測定方法，是結(jié)合仿真與試驗曲線試湊出來的。具體的方法是運(yùn)用控制變量法，不斷調(diào)整內(nèi)聚力模型的參數(shù)，直到仿真得到的應(yīng)力-應(yīng)變曲線與試驗接近。

3.3 內(nèi)聚力參數(shù)對鎢合金拉伸性質(zhì)的影響

在以上拉伸仿真中，試樣的等效應(yīng)力可以取右邊界上所有節(jié)點沿[x]方向的應(yīng)力的平均值，等效應(yīng)變?yōu)樵嚇拥恼w應(yīng)變。等效應(yīng)力和等效應(yīng)變的計算公式為：

式中：[σ]為等效應(yīng)力;[ε]為等效應(yīng)變;[N0]為試樣右邊界上的節(jié)點數(shù)目;[σxx]為節(jié)點沿[x]方向的應(yīng)力;[Ux]為右邊界沿[x]方向的位移;[L0]為試樣沿[x]方向的初始長度;[i]為非零自然數(shù)。

3.3.1 模型尺寸比對拉伸仿真結(jié)果的影響。模型的尺寸比（[Φ]）定義為模型整體尺寸與平均顆粒半徑的比值。對于一個非均質(zhì)的模型，其尺寸比應(yīng)該足夠大才能夠包含足夠多的微結(jié)構(gòu)，模型才具有代表性，同時其尺寸比也應(yīng)該足夠小，才能節(jié)省計算成本[6]。因此，確定鎢合金兩相結(jié)構(gòu)模型的尺寸比是有必要的。本研究保證平均顆粒尺寸不變，通過改變模型的整體尺寸來得到不同尺寸比的模型，如圖10所示，三個模型的尺寸比分別為10、15和20。

三個不同尺寸比仿真模型的拉伸應(yīng)力-應(yīng)變曲線與試驗曲線的對比結(jié)果如圖11所示。從圖中可以看出，當(dāng)尺寸比為10時，仿真的應(yīng)力-應(yīng)變曲線與試驗曲線保持一致的趨勢，但是還存在一定的差異，說明該尺寸比下的模型代表性不足;而當(dāng)尺寸比為15或20時，仿真曲線與試驗曲線基本重合，因此可以說明尺寸比為15或20的仿真模型具有足夠的代表性。為了方便研究，以下的仿真都采用尺寸比為20的模型。

3.3.2 內(nèi)聚力剛度對拉伸仿真曲線的影響。對于內(nèi)聚力單元剛度的選擇，許多學(xué)者提出了不同的方法。DAUDEVILLE等根據(jù)界面厚度和彈性模量確定了內(nèi)聚力剛度[7]。ZOU等根據(jù)經(jīng)驗，提出內(nèi)聚力剛度值為單位長度界面強(qiáng)度值的104～107倍[8]。內(nèi)聚力單元剛度應(yīng)足夠大，以提供合理的剛度，但也應(yīng)足夠小，以減少數(shù)值問題帶來的風(fēng)險。TURON等指出，內(nèi)聚力剛度在大于材料彈性模量的50倍時，數(shù)值計算的結(jié)果對于大多數(shù)問題能夠保證足夠的精確性[9]。本文中，內(nèi)聚力單元的剛度從1×105 N/mm3增加到5×107 N/mm3，而最大應(yīng)力和斷裂能量分別為520 MPa、55 N/mm，保持不變。仿真曲線與試驗曲線的對比結(jié)果如圖12（a）所示。當(dāng)內(nèi)聚力單元的剛度較?。?×105 N/mm3和1×106 N/mm3）時，內(nèi)聚力單元不能為整體的模型提供足夠的剛度，在拉伸過程中迅速發(fā)生損傷，導(dǎo)致拉伸的應(yīng)力-應(yīng)變曲線與試驗曲線有很大誤差;當(dāng)內(nèi)聚力單元剛度增加到1×107 N/mm3時，內(nèi)聚力單元已經(jīng)能夠為整體的模型提供一定剛度，仿真的應(yīng)力-應(yīng)變曲線能夠出現(xiàn)彈性-屈服-塑性階段，但與試驗曲線仍存在一定誤差，表現(xiàn)為彈性階段上升的斜率偏低，屈服后的塑性階段應(yīng)力值整體偏低;當(dāng)內(nèi)聚力單元剛度繼續(xù)增加到5×107 N/mm3時，仿真曲線與試驗曲線基本重合。因此，內(nèi)聚力單元的剛度選擇5×107 N/mm3比較合理，該值也符合TURON等給定的內(nèi)聚力剛度確定方法，即大于50倍的彈性模量[9]。

3.3.3 最大應(yīng)力對仿真應(yīng)力-應(yīng)變曲線的影響。最大應(yīng)力反映了顆粒界面的連接強(qiáng)度，已經(jīng)有研究表明，界面強(qiáng)度對鎢合金的拉伸性能有很大的影響，較低的界面強(qiáng)度導(dǎo)致鎢合金在拉伸過程中沿界面斷裂，從而表現(xiàn)出較低的屈服強(qiáng)度和極限強(qiáng)度[10]。本文中，最大應(yīng)力選擇為400 MPa、520 MPa、600 MPa，內(nèi)聚力單元剛度為5×107 N/mm3，斷裂能量為55 N/mm。仿真曲線與試驗曲線的對比結(jié)果如圖12（b）所示，當(dāng)最大應(yīng)力為520 MPa時，彈性和塑性階段模擬的應(yīng)力-應(yīng)變曲線與試驗結(jié)果吻合較好。當(dāng)最大應(yīng)力增大或減小時，仿真的應(yīng)力-應(yīng)變曲線彈性階段上升的斜率幾乎沒有變化，而屈服強(qiáng)度相應(yīng)增大或減小，這與已有的結(jié)論保持良好的一致性。

3.3.4 斷裂能量對拉伸仿真應(yīng)力-應(yīng)變曲線的影響。斷裂能量從0.1 N/mm增大到55 N/mm，內(nèi)聚力剛度為5×107 N/mm3，最大應(yīng)力為520 MPa。仿真結(jié)果如圖12（c）所示，從曲線中可以看出，當(dāng)斷裂能量足夠大時（5N/mm、20 N/mm、55 N/mm），該值對仿真的應(yīng)力-應(yīng)變曲線影響很小，幾乎不影響彈性階段上升的斜率，稍微影響屈服后的塑性階段，斷裂能量越小，塑性階段下降越快。而當(dāng)斷裂能量減小到一定程度（0.1 N/mm）時，仿真的應(yīng)力-應(yīng)變曲線發(fā)生很大的變化，在經(jīng)歷了一定的上升階段后，曲線迅速下降，這是由于給定的斷裂能量太小導(dǎo)致內(nèi)聚力單元在拉伸過程中完全失效，從而使拉伸應(yīng)力-應(yīng)變曲線出現(xiàn)驟降現(xiàn)象。

4 結(jié)論

根據(jù)95W-3.5Ni-1.5Fe合金粉末冶金的加工工藝，本文提出了一種基于顆粒擠壓的兩相結(jié)構(gòu)建模方法。首先根據(jù)鎢合金的實際顆粒分布特征，利用RSA算法在MATLAB軟件中生成了初始顆粒構(gòu)型，然后將所有顆粒信息輸出到ABAQUS軟件中建立一個顆粒擠壓模型。隨著擠壓過程的不斷進(jìn)行，模型模擬的顆粒體積比逐漸上升，當(dāng)其達(dá)到真實的顆粒體積比時，將所有顆粒輪廓輸出，并編寫Python腳本程序?qū)㈩w粒輪廓輸入到ABAQUS軟件中完成輪廓重建，隨后運(yùn)用布爾運(yùn)算生成粘結(jié)相，最后在顆粒與粘結(jié)相相交的區(qū)域插入內(nèi)聚力單元，以模擬顆粒界面的作用。上述模型被應(yīng)用于拉伸仿真，本研究結(jié)合準(zhǔn)靜態(tài)下95W-3.5Ni-1.5Fe合金的拉伸曲線確定了內(nèi)聚力單元的材料參數(shù)，最后分析了內(nèi)聚力單元參數(shù)對拉伸曲線的影響。

研究表明，本文基于顆粒擠壓的方法建立的95W-3.5Ni-1.5Fe合金兩相細(xì)觀結(jié)構(gòu)仿真模型是可行的，該模型考慮了95W-3.5Ni-1.5Fe合金中鎢顆粒、粘結(jié)相以及顆粒界面，近似還原了95W-3.5Ni-1.5Fe合金的真實結(jié)構(gòu)。根據(jù)仿真與試驗，本文得到了95W-3.5Ni-1.5Fe合金顆粒界面的內(nèi)聚力單元參數(shù)。其中，內(nèi)聚力單元的剛度主要影響拉伸應(yīng)力-應(yīng)變曲線彈性段的斜率，單元剛度越低，曲線彈性段的斜率就越低;最大應(yīng)力主要影響拉伸應(yīng)力-應(yīng)變曲線的屈服點，對彈性段的斜率幾乎沒有影響;斷裂能量在一定范圍內(nèi)對拉伸-應(yīng)力-應(yīng)變曲線的影響很小，且主要影響塑性階段。

參考文獻(xiàn)：

[1]羅榮梅.一種95W細(xì)晶鎢合金動態(tài)力學(xué)性能[J].遼寧化工，2015（7）：776-778.

[2]KIRAN U R，KHAPLE S K，SANKARANARAYANA M，et al.Effect of swaging and aging heat treatment on microstructure and mechanical properties of tungsten heavy alloy[J].Materials Today：Proceedings，2018（2）：3914-3918.

[3]SATYANARAYANA P V，BLESSTO B，SOKKALINGAM R，et al.Effect of Fe Addition to Binder Phase on Mechanical Properties of Tungsten Heavy Alloy[J].Transactions of the Indian Institute of Metals，2019（4）：1-9.

[4]王鷹宇.Abaqus分析用戶手冊[M].北京：機(jī)械工業(yè)出版社，2018：60-62.

[5]劉筱玲.鎢合金材料宏微觀性能分析及動態(tài)性能測試研究[D].成都：西南交通大學(xué)， 2008：20-22.

[6]SAVVAS D，STEFANOU G，PAPADRAKAKIS M.Determination of RVE size for random composites with local volume fraction variation[J].Computer Methods in Applied Mechanics and Engineering，2016（305）：340-358.

[7]DAUDEVILLE L，ALLIX O，LADEVEZE P.Delamination analysis by damage mechanics：Some applications[J].Composites Engineering，1995（1）：17-24.

[8]ZOU Z，REID S R，LI S.Modelling Interlaminar and Intralaminar Damage in Filament-Wound Pipes under Quasi-Static Indentation[J].Journal of Composite Materials，2002（4）：477-499.

[9]TURON A，DAVILA C G，CAMANHO P P.An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models[J].Engineering Fracture Mechanics，2006（10）：1665-1682.

[10]SENTHILNATHAN N，ANNAMALAI A R，VENKATACHALAM G.Microstructure and mechanical properties of spark plasma sintered tungsten heavy alloys [J].Materials Science & Engineering A，2018（710）：66-73.