基于熱-彈塑性-蠕變-相對密度耦合本構模型的熱等靜壓數值模擬

李勇進　馬迎松　屈浩宇　樊力中

摘 要：粉末冶金熱等靜壓近終（凈）成形（powder metallurgy hot isostatic pressing to near-not-shape， PM-HIP-NNS）技術用于制備高性能的復雜零部件，在粉末冶金行業具有廣闊的應用前景，其致密化過程有特殊的力學行為，需要與數值模擬相結合，以便減少實驗過程中試錯法造成的物力、財力浪費。本文基于熱-彈塑性-蠕變-相對密度耦合本構模型，對SS316L粉體材料在SS304鋼包套材料中的壓制過程進行了有限元數值模擬，探究熱等靜壓（hot isostatic pressing， HIP）過程中粉末密度場、溫度場以及應力場的變化。結果表明：高溫高壓下金屬粉末成形規律的數值模擬與試驗相吻合，證實模型可以合理預測金屬粉末在高溫高壓下的變形行為，為粉末HIP工藝優化提供了重要依據。

關鍵詞：有限元；熱等靜壓；相對密度；本構模型；包套

中圖分類號：TK123；TG306 DOI：10.16375/j.cnki.cn45-1395/t.2024.02.017

0 引言

粉末冶金熱等靜壓近終（凈）成形（powder metallurgy hot isostatic pressing to near-not-shape， PM-HIP-NNS）技術是粉末冶金領域的高端技術，在歐美等發達國家已經廣泛應用于核電、海工、航空航天等重要領域。國內在此技術領域的發展遲滯于歐美國家，主要原因是，熱等靜壓近終（凈）成形技術需要與數值模擬相結合，以便減少實驗過程中大量物力、財力的浪費，國內學者在這方面的技術研究尚未完全成熟[1]。熱等靜壓（hot isostatic pressing， HIP）粉末冶金致密化過程是一個熱力耦合過程，密度和體積均發生變化，對其致密化機理至今沒有統一的認識。熱等靜壓粉末冶金致密化過程的本構模型對于模擬結果有重要的影響，因此，此項研究工作具有非常重要的理論和實際意義。

有限元模擬仿真[2]的準確性很大程度上取決于本構模型的選擇。Teraoku[3]將研究的重點放在粉末體材料屈服準則的修正上，基于不同致密度零件的壓縮實驗對經典Shima屈服準則進行修正，獲得能準確描述Ti6Al4V合金粉末屈服行為的數值模型，并采用相關的實驗結果驗證了所建立模型的可靠性。Gillia等[4]的研究重點是將粉末熱等靜壓過程中粉末顆粒的應變強化引入屈服準則的修正中。劉國承[5]針對金屬粉末熱等靜壓致密化的數值模擬理論、致密化機制以及熱等靜壓成形件的力學性能演變等問題開展了大量的研究，建立了Inconel625粉末熱等靜壓成形的高溫本構模型。由于Inconel 718合金的顯微組織和綜合力學性能對熱響應敏感，通過掌握Inconel 718合金中相的溶解和析出規律及組織與性能之間的關系，設計不同的成形工藝，滿足不同服役環境對材料的綜合性能需求[6-8]。與致密化金屬成形相比，金屬粉末致密化過程有其特殊的力學行為[9]，粉末體發生密度和體積的變化，且經歷塑性變形、擴散及蠕變等過程，常規的Mises屈服準則不能滿足要求。為了掌握高溫高壓下金屬粉末的成形規律，明確其致密化機理，達到工業生產上節約成本、提高生產效率的目的，選取合適的本構模型，合理設計邊界條件，通過與實驗結果對比，修正有限元模型，對于預測金屬粉末在高溫高壓下的變形行為至關重要。

本文基于Van Nguyen等[10-11]提出的熱、蠕變、塑性和密度變化耦合的本構模型，采用Peirce等[12]的速率切線方法，對蠕變引起的速率依賴的黏塑性進行求導，并與單獨密度演化的速率依賴的黏塑性相耦合，建立了一個熱-彈塑性-蠕變-相對密度耦合本構模型，并使用FORTRAN語言對本構模型進行編譯，在ABAQUS大型通用有限元軟件中實現為用戶自定義材料子程序（UMAT），以SS316L粉體材料和SS304鋼包套材料為研究對象進行熱等靜壓工藝的有限元模擬，探究HIP過程中粉末密度場、溫度場以及應力場[13]的變化。

