基于內聚力界面的鎢合金兩相細觀結構拉伸模擬

楊梓清

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

關鍵詞：鎢合金;細觀結構;擠壓模型;內聚力單元;有限元法;拉伸仿真

中圖分類號：TB33;TB302.3文獻標識碼：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

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

1 鎢合金微觀結構

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

2 鎢合金微觀結構模型建立過程

鎢合金兩相微觀結構模型的建立過程主要分三個步驟，分別是顆粒擠壓模型的建立、輪廓重建與幾何修復、內聚力單元的插入。

2.1 顆粒擠壓模型

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

隨著上邊界的連續向下移動，顆粒逐漸被擠壓變形，顆粒的體積分數逐漸上升。當所有顆粒的體積分數達到鎢合金中鎢顆粒真實的體積分數時，提取所有顆粒的輪廓，輸出每個顆粒邊界上所有的節點坐標。為了保證邊界上的節點按照順序逐個輸出，在擠壓模型運算之前，要將每個顆粒邊界上所有的節點從小到大重新編號。

2.2 輪廓重建與幾何修復

完成擠壓模型后，依次輸出變形后所有顆粒邊界上的節點坐標，再編寫Python腳本程序將其按順序輸入ABAQUS軟件中，以完成建模。輪廓重建的模型如圖3（a）所示。隨后，創建一個可以包含所有變形顆粒的矩形實體，并對所有顆粒執行布爾運算，以獲得粘結相的模型，如圖3（b）所示。至此，鎢合金的兩相模型已初步建立。但是，該方法中顆粒的邊界是由輸入的節點依次連接形成的，由于計算機本身的數值誤差，相鄰顆粒之間可能會存在微小的間隙，這不利于后續的網格劃分，并且會嚴重影響網格的質量。因此，在劃分網格之前需要進行幾何修復。

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

2.3 內聚力單元的插入

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

3 鎢合金兩相結構拉伸模擬

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

3.1 鎢合金拉伸試驗

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

3.2 鎢合金兩相結構拉伸仿真模型

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

3.3 內聚力參數對鎢合金拉伸性質的影響

在以上拉伸仿真中，試樣的等效應力可以取右邊界上所有節點沿[x]方向的應力的平均值，等效應變為試樣的整體應變。等效應力和等效應變的計算公式為：

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

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

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

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

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

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

4 結論

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

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

參考文獻：

[1]羅榮梅.一種95W細晶鎢合金動態力學性能[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].北京：機械工業出版社，2018：60-62.

[5]劉筱玲.鎢合金材料宏微觀性能分析及動態性能測試研究[D].成都：西南交通大學， 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.