孪生诱发塑性(twinning induced plasticity, TWIP)钢不仅具有高抗拉强度和高硬化率, 而且具有优良的塑性、韧性和成形性能^[1], 为现代汽车发展轻量化和提高安全性等方面提供了新的选择, 引起国内外研究者的广泛关注^[2]. 一般来说, 经过变形后的多晶体金属材料通常出现各向异性行为, 这主要由于多晶体在塑性变形过程中形成变形织构^[3,4]. TWIP钢塑性变形是由位错滑移和孪生变形2种机制共同作用的结果, 当晶体内孪晶体积达到其饱和值时, 引起晶格发生转向, 转向后位错滑移在新的滑移面上继续进行. 因此, TWIP钢织构的形成是位错滑移和孪生共同作用产生的^[5]. 材料变形织构的形成与加载条件相关, 通过对不同变形条件下织构的预测, 控制材料的成形工艺以实现材料成形性能最优化.

研究人员^[2,6-8]采用现代显微测试技术和理论建模对TWIP钢进行了研究, 解释了其塑性变形是由孪生和滑移共同作用引起的; 并进一步分析了滑移和孪生的作用规律, 为定性分析滑移和孪生对宏观塑性变形的贡献提供了参考, 但没有给出各孪生系和滑移系的演化规律. 本文作者前期工作^[9-11]对单晶体塑性模型研究得到了TWIP钢单晶体塑性变形过程中各孪生系、滑移系的演化过程及孪生对位错滑移的影响, 但不能用于描述多晶体织构的演化. 此外, 研究人员以单晶体塑性模型为基础, 建立了Sachs模型^[12]、Taylor模型^[13]、自洽(self consistent, SC)模型^[14]和有限元模型^[15-19]等多晶体本构模型, 对多晶体进行了变形织构的实验和模拟研究. 目前, 由于Sachs模型假设多晶变形时各晶粒所承受的应力状态与宏观张量一致, 因此在各晶粒晶界之间产生的应变不连续, 与实验结果不太吻合; Taylor模型满足了晶界两侧的应变协调条件, 与实验结果吻合较好^[20,21], 但在晶界处的应力平衡却不能得到满足; 自洽模型在应用复杂本构关系时, 需要反复迭代计算, 应用范围也比较有限. 之后, 研究人员^[15-18]将晶体塑性理论与大型商用有限元软件结合, 在介观尺度上形成了一种新的织构模拟方法, 即晶体塑性有限元法(crystal plasticity finite element method, CPFEM), 该方法不仅能够模拟变形过程中织构的演化规律, 而且能够用于与织构相关的各种性能响应. 同时, 结合电子背散射衍射(electron back-scattered diffraction, EBSD)分析技术, 能够实时、快速和统计性地获取多晶体中各个晶粒的取向信息. 研究人员^[15,22,23]采用晶体塑性理论, 结合有限元模拟不同变形条件下的织构演化, 得到了与实验一致的良好结果, 但并未对孪生和滑移对织构演化的影响做出定量分析; Li等^[5]对TWIP钢滑移和孪生对应变硬化和织构演化的贡献做了定量分析, 认为在初始变形阶段主要依靠孪生变形, 随着应变的增加, 孪生对塑性变形的贡献很快被滑移取代, 最终孪生和滑移的贡献达到平衡, 但并未对不同加载条件对织构形成的影响进行深入分析.

本工作对退火处理后的Fe-22Mn-0.6C钢(质量分数, %)进行单轴加载实验, 获得拉伸条件下宏观力学响应及织构演化特点; 基于已建立的TWIP钢单晶体塑性本构模型^[9], 采用多晶均匀化处理建立联系单晶体与多晶体状态变量的多晶体塑性本构模型, 对该模型在ABAQUS/UMAT平台上进行二次开发并进行数值模拟, 结合实验确定模型中相关参数并验证模型可靠性; 进一步将上述本构模型应用于其他简单加载条件(压缩和扭转)下的应力应变关系及织构演化, 分析不同加载条件对应变硬化的影响及织构形成的贡献.