1 材料本構模型

2 材料參數

本文的數值模擬所用的材料是SS316L粉體材料和SS304鋼包套材料。這些材料參數提取的溫度范圍為（20， 1 125）℃[10-11]。SS316L粉體材料的初始相對密度[ρ0=0.75]。表1和表2分別給出了SS316L粉體材料和SS304鋼包套材料的比熱容系數。表3列出了這2種材料的熱膨脹系數和導熱系數。所有溫度的測量都是以攝氏溫度為單位。2種材料的密度系數、彈性模量系數和泊松比系數見表4。表5列出了每個材料的流動應力特性。表6列出了確定蠕變等效應力函數的敏感性系數。表7列出了每種材料的蠕變模型的系數。

3 致密化的有限元模型

使用FORTRAN語言對本構方程進行編譯鏈接生成自定義的材料子程序（UMAT），使用簡單模型對致密化程序進行驗證，采用ABAQUS軟件作為有限元分析模擬平臺，選用合適形狀的工業構件進行建模分析。構件模型有粉末和包套兩部分，為圓柱形構件，內部為粉末，外部為包套。為方便觀察到構件的變化，選取圓柱形構件的1/4部分進行建模，在完成建模后，采取共節點法進行網格劃分，模型構件形狀如圖1所示，包套厚度3 mm，粉末體為30 mm半徑和70 mm高的1/4圓柱，包套兩端高出包套底面12 mm。

在施加邊界條件時，由于施加的是隨時間線性變化的壓力和隨時間線性變化的溫度，所以載荷會以表格的形式進行設定，并施加在包套的外表面。模型是1/4部分的圓柱體，所以模型會向著中心線壓縮，而中心線不會在XZ平面移動，只會在豎直Y方向上移動。為了便于觀察，使模型底部的單元節點固定不動，在施加邊界條件時會對中心線上的節點賦予X、Z方向上的固定約束，對X截面上所有節點賦予X方向上的固定約束，對Z截面上所有節點賦予Z方向上的固定約束，之后在包套底部圓周線上的所有節點施加Y方向上的固定約束。在發生HIP過程之前模型存在自身的溫度，所以需要在施加載荷之前設置初始溫度為25 ℃。設置分析類型為瞬態分析，載荷大變形，設置時間步長。

本文所使用的HIP工藝為：在12 000 s內，溫度從25 ℃線性上升到1 125 ℃，壓力線性上升到130 MPa，將溫度和壓力保持2 h，然后在8 h內降至200 ℃和0。

4 對致密化模型進行分析

4.1 密度場

在經過HIP工藝工況施加、對其進行求解完成后，可以得到不同時刻構件的相對密度分布，選取4個比較典型時刻的相對密度分布（見圖2）：第44增量步，即從升溫升壓到保溫保壓的過渡時刻；第68增量步，此時為保溫保壓階段的中間時刻；第73增量步，這個時刻正是保溫保壓到降溫降壓的過渡時刻；第146增量步，即為整個HIP過程的完成時刻。

由圖2可知，包套里的粉末體在HIP過程中存在著顯著的相對密度差異，粉末體為一圓柱體。在HIP整個過程中，粉末圓柱體上下兩端的相對密度值與粉末圓柱體中間部分的相對密度值的高低呈相對變化，在HIP初期，粉末圓柱體的上底面、下底面以及圓柱側面受到包套傳遞的應力作用，粉末圓柱體的上下底面以及圓柱側面的相對密度率先發生變化。由于圓柱上下底面為平面，圓柱體側面為曲面且結構為規則的圓結構，粉末圓柱體受到的應力由包套傳遞到內部，在傳遞過程中，包套傳遞到粉末圓柱體上下底面處的應力損耗比傳遞到粉末圓柱體側面的應力損耗要小得多，這導致HIP初期粉末體上下底面的相對密度比粉末圓柱體側面的相對密度變化明顯的多，且這種狀態會持續一段時間。

