DOI: 10.11817/j.ysxb.1004.0609.2021-36677

钨动态回复过程的相场模型建立及其仿真模拟

李作胜，唐  赛，梁超平，刘文胜，马运柱

(中南大学 轻质高强结构材料国家级重点实验室，长沙 410083)

摘  要：本文通过耦合Estrin-Mecking(E-M)位错密度模型，建立了金属动态回复微观组织演变的多相场模型，并研究多晶钨在温度1523~1723 K和应变率0.001~1 s^-1范围内发生动态回复过程的力学响应行为，分析其应力-应变曲线、平均晶粒尺寸和动态回复体积分数等变化规律。结果表明：在不同温度和应变率条件下模拟获得的应力-应变曲线和位错密度分布与实验研究相吻合；随着温度升高，动态回复过程中晶界存在细微移动，平均晶粒尺寸相应地增大，且温度越高愈显著。通过量化分析动态回复体积分数，发现提高温度或应变率均可加速动态回复进程，这与已有理论和实验结果一致。

关键词：动态回复；相场法；位错密度；钨；微观组织模拟

文章编号：1004-0609(2021)-07-1767-07　　     中图分类号：TG146.4　　     文献标志码：A

引文格式：李作胜, 唐  赛, 梁超平, 等. 钨动态回复过程的相场模型建立及其仿真模拟[J]. 中国有色金属学报, 2021, 31(7): 1767-1773. DOI: 10.11817/j.ysxb.1004.0609.2021-36677

LI Zuo-sheng, TANG Sai, LIANG Chao-ping, et al. Establishment and simulation in dynamic recovery process of tungsten using phase field model[J]. The Chinese Journal of Nonferrous Metals, 2021, 31(7): 1767-1773. DOI: 10.11817/j.ysxb.1004.0609.2021-36677

钨因具有高熔点、高强度、高密度等优异性能，广泛应用于航空航天、武器制备等国防工业领域^[1-3]。由于纯钨的再结晶温度和层错能较高，在较低温变形过程中，扩展位错很窄，位错的攀移、交滑移等回复过程较容易进行，因而其在变形过程中动态回复占主导作用。动态回复和再结晶均为金属变形过程中两类软化效应，动态回复过程伴随着位错密度的变化，强烈影响动态再结晶的形核。若动态回复较强，位错密度大幅度减低，缺乏足够的能量驱动再结晶形核，使得动态再结晶难以发生。因此，研究动态再结晶的前一阶段—动态回复，对认识材料整个热变形过程的微观组织和力学响应行为变化具有重要的意义。

针对金属动态回复现象，大量研究者基于不同角度提出了动态回复过程中应力-应变曲线的数学物理等模型^[4-6]。其中E-M模型^[5]被广泛应用于描述高层错能金属变形回复过程的流变行为。CHALIMBA等^[7]采用这一模型研究VN和Nb-Ti微合金钢在热压缩变形过程中的高温流变行为，分析其动态回复软化动力学和临界值等，所得结果和实验研究相吻合。LI等^[8]通过建立基于E-M位错密度理论、加工硬化和动态回复理论的本构方程，输入实验数据，定量研究了贝氏体钢回复过程中的位错密度变化。MOSTAFAEI等^[9]基于这一模型提出一种描述动态回复过程动力学的新方法，研究了Al-Mg合金高温流变动态回复行为动力学。然而，这些数学模型均未与材料具体的微观组织结合起来，只处于数学解析层面，无法同步再现动态回复过程中的微观组织演变完整过程。

近年来，随着各类微观组织演变的数值模拟方法出现，例如蒙特卡罗法(MC)^[10-11]、元胞自动机法(CA)^[12-13]以及相场法(PF)^[14-16]等，为研究微观组织演变过程提供了新的思路。GOETZ等^[17]采用CA模拟金属动态回复过程，考虑每个时间步长随机选取一定数量的元胞，将其位错密度降低二分之一。但这一方法与材料实际位错变化存在一定差别，计算结果误差较大。而相场法因模型更准确、对界面处理更方便、参数物理意义更清晰等优点，在模拟材料凝固^[18]、相变^[19]、变形等方面得到了广泛的应用，成为材料微观组织演变研究强有力的手段之一。TAKAKI等^[20]通过建立多相场和位错密度耦合的模型模拟了纯铜在不同热变形参数下动态再结晶过程，证明了基于位错密度建立的多相场模型研究动态回复、动态再结晶过程是可行的。然而，由于目前还未有研究动态回复的相场模型，使得当前金属动态回复的相场模拟研究仍非常少见。