1 实验方法

本工作实验用TWIP钢多晶体材料的主要化学成分为Fe-22Mn-0.6C. 原材料采用Ar气氛保护, 电磁感应炉真空熔炼, 浇铸成板坯, 板坯经热轧至3.5 mm, 然后冷轧至1.5 mm, 在SRJX-8-13箱式电阻炉中进行退火处理, 温度为800 ℃, 保温时间75 min. 在二辊冷轧试验机上经过多道次轧制后厚度为1 mm, 再在四辊冷轧试验机上经过多道次轧制最终获得厚度为0.5 mm的TWIP钢钢板. 为了消除材料的初始轧制织构, 对TWIP钢钢板进行均匀化退火处理, 温度为800 ℃, 保温时间240 s^[24].

为了掌握TWIP钢多晶体在不同变形条件下的应力应变特点及织构演化规律, 对经退火处理后的TWIP钢钢板沿轧制方向取样, 试样在MTS810材料试验机上进行2组不同应变量的拉伸实验, 拉伸速率均为8 mm/min. 第一组拉伸后测得真应变为0.27; 第二组拉伸至试样断裂, 断后测得真应变为0.60. 将未拉伸试样和拉伸后的2组试样再次取样进行金相研磨并电解抛光, 采用EBSD技术观察织构状态.

2 多晶体塑性均匀化方法和数值模拟

2.1 多晶体塑性的均匀化方法

文献[9]中给出了单晶体塑性本构关系, 分析了单晶体塑性变形过程中滑移系和孪生系的演化规律, 并引入孪晶体积分数饱和值, 揭示了单晶体内孪晶体积的演化规律及孪生对位错滑移的影响, 滑移和孪生的硬化模量可用如下形式表征^[9]:

(1)sn+1,α=sn,α+hs,α1-sn,αss,α∑γ˙α∆t

式中, sn+1,α和 sn,α分别为第n+1步和第n步计算时滑移阻力, n为计算步数; γ˙α为第a滑移系的滑移速率; Dt为时间增量; hs,α和 ss,α分别为第a滑移系的硬化率和滑移阻力饱和值, 可以表示为:

(2)hs,α=hs1+c∑fβb

(3)ss,α=ss0+spr∑fβ0.5

式中, hs和 ss0分别为未发生孪生时的滑移系的硬化率和滑移阻力饱和值; spr, c, b都是材料硬化参数, fβ为第b孪生系的孪晶体积分数.

上述模型描述了TWIP钢单晶体塑性变形过程中各孪生和滑移系的应变硬化及其演化特点, 而实际金属材料多为多晶体结构, 因而需要通过多晶均匀化处理, 建立单晶体状态变量与多晶体状态变量间的关系, 使2个相邻的晶体在物理条件上满足几何协调和应力平衡条件.

Roters等^[25]基于连续介质理论提出了有限元均匀化模型, 假设一个单元代表一个晶体, 材料质点的应力响应由单个晶体的应力响应决定. 为了将文献[9]中的单晶变量参数与多晶变量参数进行区别, 在未进行特殊说明时, 将多晶变量参数上方加“—”.

根据连续介质大变形理论, 未变形时多晶体质点的参考构形 β-由材料的初始状态决定, 施加外载荷变形后多晶体的宏观塑性变形梯度^[2][5]为:

(4)F-=∇x-y-

式中, F-为宏观塑性变形梯度, 表示多晶体弹性和塑性变形梯度的总和; x-为多晶体材料质点矢径, 表示从坐标原点到物质点的矢量; y-为多晶体材料质点的非线性变形量, ∇x-y-为 y-对 x-求偏导.

F-可以分解为^[26]:

(5)F-=F-e⋅F-p

式中, F-e为多晶体弹性变形(晶格畸变)和刚性转动的合成变形梯度, F-p为多晶体滑移和孪生变形引起的塑性变形梯度.

多晶体弹性变形(晶格畸变)和刚性转动的合成变形梯度、滑移和孪生变形引起的塑性变形梯度均可表示为与所有晶体相对应变形梯度的体积平均值:

