实验-CPFEM方法在分析先进高强钢应力应变配分中的应用

韩 东， 丁 桦， 赵文娟

(1.东北大学 材料科学与工程学院，沈阳 110819; 2.辽宁省轻量化用关键金属结构材料重点实验室，沈阳 110819; 3.湘潭大学 机械工程学院，湖南 湘潭 411105)

先进高强钢(Advanced high-strength steels, AHSS) 由于具有广泛的应用前景，近年来一直是各国材料研究的重点.目前所研究的第三代AHSS实验钢多为双相或复相组织，包含高强度相(超细晶铁素体、马氏体、贝氏体等)与具有较好塑性及应变硬化能力的相(奥氏体等).在外加载荷作用下，各相间应力应变配分在变形过程中的动态变化对于材料整体变形至关重要[1].因此，分析微观应力应变配分有助于更加深入地理解现有高强钢的动态变形行为，并为通过设计新成分、相组成、相界面等途径来进一步提高高强钢性能提供理论指导[2].

受实验手段限制，采用传统实验方法研究材料应力应变行为存在以下缺点：无法测量应力场，特别是微观应力分布；难以获得高应变条件下的应变分布；无法将材料应力应变同微观组织特征及演化相结合.早期的研究工作大多仅关注微观组织演变[3]或宏观应变场[4,5].微观数字图像相关技术(Microscopic digital image correlation, μDIC)的出现和发展，为微观应变场和微观组织演化的协同分析提供了有效手段.Martin[6]等将电子背散射衍射(Electron backscatter diffraction, EBSD)技术与μDIC技术相结合，可在获得微区应变场的同时对其组织进行分析，为分析材料的不均匀变形提供了有效方法.与此同时，基于晶体塑性理论的有限元方法也应用于模拟材料的应力应变行为.早期的工作大多基于形貌简化体胞模型[7-14].为了能更好地描述微观组织特征，获得更准确的结果，基于真实微观组织的晶体塑性有限元(Crystal plasticity finite element method, CPFEM)模拟近年来被广泛使用[11,15-23].随着原位实验技术[3,5,24-28]的发展和数值模拟理论[7,8,10,15,29]的不断完善，研究者通常配合使用μDIC技术与有限元方法以分析材料的变形行为.通过μDIC技术可获得试样变形过程中的应变分布及演变规律，可用于验证有限元模型的可靠性.同时，采用有限元方法可模拟出变形过程中的应力配分，这是传统实验方法难以企及的.本文对研究应力应变配分的传统实验方法与基于μDIC的原位实验方法进行了综述，并介绍了晶体塑性有限元方法的常用几何模型；归纳了应用于不同微观组织的几何模型建立方法；总结了先进高强钢变形中应力应变配分的最新研究进展.同时本文展望了晶体塑性有限元模拟与实验手段相结合的分析方法的发展方向与研究前景.

1 分析应变配分的实验方法

1.1 准原位实验

准原位实验通常采用EBSD与常规拉伸相结合.首先对拉伸试样表面进行抛光处理，采用聚焦离子束(Focused ion beam, FIB)、显微硬度或纳米压痕等方法标出感兴趣区域(Region of interest, ROI).在拉伸实验前，先通过EBSD获得原始组织的图像和取向信息.随后将试样拉伸至较小应变量，对试样上标定ROI进行EBSD 分析，获得该区域变形后的晶体学信息.随后再对试样进行进一步拉伸.重复此过程即可获得同一区域不同变形量下的晶粒、组织、织构等演化信息.此方法的优点在于可以获得近似原位的实验结果，EBSD对试样损伤较小，可提供晶粒尺寸、应变分布、织构演化等信息.

1.2 基于μDIC的原位实验

常用原位分析技术包括原位高能XRD[30]、原位中子衍射[31-33]和数字图像相关技术[27,34-36].原位XRD和原位中子衍射实验无法获得微观尺度特征(如晶粒尺寸、形貌、各相属性等)对协调应变的作用[34]，而μDIC技术可有效解决这一问题.数字图像相关技术是一种非接触式现代光学测量实验技术，其基本原理是识别试样表面离散分布的特征点，通过比较变形前后的数字散斑图像，获得区域的位移场，进而通过计算可获得应变分布图.

特征点的选取直接决定实验结果的准确性.常用的方法有以下两种：