因此，针对目前实验上在研究不同热变形参数对回复过程的影响规律和微观组织演化全过程所面临的技术困难，本文在已有多相场模型的基础上，耦合描述回复过程位错密度变化的E-M模型，建立的相场模型将位错密度的变化和微观动态组织演化全过程相结合，可将实验数据输入相场模型，实时考察不同热变形参数(不同温度和应变率等)对动态回复过程微观组织演变的影响，并分析其动态回复动力学，模拟所得到的结果和实验一致。

1  模型与方法

1.1  相场模型耦合塑性变形E-M模型

本文采用STEINBACH等^[21]提出的多相场模型耦合描述动态回复过程位错密度变化的E-M模型。通过引入一系列连续场变量来描述多晶体中的微观组织，界面设置为具有一定厚度的弥散界面。在含有N个晶粒的体系当中，代表第i号晶粒的相场值，在晶粒的内部=1，晶粒的外部=0，界面处的相场变量值在0到1之间连续变化。并且在计算区域的每个格点上的所有相场变量总和均满足：

                                  (1)

使用的自由能泛函表达式为

                   (2)

式中：、、分别代表梯度能、势能和晶粒内部的自由能密度；V是体积。根据Ginzburg- Landau动力学方程：

, (i=1, 2, 3, …, N)            (3)

求解得到

                   (4)

将式(2)代入式(4)中，化简得到如下相场演化方程：

          (5)

式中：a_ij是晶粒i和晶粒j之间的界面梯度系数；W_jk 为能垒值大小；为相场迁移率。、、与界面厚度、界面能和晶界迁移率有关，分别满足

, ,       (6a)

晶界迁移率为

                            (6b)

式中：和D分别为晶界厚度和晶界自扩散系数，k是玻尔兹曼常数，Q_b是晶界扩散激活能。

晶界迁移驱动力等于两个相邻晶粒之间的储存能之差，即 。

本文采用的E-M模型具体如下：

                (7)

假设初始位错密度为，对式(7)两边同时积分，则

             (8)

式中：为位错密度；等式的右边第一项代表加工硬化部分，第二项代表动态回复过程的贡献，h、r分别是加工硬化速率和给定条件(温度和应变率)下的动态回复速率。通过式(8)即可求得任意时刻下每个格点的位错密度值。变形过程中流变应力与平均位错密度满足。

同时在本文开发的程序中考虑了KIM^[22]提出的优化算法，将在含有N个晶粒体系中每个格点需要存储N个相场变量简化为只需要存储5个相场变量值，这样在不影响计算结果精度的前提下，既节约资源，也提高了计算效率。

1.2  计算参数的获取设置

本文采用上述的多相场模型与E-M模型相耦合技术，以多晶钨为研究对象，模拟其在不同的热变形参数条件下的动态回复过程。计算区域大小为800×800个格子，网格尺寸为dx=dy=0.5 μm。初始微观组织通过传统理想晶粒长大获得，整个区域内每个格点赋予相同大小的初始位错密度1×10⁹/m²。应变随着时间的不断增加而增加即。对应体系的位错密度也随之改变。计算时间步长，其余参数设置详见表1。

表1  模型中的参数取值

Table 1  Value of simulation parameters

在E-M模型中h、r两个参数的数值均可由实验数据求得，具体获取方法详见文献[23]。通过计算拟合实验数据曲线，得到斜率-0.5r和截距0.5r，将不同条件下计算得到的r和代入式(9)，即可得到加工硬化速率h：

                            (9)

此外，为了计算不同条件下的h、r值，考虑剪切模量是与温度相关的函数，不同温度(K)条件下的剪切模量由计算得到，其中弹性模量E为

     (10)

2  结果与讨论

2.1  不同热加工条件下的力学响应行为

通常在实际的变形过程中，材料的变形速率和温度均对材料的性能和微观组织产生较大的影响。图1所示分别为固定温度1523 K在不同应变率和固定应变率0.01 s^-1在不同温度条件下模拟所得的应力-应变曲线，模拟结果与实验数据^[28]相比较整体误差小于5%，吻合较好，这同时也验证了模型的准确性。可以看出，在回复过程中主要是加工硬化占主导作用，随着应变的不断增加，应力先快速增长，最后逐渐趋于稳态值，且最终稳态应力随着变形温度升高或应变率减小而减小。这可能是由于温度越高，使得变形材料中存在的各类点缺陷、位错等较快地发生移动和湮灭，而应变率越小则在变形过程中有充足的时间进行动态回复，内部能量降低更多，最终稳态应力越小。