(6)F-e=1V∫βFedV

(7)F-p=1V∫βFpdV

(8)V=∫βdV

式中, Fe为单晶体弹性变形(晶格畸变)和刚性转动的合成变形梯度, Fp为单晶体滑移和孪生变形引起的塑性变形梯度, b为变形后多晶体的构形, V为多晶体体积.

由文献[9]可知, 多晶体第二Piola-Kirchhoff平均应力与 F-e的关系为:

(9)T-e=F-e-1⋅detF-eT-⋅F-e-T

式中, T-e为多晶体第二Piola-Kirchhoff平均应力, 表示所有晶体第二Piola-Kirchhoff应力的体积平均值; F-e-1为 F-e的逆矩阵; detF-e为 F-e的行列式; T-为多晶体Cauchy平均应力, 表示所有晶体Cauchy应力的体积平均值; F-e-T为 F-e转置后的逆矩阵.

而多晶体第一Piola-Kirchhoff平均应力与 T-e关系为:

(10)P-=F-e⋅T-e

(11)P-=1V∫βPdV

式中, P-为多晶体第一Piola-Kirchhoff平均应力, 表示所有晶体第一Piola-Kirchhoff应力的体积平均值; P为单晶体第一Piola-Kirchhoff应力.

上述均匀化模型中, 每个有限单元代表一个晶体, 单元的边界正好是晶体的边界. 而有限元中绝大多数单元是协调元, 在单元的边界上满足几何协调条件, 同时也在弱形式上满足应力平衡条件. 那么相邻晶体之间的几何协调和应力平衡条件就可以得到满足.

显示原图| 下载原图ZIP| 生成PPT

图1   多晶体塑性本构模型数值模拟的流程图

Fig.1   Schematic representation of polycrystal constitutive model

2.2 多晶体塑性的数值模拟过程

本工作中多晶体塑性的数值模拟过程如图1所示. 利用ABAQUS/UMAT将初始时刻所有晶体进行随机取向分布, 并定义所有晶体的材料参数和计算参数; 利用单晶体塑性本构模型计算每个晶体在变形过程中的状态变量, 同时根据多晶均匀化处理方法调节相邻晶体间应力平衡条件和几何协调关系; 计算完成后输出所需状态变量.

3 耦合孪生的多晶体塑性模型参数确定及实验验证

3.1 多晶体塑性模型的参数确定

结合文献[9,10]中TWIP钢单晶体塑性本构模型中的参数获取方法, 并取初始试探值经过不断修正可得到本构模型中的材料参数. 四阶弹性张量C中的3个材料常数分别取为C₁₁=198 GPa, C₁₂=125 GPa, C₄₄=122 GPa; 孪生体积分数阈值f₀取为0.2; 在准静态变形时参考塑性剪切率 γ˙0取为0.001 s^-1; 率敏感系数m=0.0074; 孪生和滑移阻力比值d=1.16, 其余参数见文献^[9].

3.2 多晶体拉伸过程的有限元模型

本工作通过ABAQUS/UMAT二次开发对TWIP钢多晶体塑性本构模型进行数值模拟, 并进一步建立TWIP钢多晶体拉伸过程的有限元模型, 以再现TWIP钢多晶试样单向拉伸条件下宏观应力应变关系及织构演化状态, 如图2所示. 该模型考虑实际拉伸过程中试样有效长度范围内加载、变形和初始织构分布特点, 可在其有效长度范围内取一特征段进行建模, 如图2a所示的红色区域, 特征区域尺寸为2.5 mm×0.25 mm×1 mm; 考虑到特征段的对称特点, 左端和后面分别关于xz面和zy面对称约束, 下表面z方向约束, 上表面给以z方向正向位移, 为保证计算的收敛性, 网格单元采用C3D8单元(三维8节点实体单元), 均匀划分单元数量为320, 如图2b所示; 所得某一晶体初始取向信息如图2c所示.