1) 追踪实验材料微观组织特征.例如滑移迹线的背散射电子(Backscattered electrons, BSE)[37]和二次电子(Secondary electrons, SE)[38]图像，腐蚀出的晶界[39]、EBSD图像质量图(Image quality, IQ)[40]等.该方法的优点在于不需要额外的喷涂步骤.但由于分辨率受形貌影响，较难模拟晶内变形.

2) 人工添加特征点.可采用FIB[41]、重元素沉积法(Au、Pt等)[42,43]和紫外光刻法[44]等方法.其缺点为会对试样表面产生损伤或污染，甚至使组织发生改变，从而影响实验结果.Yan等[45]提出了使纳米尺度的SiO2颗粒附着在样品表面，并通过in-lens探头进行图像获取，降低形貌、形变等引起的衬度变化，减小SiO2颗粒对图像的影响(如图1所示).

图1 两种特征点获取方法示意图[45]Fig.1 Schematic representation of two pattern imaging methodology[45]

2 晶体塑性有限元方法

2.1 晶体塑性有限元理论发展

1938年，Taylor[48]在晶体学剪切变形的基础上提出了单晶塑性运动学方程和率无关本构关系.Rice[49]和Hill[50]随后对晶体塑性变形的几何学和运动学进行了深入的研究，并给出了相应的数学描述.Asaro[51]和Peirce[52]等进一步完善了晶体塑性理论，引入了自硬化和潜硬化的概念，推广了率相关模型.率无关与率相关两种模型的区别在于是否考虑应变速率对硬化的影响.在实际情况下，应变速率对变形后的组织性能有较大的影响，率相关模型更能反映材料变形的真实状况，因此得到了更为广泛的使用.

自上世纪40年代Courant[53]等提出有限元基本思想以来，有限元方法得到了巨大的发展，被广泛应用于各行各业.晶体塑性有限元法是以晶体塑性理论为基础，结合有限元方法形成的一种全新的分析材料变形的方法.晶体塑性有限元法的优势在于可以求解复杂内部或外部边界条件下的晶体力学行为，可灵活使用各种形式的本构方程[54]，可模拟不同尺度的变形.此外还可以在模型中加入更多的微观信息特征变量，如晶粒形状、取向等，从而更准确地模拟材料的变形，预测织构的演化.

2.2 常用几何模型

2.2.1 体胞模型

体胞模型[8-10,46,47]建立在多相材料的周期性假设基础上，将多相材料视为由周期性排列的体胞构成.图2为两种不同尺寸的球形马氏体颗粒和六边形立体结构的铁素体基体组成的体胞模型.通过对多相组织形态、含量、分布情况等进行分析，建立代表性体积元模型，对代表性体积单元(RVE)模型进行求解，可获得应力应变分布情况.该方法被用于分析马氏体体积分数和马氏体尺寸分布对宏观应力应变行为的影响[8-10]以及铁素体取向对DP钢取向梯度的影响[47].

图2 含两种不同尺寸颗粒的轴对称体胞模型[8]Fig 2 Axisymmetric two-particle model RVE idealization[8]

2.2.2 基于材料真实组织模型

图3 DP980微观组织与相应有限元模型[11]Fig.3 An actual microstructure and the corresponding finite element model of DP 980[11]

双相/复相钢的变形行为与各相形貌、分布情况等紧密相关.体胞模型难以反映真实材料中的微观组织信息.因此，近年来，基于微观组织的有限元模拟成为了研究双相/复相钢变形行为的有效方法.通常使用金相显微镜、扫描电子显微镜(Scanning electron microscope，SEM)以及EBSD技术获得实验钢的微观组织图像，通过对图像进行去噪、矢量化等处理，构建模型.如图3(a)所示，使用SEM获得了DP980实验钢的微观组织图像，图中黑色和白色区域分别为铁素体和马氏体.随后使用Photoshop等软件对图像进行处理，通过调整对比度等参数将马氏体与铁素体进行区分，并生成晶界.接着将图像导入ArcMap将其转换为矢量图，利用Gridgen软件对获得的矢量图进行网格划分，并对单元所属物相进行定义，最终结果如图3 (b)所示.