图1  温度1523 K时不同应变率下和应变率0.01 s^-1时不同温度下多晶钨的应力-应变曲线实验值和模拟结果对比

Fig. 1  Comparisons between simulated and experimental stress-strain curves of polycrystalline tungsten under different deformation conditions

2.2  晶粒组织

在温度和变形量的共同作用下，外界力学响应变化的同时，势必造成内部微观组织发生变化。因此，为探究动态回复过程中晶粒大小和微观组织的变化情况，分别固定应变率和温度的条件下，系统观察平均晶粒尺寸的变化。图2(a)所示为应变率0.01 s^-1时在不同温度条件下平均晶粒尺寸的变化。由图2(a)可知，在动态回复阶段平均晶粒尺寸略有增大，且随着变形温度升高，变化越明显。图2(b)所示为在温度1623 K时应变率为0.01 s^-1不同时刻下的晶粒尺寸分布高斯拟合曲线，横坐标为每个晶粒尺寸大小和平均晶粒尺寸大小的比值，纵坐标为对应的频率；可以看出，随着晶粒不断演化，晶粒尺寸分布曲线的半高宽(FWHM)呈不断增大的趋势，对应材料内部显微组织具体表现为小晶粒变小，大晶粒长大，存在大晶粒吞噬小晶粒的现象。

同时，结合图3所示在不同时刻下相场模拟所得的晶粒微观组织的演化过程，从黄色圈中观察可知，在动态回复阶段存在一些小角度晶界细微移动，但变化较小且缓慢。再者由于本文在设置位错密度时晶粒内部和外部给予不同的权重，使得晶界上的位错密度普遍高于晶粒内部，由图4所示在2000 dt时刻下的晶粒位错密度分布图发现，模拟结果较为理想，晶界附近的位错密度高于晶内，多个晶粒的晶界交点处位错密度最高，这和已有经典动态回复过程中的晶粒形态、大小和初始状态基本相同的现象一致。

2.3  动态回复动力学

在动态回复过程中，应变率为

                            (11)

式中：表示动态回复开始时对应的应变值，t为应变由变化至所需的时间。

图2  应变速率0.01 s^-1时不同温度下平均晶粒尺寸的变化和温度1623 K、应变率0.01 s^-1时不同时刻的晶粒尺寸分布规律

Fig. 2  Average grain size versus true strain when strain rate =0.01 s^-1 at different temperatures(a) and grain size distribution curves at different times with T=1623 K and =0.01 s^-1(b)

动态回复体积分数可通过式(12)求得

                  (12)

动态回复体积分数表示为任一时刻应力的增量和总增量值之间的比值，动态回复体积分数可直接通过模拟所得的应力-应变曲线得到。图5所示为应变率0.01 s^-1时在不同温度下和温度1573 K时不同应变率下动态回复体积分数分别随时间变化的情况。由图5可知，在一定温度或者应变率下，初期的回复速率很大，之后随着时间增加，回复速率逐渐趋于零，但总动态回复体积分数则不断增大，最终接近于1，证明此时回复处于稳定的动态阶段。温度为1523 K的回复速率明显低于其他温度，以及应变率为0.001 s^-1的回复速率比其他更高速率的回复缓慢很多，因此，通过提高温度或者增大变形应变率，将有利于动态回复进程。这可为在实际生产中控制金属材料在变形过程中发生动态回复的工艺参数选择提供理论参考。

图3  温度1623 K、应变率0.01 s^-1条件下微观组织演变过程

Fig. 3  Time evolution of microstructure when T=1623 K and =0.01 s^-1

图4  温度1623 K、应变率0.01 s^-1条件2000 dt时刻的晶粒位错密度分布图

Fig. 4  Dislocation density distribution map(a) and its partial map(b) obtained when T=1623 K and =0.01 s^-1 at 2000 dt

图5  应变率0.01 s^-1时不同温度下和温度1573 K时不同应变率下的动态回复体积分数曲线

Fig. 5  Dynamic recovery fraction curves under different conditions

3  结论

1) 通过耦合Estrin-Mecking位错密度模型，建立了模拟金属变形动态回复过程的多相场模型。