图3a为TWIP钢多晶体经退火处理后未变形时的极图. 可以看出, 经过退火后的TWIP钢在{100}, {110}和{111}面上没有出现明显的退火织构, 3个极图中的极密度比较均匀且较低, 晶粒取向为任意分布. Barbier等^[27]所研究的TWIP钢经过退火处理后仍有明显织构, 因此考虑了初始织构对TWIP钢拉伸变形过程中织构演化的影响. 本工作忽略初始退火织构对TWIP钢多晶材料影响, 认为经过退火处理后所有晶粒取向均匀分布. 结合初始条件下的实验结果, 对所有晶粒进行随机取向分布, 所有晶粒随机取向分布后织构如图3b所示. 可以看出, 晶粒随机取向分布后的TWIP钢在{100}, {110}和{111}面上没有出现明显的织构, 晶粒取向为任意分布, 与退火条件下的晶体织构有良好的一致性.

显示原图| 下载原图ZIP| 生成PPT

图2   孪生诱发塑性(TWIP)钢多晶体拉伸过程的有限元模型

Fig.2   Finite element model of macroscopic tensile process for twinning induced plasticity (TWIP) steel by using CPFE model (unit: mm)

显示原图| 下载原图ZIP| 生成PPT

图3   退火处理和晶粒取向随机分布后的TWIP钢在{100}, {110}和{111}面织构

Fig.3   Experimental textures after annealing treatments (a) and simulated textures after grains random orientation (b) in {100}, {110} and {111} planes for TWIP steel (RD—rolling direction; TD—transverse direction)

3.3 宏观应力应变和晶体取向的验证

本工作结合TWIP钢多晶体拉伸实验和所建立的拉伸过程有限元模型, 获得拉伸实验的应力、应变硬化率和孪晶体积分数随应变的演化与其模拟结果, 如图4所示. 可以看出, 在应变为0.27和0.60条件下模拟获得的应力应变曲线与实验获得的结果有较好的一致性; TWIP钢有较高的强度和较好的塑性, 抗拉强度达到1200 MPa, 断后延伸率达到80%, 随着应变的不断增加, 应变硬化增强, 应变硬化率逐渐降低; 研究人员^[28,29]采用TEM 观察了TWIP钢在不同拉伸变形量时的孪晶分布, 结果表明, 随着应变增加孪晶体积不断增加, 从图4模拟结果中也得到了相一致的结论, 在断裂时孪晶体积分数达到15%.

显示原图| 下载原图ZIP| 生成PPT

图4   实验和模拟获得的应力、应变硬化率及孪晶体积分数的演化对比

Fig.4   Simulated and experimental results of stress, strain hardening rate and volume fraction of twin (F_T —tension)

图5为真应变为0.27和0.6时的织构图和取向分布函数(orientation distribution function, ODF)图的实验和模拟结果. 可以看出, 随应变的增加TWIP钢多晶体逐渐形成织构, 且织构强度不断增强. 从ODF图实验和模拟结果可以看出, 单向拉伸后的TWIP钢多晶体主要是产生黄铜R型织构, 包括{111}<112>织构和{111}<110>织构. 其中, 应变为0.27时, 主要为{111}<112>织构, 其密度水平在4左右. 当应变为0.6时, {111}<112>织构增强, 同时产生新的{111}<110>织构, 其密度水平均在8左右. 同时随着拉伸应变的增加, 拉伸织构的锋锐程度增加. 通过对比结果可知, 该模型能较好地预测织构的演化, 具有一定可靠性, 可在实际的TWIP钢多晶材料塑性加工过程中进行应用.

4 不同加载条件下TWIP钢多晶体塑性变形行为

为了进一步分析不同加载条件下TWIP钢多晶体塑性变形特点及织构的形成和演化规律, 基于所建立的多晶体塑性有限元模型, 分别应用于TWIP钢多晶体压缩和扭转变形过程的数值模拟. 与上述的拉伸变形有限元模型相似, 所有网格单元采用C3D8单元, 压缩变形过程数值建模采用与上述拉伸变形模型相同的特征段, 尺寸均为2.5 mm×0.25 mm×1 mm, 均匀划分为320个有限单元, 位移边界条件与拉伸变形数值模型相同, 位移加载与拉伸变形方向相反, 通过加载不同位移量来控制变形量的大小. 扭转变形有限元模型采用直径1 mm×1 mm圆柱体模型, 自由划分为624个有限元网格, 在圆柱的一端进行端面固定, 在另一端面施加扭转变形, 通过施加不同扭转角度来控制变形量的大小.

