DOI：10.19476/j.ysxb.1004.0609.2019.08.11
基于滑移与孪生镁合金晶体塑性本构模型及微结构关联分析
蓝永庭1，陈 渊3，任一方2，张克实4，王 帅2
(1. 广西科技大学 职业技术教育学院，柳州 545006；
2. 广西科技大学 机械工程学院，柳州 545006；
3. 广西大学 机械工程学院，南宁 530004；
4. 广西大学 土木建筑工程学院，南宁 530004)
摘 要：为了揭示镁合金晶体微结构与孪生变形之间的关联性，采用伪滑移模式描述孪生变形，建立包含滑移与孪生变形机制的晶体塑性本构关系，推导以剪切应变率作为自变量的本构方程牛顿-拉普森迭代式。采用基于体素方法建立的多晶三维微结构有限元模型，对挤压镁合金材料沿挤压方向加载的单轴拉伸和压缩变形行为进行模拟。模拟结果表明，单轴加载下模拟与试验的材料宏观硬化行为基本吻合，及压缩变形的模拟织构演化与试验织构演化也基本趋同。对多晶体内基本物理量的统计分析表明，该模型能够展现多晶内的应力应变、孪晶体分数在空间的不均匀分布特性，而孪晶体分数分布与晶粒尺寸、晶界位向角及晶界倾角等微结构之间缺乏确切的关联性，表明非均匀微结构所导致的独特应力非均匀状态并不是出现孪生变形与微结构强关联性的原因。
关键词： 镁合金；孪生变形；微结构；塑性本构模型；织构
文章编号：1004-0609(2019)-08-1660-16 中图分类号：TG1；O34 文献标志码：A
航空及运输工业对轻质合金的需求促使人们对比强度高的镁及其合金的变形机制开展了大量的研究工作。由于室温下仅有有限的独立易滑滑移系，镁合金的成形性通常比其他金属材料的低。为了提高金属成形性需要更多地了解多尺度下材料的变形机理。镁合金的晶体结构为密排六方(HCP)，在镁单晶晶格中，原子排列最紧密的是{0001}基面；而在基面上原子最密排方向为
或
，镁合金最常见的滑移系有基面滑移、柱面滑移和锥面滑移，然而基面滑移系和柱面滑移系都只能提供2个独立滑移系，而锥面滑移系可以提供4个独立滑移系。由于位错的基面、柱面和锥面滑移的滑移方向均为垂直于c轴的方向，因此无法协调沿c轴方向的应变。尽管
位错滑移沿
方向可以协调沿c轴方向的应变，但在温室环境下锥面滑移系的启动需要相对比较高的临界分解剪切力(CRSS)，一般情况下难以激发启动[1]，且位错形成的主要区域往是在诸如近表面、晶粒之间的界面、残余孪晶界，以及晶粒内部晶体结构突变和位错纠结积塞等应力集中之处[2]。故常温下，HCP晶体缺少足够独立活动的滑移系数量，而孪生变形可提供额外的独立变形模式来满足Von Mises准则，因而，孪生变形在镁合金塑性变形中扮演着重要角色，成为发挥促进塑性能力作用的一种变形机制[3-5]。
镁合金材料在单轴加载条件下的塑性变形行为具有明显的各向异性性质，不同加载路径下所表现出的应力应变曲线、屈服强度以及硬化率演化等塑性行为特征都存在显著差异[6-8]。与滑移变形相比，孪生变形具有单极性质，且改变晶体取向是一个突发过程，这种特性对镁合金材料的塑性变形行为有着显著影响。因此，从细观的角度研究复杂条件下的镁合金塑性行为及破坏机理是非常有必要的。
早期，研究者们对Cd、Zn晶体开展的试验中发现，孪生应力有一个比较宽的范围，并且样品在轻微弯曲后出现微小孪晶尖端，或存在凹痕的晶体发生孪生变形所需要的应力低于没有弯曲或小凹痕的晶体。这表明孪生形核所需的应力高于孪晶长大所需的应力。研究[9]发现，某种孪晶变体占主导地位取决于由晶体取向确定的Schmid因子(SF)。GODET等[10]针对拉伸孪晶采用SF预测镁合金管的织构变化。BARNETT等[11]的研究表明，在AZ31中孪生系SF对孪生形核比孪晶变厚的影响更明显。这些说明孪生系SF对孪生变形具有重要的影响作用。然而，BEYERLEIN等[12-13]中在变形镁合金晶体中可观测到低SF孪晶出现和高SF孪晶缺失的情况，从另一方面说明了孪生形核及长大不仅仅受到SF的控制。JONAS等[14]指出，孪晶变体的选择可以受到相邻晶粒的应变协调影响，当孪晶在晶界处形成时，相邻晶粒必须发生变形以协调形状变化，即挤出效应。然后相邻晶粒的变形在孪晶和晶体之间产生了驱动应力，同时孪晶变体的选择受到必需变形类型(如柱面滑移或锥面滑移等)和数量大小的影响。在HCP晶体中，晶粒尺寸效应在Zn[15]、Zr[16]、Ti[17]、Mg及变形镁合金[18]中得到了研究。研究结果表明，孪生变形的晶粒尺寸效应比滑移变形的更加强烈，孪生变形与晶粒尺寸存在递增变化关系，大晶粒更容易发生孪生变形。在大量的EBSD研究中，观察到大多数孪晶形核发生在晶界处[12]。KNEZEVIC等[19]的研究揭示了孪生变形机制以及孪晶形核多数发生在晶界区域，观察到孪晶形核始于对称倾斜晶界，形核路径依赖于倾斜角度及晶界的位错状态，特别是孪晶形核更倾向于小角度晶界。综合以上研究表明，孪生变形与微结构有很强的关联性，这种关联可能由两方面的原因造成的：一是微结构形态决定位错组态分布，而位错组态则决定孪生位错的运动结果，即反映孪生形核的难易程度；二是由微结构在空间分布的不均匀性形成独特的应力分布状态而促成孪生变形与微结构的联系。第一种原因是由孪生变形本身性质决定的，使孪生形核的启动强度与微结构性质联系起来；第二种原因是孪生变形外部应力因素促成的，使孪生变形的孪生系分解剪切应力与微结构的不均匀分布联系起来。然而，目前仍没有足够的证据说明哪一种原因在孪生变形与微结构关系中起关键作用，许多关于孪生变形的基本问题仍然没有得到解决。采用各种数值方法有助于解释试验现象，在细观尺度上，晶体塑性模型在理解孪生变形的基本性质方面起了很重要的作用。在晶体模型中，孪生变形通常被假定为按滑移模式来处理，即伪滑移模式。这种伪滑移模式已应用于各种数值分析中来分析孪生形核[12, 20]、孪晶与母晶的相互作用[21-23]、孪生变形的晶格应变演化[24-26]、织构演化[27-29]、各种不同的硬化过程[30]，以及孪生与去孪生变形[31-33]等。然而，采用伪滑移方法模拟孪生变形未能描述孪生变形的动态特性，特别是在孪生变形的开端和形核阶段[18, 34]。这种模式在描述孪生变形与微结构的相关性方面仍缺乏系统分析，但许多学者开始关注这方面的研究[35-36]。本文作者尝试采用伪滑移方法，结合3D建模技术和晶体塑性有限元模型研究镁合金细观变形机制对宏观塑性行为的影响以及孪生变形与微结构之间的关系。探索伪滑移模式能够在多大程度上可以描述孪生变形与微结构的相关性，从而评价伪滑移方法对孪生变形的描述能力，并辨明导致孪生变形与微结构之间强关联性的原因。
1 试样制备及力学实验
试样采用直径为d 20 mm的镁合金AZ31B挤压棒，试样材料的初始平均晶粒尺寸约为20 μm。在500 ℃下进行1 h的时效退火热处理，并氩气保护中空冷以减少试样氧化。随着晶粒再结晶并长大，最终获得平均尺寸约为40 μm等轴晶粒，其中最小晶粒尺寸约为10 μm，最大超过80 μm，如图1所示。
图1 AZ31B试样材料显微组织
Fig. 1 Optical micrograph of AZ31B samples
应用电火花加工(EDM)设备将材料加工成拉伸和压缩试样，拉伸试样总长度为120 mm，其横截面为4 mm×4 mm，如图2(a)所示。压缩试样的截面直径为d 6 mm，长度为12 mm，如图2(b)所示。沿挤压方向(ED)压缩加载能够直接导致拉伸孪生的启动。为了消除电火花加工表面的损伤，在溶液中进行化学腐蚀。在室温下，应用MTS908测试系统以常应变率为1.0×10-4 s-1开展单轴拉伸和压缩测试。
图2 拉伸试样和压缩试样的尺寸示意图
Fig. 2 Schematic diagram of tensile specimen(a) and compression specimen(b) sizes (Unit: mm)
2 晶体塑性本构模型
在有限变形运动中将变形分解为弹性部分和塑性部分，则物质点的总变形梯度按乘法分解可表示为
＞0 (1)
式中：F表示总的变形梯度张量；
表示弹性变形梯度张量，其可描述晶格伸长与刚体转动等部分的弹性变形；
表示不可压缩的塑性变形梯度张量。
同样，速度梯度L由弹性速度梯度Le和塑性速度梯度Lp所组成，即
(2)
其中，
式中：上标“
”表示对时间的导数，上标“-1”表示矩阵求逆。
考虑到镁合金塑性变形主要是由滑移系上的位错滑移和孪生变形所主导。根据塑性速度梯度Lp的定义，塑性流动律的表达式为[37]
(3)
式中：
表示晶粒母体滑移系滑移剪切率；
表示孪晶体滑移系滑移剪切率；
表示孪生剪切应变率；
表示孪生系
的孪晶体积分数；
表示晶粒母体滑移系的启动数量；
表示孪生系的启动数量；
表示孪晶区启动滑移系数量；
表示晶体母体滑移系张量；
表示孪晶区孪生系张量；
表示孪晶区滑移系张量。
于是塑性变形率张量
及塑性旋转张量
分别表示为
(4)
且
(5)
式中：
，
，
为变形后滑移面或孪生面法向向量，
为变形后滑移或孪生方向向量。
2.1 滑移与孪生变形模式
镁合金晶体塑性变形模式主要有滑移变形和孪生变形，通过一个与分解剪切应力有关的指数函数来描述滑移变形率和孪生变形率，根据HUTCHINSON[38]建议的粘塑性关系，滑移系与孪生系分解剪切应变率可分别表示为
(6)
式中：
称为参考剪应变率，是待定材料常数；k为反映材料的率敏感性质的材料常数；、分别表示单晶体第
滑移系、第孪生系的分解剪应变率；
、
分别表示第滑移系、第孪生系的分解剪切应力，且由于孪生变形的单极性质，只有当分解剪切应力与孪生方向一致时才能启动；
、
分别表示第滑移系、第孪生系的应变硬化函数。
晶体发生孪生变形后产生的孪晶体分数的演化与剪切应变率的关系通过孪生系的理论剪切应变联系起来。因此，在每个晶粒中孪生系中孪晶体积分数
的演化可定义为
(7)
式中：表示孪生系的孪晶体积分数变化率；表示孪生系剪切应变率；
表示孪生系的理论剪切应变。
2.2 滑移系与孪生系的强度演化
根据滑移系与孪生系变形特点，其硬化率
演化规律采取以下形式[39]，
(8)
表示滑移或孪生硬化模量，并建议为
(9)
式中：q为常数，通常取q=0.5；
为克罗内克函数。参考文献[40]，滑移系与孪生系的硬化参数
分别采取以下形式：
(10)
式中：各变量下标“S”、“T”分别表示对应滑移系和孪生系的变量；
、
分别表示滑移系与孪生系的累积剪切应变；h0为滑移系和孪生系初始硬化模量；
为滑移系和孪生系初始屈服剪切应力；
为滑移系和孪生系饱和剪切应力，计算时将作为g的初始值；f表示总的孪晶体积分数；r表示孪晶对滑移硬化的耦合参数，通过试验数据拟合确定。
3 晶体本构模型的数值实现
晶体塑性的数值计算是十分复杂且面临很多困难，研究者们提出一种在切线法架构下的显性数值算法，采用Jaumann应力率和变形率来描述本构模型。STEINMANN等[41]提出完全隐性算法，把滑移剪切率看作晶体变形数值模拟的基本变量。ZHANG等[42]将应变分量作为基本变量发展了Newton–Raphson迭代算法。这里，文中结合孪生变形机制，将滑移应变率和孪生应变率作为Newton–Raphson迭代算法的基本变量推导本构模型及求解算法。
3.1 晶体塑性本构方程的时间积分算法
根据剪切应变率是剪切应力和硬化强度的函数，那么应变率的泰勒一阶近似表达式可以写为
(11)
式中：
、
分别表示第变形系在时刻t、
的剪切应变率；
、
分别表示剪切应变率对分解剪切应力、变形阻力的偏微分；
、
分别表示
时间间隔内的应力增量和变形阻力增量。
式(11)中各项由以下方法给出，对式(6)求偏微分可得
(12)
(13)
根据文献[43]，分解剪切力增量可表示为
(14)
式中：
，其中，P为Schmid张量，C为弹性模量张量，W为旋转张量，
为Cauchy应力张量。弹性应变张量
，且考虑到晶粒中的孪晶片层结构交错复杂，孪晶中的滑移变形很难启动，故忽略小量的孪晶区滑移变形，故
(15)
根据变形系阻力的定义，变形阻力增量可表示为
(16)
将以上各式(12)~(16)代入式(11)，整理得
(17)
将对应变形系的剪切应变增量按线性叠加计算，即
(18)
这里，
为控制参数(一般取0＜＜1)，并代入式(17)整理得
(19)
对上式进行进一步整理，可得
(20)
3.2 基于应变变量的牛顿-拉普森迭代法
在数值模拟中，将式(20)改写为(21)的形式，主要是考虑到算法植入ABAQUS/Standard有限元软件的方便性，于是有
(21)
式中：
表示各变形系的剪切应变率。
以剪切应变率为变量的牛顿-拉普森迭代法可表示为
(22)
式中，
，
式中：
。
3.3 一致切线刚度矩阵
由于有限元计算采用Newton–Raphson迭代方法需确定当前刚度矩阵，因此，对一致切线刚度矩阵的估计严重影响到计算收敛速度和计算效率。因此，在本构方程的迭代计算过程中需要不断地更新一致切线刚度矩阵。
分解剪切应力可以由以下公式表示
(23)
对式(23)两边求其物质导数，可推导得[43]
(24)
式中：
及
，
表示弹性变形速度张量。
对于临界滑移系或孪生系，处于临界状态时，有
(25)
式中：D为总的变形率张量。于是，
(26)
解得
(27)
式中：
表示
的逆矩阵的第行向量。
以现时构型为基准，应力增量可表示为
(28)
式中：
为Cauchy应力增量张量，C为弹性模量张量。将式(27)代入式(28)，可得
(29)
令
即为一致切线刚度矩阵。
3.4 弹性变形梯度的计算
由式(2)定义了弹性速度梯度张量，考虑在时间增量内，可得
(30)
由弹性速度梯度、弹性变形率张量和弹性旋转张量的关系，可以得到
(31)
将
及
代入式(31)，可得
(32)
又因为及
，则
(33)
根据文献[44]，有
，所以
(34)
式中：
为时间内刚体旋转矩阵。因此，从t到时刻弹性变形梯度可由t时刻的累积弹性变形梯度表示如下：
(35)
通过
就可以将t时刻当前构型的张量转换为时刻当前构型的张量。
4 多晶镁合金的有限元模拟
4.1 3D微结构模型的构建
采用重构镁合金等效微结构体作为晶体塑性模拟的有限元模型。这里3D多晶微结构数据是由相应的2D扫描成像图中采集而来。采用EBSD技术，从一系列试验中统计晶体形态结构数据，其中统计生成微结构特征的参数包括晶体形态参数的分布函数，如晶粒体积、相邻晶粒数，以及晶体(结晶)结构参数，如晶体取向、晶界角度和微织构。采用分布函数对微结构参数进行描述，以便生成统计等价的3D重构晶体结构模型。以真实晶体形态结构的统计数据为基础，系列生成3D微结构模型的过程包括：1) 生成晶核，通过控制相邻晶核最大距离和最小距离，匹配晶粒尺寸统计分布；2) 采用Voronoi方法，生成等轴晶的晶粒，并计算统计晶粒等效直径、每个晶粒的晶界面的顶点(晶界面为多变形)、晶界面法向，相邻晶粒数以及相邻晶粒序号等信息；3) 根据统计结果，赋予每一个晶粒给定的晶体取向。获得等效织构和晶界角、Schmid因子等数据。
重构3D微结构模型是对真实晶体形态结构的统计数据进行匹配，微结构重构包含的结构特征与真实微结构在统计上相一致。文中所建的镁合金等效微结构体模型，如图3所示。
图3 Voronoi微结构模型
Fig. 3 Voronoi microstructural models
整个微结构模型是边长为0.25 mm的正方体，包含373个晶粒。其晶粒直径、孪晶Schmid因子、晶界倾角的分布对比测量数据在统计上保持相一致，如图4~6所示，分别表示晶粒直径、孪晶Schmid因子、晶界倾角的统计分布。
4.2 加载与模型边界条件
如图3所示，分别将挤出棒的RD方向、TD方向及ED方向定义为坐标系的x轴、y轴及z轴。计算模型施加的边界条件为：1) 约束x =0、y =0及z =0三个平面的法向位移为0；2) 保持平面x =1与y =1为自由平面状态，平面z=1施加z方向的均匀拉伸或压缩应变，最大加载变形应变量为20%。
图4 等效晶粒直径(尺寸)的分布
Fig. 4 Distribution of equivalent grain diameters (size)
图5 晶粒孪生系Schmid因子的分布
Fig. 5 Distribution of Schmid factor of twinning system
图6 晶界倾角的分布
Fig. 6 Distribution of grain boundary misorientation angles
4.3 本构模型参数
将上述的本构模型作为材料子程序(UMAT)植入ABAQUS/Standard有限元分析软件中，采用逐渐逼近法，将应力应变数值计算结果与试验数据不断比较逼近，筛选出一组与应力应变试验曲线相匹配的材料参数，如表1所示。
表1 AZ31单晶硬化参数
Table 1 Single crystal hardening parameters of AZ31
值得指出的是，利用这组硬化参数进行数值模拟，我们可以得到类似于图7所示的应力-应变曲线，表明由这组参数模拟得到的单调拉伸和压缩的应力应变曲线与试验结果比较吻合。
5 实验结果
5.1 镁合金拉伸和压缩加载下的宏观硬化行为
图7所示为挤压镁合金沿挤压方向拉伸和压缩加载下的宏观应力应变试验和模拟曲线。比较图中的试验与模拟曲线可知，两者的宏观硬化行为基本吻合。对拉伸与压缩曲线进行比较分析，可以看出拉伸和压缩宏观硬化行为具有显著的差异，表明镁合金材料力学行为具有明显各向异性特征。从细观的角度来看，拉伸变形的硬化行为主要是由滑移变形来主导，所以表现出上凸的硬化曲线。而压缩变形的硬化行为有孪生变形的参与，具体来说，在压缩开始屈服阶段，塑性变形主要是孪生变形主导，孪生变形开始呈现低硬化率，甚至是软化现象，曲线表现为一段接近水平的线段；随着材料变形的不断增加，孪生应变硬化也随之增加，孪生变形逐渐减少直至耗尽。随应力水平的提高，滑移变形启动并主导后续的变形，因此，整个压缩变形硬化行为表现出不寻常的“S”形曲线。可见，拉伸与压缩宏观硬化行为的差异主要是由变形机制活动情况不同造成的。因此，深入研究细观变形机制是深化理解镁合金各向异性性质的基础。
图7 挤压AZ31单轴拉伸与压缩下的应力-应变试验与模拟结果比较
Fig. 7 Comparison of simulation and experimental results for loading direction stress–strain response of polycrystalline AZ31 in uniaxial tension and compression
5.2 模拟织构及演化分析
图8 挤压态镁合金AZ31的取向极图[45]
Fig. 8 Pole figures of as-extruded AZ31 magnesium alloy[45]
镁合金挤压棒具有典型的挤压织构，即c轴方向多数近似垂直于ED方向(挤压方向)，如图8(a)所示。图8中ED、RD、TD分别代表挤压(轴线)方向、半径方向和切向方向。如图8所示，多晶体模拟初始织构与文献[45]测量的晶体织构基本吻合。
图9所示为微结构数值模型中初始状态下各晶粒的c轴晶向示意图。彩色圆点表示晶粒的晶核及其空间位置，箭头指向为晶粒的c轴方向。从图9可以看出，微结构模型中的晶粒晶格c轴方向大部分垂直或接近垂直ED方向(z轴)，这样的分布符合实际挤压棒材的晶粒取向分布统计。图9形象地呈现晶粒的晶核及取向在空间的位向分布。
图9 微结构模型中晶粒的晶核位置及晶粒c轴方向(其中圆球表示晶核，箭头表示c轴方向)
Fig. 9 Crystal nucleus position and c-axis of grains in microstructure model (Each sphere represents crystal nucleus, arrows represent c-axis)
图10 在ED方向上不同压缩应变的模拟和测量极图[45]
Fig. 10 Simulated and measured pole figures[45] during compressive deformation in ED at different strains
图10所示为镁合金微结构模型在压缩加载下应变分别达到5.5%和10%时各晶粒的取向分布图像。图中分别表示(0001)方向的取向图，以及
方向的取向图。从图10(a)可以看出，(0001)方向已由初始状态时分布于图像边缘变为集中向图像中心位置偏转，而方向则由比较靠近图像中心向图像边缘扩散。这样的取向变化主要是因为孪生变形引起取向的大角度转动造成的。文中微结构织构模拟结果与文献[45]的试验测量结果基本吻合。图10(b)所示应变量为10%时的晶向图像，与图10(a)比较，其图像的变化并不大。这是因为材料压缩变形达到5.5%时，大部分晶粒的孪生变形基本完成，接下来进一步的变形主要是由滑移变形主导的变形过程。滑移变形引起的晶向转动主要是弹性变形梯度导致的刚性转动，其累积转动角度比较小，反映在织构图像的变化不够明显，这与宏观硬化行为的分析结果相一致。
5.3 应力、应变及孪晶的空间分布非均匀性
图11所示为压缩变形量为5.5%时加载方向上的应力S33、应变L33及孪晶体分数的空间分布云图。从图11可以看出，应力S33、应变L33及孪晶体分数在空间的分布极不均匀。在图11(a)中，压缩应力的变化范围从-380~17 MPa，平均应力约为-190 MPa。应力变化从正应力到负应力在空间连续分布，说明微结构内部存在明显的空间应力梯度，这是导致局部损伤或裂纹萌生的基本驱动力。与应力空间分布相似，在图11(b)中应变在空间的分布也是极不均匀，应变变化幅度范围从-0.2%~-15.6%。与应力梯度分布对应，多晶内部存在明显的应变梯度。图11(c)表示孪晶体分数的空间分布，不均匀特征也是明显的，说明孪生变形在各晶粒的变形程度存在较大的差异。从图11(c)可以看到，有些晶粒的孪生体积分数已经达到接近100%，而其中一部分晶粒孪晶体积分数还不到10%。这表明微结构体内各晶粒产生的孪生变形很不均匀，造成孪晶在空间的分布不均匀，从而表现出多晶体内晶粒结构形态的复杂特性，对多晶材料力学性能、损伤或裂纹产生的预测因此会变得更加困难。
以上不均匀现象通常与晶体内产生的变形机制不同有很大关系。晶粒的大小、晶向、晶界倾角等微结构都可能影响到孪生变形和滑移变形的启动，而孪生变形和滑移变形启动的不同组合决定了多晶材料的微结构演化以及宏观硬化行为。
图11 加载方向上的应力、应变空间分布云图和孪晶体分数空间分布云图
Fig. 11 Spatial distribution contours of loading direction stress(a), loading direction strain(b) and volume fraction of twinning(c)
5.4 晶粒尺寸与晶粒应力、应变、孪晶分数的关联性分析
图12~14所示分别为微结构模型中晶粒尺寸与压缩应力、压缩应变和孪晶体积分数的统计关系。图12反映了加载方向的压缩应力S33随晶粒尺寸的变化而变化的关系图。从图12中可以看出，两者之间没有形成线性比例关系，而是在一定范围内无规律的变化。然而，模拟统计数据显示，小尺寸晶粒的应力分布比大尺寸晶粒的应力分布显得更集中一些。WANG等[46]采用三维X射线技术测量了1000多个晶粒的等效应力与晶粒半径的关系分布，如图13所示，等效应力在一定范围内非均匀分布，且与晶粒半径没有确切的函数关系。
图14所示为加载方向的压缩应变L33随晶粒尺寸的变化而变化的关系图。从图14可以看出，与应力分布相似，两者之间也没有形成线性比例关系，而是在一定范围内无规律的变化，在小尺寸晶粒的应变分布比大尺寸晶粒的应变分布也更集中一些。图15所示为晶粒孪晶体分数与晶粒直径关系。从图15中可以看出，在各晶粒的孪晶体分数与晶粒尺寸并没有形成简单的线性关系，其孪晶体分数随晶粒尺寸的变化而变化，且变化幅度比较相似，变化范围在0.1~0.8之间；图中并没有发现在某一尺寸范围内孪晶体分数的变化特别突出。
图12 应变为5.5%时加载方向的应力S33与晶粒直径关系图
Fig. 12 Relationship between calculated stress S33 in loading direction for average strain of 5.5% and grain diameter
图13 等效应力σVM随晶粒半径的变化[46]
Fig. 13 Variation of σVM with grain radius at different strains[46]
图14 加载方向的应变L33与晶粒直径关系
Fig. 14 Relationship between calculated strain in loading direction for average strain of 5.5% and grain diameter
图15 晶粒孪晶体分数与晶粒直径关系
Fig. 15 Relationship between volume fraction of twinning in each grain for average strain of 5.5% and grain diameter
以上对模拟结果的统计分析表明，晶粒尺寸对应力、应变及孪晶体分数的分布不存在确切的函数关系，说明采用伪滑移方法描述孪生变形进行数值模拟时，虽然能够描述应力应变在多晶微结构内部产生的不均匀特征，但在描述孪生变形与晶粒尺寸的关系上与试验现象仍存在较大的偏差，关于这方面的问题将在讨论部分进一步阐述。
5.5 晶界位向角与晶界邻域的平均应力、应变、孪晶体分数的关联性分析
定义晶粒中晶界面的法线方向与加载轴的角度，称为晶界位向角。图16~18所示分别为结构模型中晶粒的晶界位向角与晶界邻域的应力、应变和孪晶体积分数的统计关系。图16所示为加载方向的压缩应力S33与晶界位向角的关系图。从图16可以看出，约90°的晶界位向角对应的应力在125~260 MPa之间分布，而约20°或160°晶界位向角对应的应力变化约在160~225 MPa之间分布，整个分布范围呈现椭圆形区域，即以90°的晶界位向角为中心，应力变化范围随晶界位向角向两端变化而变小。
图16 晶界邻域加载方向的应力S33与晶界位向角的关系
Fig. 16 Relationship between calculated stress in loading direction for average strain of 5.5% and misorientation between grain boundary plane normal and loading direction
图17 晶界邻域加载方向的应变L33与晶界位向角的关系
Fig. 17 Relationship between calculated strain in loading direction for average strain of 5.5% and misorientation between grain boundary plane normal and loading direction
图18 晶界邻域晶粒孪晶体分数与晶界位向角的关系
Fig. 18 Relationship between volume fraction of twinning in each grain for average strain of 5.5% and misorientation between grain boundary plane normal and loading direction
图17所示为加载方向的压缩应变L33与晶界位向角的关系图。从图17可以看出，所有晶界位向角对应的压缩应变变化约为0.025~0.075，平均压缩应变为0.055。与压缩应力的分布相似，应变分布对晶界位向角没有特殊偏好，因此整个分布范围呈现长方形区域。
图18所示为加载方向的孪晶体分数与晶界位向角的关系图。从图18可以看出，所有晶界位向角对应的孪晶体分数变化约为0.2~0.7，平均值约为0.5。与应变的分布相类似，孪晶体分数分布对晶界位向角没有特殊偏好，因此整个分布范围呈现近似的长方形区域。
5.6 晶界倾角与晶界邻域的平均应力、应变、孪晶体分数的关联性分析
定义两相邻晶粒的晶界倾角为两个相邻晶粒取向之间的角度，根据式(36)进行计算两相邻晶粒A、B之间的晶界倾角。
(36)
图19~21所示为微结构模型中晶粒的晶界倾角与晶界邻域的压缩应力、压缩应变和孪晶体积分数的统计关系。图19所示为加载方向的压缩应力S33与晶界倾角的关系图。从图19可以看出，晶界倾角在0°~90°之间变化，晶界倾角对应的压缩应力在125~260 MPa之间分布，整个分布范围呈现长方形。但压缩应力的分布对晶界倾角的大小没有出现特殊倾向。
图20所示为晶界倾角与加载方向的压缩应变L33的关系图。从图20可以看出，所有晶界倾角对应的压缩应变在0.03~0.07之间变化，平均应变为0.055。与晶界倾角对应的应变分布没有出现特殊倾向，因此，整个分布范围呈现长方形。
图21所示为晶界倾角与晶粒的孪晶体积分数的关系图。从图21可以看出，与图17的统计结果相似，所有晶界倾角对应的孪晶体分数变化在0.02~0.08之间分布，平均值为0.5。孪晶体积分数分布对晶界倾角没有特殊倾向，因此整个分布范围呈现近似的长方形区域。
综合以上分析，针对晶粒尺寸、晶向、晶界等微结构与应力、应变、孪晶体积分数的模拟结果进行了统计关联分析。分析结果表明该模型能够再现应力、应变在空间的不均匀分布，在统计上与实验结果基本吻合，但在统计微结构与孪晶体分数的分布关系上与试验结果仍存在不相符的情况。
图19 晶界邻域加载方向上应力S33与晶界倾角的关系
Fig. 19 Relationship between calculated stress in loading direction for average strain of 5.5% and grain boundary misorientation angles
图20 晶界邻域加载方向上应变L33与晶界倾角的关系
Fig. 20 Relationship between calculated strain in loading direction for average strain of 5.5% and grain boundary misorientation angles
图21 晶界邻域晶粒孪晶体分数与晶界倾角的关系
Fig. 21 Relationship between volume fraction of twinning in each grain for average strain of 5.5% and grain boundary misorientation angles
6 分析与讨论
前面提到的文献及本文所开展的模拟分析结果表明，从按体积平均的物理量演化来看，伪滑移模式在一定程度上能够有效描述由孪生变形引起的非均匀应力应变演化、应变硬化率演化、多晶织构演化、孪生与去孪生变形特征等。正如文中模拟的结果一样，宏观硬化行为与织构演化的模拟结果与试验比较吻合，然而，进一步地统计分析发现，微结构与孪晶分布之间缺乏确切的关联性，这与文献[12]中呈现的试验现象并不一致。究其原因，可以从孪生变形过程去分析。原子尺度上的研究表明[47]，稳定的孪晶形核由多重原子层构成。分子动力学研究[48]也表明，孪晶形核由孪晶偏位错合并而成，孪晶偏位错来源于大错配晶粒晶界，位错均匀排列在小角倾斜晶界处。孪生形核是位错堆积导致孪晶芽尖形成，在此基础上，孪晶薄片快速形成，且沿晶粒内部呈贯穿性增殖。随后孪晶进入逐渐增厚过程。因此，孪晶的形成过程可分为细小孪晶的形成，然后长大的过程。
在某些金属材料中，如镁、钛，孪晶变厚过程持续进行，直到整个晶粒取向发生改变，即晶向重定向[49-50]。从以上模拟结果看，伪滑移模式对孪生变厚过程的描述是合理有效的。然而，由于孪晶形核的动力机制特性，在孪生开端和增殖阶段描述孪生形核与周围环境间相互作用的关系仍然比较困难。伪滑移模式正是由于缺乏对孪晶形核的动态特性及形核位置邻域微结构演化的描述能力，结果导致难以对孪晶微结构特征进行识别和对离散孪晶带形态的非均匀空间分布进行预测，致使模拟统计结果与试验现象不一致。这也表明伪滑移模式缺乏描述镁合金多晶体变形细节的能力，且非均匀微结构所导致的独特应力非均匀状态并不是孪生变形与微结构强关联性的原因。显然，孪生变形所导致的晶格取向快速变向以及塑性流动的不均匀性，与微结构形态和结构特征之间存在错综复杂的关系，因此深入研究孪生形核机制及其与微结构之间的关系仍是一项很重要的工作。
7 结论
1) 文中构建的本构模型基本能够描述单调加载下镁合金材料的宏观硬化行为，并能够从细观的变形机制解释镁合金力学性能的各向异性特性。
2) 对于拉伸与压缩宏观硬化行为的差异，其原因是不同加载情形导致滑移与孪生变形的不同组合造成的。
3) 镁合金多晶织构的变化主要是孪生变形引起晶粒取向的大角度、快速转动造成的，而滑移变形引起的晶向变化是连续的、小角度转动的。
4) 针对晶粒尺寸、晶向、晶界等微结构与应力、应变、孪晶体分数的模拟结果进行了统计分析。分析结果表明，该模型能够展现多晶内的应力、应变在空间的不均匀分布特性，在统计上与试验结果基本吻合，因而伪滑移模式对孪生变厚长大过程的描述是合理有效的。
5) 统计分析表明，微结构与孪晶体分数的关联分布与试验现象存在不相符的情况，说明采用伪滑移模型描述孪生变形的方法仍存在不足，无法进一步描述更多的变形细节；这也表明非均匀微结构所导致的独特应力非均匀状态并不是造成孪生变形与微结构强关联性的原因。
REFERENCES
[1] AGNEW S R, YOO M H, TOME C N. Application of texture simulation to understanding mechanical behavior of Mg and solid solution alloys containing Li or Y[J]. Acta Mater, 2001, 49(20): 4277-4289.
[2] YOO M H, AGNEW S R, MORRIS J R, HO K M. Non-basal slip systems in HCP metals and alloys: Source mechanisms[J]. Mater Sci Eng, 2001, A319/321: 87-92.
[3] XIN Y C, WANG M Y, ZENG Z, NIE M G, LIU Q. Strengthening and toughening of magnesium alloy by 
extension twins[J]. Scripta Mater, 2012, 66(1): 25-28.
[4] LIU G D, XIN R L, LIU F Y, LIU Q. Twinning characteristic in tension of magnesium alloys and its effect on mechanical properties[J]. Materials and Design, 2016, 107: 503-510.
[5] ZHOU G W, JAIN M K, WU P D, SHAO Y C, LI D Y, PENG Y H. Experiment and crystal plasticity analysis on plastic deformation of AZ31B Mg alloy sheet under intermediate temperatures: How deformation mechanisms evolve[J]. Int J Plast, 2016, 79: 19-47.
[6] 吴国华, 肖 寒, 周慧子, 王瑞雪, 程 明, 张士宏. 挤压态AZ31 镁合金温热拉伸性能的各向异性[J]. 中国有色金属学报, 2017, 27(1): 57-63.
WU Guo-hua, XIAO Han, ZHOU Hui-zi, WANG Rui-xue, CHENG Ming, ZHANG Shi-hong. Anisotropy of warm- temperature tensile properties of extruded AZ31 magnesium alloy[J]. The Chinese Journal of Nonferrous Metals, 2017, 27(1): 57-63.
[7] PANDEY A, KABIRIAN F, HWANG J H, CHOI S H, KHAN A S. Mechanical responses and deformation mechanisms of an AZ31 Mg alloy sheet under dynamic and simple shear deformations[J]. Int J Plast, 2015, 68: 111-131.
[8] DONG S, YU Q, JIANG Y Y, DONG J, WANG F H, LI J, DING W J. Characteristic cyclic plastic deformation in ZK60 magnesium alloy[J]. Int J Plast, 2017, 91: 25-47.
[9] PARK S H, HONG S G, LEE C S. Activation mode dependent
twinning characteristics in a polycrystalline magnesium alloy[J]. Scr Mater, 2010, 62: 202-205.
[10] GODET S, JIANG L, LUO A A, JONAS J J. Use of Schmid factors to select extension twin variants in extruded magnesium alloy tubes[J]. Scr Mater, 2006, 55: 1055-1058.
[11] BARNETT M R, GHADERI A, QUINTA DA FONSECA J, ROBSON J D. Influence of orientation on twin nucleation and growth at low strains in a magnesium alloy[J]. Acta Mater, 2014, 80: 380-391.
[12] BEYERLEIN I J, TOME C N. A probabilistic twin nucleation model for HCP polycrystalline metals [J]. Proc R Soc A Math Phys Eng Sci, 2010, 466: 2517-2544
[13] KNEZEVIC M, BEYERLEIN I J, LOVATO M L, TOME C N, RICHARDS A W, MCCABE R J. A strain-rate and temperature dependent constitutive model for BCC metals incorporating non-Schmid effects: Application to tantalum- tungsten alloys[J]. Int J Plast, 2014, 62: 93-104.
[14] JONAS J J, MU S, AL-SAMMAN T, GOTTSTEIN G, JIANG L, MARTIN E. The role of strain accommodation during the variant selection of primary twins in magnesium[J]. Acta Mater, 2011, 59: 2046-2056.
[15] KNEZEVIC M, CAPOLUNGO L, TOME C N, LEBENSOHN R A, ALEXANDER D J, MIHAILA B, MCCABE R J. Anisotropic stress-strain response and microstructure evolution of textured a-uranium[J]. Acta Mater, 2012, 60: 702-715.
[16] GHADERI A, BARNETT M R. Sensitivity of deformation twinning to grain size in titanium and magnesium[J]. Acta Mater, 2011, 59: 7824-7839.
[17] YU Q, MISHRA R, MINOR A. The effect of size on the deformation twinning behavior in hexagonal close-packed Ti and Mg[J]. JOM, 2012, 64: 1235-1240.
[18] ABDOLVAND H, DAYMOND M R. Multi-scale modeling and experimental study of twin inception and propagation in hexagonal close-packed materials using a crystal plasticity finite element approach; part II: Local behavior[J]. J Mech Phys Solids, 2013, 61: 803-818.
[19] KNEZEVIC M, LEVINSON A, HARRIS R, MISHRA R K, DOHERTY R D, KALIDINDI S R. Deformation twinning in AZ31: Influence on strain hardening[J]. Acta Mater, 2010, 58: 6230-6242.
[20] NIEZGODA S R, KANJARLA A K, BEYERLEIN I J, TOME C N. Stochastic modeling of twin nucleation in polycrystals: An application in hexagonal close-packed metals[J]. Int J Plast, 2014, 56: 119-138.
[21] FERNANDEZ A, JERUSALEM A, GUTIERREZ- URRUTIA I, PEREZ-PRADO, M T. Three-dimensional investigation of grain boundary-twin interactions in a Mg AZ31 alloy by electron backscatter diffraction and continuum modeling[J]. Acta Mater, 2013, 61: 7679-7692.
[22] MAREAU C, DAYMOND M R. Micromechanical modeling of twinning in polycrystalline materials application to magnesium[J]. Int J Plast, 2016, 85(3): 803-818.
[23] JUAN P A, BERBENNI S, BARNETT M R, TOME C N, CAPOLUNGO L. A double inclusion homogenization scheme for polycrystals with hierarchal topologies: application to twinning in Mg alloys[J]. Int J Plast, 2014, 60: 182-196.
[24] MAREAU C, DAYMOND M R. Study of internal strain evolution in Zircaloy-2 using polycrystalline models: comparison between a rate-dependent and a rate-independent formulation[J]. Acta Mater, 2010, 58: 3313-3325.
[25] WARWICK J L W, JONES N G, RAHMAN K M, DYE D. Lattice strain evolution during tensile and compressive loading of CP Ti[J]. Acta Mater, 2012, 60: 6720-6731.
[26] XI B L, FANG G, XU S W. Multiscale mechanical behavior and microstructure evolution of extruded magnesium alloy sheets: Experimental and crystal plasticity analysis[J]. Materials Characterization, 2018, 135: 115-123.
[27] COGHE F, TIRRY W, RABET L, SCHRYVERS D, VAN HOUTTE P. Importance of twinning in static and dynamic compression of a Tie6Ale4V titanium alloy with an equiaxed microstructure[J]. Mater Sci Eng A, 2012, 537: 1-10.
[28] KNEZEVIC M, NIZOLEK T, ARDELJAN M, BEYERLEIN I J, MARA N A, POLLOCK T M. Texture evolution in two-phase Zr/Nb lamellar composites during accumulative roll bonding[J]. Int J Plast, 2014, 57: 16-28.
[29] ARDELJAN M, BEYERLEIN I J, MCWILLIAMS B A, KNEZEVIC M. Strain rate and temperature sensitive multi-level crystal plasticity model for large plastic deformation behavior: Application to AZ31 magnesium alloy[J]. Int J Plast, 2016, 83: 90-109.
[30] KNEZEVIC M, MCCABE R J, LEBENSOHN R A, TOME C N, LIU C, LOVATO M L, MIHAILA B. Integration of self-consistent polycrystal plasticity with dislocation density based hardening laws within an implicit finite element framework: Application to low-symmetry metals[J]. J Mech Phys Solids, 2013, 61: 2034-2046.
[31] KONDO R, TADANO Y, SHIZAWA K. A phase-field model of twinning and detwinning coupled with dislocation-based crystal plasticity for HCP metals[J]. Comput Mater Sci, 2014, 95: 672-683.
[32] WANG H, WU P D, WANG J, TOME C N. A crystal plasticity model for hexagonal close packed (HCP) crystals including twinning and de-twinning mechanisms[J]. Int J Plast, 2013, 49: 36-52.
[33] QIAO H, AGNEW S R, WU P D. Modeling twinning and detwinning behavior of Mg alloy ZK60A during monotonic and cyclic loading[J]. Int J Plast, 2015, 65: 61-68.
[34] BIELER T, WANG L, BEAUDOIN A, KENESEI P, LIENERT U. In situ characterization of twin nucleation in pure Ti using 3D-XRD[J]. Metall Mater Trans A, 2014, 45: 109-122.
[35] JUAN P A, PRADALIER C, BERBENNI S, MCCABE R J, TOME C N, CAPOLUNGO L. A statistical analysis of the influence of microstructure and twin-twin junctions on twin nucleation and twin growth in Zr[J]. Acta Mater, 2015, 95: 399-410.
[36] CHENG J H, GHOSH S. Crystal plasticity finite element modeling of discrete twin evolution in polycrystalline magnesium[J]. J Mech Phys Solids, 2017, 99: 512-538.
[37] KALIDINDI S R. Incorporation of deformation twinning in crystal plasticity models[J]. J Mech Phys Solids, 1998, 46(2): 267-290.
[38] HUTCHINSON J W. Bounds and self-consistent estimates for creep of polycrystalline materials[J]. Proc Roy Soc Lond, 1976, A348: 101-127.
[39] PAN J, RICE J R. Rate sensitively of plastic flow and implications for yield surface vertices[J]. Int J Solids Struct, 1983, 19: 973-987.
[40] LAN Yong-ting, ZHONG Xian-ci, QUAN Gao-feng, LIN Ruo-cheng, ZHANG Ke-shi. Crystal anisotropy of AZ31 magnesium alloy under uniaxial tension and compression[J]. Transactions of Nonferrous Metals Society of China, 2015, 25: 249-260.
[41] STEINMANN P, STEIN E. On the numerical treatment and analysis of finite deformation ductile single crystal plasticity[J]. Comput Methods Appl Mech Engrg, 1996, 129: 235-254.
[42] ZHANG K S, WU M S, FENG R. Simulation of microplasticity-induced deformation in uniaxially strained ceramics by 3-D Voronoi polycrystal modeling[J]. Int J Plast, 2005, 21: 801-834.
[43] 王自强, 段祝平. 细观塑性力学[M]. 北京: 科学出版社, 1995: 120.
WANG Zi-qiang, DUAN Zhu-ping. Meso-mechanics of plasticity[M]. Beijing: Science Press, 1995: 120.
[44] HUGHES T J R, WINGET J. Finite rotation effects in numerical integration of rate constitutive equations arising in large-deformation analysis[J]. Int J Numer Meth Eng, 1980, 15: 1862-1867.
[45] KABIRIAN F, KHAN A S, GNAUPEL-HERLOD T. Visco-plastic modeling of mechanical responses and texture evolution in extruded AZ31 magnesium alloy for various loading conditions[J]. Int J Plast, 2015, 68: 1-20.
[46] WANG L Y, HUANG Z H, WANG H M, MALDAR A, YI S B, PARK J S, KENESEI P, LILLEODDEN E, ZENG X Q. Study of slip activity in a Mg-Y alloy by in situ high energy X-ray diffraction microscopy and elastic viscoplastic self-consistent modeling[J]. Acta Mater, 2018, 155: 138-152.
[47] WANG J, HOAGLAND R G, HIRTH J P, CAPOLUNGO L, BEYERLEIN I J, TOME C N. Nucleation of a (1-012) twin in hexagonal close-packed crystals[J]. Scr Mater, 2009, 61(9): 903-906.
[48] WANG J, BEYERLEIN I J, TOME C N. An atomic and probabilistic perspective on twin nucleation in Mg[J]. Scr Mater, 2010, 63(7): 741-746.
[49] CHUN Y B, DAVIES C H J. Twinning-induced negative strain rate sensitivity in wrought Mg alloy AZ31[J]. Mater Sci Eng A, 2011, 528(18): 5713-5722.
[50] PRAKASH D G L, DING R, MOAT R J, JONES I, WITHERS P J, FONSECA J Q D. Deformation twinning in Ti-6Al-4V during low strain rate deformation to moderate strains at room temperature[J]. Mater Sci Eng A, 2010, 527(21/22): 5734-5744.
Constitutive model based on slip and twinning of AZ31 Mg alloy and analysis of microstructural relatedness
LAN Yong-ting1, CHEN Yuan3, REN Yi-fang2, ZHANG Ke-shi4, WANG Shuai2
(1. School of Vocational and Technical Education, Guangxi University of Science and Technology, Liuzhou 545006, China;
2. School of Mechanical Engineering, Guangxi University of Science and Technology, Liuzhou 545006, China;
3. School of Mechanical Engineering, Guangxi University, Nanning 530004, China;
4. College of Civil Engineering and Architecture, Guangxi University, Nanning 530004, China)
Abstract: In order to discover the relationship between the twinning deformation and microstructure in Mg alloy, describing twinning by pseudo-slip mode, the crystal plastic constitutive relation including slipping and twinning was established. And the Newton-Raphson iteration method was developed with the shear strain rate directly as the basic variables of iteration. The deformation behaviors of extruded Mg alloy loaded by uniaxial tension and compression along the extrusion direction were simulated using the finite element model of polycrystal based on voxel algorithm. The simulation results show that the macroscopic hardening behaviors predicted using the proposed model are in good agreement with the experimental results, and there are the same evolutions of polycrystalline texture during testing and simulating for Mg alloy. The statistical analyses of principal physical variables show that the proposed constitutive model is available to demonstrate the inhomogeneous characteristic of spatial distribution of stress-strain, the fraction of twinning in polycrystal, but there is a short age of the explicit relatedness of the volume fraction of twinning with respect to grain size, the misorientation between grain boundary plane normal and loading direction, the grain boundary misorientation angles in polycrystal, respectively. It is shown that the inhomogeneous state of unique stress resulted from the inhomogeneity of microstructure is not the reason for the strong relatedness between twinning deformation and microstructure.
Key words: magnesium alloy; twinning deformation; microstructure; plastic constitutive model; texture
Foundation item: Projects(11462002, 11862002) supported by the National Natural Science Foundation of China; Project(2016GXNSFAA380218) supported by the Guangxi Natural Science Foundation, China; Project(Z01) supported by the Science Foundation for Doctorate Research of Guangxi University of Science and Technology, China
Received date: 2018-07-12; Accepted date: 2018-10-16
Corresponding author: LAN Yong-ting; Tel: +86-772-2686126; E-mail: 100000877@gxust.edu.cn
(编辑 李艳红)
基金项目：国家自然科学基金资助项目(11462002，11862002)；广西自然科学基金资助项目(2016GXNSFAA380218)；广西科技大学博士科学基金资助项目(校科博Z01)
收稿日期：2018-07-12；修订日期：2018-10-16
通信作者：蓝永庭，副教授，博士；电话：0772-2686126；E-mail：100000877@gxust.edu.cn