在HIP進行一段時間后，應力不斷被傳遞進入粉末圓柱體內部，圓柱體中心處的密度開始變大。由于圓柱體側面的受力面積比上下底面的受力面積大得多，側面隨著上下底面間距的減小，圓柱粉末體受到的軸向力逐漸低于徑向力，導致粉末圓柱體上下底面的相對密度增長量逐漸小于粉末圓柱體側面的相對密度增長量。在圓柱體上下底面的相對密度值增大到一定數值后，粉末圓柱體側面的相對密度值在一定時間后便與上下底面相對密度值一致。經過一段時間后，圓柱粉末體受到的軸向力低于徑向力，此時粉末圓柱體側面的相對密度值大于圓柱體上下底面處的相對密度值。由于應力和溫度不斷傳遞進入粉末圓柱體內部，粉末體的間隙不斷被擠壓，粉末進行流動和顆粒重排，并主要進行蠕變機制，此時粉末圓柱體的內部的相對密度逐漸由外至內地呈階梯狀有序增大，直至粉末圓柱體內部的階梯形相對密度梯度消失，轉變為內部相對密度均勻，圓柱體的相對密度呈階梯狀從大到小由中間處向上下底面方向傳遞，此時粉末圓柱體的上下底面處的相對密度值小于粉末體中心處的相對密度值。粉末圓柱體上下底面與其側面共線處位置的相對密度值在整個HIP過程中，除了HIP初期粉末體存在壓力屏蔽效應導致它的值大于粉末體中心處的相對密度值外，在之后的時間里，它的相對密度值一直是最低的。

4.2 溫度場

在HIP過程中，構件內部的溫度場對粉末HIP致密化起著重要的作用。圖3為4個典型時刻的溫度云圖，圖3（a）為處于升溫升壓到加溫加壓過渡時刻的溫度云圖，圖3（b）為即將保溫保壓完成結束時刻的溫度云圖，圖3（c）為卸載階段中間時刻的溫度云圖，圖3（d）為HIP過程的完成時刻。從圖3中可知，整個HIP過程中，構件的溫度場分布都是不均勻的，存在較大的溫度梯度。在升溫到即將進行保溫保壓的時刻，內外溫差約2 ℃。經過長時間的保溫保壓，構件內部溫差為0，長時間的保溫保壓使得構件內部的溫度達到平衡。卸載后，隨著相對密度的提高，粉末體內部的空隙總體積減少，阻礙傳熱的介質變少了，構件內外還會存在溫差，只是溫差會變小。載荷卸載到一半時，內外溫差約為1 ℃，到了完全卸載完成時，溫差約為1 ℃。

內外存在明顯區別的溫度梯度，在升溫階段，包套溫度比粉末的溫度高，導致粉末體的外層溫度比內部溫度高，外部粉末體的致密化速度要比粉末體內部的致密化速度快，使得粉末體形成一個另類的空腔包套，所以在對構件施加壓力時，粉末體的外部會承受更大的壓力，影響了粉末體內部的受力狀態。然后隨著升溫升壓的進行，粉末體吸收更多的熱量，溫度梯度不斷拉大，相對密度差異不斷變大。直到保溫保壓的進行，粉末的蠕變和擴散效應變得更強，粉末體的內外溫度逐漸均勻化，溫度梯度減小，相對密度也逐漸變得均勻。到了卸載階段，構件的外部溫度逐漸降低，形成內外溫差，引起構件收縮的熱應變，粉末體得到更大的熱應力，在此時相對較高的溫度壓力作用下，促進了構件相對密度進一步提高。

4.3 Mises等效應力場

圖4是HIP完成后構件的Mises等效應力云圖。由圖4可知，粉末體外部的等效應力值比粉末體內部的等效應力值高，且等效應力值從粉末體內部到粉末體外部由小到大呈階梯狀分層次排列。由于粉末體受到來自包套側面的徑向擠壓作用，以及包套上下底面對粉末體的軸向擠壓作用，使得粉末圓柱體的上下底面與側面相接觸的直角區域的等效應力值為粉末體的等效應力值分布的最大處，沿粉末圓柱體的4個直角的對角線會產生較大的剪切應力。所以，在粉末圓柱體中，沿著包套拐角對角方向的Mises等效應力最大，這種變化很可能對構件的致密化產生不利影響。