为模拟宏观试样的变形行为或分析不同相的平均应力应变，常通过微观组织图像建立代表性体积单元.通常代表性体积单元的建立基于以下假设：微观组织在所研究宏观区域内均匀分布，借助电子显微镜或金相显微镜获得单个或多个区域的图像作为代表性体积单元.Zhou[48]研究了DP980实验钢的应力应变行为，随机选取了20个尺寸为110 μm ×110 μm的区域进行模拟.结果显示，各区域模拟得出的宏观应力-应变行为并不相同.但是，若将这些区域随机分为两组(每组10个视场)，分别对每组结果取平均值，则两组结果极为相近.因此，对于组织不均匀或变形不均匀的材料，选取代表性体积单元时，应更为谨慎.选取数量较多或面积较大的区域进行模拟可以提高模拟结果的准确度.

目前对双相/复相钢的应力应变行为的模拟研究大多基于二维模型，对三维试样不同位置、不同截面的变形模拟较少.由于在不同位置所受载荷的方向不同，边界条件也不同，会对应力应变分布的模拟结果产生较大影响.Paul[49,50]使用2D RVE模型对DP590和DP780两种实验钢的力学行为和微观组织演化进行了研究.如图4所示，作者选取了三个不同位置进行模拟，探究不同载荷和边界条件的影响.其中A区域为平面应力状态，由于位于试样中心位置，其侧面在Y轴方向没有位移.B区域为平面应力状态，但由于接近试样表面，其侧面被认为没有约束.C区域为平面应变状态.模拟结果显示，在A区域中，由于侧面约束的存在，两种实验钢试样内均出现垂直延伸的高应变带.应变大多集中在铁素体中.B区域中，高应变带出现在最大剪应力方向上.C区域中，在加载初期即出现局部剪切集中的现象，且应变水平高于A、B两个区域.可见，二维模型难以全面反映出试样不同位置的应力应变情况.

图4 拉伸试样上选取不同位置微元[49]Fig.4 Microelements in three different locations of a sheet specimen under tensile loading[49]

2.3 本构模型

2.3.1 唯象晶体塑性本构模型

在晶体塑性理论中，流动法则反映了变形过程中材料的分切应变率同分切应力之间的关系.率相关流动法则采用幂指数形式描述如下：

(1)

实际晶体的硬化函数非常复杂，为了便于计算，将其简化为线性模型：

(2)

式中，hαβ为滑移的硬化模量，表示由第β个滑移系的塑性剪切引起的第α个滑移系的滑移阻力的增加，既包含了自硬化的影响，也包含了不同滑移系之间的潜硬化影响.

自上世纪30年代Taylor提出各向同性硬化假设以来，学者们提出众多的硬化模型.目前较为常见的为Asaro硬化模型[51]以及Bassani和Wu[52,53]提出的硬化模型.

Asaro硬化模型如下：

(3)

hαβ=qhαα,α≠β

(4)

上述模型的不足之处是无法描述潜硬化的影响.上世纪90年代，Bassani和Wu提出了新的硬化模型，通过对不同滑移系间的硬化使用不同的潜硬化系数，能更准确地反映实际硬化情况.其表达式为：

hαα=F(γα)G(γβ)

(5)

hαβ=qhαα,α≠β

(6)

(7)

(8)

式中，F(γα)和G(γβ)分别为滑移系单滑移的硬化模量和具有交互作用的潜在硬化模量，hs为硬化第一阶段(易滑移阶段)的硬化模量，fαβ为滑移系α与滑移系β间的相互作用系数.

2.3.2 基于位错的晶体塑性本构模型

(9)

(10)

(11)

(12)

(13)

式中：kB为玻尔兹曼常数，G为剪切模量.假设林位错切过过程由速率决定，可动位错速率可表示为：

(14)

式中：λα为移动宽度，与林位错间距成反比，νattack为相互作用频率，Qslip为位错滑移有效激活能，Vα为激活体积.有效剪切应力可表示为：

(15)

此模型中的硬化法则为位错密度的演化.文献中考虑了位错锁、位错偶极子、热湮灭和非热湮灭四种过程.其中前两种增加位错密度，后两种降低位错密度[54].

3 研究现状

3.1 DP钢

DP钢通常为铁素体-马氏体双相组织，马氏体呈颗粒状分布于铁素体基体上.马氏体硬度较高，不容易变形，因此应变主要集中于铁素体中.由于马氏体与铁素体晶格结构相似，因此无法通过EBSD将两者标定出.但由于马氏体与铁素体IQ(Image quality，衍射花样质量)相差较大，所以可以在IQ图中由不同衬度区分.

