聂君锋 , 汤镇睿 , 张海泉 , 李红克 , 王鑫
清华大学 核能与新能源技术研究院, 北京 100084

收稿日期：2016-03-24
基金项目：国家自然科学基金资助项目（11202114）；北京高等学校青年英才计划项目（YETP0156）
作者简介：聂君锋(1982—), 男, 副研究员。E-mail:niejf@tsinghua.edu.cn


摘要：晶体塑性理论是将晶体微观尺度的位错运动与宏观尺度的塑性形变相结合的重要理论，提供了在细观尺度内研究材料力学行为的有效方法。位错的密度变化对金属晶体的硬化行为有着重要的影响。该文在晶体塑性理论的基础上引入位错运动理论，建立基于位错密度的体心立方晶体（body center cubic，BCC）塑性本构模型，研究BCC的力学行为；并借助ABAQUS有限元软件，编写UMAT子程序，实现对BCC结构的铁单晶及多晶单轴拉伸试验的数值模拟。结果表明：该本构模型能有效地模拟铁单晶及多晶单轴拉伸的力学行为。
关键词：晶体塑性位错运动体心立方晶体（BCC）单轴拉伸有限元
Crystal plasticity constitutive model for BCC based on the dislocation density
NIE Junfeng, TANG Zhenrui, ZHANG Haiquan, LI Hongke, WANG Xin
Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China


Abstract: Crystal plasticity theory is a fundamental theory that combines the crystal microscopic slip mechanism with macroscopic plastic deformation to predict meso-scale plastic deformation. The dislocation density has an important influence on the hardening behavior of metal crystals. This paper presents a constitutive model based on crystal plasticity theory and dislocation motion theory for the BCC crystal structure. The model is used to study the mechanical behavior of a BCC lattice. Using the UMAT subroutine in ABAQUS for numerical simulations of a uniaxial tensile tests of single crystal and polycrystal iron. The results show that the constitutive model effectively simulates the mechanical behavior of the uniaxial tensile test for single crystal and polycrystal iron.
Key words: crystal plasticitydislocation motionbody center cubic (BCC)uniaxial tensilefinite element method (FEM)
金属晶体材料变形是一个复杂变化的过程。影响金属材料力学行为的因素有很多，例如：环境温度、加载方向、加载速度、材料处理方式、滑移方向等。金属晶体塑性理论从微观的运动机制出发，与宏观的塑性变形行为相结合，从而揭示晶体材料的变形规律。在已有研究基础上，Hill和Rice^[1-3]对晶体塑性变形几何学和运动学进行了严格的数学描述，并把模型应用到率无关的弹塑性有限变形分析。Asaro^[4-5]和Peirce等^[6]进一步发展了晶体塑性本构理论，并应用到率相关的黏塑性分析，使得金属晶体塑性理论模型逐步完善。
自20世纪40年代起，由于航空事业的飞速发展，人们对结构设计和计算精度提出了更高的要求，从而逐渐产生了矩阵力学分析法。1960年Clough^[7]在处理平面弹性问题时，首次提出并使用有限元方法，其基本思想是将结构离散化，将复杂的整体结构分解为一系列单元，再对标准的单元进行处理，然后将处理后的单元进行组装，最后进行数值求解。这种思想在诞生之初受到计算能力的影响，并未达到广泛使用。但随着计算机性能的提高，如今已使得大规模计算成为现实。
晶体塑性有限元法是将晶体塑性理论与有限元理论相结合得到了一种新的计算细观尺度内材料塑性变形的方法。这种方法与传统有限元方法相比，更加接近金属晶体变形的本质，使得模拟变形过程更加真实可靠。
本文在经典的晶体塑性理论框架上，结合位错运动理论，以BCC(body center cubic)晶体结构为研究对象，建立一种模拟BCC结构金属材料力学性能的本构模型。
通过ABAQUS有限元平台编写了用户子程序UMAT，对该本构模型进行数值实现。对BCC结构铁单晶的力学行为进行数值模拟，该本构模型能够很好地模拟铁单晶单轴拉伸过程中的力学行为。同时，通过采用Voronoi生成铁多晶结构，模拟铁多晶平面应变过程中的力学行为。
1 基于位错运动的BCC晶体塑性本构模型1.1 晶体塑性理论在晶体塑性理论^[8]中，总的变形梯度F可以分为2部分：晶体沿着滑移系滑移所对应的变形梯度F^p；晶格畸变和刚性转动所对应的变形梯度F^*。

$F = {F^ * }{F^P}.$

(1)

令L表示速度梯度张量，将其分解为弹性变形部分L^*和塑性变形部分L^p。

$L = \dot F{F^{ - 1}} = {L^ * } + {L^P},$

(2)

${L^ * } = {{\dot F}^ * }{F^{ * - 1}},$

(3)