由圖4的包套外側與內側圖可知，在包套上下底面與包套體側面的相接處表現出高應力帶，如果發生開裂，HIP工藝可能會發生失敗。在上底面與側面相接觸的直角位置處，包套受力變形會發生大的軸向位移，使得粉末體被壓縮，包套對粉末體施加了一個軸向力，粉末體會對包套產生一個軸向力抵抗力，所以在包套上底面的直角部位會產生應力集中。同時，包套在變形過程中，包套側面也會對粉末體施加壓力，使得粉末上底面處產生了抵抗側向包套的徑向力，在這里同時產生了軸向力與徑向力，使得此部位有著很大的應力抵抗，集中形成了高應力帶。下底面與側面相接觸的直角位置處的帶狀區域與上底面與側面相接觸的直角位置處的受力狀態及引發原因都基本相同，也形成了高應力區，這2個位置是包套容易破裂的地方。包套的側面也是同時承受著軸向力與徑向力，在這里，包套受到的作用力非常大，形成了一條高壓力帶，也是包套容易破裂的地方。通過對粉末體與包套的Mises等效應力分析，可以更合理地優化構件和包套的結構設計，為HIP工藝優化提供重要的依據。

為研究不同位置等效應力的變化規律，選取幾個關鍵節點進行研究（見圖5）。由圖6可知，幾個節點在保溫保壓階段的等效應力最小，在升溫升壓階段和卸載階段，等效應力開始產生，卸載階段的等效應力最大值比升溫升壓階段的等效應力最大值要大得多。在升溫升壓階段，應力先是快速增大到一定數值，然后又快速減小為幾乎是整個HIP過程的最小值。而在進入卸載階段后，隨著溫度的降低，等效應力快速增大，產生遠遠超出升溫升壓階段的殘余應力。綜合來看，卸載階段粉末體內部產生高應力的原因是溫度下降引起熱應力，同時也是粉末體在此階段致密化的主要動力。

5 結論

本文基于熱-彈塑性-蠕變-相對密度耦合本構模型，對SS316L粉體材料圓柱試樣在SS304鋼包套材料中的壓制過程進行了有限元數值模擬，HIP過程中粉末密度場、溫度場以及應力場的變化主要表現為：

1）自HIP過程開始，粉末體上下底面位置處的密度先發生變化，在HIP進行一段時間后，粉末體側面位置處的密度變得越來越大，超越了粉末體上下底面處的密度值。粉末體中心的密度也逐漸變大，直至比粉末體上下底面位置的密度值大，整個HIP過程中粉末體的密度值大小以粉末體中心處為中心呈階梯狀分層次分布，HIP初期密度值自粉末體中心處從小到大向四周排列分布，隨著HIP的進行，密度值又自粉末體中心處從大到小分層次向兩底面排列分布。

2）整個HIP過程中，粉末體存在著較大的溫度梯度。從HIP進行到保溫階段結束，粉末體由內到外、溫度由小到大呈階梯狀分層次分布。到了卸載階段，粉末體由內到外、溫度由大到小呈階梯狀分層次分布。

3）粉末體外部的等效柯西應力值比粉末體內部的等效柯西值應力高，且等效柯西應力值從粉末體內部到粉末體外部由小到大呈階梯狀分層次排列。粉末體的上下底面與側面相接觸的直角區域處的等效柯西應力值為粉末體的等效柯西應力最大值，沿粉末圓柱體的直角的對角線周圍有著較大的剪切應力。

4）包套的上下底面與側面相接觸的直角位置處和側面的中心位置處形成了高應力區，形成了高壓力帶，這是包套容易破裂的地方。

本研究為HIP工藝能更合理地優化構件和包套的結構設計奠定基礎。

參考文獻

[1] 馬宏昌，何力軍，陳磊，等.鈹粉熱等靜壓成型的數值模擬[J].寧夏大學學報（自然科學版），2020，41（1）：45-49.

[2] 李勇進，馬迎松，屈浩宇，等.基于ANSYS的鋁合金球淬火熱分析[J].廣西科技大學學報，2023，34（1）：105-111.

[3] TERAOKU T.Hot isostatic pressing simulation for titanium alloys[J].International Journal of Powder Metallurgy，2008，44（5）：57-61.