Tasan等[1]研究了DP钢的应变硬化行为，提出了实验与数值模拟相全方位结合的新方法.在原位变形过程中，通过μDIC技术，可以获得更高质量的应变图.同时，借助EBSD、电子通道衬度成像(Electron channeling contrast imaging, ECCI)、SE图像对相应的微观组织进行分析.通过EBSD图像建立CPFEM几何模型.铁素体属性由纳米压痕实验配合CPFEM模拟获得，对于晶粒大小与压痕接近的马氏体，其属性由拟合宏观应力应变曲线得到.针对该实验钢较高的相衬度和非线性应力-应变响应的特点，采用了新开发的基于FFT 的谱求解器.通过将模拟结果与实验应变场进行对比，可以获得局部应力场.图5为不同应变量下应力应变模拟结果.图5(a)显示，铁素体和马氏体存在显著的应变配分，铁素体承担了绝大多数的变形，且内部存在应变分布的不均匀，应变集中于与载荷方向成40～50°角的条带中.如图5(b)所示，应力集中于马氏体中.马氏体形貌和强度直接影响应力在马氏体区域分布的均匀性.变形过程中，马氏体较早达到应变峰值，且马氏体形状和与载荷方向所成角度对材料变形有着重要影响.由于所使用2D模型的局限性，无法描述马氏体在三维空间的分布情况.而实际材料中存在马氏体与铁素体互相重叠的情形，导致模拟结果与实际情况的偏差.

图5 不同应变量(左：0.03，中：0.05，右：0.08)的数值模拟结果[1]Fig.5 Numerical results obtained from the CP simulations at an average strain of ε=0.03 (left), ε=0.05 (center) and ε=0.08 (right)[1]

回火处理不但可以影响DP钢中马氏体的体积分数，也会影响两相性能.马氏体性能的变化对DP钢应力应变配分有较大影响.Han[57]对不同温度回火的DP双相钢进行了研究.模拟结果显示，随着铁素体与马氏体间强度差的增大，铁素体所承担的应变也增加，并且铁素体内应变分布更加不均匀，出现明显的高应变区和低应变区.当铁素体与马氏体属性相近时，应变配分现象更明显，但是铁素体内部应变分布更为均匀.

3.2 贝氏体钢

Fujita等[58]研究了含有马氏体-奥氏体的贝氏体钢.其基体为贝氏体，马氏体和残余奥氏体形成马奥岛.通过组织图像构建有限元模型时，EBSD可将奥氏体标定出.但由于马氏体与贝氏体具有相似的晶体结构，难以通过EBSD进行标定区分.因此作者通过置信度因子(Confidence index, CI)对其进行区分.置信度因子越高，表明晶体学取向标定更准确.马氏体和晶界处置信度因子较低，贝氏体置信度因子较高，可据此将各相识别出.后续模拟过程同DP钢相似.

部分研究者对铁素体-马氏体钢[1,45,59]和铁素体-贝氏体钢[46]的研究显示，对于含有较高体积分数第二相的双相钢而言，应变集中往往出现在粗大的岛状第二相附近.而Fujita等研究结果显示，在变形的早期阶段，应变集中出现在贝氏体晶界附近而非马奥岛附近.这是由于该贝氏体钢中马奥岛体积分数和尺寸都较小，且马氏体-铁素体在变形初期有着与贝氏体相近的硬化能力.随着应变程度的增加，粗大的马奥岛附近出现较高应力梯度，马奥岛晶粒尺寸和分布都影响着应变集中.细小均匀分布的马奥岛可以防止严重的局部应变集中，有助于提高贝氏体-马氏体-奥氏体复相钢的应变容纳能力.

3.3 双相不锈钢

双相不锈钢组织由体积分数近似相同的铁素体和奥氏体组成.Tao[60]研究了2205双相不锈钢的应力应变配分.实验结果表明，高应力区域大多出现于奥氏体中，高应变区域则大多位于铁素体中.这是由于相比奥氏体，铁素体具有较低的屈服强度，先于奥氏体发生变形.铁素体中最大应变出现在岛状奥氏体之间的锥形区域，其方向近似垂直于施加载荷方向.基于幂率硬化法则模拟的力学曲线在真应变小于0.15时与实验结果符合较好.随着真应变进增大，模拟应力值与实际应力值偏差增大，反映出所用幂率硬化法则不能很好吻合全局变形行为.