图6给出了压缩和扭转变形过程中应力、硬化率和孪晶体积分数随应变演化的模拟结果. 可以看出, 压缩和扭转过程中随应变的增加应力和孪晶体积分数逐渐增加, 应变硬化率逐渐降低. 对比图4和6可以看出, 在拉伸和压缩变形过程中, 在相同变形量条件下, 由于材料的拉伸应变硬化模量小于压缩应变硬化模量, 致使压缩应力始终大于拉伸应力^[5].

图7a~d给出了沿轧制方向压缩应变为0.8和1.52时的{100}, {110}和{111}面上的极图和ODF图. 可以看出, 单向压缩后的TWIP钢主要产生黄铜和A织构, 包括{110}<112>织构和{110}<111>织构. 其中, 应变为0.8时, 主要为{110}<112>织构, 其密度水平在3左右. 当应变为1.52时, {110}<111>织构增强, 其密度水平均在3左右, 随着应变增加, 虽然织构类型有变化, 但是织构密度水平均较低.

显示原图| 下载原图ZIP| 生成PPT

图5   真应变为0.27和0.6时的织构图和ODF图的实验和模拟结果

Fig.5   Experimental (a, e) and simulated (c, g) textures for TWIP steel and corresponding ODFs (b, d, f, h) while strain e=0.27 (a~d) and e=0.6 (e~h) along rolling direction

显示原图| 下载原图ZIP| 生成PPT

图6   压缩和扭转过程应力、硬化率和孪晶体积分数随应变的演化

Fig.6   Evolution curves of stress, strain hardening rate and volume fraction of twin during compression and torsion processes (F_C —compression; T—torsion)

图7e~h给出了绕轧制方向所在轴扭转后的{100}, {110}和{111}面上的极图和ODF图. 可以看出, 在扭转量较小时, TWIP钢多晶体几乎没有织构形成, 主要因为圆柱体扭转变形时材料变形量沿径向方向逐渐增加, 变形较小时, 材料内部几乎没有变形, 多晶体内部接近初始状态, 故织构状态不明显; 随着变形量的增加, 多晶材料内部变形增加, 织构逐渐明显. 而在拉伸和压缩变形条件下, 所有晶体都将受到外力而进行变形, 故即使在小变形条件下, 多晶材料形成织构也较明显.

显示原图| 下载原图ZIP| 生成PPT

图7   TWIP钢压缩应变分别为0.8和1.52以及扭转应变分别为0.72和1.57时模拟的极图和ODF图

Fig.7   Simulated compression with e=0.8 (a) and e=1.52 (c), and torsion with e=0.72 (e) and e=1.57 (g) textures for TWIP steel along rolling direction, respectively, and corresponding ODFs (b, d, f, h)

5 结论

(1) 成功预测了拉伸过程中黄铜R型织构的生成, 同时发现随着拉伸应变量的增加, 拉伸织构的锋锐程度增加, 与实验结果有较好的一致性, 验证了模型的可靠性.

(2) 孪生诱发塑性钢多晶材料的拉伸应变硬化模量小于压缩应变硬化模量, 所以在相同变形条件下材料的压缩应力始终大于拉伸应力; 同时, 拉伸过程中, 随着应变增加, 主要织构类型发生变化, 平均密度水平增加; 压缩过程中, 随应变的增加, 虽然织构类型变化, 但是平均密度水平基本不变.

(3) 圆柱体扭转变形时材料变形量沿径向方向逐渐增加, 当扭转变形较小时, 由于材料内部变形较小, 织构形成不明显, 随着应变增加, 织构逐渐明显.