2) 在固定温度为1523 K不同应变率和固定应变率为0.01 s^-1不同温度两种情况下得到的应力-应变曲线和实验结果吻合较好。

3) 在不同温度下，动态回复过程中晶界存在细微移动，随着温度升高，平均晶粒尺寸相应增大，且温度越高愈显著。位错密度大小分布满足晶粒内部＜两晶界处＜多晶界交点的规律。

4) 变形温度和应变率的增大加速动态回复的进程。

REFERENCES

[1] DUMMER T, LASALVIA J C, RAVICHANDRAN G. Effect of strain rate on plastic flow and failure in polycrystalline tungsten[J]. Acta Materialia, 1998, 46(17): 6267-6290.

[2] LENNON A M, RAMESH K T. The thermoviscoplastic response of polycrystalline tungsten in compression[J]. Materials Science and Engineering A, 2000, 276(1/2): 9-21.

[3] ZHANG X, YAN Q, YANG C, et al. Microstructure, mechanical properties and bonding characteristic of deformed tungsten[J]. International Journal of Refractory Metals & Hard Materials, 2014, 43: 302-308.

[4] LIN Y C, WEN D X, HUANG Y C, et al. A unified physically based constitutive model for describing strain hardening effect and dynamic recovery behavior of a Ni-based superalloy[J]. Journal of Materials Research. 2015, 30(24): 3784-3794.

[5] ESTRIN Y, MECHING H. A unified phenomenological description of work hardening and creep based on one-parameter models[J]. Acta Metallurgica, 1984, 32(1): 57-70.

[6] HUANG C Q, DENG J, WANG S X, et al. A physical-based constitutive model to describe the strain-hardening and dynamic recovery behaviors of 5754 aluminum alloy[J]. Materials Science and Engineering A, 2017, 699: 106-113.

[7] CHALIMBA S A J, MOSTERT R J, STUMPF W E, et al. Modeling of high-temperature flow stress of VN and Nb-Ti microalloyed steels during hot compressive deformation[J]. Journal of Materials Engineering and Performance, 2020, 29(1): 330-341.

[8] LI L, YE B, LIU S, et al. Inverse analysis of the stress–strain curve to determine the materials models of work hardening and dynamic recovery[J]. Materials Science and Engineering A, 2015, 636: 243-248.

[9] MOSTAFAEI M A, KAZEMINEZHAD M. A novel approach to find the kinetics of dynamic recovery based on hot flow curves[J]. Materials Science and Engineering A, 2012, 544: 88-91.

[10] IVASISHIN O M, SHEVCHENKO S V, SEMIATIN S L. Implementation of exact grain-boundary geometry into a 3-D Monte-Carlo (Potts) model for microstructure evolution[J]. Acta Mater, 2009, 57(9): 2834-2844.

[11] RADHAKRISHNAN B, SARMA G B, ZACHARIA T. Modeling the kinetics and microstructural evolution during static recrystallization—Monte Carlo simulation of recrystallization[J]. Acta Mater, 1998, 46(12): 4415-4433.

[12] KROC J. Application of cellular automata simulations to modeling of dynamic recrystallization[J]. Lecture Notes in Computer Science, 2002, 2329: 773-782.

[13] PECZAK P, LUTON M J. A Monte Carlo study of the influence of dynamic recovery on dynamic recrystallization[J]. Acta Metallurgica et Materialia, 1993, 41(1): 59-71.

[14] GAO Y, HUANG L, DENG Q, et al. Phase field crystal simulation of dislocation configuration evolution in dynamic recovery in two dimensions[J]. Acta Mater, 2016, 117: 238-251.

[15] 徐树杰, 师春生, 赵乃勤. 热加工过程中动态再结晶现象的多相场研究[J]. 物理学报, 2012, 61(11): 399-408.

XU Shu-jie, SHI Chun-sheng, ZHAO Nai-qin. Dynamic recrystallization phenomenon in hot-working process by multi-phase-field model[J]. Acta Phys Sin, 2012, 61(11): 399-408.

[16] 罗志荣, 高英俊, 朱甜霞. 相场方法研究变形镁合金的晶粒分布[J]. 中国有色金属学报, 2013, 23(6): 1496-1506.

LUO Zhi-rong, GAO Ying-jun, ZHU Tian-xia. Phase field model for grain distribution of deformation magnesium alloy[J]. The Chinese Journal of Nonferrous Metals, 2013, 23(6): 1496-1506.