Jeong[33]使用CPFEM与EBSD相结合的手段研究了原始组织对双相不锈钢细观力学行为的影响.该实验采用了准原位拉伸-EBSD分析的方法获得试样在不同变形量下的应变分布，并基于真实微观组织建立了CPFEM模型.作者使用了唯象本构模型.双相的微观硬化参数由测得的宏观应力与沿载荷方向的晶格畸变拟合得到，应力-应变关系通过原位中子衍射与基于简化代表体积元模型的晶体塑性有限元模拟得到.Kernel取向平均差(Kernel average misorientation, KAM)图显示(图6)，在变形的各个阶段，奥氏体中的KAM值均大于铁素体KAM值.这主要是由于奥氏体具有较低的流变应力，且晶内存在高密度的大角晶界.晶体塑性有限元模拟较好地反映出了材料微观织构的演化以及双相的KAM分配情况.但是，模拟得出的铁素体中的KAM值低于实验结果.作者使用TEM对该差异进行了分析，认为其原因在于铁素体与奥氏体不同的变形机制：奥氏体发生平面滑移，形成层错；铁素体内形成位错胞.由于位错胞结构会在局部形成晶粒取向梯度，因此会导致模拟结果出现偏差.

图6 DSS钢中铁素体与奥氏体KAM图[33]Fig.6 Evolution of the KAM in the constituent phases of DSS during uniaxial tension simulation[33]

图7 基于Abaqus/Standard有限元实现流程图[62]Fig.7 The flow chart of numerical algorithm into implicit FEM[62]

图8 不同应变量von Mises 应变场有限元模拟结果[63]Fig.8 Results of microstructure-based FE simulations of plane-strain tensile deformation[63]

图9 不同区域的平均局部应力应变曲线和von Mises应变场[63]Fig.9 Averaged local stress-strain curves and von Mises strain fields obtained for plane-strain tension[63]

3.4 TWIP/TRIP钢

在TWIP钢和TRIP钢中，由于TWIP和TRIP现象的出现，使得材料具有优异的力学性能.为了模拟其变形行为，在经典晶体塑性理论中，除滑移外还需加入孪生和马氏体相变的描述.

Sun[61]和郭宁[62]提出了耦合TRIP和TWIP效应的晶体塑性有限元模型.采用率相关方法描述塑性流动，引入滑移和孪生阻力以及相应临界剪切应力之比作为屈服准则.用易于数值实现的方式建立了马氏体体积演化分数的率相关唯象描述，并以此作为相变流动法则.基于传统单晶模型，考虑较大塑性变形条件和TWIP效应的影响，建立了阶段式应变硬化拓展法则，并忽略马氏体中的塑性变形.采用Abaqus有限元软件，将耦合孪生和马氏体相变的晶体塑性本构模型以UMAT的方式嵌入，具体流程如图7所示.

Latypov[63]对Fe-6Mn-0.15C-1.5Si-3Al实验钢进行了研究.作者采用了中锰TRIP+TWIP钢的唯象本构模型，考虑了位错密度演化的硬化作用.孪生和相变的动力学分别用唯象法则和考虑孪生的Olson-Cohen模型来描述.图8(a) 和图8(b)分别模拟了试样在实际条件下(存在TWIP和TRIP效应)和忽略TWIP+TRIP机制条件下的变形行为.在TWIP和TRIP机制激活条件下，当应变为0.005时，奥氏体区域之间出现应力集中，呈条带状分布，与载荷方向夹角为45°.大部分铁素体未发生变形，仅分布在应力集中区域附近的铁素体发生变形.随着宏观应变量的增加，铁素体开始变形，应力分布趋于均匀.当应变增加到0.2时，条带状应力集中区再次出现，分布在铁素体区域内，并随应变量继续增加更加明显.此外，作者分析了不同应变量下的应变配分情况.如图9所示，变形初期，应变集中在铁素体内部的窄带，随着TWIP、TRIP效应的发生，高应变带转移至奥氏体-马氏体区域.该结果表明，TWIP和TRIP效应的激活使得部分奥氏体转变为马氏体，形成奥氏体-马氏体混合区域，使其流变应力超过铁素体，从而发生应力应变的转移.

