18MnNiMo钢核电厚壁锻件淬火冷却温度场数值模拟研究

杜东旭 任利国 赵德利 李家驹

(天津重型装备工程研究有限公司，天津300457)

18MnNiMo钢核电厚壁锻件淬火冷却温度场数值模拟研究

杜东旭 任利国 赵德利 李家驹

(天津重型装备工程研究有限公司，天津300457)

使用有限元模拟软件对300 mm壁厚锻件冷却过程进行数值模拟。通过与实测冷却速率对比，优化了数值模拟精度。

18MnNiMo；锻件；淬火；温度场；数值模拟

核电是一种安全、清洁、经济、高效的能源，积极发展核电是我国解决能源供需矛盾的主要选择之一。淬火是核电大型锻件生产过程中的关键一步。随着核电设计的日益创新，核电装机容量不断扩大，锻件壁厚也随之增加。锻件厚度的增加必然导致淬火时冷却能力下降，进而导致产品性能异常。核电大型锻件在强度、韧性方面有着极高的要求，因此预测锻件淬火冷却时的温度场对核电产品的生产具有重要指导意义。尽管目前的淬火过程计算机模拟计算还不是成熟的技术，但是它与必要的试验测试相结合，可以成为研究淬火方法和制订工艺参数的辅助决策工具。本文采用ABAQUS有限元软件，对18MnNiMo钢材质的核电蒸发器328 mm厚测温件进行淬火冷却温度场模拟。通过对比实测温度曲线，优化提升模拟精度，得到近似准确的淬火冷却温度场。

1 大锻件冷却过程有限元模拟基础

1.1 导热方程

锻件在冷却时的传热现象可以用Fourier导热微分方程和反映锻件与环境交互关系的边界条件描述，若大锻件用钢的物性参数ρ、cp和λ在x、y、z三个方向上各向同性，则直角坐标系下的导热方程为：

式中，λ为材料的热导率，单位为W/(m·K)；q′为内热源的热流密度，单位为W/m2；ρ为材料的密度，单位为kg/m3；cp为材料的定压比热，单位为J/(kg·K)；T为温度，单位为K；t为时间，单位为s。

相应的初始条件为：

式中，Tinitial为锻件的初始温度。

大锻件冷却过程最常用的是对流和辐射换热边界，可用第三类边界条件表示：

式中，H为换热系数，单位为W/(m2·K)；Γ为换热边界；∂T/∂n为延向的温度梯度；Ts为零件表面温度，单位为K；Tf为介质温度，单位为K。

1.2 大锻件淬火时温度场、组织场、应力场三场耦合模拟

热处理过程是温度、组织转变、应力三方面相互作用的复杂过程，其关系如图1所示。①表示温度对应力影响，热处理时锻件内外冷速不均匀导致产生热应力。②表示温度对组织转变的作用。③表示组织转变对应力的影响。④表示组织转变对温度场的影响。主要是组织转变时产生潜热，影响温度场分布。⑤表示应力对组织转变的影响，主要是指应力作用改变等温转变的初始和终了时间，影响组织转变。⑥表示应力对温度的影响。主要是指大应力作用使产品发生变形，变形功转化为热能进而影响温度场。

图1 温度场、组织场、应力场的三场耦合关系Figure 1 Three coupling relations of temperature field, microstructure field and stress field

大型核电锻件淬火时变形量基本可以忽略不计，因此，应力场可以近似看做对温度场没有影响。从三场的耦合关系可以看出，温度场的模拟是三场耦合模拟的基础。因此，我们应选取正确的热物参数、相变潜热参数、换热系数，以保证模拟结果的准确性。

2 有限元模型的建立

运用ABAQUS有限元软件对材质为18MnNiMo钢的328 mm厚测温件的淬火过程进行模拟。测温件实体尺寸及热电偶位置如图2所示。其中A、B、C、D四处测温点距离其最近表面分别为23 mm、 88 mm、155 mm、75 mm。淬火设备为5 m×3 m×5 m浸水淬火槽。模拟淬火起始

温度为890℃，模拟淬火时间为2 h。

根据328 mm厚测温件的实际尺寸建立了有限元几何模型，建模时选用了精度较高的六面体网格。为节省计算时间，并且尽量在软件允许的条件下细化网格，所以在保留完整厚度方向(328 mm)的前提下，取测温件的1/4部分做建模处理。模型如图3所示。最终，模型共有21 252个单元，92 371个节点，单元类型为DC3D20。

模拟程序采用本构模型UMAT编辑328 mm测温件的材料18MnNiMo钢的热物性能参数。具体的热物参数如表1所示。热物参数随温度变化而变化，在温度场模拟计算时，ABAQUS软件对两个温度点间的参数自动进行线性插值拟合。

