基于热-弹塑性-蠕变-相对密度耦合本构模型的热等静压数值模拟

李勇进　马迎松　屈浩宇　樊力中

摘 要：粉末冶金热等静压近终（净）成形（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