4 尚存在的问题

准原位实验方法的主要缺点为： ① 在较高的应变量下，EBSD的标定率会大幅下降，影响结果的准确性； ② 使用FIB和压痕等方法对试样表面进行标记时，会损伤样品，并对标记周围晶粒变形产生一定影响； ③ 实验流程较为繁琐，试样易受污染，影响实验结果.目前普遍采用的μDIC结合原位试验的方法可以较好地获得材料的微观区域应变分布，避免准原位实验的缺点.但由于μDIC是根据特征点的位移计算应力场，因此其精度与特征点密度正相关.即使采用不同探头，采用追踪微观组织特征点的方法通常都难以获得较好的特征点衬度，特别是当晶粒尺寸较小时，在较高放大倍数时获得清晰的组织细节较为困难.相较之下人工添加的特征点可以获得更高的衬度，特别是在选用尺寸为纳米级的喷涂颗粒时，能对细晶组织进行分析.但是目前的方法无法精确控制附着颗粒密度.

采用晶体塑性有限元模拟双相或多相钢应力应变配分存在以下不足：一是现有模型未包含对材料强化机制及变形组织的描述.材料在变形时会出现不同的亚结构，如位错胞、位错墙、微带等，其演变过程直接影响材料的硬化行为.此外，部分材料存在固溶强化或析出强化机制，但现有模型并未引入相应变量[64].二是基于二维模型的模拟难以反映实际材料变形情况.目前3D-EBSD技术可以重构出材料的三维真实组织模型，但由于3D-EBSD采用FIB对样品表面进行连续切割，实验结束后研究区域完全遭到破坏，无法通过原位拉伸实验获得材料的实际变形行为.因此目前基于真实组织的模型大多为二维模型.虽然模拟结果能较好体现不同相之间的应力应变配分情况，但由于二维模型无法全面描述晶粒的形状、尺寸特点(例如板条状晶粒)，也未考虑实际三维空间中晶粒与周围晶粒之间的相对位置，因此模拟结果具有局限性.

此外，使用晶体塑性有限元模拟材料的应力应变配分行为需定义材料中各物相属性，但目前仍没有准确获得材料中不同相的属性的有效手段.由于材料属性受成分、热处理工艺等影响较大，难以通过查阅已有文献获得准确数据.部分学者通过纳米压痕实验获得相关参数[1,57,58]，但目前对纳米压痕硬度与屈服和抗拉强度之间换算关系的研究仍处于探索阶段.此外，载荷选取对纳米压痕硬度值有较大影响，进而影响结果的准确性.

5 结语与展望

基于微观组织的实验-晶体塑性有限元模拟方法虽存在上述不足之处，但成熟的晶体塑性理论、丰富的本构方程以及日益发展的微观组织观测手段，结合高效灵活的有限元方法，使得该方法在研究双相/复相钢应力应变配分上具有传统方法无法比拟的优越性.目前已经可以较好地模拟DP钢、双相不锈钢以及含有TRIP和TWIP效应的多相钢的微观变形行为.

未来应继续提高μDIC技术的分析精度.一方面，应继续开发和改进喷涂点颗粒制备及喷涂工艺，使颗粒尺寸及喷涂密度可控，从而满足不同尺度、不同精度的分析需求.另一方面，将材料本身组织形貌作为特征，可避免喷涂颗粒对结果的影响.因此，应结合表征手段的发展，提高获取图像质量及特征点衬度，优化算法，从而扩大该方法的应用范围.

此外，应继续探索获得材料准确属性的实验方法；完善晶体塑性有限元理论，建立基于真实组织的三维模型，完善材料本构方程，使其可以模拟不同成分、组织、变形机制的高强钢变形行为.特别是对于第三代AHSS钢中的高锰钢和中锰钢，不同的热处理工艺会使实验钢中出现多种碳化物与金属间化合物，显著影响力学性能.同时，随着合金元素含量的改变，中锰钢中可存在TRIP、TWIP及微带诱发塑性等不同变形机制.将其考虑在内有助于深入分析微观组织对变形行为和力学性能的影响，为合金设计和热处理工艺优化提供理论依据.