3 模拟结果与分析讨论

3.1 恒定换热系数无相变潜热条件下的模拟结果

通常认为，界面换热要经历三个机理不同的阶段[1]：膜沸腾、核沸腾及对流换热阶段，相应的换热系数取值范围也有所不同。例如，水的膜态沸腾阶段为450℃以上，这个阶段由于在工件表面形成一层过热的蒸汽膜，将工件表面和液体隔开，而蒸汽的导热性很差，工件的冷却主要靠辐射和蒸汽导热来实现，所以此时换热系数较小。当蒸汽膜破裂以后，淬火介质与工件直接接触，不断产生强烈的沸腾，由于需要大量汽化潜热，冷却速度达到最快，此时水的换热系数最大。当淬火工件表面的温度降至淬火介质的沸点温度以下，沸腾停止，开始对流冷却阶段，此时冷却速度又变小。因此，换热系数通常情况下不是一个定值。

图2 328 mm厚测温件尺寸和敷偶位置图Figure 2 Dimension of workpiece with thickness 328 mm used for temperature test and the thermocouples location map

图3 328 mm测温件有限元模型Figure 3 Finite element model of workpiece with thickness 328 mm used for temperature test

温度K弹性模量GPa泊松比导热系数W/(m·K)比热容J/(kg·K)密度kg/m3293373473573673773873973112311731223209.8203.8198.8191.4183.3175164151-1221170.310.30.310.310.320.320.320.31-0.30.341.243.243.242.54240.438.131.929.728.329.6454478517574647700786955---7922----------

但是，大锻件淬火时冷却介质都是不断搅动的，通过鼓风、循环、喷淬的方法破坏蒸汽膜，缩短甚至消除膜沸腾阶段。因此，如果破坏效果好，可以近似认为换热系数在绝大多数时间内是稳定状态下的定值。ABAQUS软件中反映淬火过程换热的参数是锻件的表面换热系数。水的稳定换热系数约为2 324 W/(m2·K)，这个系数多次应用于有限元淬火模拟过程[2]。我们首先选用这个数值作为328 mm厚测温件的表面换热系数进行淬火过程温度场模拟。

淬火结束时的温度场模拟结果如图4所示。模拟结果显示，经2 h淬火后，锻件中心部位温度在40℃左右，表面温度在10℃左右，此结果与实际淬火的最终结果相符。但是，锻件实际的温度场是时间的函数，最终结果的正确并不能说明模拟结果的正确性。因此需要对比模拟冷却温度曲线和实测冷却温度曲线，我们以B测温点的模拟结果为例。

B测温点的模拟温度曲线同实测温度曲线的对比结果如图5所示。可以看到，表面换热系数取2 324 W/(m2·K)所得到的模拟结果同实测结果差异巨大。高温段840～400℃的B点实测平均冷速为8.2℃/min，而模拟出的高温段B点冷速高达586℃/min。究其原因，从B点实测冷却曲线可以看出，5 m×3 m×5 m淬火水槽并没有起到缩短消除膜沸腾阶段的作用。高温起始阶段即膜沸腾阶段测温件冷却速度极慢。所以，用换热系数取定值来模拟大型锻件淬火过程是不合理的，根据文献取的表面换热系数2 324 W/(m2·K)明显偏大[2]。另外，本次模拟没有考虑到相变潜热对冷却温度场的影响，这也是造成这次模拟准确度较低的原因。

图4 328 mm厚测温件温度场模拟结果Figure 4 Temperature field simulation results of workpiece with thickness 328 mm used for temperature test

图5 B点模拟冷却曲线Figure 5 Simulated cooling curve of B point

3.2 考虑相变潜热后的模拟结果

热处理过程中发生组织转变时会吸收或释放潜热。固态组织转变的潜热虽不像熔化或凝固时那么大，但也是不可忽略的一个因素。ABAQUS模拟软件中没有相变潜热的概念，因此需要自行定义。模拟计算中，处理潜热通常有三种方法[2]，即等效热量法(温度回升法)、等效热熔法和比焓法。以等效热熔法为例，即将相变潜热对温度场的影响换算成等效的热熔，通过热熔的变化反应潜热的变化。通过公式(1)可以得到等效热熔的公式：

式中，ΔV为Δt时间内组织变化的增量；ΔT为Δt时间内温度变化的增量；L为相变潜热。