[4] GILLIA O，BOIREAU B，BOUDOT C，et al.Modelling and computer simulation for the manufacture by powder HIPing of blanket shield components for ITER[J].Fusion Engineering and Design，2007，82（15）：2001-2007.

[5] 劉國承.金屬粉末熱等靜壓致密化數值模擬與試驗研究[D].武漢：華中科技大學，2011.

[6] 徐文才，郎利輝，黃西娜，等.鋁合金復雜薄壁件熱等靜壓成形數值模擬[J].鍛壓技術，2019，44（6）：65-72.

[7] 成雅徽. GH4169合金粉末選區激光熔化成形數值模擬及試驗研究[D].太原：中北大學，2016.

[8] 姚草根，孟爍，李秀林，等.粉末氧含量對熱等靜壓FGH4169合金力學性能與組織的影響[J].粉末冶金材料科學與工程，2017，22（1）：33-40.

[9] 曾旺，顧連勝，陳松靖，等.纖維對超高性能混凝土力學性能及微觀結構影響研究[J].廣西科技大學學報，2023，34（1）：20-27.

[10] VAN NGUYEN C，BEZOLD A，BROECKMANN C.Anisotropic shrinkage during hip of encapsulated powder[J].Journal of Materials Processing Technology，2015，226：134-145.

[11] VAN NGUYEN C，DENG Y，BEZOLD A，et al.A combined model to simulate the powder densification and shape changes during hot isostatic pressing[J].Computer Methods in Applied Mechanics and Engineering，2017，315：302-315.

[12] PEIRCE D，SHIH C F，NEEDLEMAN A.A tangent modulus method for rate dependent solids[J].Computers and Structures，1984，18（5）：875-887.

[13] 郝子龍，石光林.基于ANSYS的平板堆焊凝固過程的應力分析[J].廣西科技大學學報，2016，27（1）：54-57.

[14] KUHN H，DOWNEY C.Deformation characteristics and plasticity theory of sintered powder metals[J].International Journal of Powder Metallurgy，1971，7（1）：15-25.

[15] ABOUAF M，CHENOT J L，RAISSON G，et al.Finite element simulation of hot isostatic pressing of metal powders[J].International Journal for Numerical Methods in Engineering，1988，25（1）：191-212.

[16] KOHAR C P，MARTIN E，CONNOLLYD S，et al.A new and efficient thermo-elasto-viscoplastic numerical implementation for implicit finite element simulations of powder metals：an application to hot isostatic pressing[J].International Journal of Mechanical Sciences，2019，155：222-234.

Hot isostatic pressure numerical simulation based on thermo-elastoplastic-creep-relative density coupled constitutive model

LI Yongjin， MA Yingsong*， QU Haoyu， FAN Lizhong

（School of Mechanical and Automotive Engineering， Guangxi University of Science and Technology， Liuzhou 545616， China）

Abstract： Powder metallurgy hot lsostatic pressing to near-not-shape （PM-HIP-NNS） technology is used to prepare complex components with high performance， and has a good application prospect in powder metallurgy industry. The densification process has special mechanical behavior， which needs to be combined with numerical simulation to reduce the waste of material and financial resources caused by trial and error in the experimental process.? Based on the thermo-elastoplastic-creep-relative density coupling constitutive model， the finite element numerical simulation of the compression process of SS316L powder material in SS304 ladle capsule material was carried out in this paper， and the changes of powder density field， temperature field and stress field during HIP process were investigated.? The results show that the numerical simulation of metal powder forming law under high temperature and high pressure is consistent with the experiment， which confirms that the model can reasonably predict the deformation behavior of metal powder under high temperature and high pressure， and provides an important basis for the optimization of powder HIP process.

Keywords： finite element; hot isostatic pressure; relative density; constitutive model; capsule

（責任編輯：黎 婭）

收稿日期：2023-03-14；修回日期：2023-06-06

基金項目：廣西高校中青年教師科研基礎能力提升項目（2019KY0379）；廣西科技基地和人才專項項目（桂科AD19110111）資助

第一作者：李勇進，在讀碩士研究生

*通信作者：馬迎松，博士，講師，研究方向：材料與結構損傷破壞的宏細管研究，E-mail：mys@gxust.edu.cn