[17] GOETZ R L, SEETHARAMAN V. Modeling dynamic recrystallization using cellular automata[J]. Scripta Materialia, 1998, 38(3): 405-413.

[18] BOETTINGER W J, WARREN J A, BECKERMANN C, et al. Phase-field simulation of solidification[J]. Annu Rev Mater Res, 2002, 32:163-194.

[19] YAMANAKA A, TAKAKI T, TOMITA Y. Multi-phase- field modeling of diffusive solid phase transition in carbon steel during continuous cooling transformation[J]. Journal of Crystal Growth, 2008, 310(7): 1337-1342.

[20] TAKAKI T, HIROUCHI T, HISAKUNI Y, et al. Multi- phase-field model to simulate microstructure evolutions during dynamic recrystallization[J]. Materials Transactions, 2008, 49(11): 2559-2565.

[21] STEINBACH I, PEZZOLLA F. A generalized field method for multiphase transformations using interface fields[J]. Physica D Nonlinear Phenomena, 1999, 134(4): 385-393.

[22] KIM S.G. Computer simulations of two-dimensional and three-dimensional ideal grain growth[J]. Phys Rev E, 2006, 74(6): 061605.

[23] MOMENI A, DEHGHANI K, EBRAHIMI G R. Modeling the initiation of dynamic recrystallization using a dynamic recovery model[J]. Journal of Alloys and Compounds, 2011, 509(39): 9387-9393.

[24] ALFONSO L A. Thermal stability of warm-rolled tungsten[D]. Danmarks: Technical University of Denmark, 2016: 7-187

[25] GIANNATTASIO A, YAO Z, TARLETON E, et al. Brittle-ductile transitions in polycrystalline tungsten[J]. Philosophical Magazine, 2010, 90(30): 3947-3959.

[26] MANNHEIM A, VAN DOMMELEN J A W, GEERS M G D. Long-term microstructural evolution of tungsten under heat and neutron loads[J]. Comput Mater Sci, 2019, 170: 109146.

[27] MASON D R, YI X, KIRK M A, et al. Elastic trapping of dislocation loops in cascades in ion-irradiated tungsten foils[J]. Journal of Physics: Condensed Matter an Institute of Physics Journal, 2014, 26(37): 375701.

[28] WANG J, ZHAO G, CHEN L, LI J. A comparative study of several constitutive models for powder metallurgy tungsten at elevated temperature[J]. Materials & Design, 2016, 90: 91-100.

Establishment and simulation in dynamic recovery process of tungsten using phase field model

LI Zuo-sheng, TANG Sai, LIANG Chao-ping, LIU Wen-sheng, MA Yun-zhu

(National Key Laboratory of Science and Technology for National Defense on High-strength Structural Materials, Central South University, Changsha 410083, China)

Abstract: Coupling with the Estrin-Mecking(E-M) dislocation density model, a multiphase field model was established to simulate the dynamic recovery process of tungsten in this work. The mechanical response behavior of polycrystalline tungsten during the dynamic recovery process was studied under the circumstances that strain rate ranging from 0.001 s^-1 to 1 s^-1, and temperature ranging from 1523 K to 1723 K. The dependences of the stress-strain curves, average grain size and dynamic recovery volume fraction on strain rate and temperature were carefully analyzed. The results reveal that stress-strain curves and dislocation density distribution at different temperatures and strain rates are in good agreement with experimental observations. As the temperature increases, the grain boundaries moves slightly, and the averaged grain size increases accordingly. Moreover, by quantifying the volume fraction of dynamic recovery, it is found that increasing the temperature or strain rate can accelerate the dynamic recovery process. This agrees well with previous theories and experimental results.

Key words: dynamic recovery; phase-field method; dislocation density; tungsten; microstructure evolution simulation

Foundation item: Project(51931012) supported by the National Natural Science Foundation of China

Received date: 2021-02-25; Accepted date: 2021-05-10

Corresponding authors: LIU Wen-sheng; Tel: +86-731-88877998; E-mail: liuwensheng@csu.edu.cn

MA Yun-zhu; Tel: +86-731-88877825; E-mail: zhuzipm@csu.edu.cn

(编辑  何学锋)

基金项目：国家自然科学基金资助项目(51931012)

收稿日期：2021-02-25；修订日期：2021-05-10

通信作者：刘文胜，教授，博士；电话：0731-88877998；E-mail：liuwensheng@csu.edu.cn

马运柱，教授，博士；电话：0731-88877825；E-mail：zhuzipm@csu.edu.cn