但是，组织的变化量的模拟计算过程极其复杂。对于扩散型相变，按照连续转变曲线即CCT曲线进行组织模拟，首先需要根据Scheil叠加原理(孕育期叠加原理)[3]将CCT连续冷却过程转变成TTT等温转变过程进行处理。根据孕育期叠加判据来定义有限元模型中各个结点或单元是否发生铁素体、珠光体和贝氏体相变。之后，根据扩散型等温转变公式[4]求得铁素体、珠光体和贝氏体的转变量， 如公式(5)所示。对于非扩散型的马氏体相变，根据K-M公式[5]求得转变量，如公式(6)所示。

V=1-exp(-btn)

(5)

V=1-exp[-α(Ms-T)]

(6)

式中，V为转变量；n、b、α为根据TTT图计算出的常数；t为转变时间；Ms为马氏体转变点。

根据大型核电锻件实际生产经验，由于核电材料18MnNiMo淬透性不佳[6]，厚壁锻件各处淬火冷却速度主要集中在(10～40)℃/min这个范围内。因此，我们可以利用锻件的这一特点，简化组织转变判定程序和组织转变量模拟计算程序。根据图6所示的18MnNiMo的CCT曲线图，可以看到，大型核电锻件在冷速(10～40)℃/min这个范围内主要发生铁素体、珠光体转变和贝氏体相变。所以我们分别设置两种极限情况的模拟程序，即只发生贝氏体转变的情况和发生珠光体+铁素体+贝氏体转变的情况。文献认为发生贝氏体相变的相变潜热约为(56 000～79 500)J/kg，珠光体和铁素体相变的相变潜热约为77 000J/kg和75 000J/kg[7-8]。我们分别选用56 000J/kg和79 500J/kg作为贝氏体相变潜热进行两次只发生贝氏体转变的模拟，然后选用79 500J/kg作为贝氏体相变潜热模拟珠光体+铁素体+贝氏体转变过程。

B点模拟结果如图7所示。可以看到，同未考虑相变潜热的模拟曲线相比，加入相变潜热参数后模拟出的三条冷却曲线更接近于实测曲线。但是，考虑相变潜热模拟出的三条曲线高度接近，高温段840～400℃冷速均为21.5℃/min。此结果说明，对于淬透性不佳的18MnNiMo钢材质的大壁厚核电锻件而言，当模拟参数相变潜热在105J/kg的量级范围内变动以及改变组织转变的类型均对模拟结果影响不大。即考虑组织场的相变潜热作用会提升温度场的模拟准确度，使得模拟结果更贴近实测温度曲线，但相变潜热参数的具体数值变化并不显著影响模拟曲线的精度。如何计算真实的能反映淬火设备和淬火介质的表面换热系数参数是进一步得到更贴近真实测温曲线的模拟结果的关键性因素。

图6 18MnNiMo钢CCT连续转变图Figure 6 CCT continuous transitiondiagram of 18MnNiMo steel

图7 加入相变潜热的模拟结果Figure 7 Simulation results after addingphase change latent heat parameter

3.3 计算表面换热系数后的模拟结果

目前计算表面换热系数H的方法主要是反传热法，即通过实际的测温曲线反推导换热系数H。反传热法又分多种计算方式，应用较广的是近表面两测点差分解法和非线性迭代估算法。前者比较直观，但误差大。后者需要大量的迭代估值计算，逐步缩小估算值同实际值的误差，但不直观且易不收敛。这里我们选用较为简便直观的近表面两测点差分解法来计算表面换热系数H。

近表面两测点差分解法计算示意图[9]，如图8所示。由此方法，第三类边界条件可以转化成以下公式：

图8 近表面两测点差分解法计算示意图Figure 8 Schematic drawing of calculating the differences oftwo measured points near surface by decomposition method

将图9中328 mm测温件a、b两点的实测温度曲线代入公式(7)中，得到表面换热系数H的计算结果，如图10所示。可以看到，初始淬火时表面换热系数很低，高温阶段900～600℃表面换热系数H值在逐渐增加。这再次说明我们的淬火设备没有明显缩减淬火过程的膜沸腾阶段，因此用固定的表面换热系数模拟大壁厚核电产品的淬火过程是不准确的。换热系数H曲线在600～400℃阶段达到稳定，约600 W/(m2·K)。200～50℃阶段换热系数H较大，达到(800～1 000)W/(m2·K)。

将图10中计算得出的表面换热系数结果添加到ABAQUS模拟程序中，得到328 mm厚测温件B点的模拟冷却曲线，如图11所示。可以看到，加入图10的表面换热系数后，模拟曲线和实测曲线的吻合度大幅度提升，尤其是450℃以上的高温阶段，吻合度较高，两条曲线几乎重叠。模拟出的高温阶段840～400℃的冷却速度为11.3℃/min，也同实测曲线的冷却速度8.2℃/min十分接近。450℃以下的模拟结果同实测曲线仍有一定的差距。这可能是此阶段换热系数计算误差大造成的，也可能是此阶段的一些热物参数信息不够准确造成的。但是，由于18MnMoNi钢的主要相变温度区间为840～400℃，所以，此模拟结果具有对实际产品生产的指导预测意义。