${L^P} = {{\dot F}^p}{F^{p - 1}} = \sum\limits_{a = 1}^n {{{\dot \gamma }^\alpha }{s^{ * \alpha }} \otimes {m^{ * \alpha }}} .$

(4)

其中s^*α和m^*α分别为第α滑移系在中间构型中的滑移方向和滑移面的法向。式(1)—(4) 将晶体的宏观变形与微观滑移剪切应变率联系起来。
定义Schmid分解剪应力τ^α为

${\tau ^\alpha } = \tau :{P^\alpha }.$

(5)

其中：τ表示Kirchhof应力张量，

${P^\alpha } = \frac{{{s^\alpha }{m^\alpha } + {m^\alpha }{s^\alpha }}}{2}.$

1.2 位错滑移运动在连续本构模型^[8]中，晶体第α个滑移系的剪切应变率${\dot \gamma ^\alpha }$通常采用如下的热激活形式^[9]：

$\left\{ \begin{array}{l}{{\dot \gamma }^\alpha } = 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left| {{\tau ^\alpha }} \right| < {g^\alpha };\\{{\dot \gamma }^\alpha } = \dot \gamma _0^\alpha \exp \left\{ { - \frac{{{Q_0}}}{{kT}}{{\left[ {1 - {{\left( {\frac{{\left| {{\tau ^\alpha }} \right| - {g^\alpha }}}{{{{\hat \tau }^\alpha }}}} \right)}^p}} \right]}^q}} \right\}{\mathop{\rm sgn}} \left( {{\tau ^\alpha }} \right),\;\;\;\;\;\;\;\;\;\;\;\left| {{\tau ^\alpha }} \right| > {g^\alpha }.\end{array} \right.$

(6)

${{\hat \tau }^\alpha } = \hat \tau _0^\alpha \frac{G}{{{G_0}}}.$

(7)

其中：Q₀为在没有外力功作用下能够克服障碍产生滑移的激活能，p和q是流动法则的相关参数，G和G₀分别表示当前温度T和绝对零度时弹性剪切模量，${{\dot \tau }^\alpha }$表示没有热激活能时位错运动的最大滑移阻力，g^α为由材料本身的位错造成的滑移阻力。
研究发现，位错密度的变化对晶体的硬化行为有着重要的影响。在晶体塑性变形过程中，随着晶体中位错密度的不断增加，位错之间的交互作用不断增强，导致位错滑移越来越困难，继续运动所需要的应力便会增强。具体来说，引起晶体应变硬化的机制主要存在以下形式：位错的塞积、交割、移动位错俘获成为固定位错、Frank-Read位错源不断消耗等。为了描述晶体的硬化行为，本文在Taylor硬化律基础上，衍生出一个基于位错密度的硬化公式：

${g^\alpha } = Gb\sqrt {{q_\rho }\sum\limits_{\beta + 1}^N {\left[ {{A^{\alpha \beta }}\left( {\rho _{\rm{M}}^\beta + \rho _{\rm{I}}^\beta } \right)} \right]} } .$

(8)

其中：q_ρ是一个统计参数，用来表示实际位错分布与假想位错的有规则空间分布之间的偏差；b为晶体的Burgers矢量的模；A^αα和A^αβ(α≠β)分别表示自硬化系数与潜硬化系数；ρ_M^β表示移动位错的密度；ρ_I^β表示螺型位错的密度。
对于体心立方晶体而言，位错演变形式包括移动位错的增殖及湮灭、移动位错俘获成为固定位错、固定位错的动态恢复、移动螺旋位错的交叉滑移以及移动位错的攀升等。本文主要考虑前3种位错演变机制，得到移动位错和固定位错的演化公式分别为

$\dot \rho _{\rm{M}}^\alpha = \left( {\frac{{{k_{{\rm{mul}}}}}}{{b{l_{\rm{d}}}}} - \frac{{2{R_{\rm{c}}}}}{b}\rho _{\rm{M}}^\alpha - \frac{1}{{b{\lambda ^\alpha }}}} \right)\left| {{{\dot \gamma }^\alpha }} \right|,$

(9)

$\dot \rho _{\rm{I}}^\alpha = \left( {\frac{1}{{b{\lambda ^\alpha }}} - {k_{{\rm{dyn}}}}\rho _{\rm{I}}^\alpha } \right)\left| {{{\dot \gamma }^\alpha }} \right|.$

(10)

其中：k_mul表示移动位错增殖系数；l_d表示移动位错片段的平均长度，l_d≈ $1/\sqrt {\sum\limits_{\beta - 1}^N {\rho _M^\beta } } $；R_c表示Burgers矢量相反的2个移动位错发生相互湮灭的临界尺寸；λ^α表示在第α滑移系上移动位错俘获为固定位错的平均自由程，λ^α= $1/{\beta _\rho }\sqrt {\rho _M^\alpha + \rho _1^\alpha } $；k_dyn表示固定位错的动态恢复系数。
2 材料参数确定本文选取Keh^[10]的铁单晶单轴拉伸实验数据进行模拟对比，采用相同的材料物性参数及实验条件。BCC晶体的滑移系统包括3组，分别为{110}〈111〉、{112}〈111〉和{123}〈111〉。在晶体学上，3组滑移系分别有12、12和24个等价的滑移系统，假设每个滑移系具有相同的临界分解滑移阻力初始值$\hat \tau _0^\alpha $，均为390 MPa^[11]。对于每个滑移系上的移动位错及固定位错密度分布，实验中并没有明确的数据统计分布，本文选取常用的假设，即初始移动位错密度和固定位错密度所占份额相同，则密度值约等于80 000 mm^-2，不超过铁单晶单轴拉伸实验^[10]中试样的位错密度上限。弹性常数及物理常数可以从文[11-14]中获得，其中弹性模量为C₁₁=263 GPa、C₂₂=134 GPa、C₃₃=119 GPa、G₀=87.6 GPa、G=82.534 GPa，硬化参数选择为q_ρ=0.06、b=0.248 nm、R_c=1.5 nm、A^αα=1.0、A^αβ=0.2、k_mul=0.034 5、β_ρ=0.074、k_dyn=275，温度设定为298 K，拉伸应变速率为实验速率3.3×10^-4 s^-1[10]。流动参数根据Kocks建议^[15]选择为${\dot \gamma }$=10⁷ s^-1、p=0.47、q=1.1、Q₀=1.32×10^-19 J、${{\hat \tau }_0}$=390 MPa。
3 计算结果与分析3.1 单晶单轴拉伸实验模拟针对单晶拉伸实验模拟，本文建立了22.5 mm×5 mm×2 mm的单晶模型。为了尽可能模拟真实的单轴拉伸实验过程，模型沿长度方向一端采用对称边界条件，另一端施加位移加载。根据金属单晶具有各向异性的特点以及真实实验中拉伸方向，本模型分别在[111]、[110]、[100] 3个晶向进行加载计算，得到相应的数值计算结果与Keh实验^[10]对比，结果如图 1所示。

图 1 晶向计算与实验单轴拉伸应力-应变曲线

图选项

从图 1可以看出，铁单晶在不同晶向上的单轴拉伸力学行为不同，这反映了单晶的各向异性特点，其中[111]方向是密排方向，屈服强度最高。从3个晶向对比中可以发现，本文的模型可以很好地预测铁单晶单轴拉伸过程中的力学行为。其中，对[110]晶向的力学行为可以精确地模拟，但对[111]晶向的模拟过程中会产生一定的偏差。这是由于在滑移系开动后期，移动螺旋位错的交叉滑移以及移动位错的攀升等因素的影响增大，使得模型预测结果产生偏移。整体而言，本模型模拟屈服强度较实验偏低，这是由于本模型中尚未考虑非Schimid效应的影响，降低了其屈服强度。
3.2 多晶单轴拉伸模拟在此基础上，利用本文建立的晶体塑性本构模型模拟了铁多晶的单轴拉伸行为。在模拟多晶平面应变过程中，采用Voroni原理生成多晶模型，在0.45 mm×0.1 mm平面中生成60个晶粒图，在Abaqus中建立二维模型，如图 2所示。在有限元模型中，每个晶粒赋予随机生成的晶向。

图 2 多晶单轴拉伸模型

图选项

采用与单晶相同的材料参数，通过UMAT子程序进行模拟计算。对于铁多晶而言，其力学性能是通过其包含的晶粒的平均响应状态来确定的。本文中铁多晶的晶向是随机的，晶粒之间符合纯几何边界条件假定。在计算中，采用体积平均的方法计算应力和应变：

$\left\{ \begin{array}{l}\bar \sigma = \sum\limits_{k = 1}^N {{v^k}{\sigma ^k}} ,\\\bar \varepsilon = \sum\limits_{k = 1}^N {{v^k}{\varepsilon ^k}} .\end{array} \right.$

(11)

其中：σ和ε分别为平均应力和平均应变，σ^k和ε^k分别为第k个积分点的平均应力和平均应变，υ^k为第k个积分点所占据的体积。
根据式(11) 得到拉伸方向的平均应力-应变曲线如图 3所示。

图 3 多晶与单晶模型单轴拉伸应力-应变曲线

图选项

从图 3可以看出，模拟多晶单轴拉伸曲线介于3个晶向之间。这是由于多晶的晶向是随机给定的，因此平均化后得到的多晶单轴拉伸的应力-应变曲线，应介于单晶的屈服强度最高和最低晶向的单轴拉伸应力-应变曲线之间，这也符合理论分析的结果^[16]。
多晶模型单轴拉伸得到的应力云图及应变云图如图 4所示。对比模型晶粒分布可以发现，在晶粒交界处应力、应变较大，而晶粒内部的应力、应变较小，这体现了晶界对位错的塞积作用，从而使得晶界处产生应力集中，发生非均匀变形。

图 4 多晶模型单轴拉伸得到的应力云图及应变云图

图选项

4 结论本文在经典的晶体塑性理论框架下，结合位错运动理论，建立了BCC晶体塑性本构模型，并借助ABAQUS有限元平台，编写了用户材料子程序UMAT进行了数值实现。对铁单晶及多晶的拉伸实验过程进行了数值模拟，结果表明：该本构模型能够很好地模拟铁单晶及多晶的单轴拉伸力学行为。对于铁单晶而言，不同晶向上铁单晶的单轴拉伸力学行为不同，反映了单晶的各向异性特点，其中[111]方向是密排方向，屈服强度最高。3个晶向的计算结果均与实验符合较好，证明了该模型的有效性。对于铁多晶而言，Voronoi模型能够较好地反映多晶的单轴拉伸行为，模拟结果与理论预测趋势一致，并且体现了晶界处的变形不均匀性。
通过以上结果可以看出，本文建立的本构模型能够较好地反映出BCC材料的力学行为，对于研究和预测材料的力学行为具有重要意义。

参考文献

[1]	Hill R. Generalized constitutive relations for incremental deformation of metal crystals by multislip[J]. Journal of the Mechanics and Physics of Solids, 1966, 14(2): 95–102. DOI:10.1016/0022-5096(66)90040-8
[2]	Hill R. The essential structure of constitutive laws for metal composites and polycrystals[J]. Journal of the Mechanics and Physics of Solids, 1967, 15(2): 79–95. DOI:10.1016/0022-5096(67)90018-X
[3]	Hill R, Rice J R. Constitutive analysis of elastic-plastic crystals at arbitrary strain[J]. Journal of the Mechanics and Physics of Solids, 1972, 20(6): 401–413. DOI:10.1016/0022-5096(72)90017-8
[4]	Asaro R J, Rice J R. Strain localization in ductile single crystals[J]. Journal of the Mechanics and Physics of Solids, 1977, 25(5): 309–338. DOI:10.1016/0022-5096(77)90001-1
[5]	Asaro R J. Micromechanics of crystals and polycrystals[J]. Advances in Applied Mechanics, 1983, 23(8): 1–115.
[6]	Peirce D, Shih C F, Needleman A. A tangent modulus method for rate dependent solids[J]. Computers & Structures, 1984, 18(5): 875–887.
[7]	Clough R W. The finite element method in plane stress analysis[C]//Proceedings of the 2nd ASCE Conference on Electronic Computation. Pittsburgh, USA: ASCE, 1960.
[8]	王自强, 段祝平. 塑性细观力学[M]. 北京: 北京大学出版社, 1995.WANG Ziqiang, DUAN Zhuping. Plastic Meso Mechanics[M]. Beijing: Peking University Press, 1995. (in Chinese)
[9]	Busso E P. Cyclic Deformation of Monocrystalline Nickel Aluminide and High Temperature Coatings[D]. Cambridge: Massachusetts Institute of Technology, 2005.
[10]	Keh A S. Work hardening and deformation sub-structure in iron single crystals deformed in tension at 298 K[J]. Philosophical Magazine, 1965, 12: 9–30. DOI:10.1080/14786436508224942
[11]	Suzuki T, Koizumi H, Kirchner H O K. Plastic flow stress of b.c.c. transition metals and the Peierls potential[J]. Acta Metallurgica Et Materialia, 1995, 43(6): 2177–2187. DOI:10.1016/0956-7151(94)00451-X
[12]	Johnson R A, Oh D J. Analytic embedded atom method model for BCC metals[J]. Journal of Materials Research, 1989, 4(5): 1195–1201. DOI:10.1557/JMR.1989.1195
[13]	Brunner D, Diehl J. Strain-rate and temperature dependence of the tensile flow stress of high-purity α-iron above 250 K (regime I) studied by means of stress-relaxation tests[J]. Physica Status Solidi, 1991, 124(1): 155–170. DOI:10.1002/(ISSN)1521-396X
[14]	Spitzig W A, Keh A S. The role of internal and effective stresses in the plastic flow of iron single crystals[J]. Metallurgical & Materials Transactions B, 1970, 1(12): 3325–3331.
[15]	Kocks U F. Thermodynamics and Kinetics of Slip[M]. Oxford: Pergamon Press, 1975.
[16]	Hutchinson J W. Plastic stress-strain relations of F.C.C polycrystalline metals hardening according to Taylor's rule[J]. Journal of the Mechanics & Physics of Solids, 1964, 12(1): 11–24.