图9 328 mm测温件A、B两点实测温度曲线Figure 9 The curve of measured temperatures of A, B points in workpiece with thickness 328 mm used for temperature test

图10 328 mm测温件表面换热系数计算结果Figure 10 Calculation results of surface heat transfer coefficient of workpiece with thickness 328 mm used for temperature test

图11 计算表面换热系数后的模拟结果Figure 11 Simulation results after calculating surface heat transfer coefficient

4 进一步提高模拟精确度的分析

可以看到，模拟冷却曲线目前在450℃以下同实测曲线仍有一定的差距，主要原因可能是我们用于计算表面换热系数的A、B两点其纵向位置并不在一条直线上，不是严格按照如图7所示的两点位置布置的。这样造成了两点不是一维传热，即这两点并非只有厚度方向上的传热差异，在横向上也有传热差异，这样易产生误差。同时，A、B两点距离表面的距离仍然较大。大锻件在冷却时厚度方向上温差很大，这样计算出的表面换热系数很不稳定，误差大。所以，在测定计算表面换热系数时，选取的点应在同一厚度方向上，并且尽可能的贴近淬火表面。

表面换热系数能够较好的反应淬火设备的冷却能力，所以应尽可能多的测定锻件不同位置的淬火表面换热系数。针对不同淬火设备和锻件产品，建立换热系数场，形成数据库。另外，准确测定材料的热物参数。材料在室温至900℃的导热系数、热扩散系数、比热容、弹性模量、泊松比等参数都应准确测定。这样可以大幅度提升今后的模拟精度。

5 结论

(1)运用ABAQUS有限元软件，充分考虑相变潜热和计算换热系数后，模拟了淬火温度场， 840～400℃冷却阶段模拟精度高，具有指导预测意义。

(2)对于18MnMoNi钢核电大壁厚锻件而言，考虑组织场的相变潜热作用会提升温度场的模拟准确度，使得模拟结果更贴近实测温度曲线，但具体相变的类型及潜热参数的数值变化并不显著影响模拟曲线的精度。

(3)对于大型核电锻件，淬火设备没有明显缩减淬火过程的膜沸腾阶段，不可以用固定的表面换热系数作为参数模拟大壁厚核电产品的淬火过程。

[1] ASM Handbook. Vol.4 Heat Treating[M]. 1991:69.

[2] 刘庄, 吴肇基, 吴景之, 等. 热处理过程的数值模拟[M]. 北京: 科技出版社, 1996.

[3] B.Hildenwall，T.Ericsson. Prediction of Residual Stresses in Case-Hardening Steels[C]. Harden ability Concepts with Applications to Steel. TMS-AIME, Warrendale, PA, 1978:579-606.

[4] W.A.Johnson and R.F.Mehl. Reaction Kinetics of Nucleation and Growth[J]. Trans. AIME, 1939,135:416-518.

[5] D.F.Koistinen and J.H.Woodhead. Kinetics of the Nucleation and Growth of Proeutectoid Ferrite in some Iron-Carbon-Chromium Alloys[J]. JISI, 209(11):883-889.

[6] 陈红宇, 杜军毅, 邓林涛, 等. 核反应堆压力容器锻件用SA508系列钢的比较和分析[J]. 大型铸锻件, 2008(1):1-3.

[7] Fletcher A J. Thermal Stress and Strain Generation in Heat Treatment[M]. Springer Netherlands, 1989.

[8] 张立文, 赵志国, 范权利. 35CrMo钢大锻件淬火过程组织分布的计算机预测[J]. 金属热处理学报, 1994, 15(4): 15-20.

[9] Kamamoto S, Nishimori T, Kinoshita S. Analysis of residual stress and distortion resulting from quenching in large low-alloy steel shafts[J]. Materials Science & Technology, 2013, 1(10):798-804.

编辑 杜青泉

Research on Numerical Simulation of Quenching Temperature Field of 18MnNiMo Steel Thick Walled Nuclear Power Forging

Du Dongxu, Ren Liguo, Zhao Deli, Li Jiaju

Numerical simulation has been carried out for cooling process of forging with wall thickness 300 mm by finite element simulation software. By comparing with the measured cooling rate, the precision of numerical simulation has been optimized.

18MnNiMo；forging; quenching; temperature field; numerical simulation

2016—09—20

TG156

B