分享：基于位错密度的Fe-22Mn-0.6C型TWIP钢物理本构模型研究

孙朝阳, 黄杰, 郭宁, 杨竞

北京科技大学机械工程学院, 北京 100083

摘要

关键词： TWIP钢 ; 位错密度 ; 孪晶诱导塑性 ; 本构模型 ; 应变速率

TWIP (twinning induced plasticity)钢具有较高的抗拉强度(600~1100 MPa)以及相当大的延伸率(60%~95%), 其强塑积能达到50000 MPa·%[1], 成为重点研究的高强钢之一. 该钢种应用于汽车车身零部件可使车身框架的抗冲击性能得到提高[2], 同时在不牺牲力学性能的情况下达到减轻车重的目的[3]. 此外, 其优异的成形性能使常温下复杂零部件的拉深成形成为可能[2]. TWIP钢的高强高塑特性源于TWIP效应[4]: (1) 形变孪晶首先在变形最大且产生局部应力集中的部位形成, 孪晶界阻碍该区域位错的滑移而导致位错塞积, 使局部强度提高, 变形难以继续, 因此变形向相邻应变较低的区域转移, 从而推迟颈缩的形成, 大大提高了延伸率, 表现出高塑性的特性; (2) 形变孪晶在奥氏体晶界处向晶内贯穿, 起到了细化晶粒的作用, 降低位错滑移的平均自由程, 阻碍了位错的运动, 从而提高了应变硬化率, 表现出高强度的特性. 目前, 实验室已开发出多种不同成分的TWIP钢[5], 并采用不同微观实验方法对TWIP钢塑性变形机制开展了较为深入的研究, 已定性地分析了孪生与滑移共同作用引起TWIP钢高强高韧的本质[6-10]. 为了更好地掌握TWIP钢塑性变形过程各机制的相关作用及其对宏观力学性能的贡献, 涉及孪生和滑移机制的物理关系模型研究成为近年来研究的热点[11-13].

但是, 单一方程式的本构关系如Johnson-Cook经验型的本构模型[14]或Zerilli-Armstrong物理模型[15]不能满足要求, 虽然它们能模拟材料发生孪生后的应变硬化曲线, 但却无法记录孪晶体积分数和位错密度的演化规律, 也不能说明孪晶对位错滑移的影响. 因而, 为了描述TWIP钢滑移及孪生机制的相互作用, Bouaziz等[5,11]采用基于位错动力学的一系列微分方程建立了TWIP钢的加工硬化模型, 认为孪晶的产生降低了位错运动的自由程, 从而提高了硬化率, 并研究了晶粒尺寸效应及Bauschinger效应对孪生诱导塑性的影响. 之后, Kim等[13]进一步引入动态应变时效(dynamic strain aging)建立了Al-TWIP钢的本构模型, 发现应变硬化取决于滑移和孪生, 孪晶为位错运动的有效壁垒. 根据上述2种模型均能分析孪生体积分数和位错密度演化之间的关系, 且模拟的应力-应变曲线与实验结果也符合较好. Bouaziz等[5,11]和Kim等[13]均假设孪晶区内不再发生滑移, 但研究[16-18]表明孪晶内的滑移是塑性变形的一种重要模式, 对材料整体塑性变形产生贡献, 因而在描述整体塑性变形时考虑孪晶内滑移的贡献是有必要的. 此外, TWIP钢冲压成形性能及其它力学性能受应变速率影响[4], 而这直接关系到其在汽车结构件中的生产和使用安全性[19]. 所以建立应变速率相关的本构模型, 更具有实际意义和应用价值.

本工作以Fe-22Mn-0.6C型TWIP钢为例, 基于位错密度的不断增殖及回复和孪生体积分数的演化, 考虑孪晶内部的滑移对整体塑性变形的贡献, 建立TWIP钢应变速率相关的应变硬化模型. 根据孪生区和基体区Taylor因子的不同[16], 提出塑性变形增量新的描述方法. 结合实验确定模型中相关参数, 分析不同应变速率加载下的TWIP钢力学行为, 重点分析应变速率、孪晶体积分数演化对TWIP钢应变硬化的影响.

1 基于位错密度的物理模型

1.1 特征微元的应变分解

当材料中存在形变孪晶时, 塑性变形由孪生和滑移2种机制共同作用[13]. 在塑性变形过程中, 孪晶平行分布于每个晶粒中, 但不是所有孪晶都贯穿整个晶粒[20]. 如果考虑到孪晶本身对整体变形的贡献, 同时假设孪晶内部不再发生滑移, 剪切变形增量dgp可采用图1所示的滑移和孪生2种机制作用的微元体描述[21]:

式中, F为孪晶体积分数, gp为整体塑性剪切变形, γs为滑移引起的塑性剪切变形, gt为孪生引起的塑性剪切变形, gt= $\begin{array}{c}\sqrt{2}/2\end{array}$, M为平均Taylor因子, e为应变.

联立式(1)和(2)得到[21]:

式中, εp为整体塑性变形, εs为滑移引起的塑性变形, εt为孪生引起的塑性变形.

余勇等[18]研究V-5Cr-5Ti的应变硬化行为时认为, 孪晶内的位错仍然对材料塑性变形产生贡献. 一方面孪晶可作为亚晶粒处理, 其中的位错滑移也引起塑性变形; 另一方面, 当孪晶形成后, 孪晶的取向更有助于位错滑移而引起塑性变形. 因此材料的塑性变形应表示为:

图1   晶粒中孪晶示意图

Fig.1   Schematic morphology of twins in the representative element

相对于式(3), 式(4)考虑了孪晶内的滑移Fdεs对整体塑性变形的贡献.

Salem等[16,17]研究钛合金时发现, 在压缩实验中孪晶区域的平均Taylor因子比基体区域小40%左右, 一方面说明孪晶内部更容易发生滑移, 另一方面, 在考虑孪晶内滑移对整体塑性变形贡献的同时还需要考虑孪晶区域和基体区域不应采用相同的Taylor因子对变形的影响. 故式(2)中的整体塑性变形可分解为基体区的滑移、孪晶区的孪生及孪晶内的滑移3部分:

式中, Ms为基体区Taylor因子, Mt为孪晶区Taylor因子.

因为孪晶体积分数较小, 基体区Taylor因子简化为Ms≈M, 而孪晶区Taylor因子根据Salem等[16,17]的实验结果, 可近似取Mt≈0.6M. 则式(5)可化简为:

此分解考虑了孪晶内的滑移对整体塑性变形的贡献, 同时引入孪晶区与基体区具有不同的Taylor因子的概念, 更准确地描述了孪晶引起的塑性变形及孪晶内的滑移引起的塑性变形对整体塑性变形的影响.

下面将会根据计算结果进一步分析上述3种塑性变形的分解方式对模型计算的影响.

1.2 流变应力分析

基于位错理论, 一般可将流变应力σ分解为热激活应力s0和非热应力sa[22]:

式中, 热激活应力s0表征位错密度为0时的流变应力, 只与温度和应变速率有关. 非热应力sa表征位错交互作用对位错滑移的阻力.

(1) 非热应力

非热应力主要由位错产生的长程应力场所决定, 表示位错相互之间的阻力, 可用位错密度表示为[11]:

式中, a表征位错间的相互作用强度, m为剪切模量, b为Burgers矢量模, r为位错密度. 因为存在位错累积和动态回复之间的竞争, 位错密度随滑移引起的塑性变形es的演化规律表示为[23,24]:

式中, l为位错的平均自由程, f为动态回复参数, 表征动态回复的几率.

如果把孪晶看作位错不可逾越的障碍, 平均自由程在受晶粒尺寸和位错密度影响的同时, 还受孪晶平均间距影响, 可表示为[11]:

式中, d为晶粒尺寸, t为孪晶平均间距, k为一个待拟合的常数, 与加工硬化极限有关. 位错密度越大, 平均自由程越小.

而孪晶平均间距t可表示为[11]:

式中, e为孪晶平均厚度, 其与应变大小无关; 上式中的数字“2”表示一个修正因子[13]. 如图1所示, 随着孪晶体积分数的增多, 孪晶平均间距减小, 平均自由程也跟着减小, 最终影响位错密度的演化.

除了位错密度的演化规律, 还需要描述孪晶体积分数随塑性变形的演化规律, 一般采用经验方程[5]:

式中, einit为形变孪晶产生时的临界应变值, F0为孪晶体积分数的饱和值, b和m为材料常数.

(2) 热激活应力

余勇等[18]]认为恒温下的热激活应力s0是 $\begin{array}{c}ln{{\displaystyle \dot{\epsilon }}}_p\end{array}$的线性函数:

式中, A0和B0为与温度相关的常数, $\begin{array}{c}{{\displaystyle \dot{\epsilon }}}_p\end{array}$为塑性应变速率.

而材料在屈服时塑性应变ep≈0.002, 根据式(12)可知, 其小于设定的临界值einit=0.05, 表明尚未发生孪生, 塑性变形完全由滑移引起, 即dep=des. 根据式(7)和(8), 应力增量可表示为:

图2   Fe-22Mn-0.6C钢在2种应变速率下的真应力-真应变曲线

Fig.2   True stress-true strain curves of Fe-22Mn-0.6C steel tested at room temperature with strain rates of 0.01 and 0.1 s-1

式中, r0为材料屈服前的位错密度. 一般r0为1010 m-2量级[12].

通过Taylor展开, 可做以下简化:

综合式(9), (14)和(15)可得:

同时, 考虑到应力s与应变e的关系:

式中, E为弹性模量.

联立式(16)和(17)得到:

对式(18)两边求时间微商, 即得塑性应变速率 $\begin{array}{c}{{\displaystyle \dot{\epsilon }}}_p\end{array}$与应变速率 $\begin{array}{c}{\displaystyle \dot{\epsilon }}\end{array}$之间的关系:

根据式(7), 屈服应力sy可表示为热激活应力与屈服前的非热应力之和:

联立式(13), (19)和(20), 可得到屈服应力sy与 $\begin{array}{c}ln{\displaystyle \dot{\epsilon }}\end{array}$之间的关系:

图3   Fe-22Mn-0.6C钢在2种应变速率下的应变硬化率

Fig.3   Variation of strain hardening rate of Fe-22Mn-0.6C steel as a function of strain

式中, A和B为与温度相关的常数.

2 实验方法与结果

实验用TWIP钢的化学成分(质量分数, %)为Fe-22Mn-0.6C, 采用电磁感应炉真空熔炼, Ar气气氛保护, 浇铸成板坯, 板坯经热轧至3.5 mm然后冷轧为1.5 mm, 在SRJX-8-13箱式电阻炉中进行退火处理, 退火温度为800 ℃, 保温时间75 min. 然后在二辊冷轧试验机上经过多道次轧制后厚度轧为1 mm, 由于材料的加工硬化, 随着变形量逐渐增大, 二辊轧机的压下强度已不能继续压下, 最后将钢板转移至四辊冷轧试验机经过多道次轧制后最终获宽度约为130 mm, 厚度为0.5 mm的TWIP钢板. 按照国标GB/T 228-2010, 加工拉伸试样. 在AG-X拉伸试验机上沿轧制方向进行拉伸实验. 实验温度为室温(20 ℃), 应变速率为0.001~0.1 s-1.

图2为实验获得Fe-22Mn-0.6C钢在应变速率为0.01和0.1 s-1的真应力-真应变曲线. 可以看出, 应变速率为0.01 s-1时, 屈服应力为120 MPa, 而发生颈缩时的应力达到1150 MPa, 真应变达到0.77, 强塑积为88550 MPa·%. 应变速率为0.1 s-1时, 屈服应力为200 MPa, 发生颈缩时的应力达到1000 MPa, 真应变达到0.65, 强塑积为65000 MPa·%. 随着应变速率的提高, 屈服应力提高80 MPa, 抗拉强度降低150 MPa, 其强塑积降低23550 MPa·%.

图3为Fe-22Mn-0.6C钢的应变硬化率随着真应变的变化曲线. 可以看出, 2种速率下的曲线没有比较明显的差别. 计算应变硬化率时, 先光滑应力应变曲线, 每10组应力采集点求平均值, 然后利用Origin软件计算其硬化率. 由图看到, 这2条曲线可分为3个阶段. 在阶段A, 硬化率急剧下降, 此时未发生孪生, 位错滑移是主要的变形机制, 其硬化率完全取决于位错的增殖及动态回复. 动态回复时, 应变硬化率持续下降, 这与Mecking和Kocks[24] 的研究成果一致. 阶段B约从e=0.05开始, 硬化率开始逐步提高, 但提高幅度不大, Gutierrez-Urrutia和Raabe[9]认为这是因为形变孪晶的形成降低了位错间的平均自由程, 促进了位错密度的增殖, 最终提高了应变硬化率. 进入阶段 C后, 应变硬化率再次下降, 直至颈缩现象发生, 这归因于孪晶体积分数达到饱和值[7], 所以孪生体积分数的增量很小.

3 模型计算与结果验证

3.1 数值计算及模型参数确定

上述考虑孪晶和滑移耦合的TWIP钢的塑性变形的物理本构模型为一组微分方程组, 无法直接求解. 为此, 本工作提出该物理模型数值计算过程. 首先根据屈服应力与应变速率之间的关系(式(21)),得到屈服应力值, 然后由式(19)计算塑性应变速率, 并得到整体塑性变形增量及塑性变形值. 对比孪生临界变形值及整体塑性变形的大小, 判断是否发生孪生. 若发生孪生, 则采用分解式(6), 得到由滑移引起的塑性变形值; 若未发生孪生, 则塑性变形完全由滑移引起, Δes=Δep. 然后根据位错密度的演化规律式(9), 得到位错密度增量. 最终由式(14)得到其流变应力增量, 并更新当前流变应力值. 时间增量步设置为10-5 s. 根据上述数值计算过程可得到具体流程如图4所示.

图4   物理本构模型的计算流程图

Fig.4   Schematic representation of physical constitutive model (t1—step time, einit—critical strain when twinning begins)

采用应变速率0.001, 0.01, 0.02, 0.04和0.1 s-1, 获得不同应变速率下的屈服应力值sy和 $\begin{array}{c}ln{\displaystyle \dot{\epsilon }}\end{array}$的关系, 如图5所示, 通过线性拟合得到其方程为:

模拟计算中涉及各类材料常数及拟合参数如表1所示, 其中硬化参数k, 动态回复参数f及参数a, b, m均通过拟合应力应变曲线获取, 而孪晶体积分数饱和值F0参考文献[13]中实验取为0.2, 孪生临界应变值einit表征开始发生孪生时的塑性变形值, 而阶段B从e=0.05开始, 硬化率开始逐步提高, 认为此时开始发生孪生, 故设定einit=0.05, 其余参数如剪切模量、弹性模量、Taylor因子及孪晶平均厚度等均为材料常数由实验检测获得, 参考文献[8,12]中对Fe-22Mn-0.6C系TWIP钢的研究工作.

图5   Fe-22Mn-0.6C钢室温下的屈服应力sy与应变速率lnε? 的关系

Fig.5   Yield stress sy of Fe-22Mn-0.6C steel at room temperature as a function of lnε? ( ε? —strain rate)

3.2 模型计算结果验证

2种应变速率下Fe-22Mn-0.6C钢实验及模拟真应力-真应变曲线对比结果如图6所示. 可以看出, 其计算和实验结果吻合非常好, 说明建立的屈服应力与应变速率的关系是合理的. 图7为不同应变速率下Fe-22Mn-0.6C钢的应变硬化率计算与实验结果的对比. 因为模拟计算的真应力-真应变曲线较光滑, 没有震荡现象, 故可对计算出的真应力-真应变曲线直接求导得到应变硬化率的模拟值. 不同应变速率下的应变硬化率计算与实验结果的趋势基本一致, 呈现出3个明显阶段, 很好地反映了孪晶在阶段B时的强化作用以及在阶段C的软化作用. 此外, 不同应变速率下的应变硬化率基本一致, 说明应变速率对应变硬化的影响较小.

前文中, 式(3), (4)和(6)分别对应3种分解形式, 图8为不同分解下的真应力-真应变曲线. 3种分解在前期未发生孪生时是等价的, 故前期计算结果曲线完全一致, 孪生体积分数饱和后, 3种曲线开始分化. 在e≈0.2时, 已经发生孪生, 曲线开始分开, 如在e≈0.65时, 不同分解下的应力值s (式(3))>s (式(6))>s (式(4)), 最大可相差100 MPa左右.

图6   Fe-22Mn-0.6C钢在2种应变速率下的实验及模拟真应力-真应变曲线对比

Fig.6   Comparison of simulated and experimental true stress-true strain curves for Fe-22Mn-0.6C steel with strain rates of 0.1 and 0.01 s-1

为了描述建立本构关系的精确程度, 引入平均相对误差d的概念[25]:

图7   Fe-22Mn-0.6C钢在2种应变速率时的实验及模拟应变硬化率曲线对比

Fig.7   Comparison of simulated and experimental strain hardening rates for Fe-22Mn-0.6C steel with strain rates of 0.01 s-1 (a) and 0.1 s-1 (b)

式中, Ei为实验流变应力值, Pi为计算应力值, N为实验数据采集点i的数目(本工作为70). 通过计算得到3种分解方式下的平均相对误差值, 如表2所示. 可以看出, 采用本工作的模型总平均误差值只有0.84%, 比未考虑孪晶内滑移贡献的模型[11]和未考虑孪晶区和基体区具有不同Taylor因子的模型[18]明显要小, 分别减小1.1%和2.9%, 表明采用式(6)分解的本构模型具有较高精度, 能很好地满足工程要求.

进一步分析发现, 3个模型在未孪生及孪生很小时, 应力值基本一致, 如表2所示, 0<e≤0.2时, 平均相对误差基本相当. 只有当大量孪生时, 即e >0.2时, 孪生开始大幅度增加, 在e =0.4时体积分数趋近饱和值, 随着孪生体积分数的增加, 未考虑孪晶内滑移的模型误差先减小后增大, 考虑孪晶内滑移的模型误差持续增大, 而本工作模型的误差却持续减小, 由此可推断, 经修正后的模型能更好地描述孪生后的TWIP应变硬化行为.

图8   应变速率为0.1 s-1时Fe-22Mn-0.6C钢的3种变形分解下的模拟真应力-真应变曲线

Fig.8   True stress-true strain curves of Fe-22Mn-0.6C steel subjected to three decompositions with strain rate of 0.1 s-1

图9   应变速率为0.1 s-1时Fe-22Mn-0.6C钢的塑性应变速率、滑移率和孪生率

Fig.9   Evolution of plastic strain rate, slip rate and twinning rate for Fe-22Mn-0.6C steel with strain rate of 0.1 s-1

3 讨论

研究[19,20]表明, TWIP钢变形时除了产生形变孪晶, 还可能发生马氏体相变, 甚至TWIP/TRIP (transformation induced plasticity)效应共存. 其塑性变形机制及力学性能都与层错能有关[26], 对于Fe-Mn-C系奥氏体钢, 层错能Γ<18 mJ/m2时, 发生马氏体相变; 12<Γ<35 mJ/m2时, 形变孪晶对变形起主导作用. 本工作研究的TWIP钢Fe-22Mn-0.6C在293 K 变形时, 仅有孪生, 而无马氏体相变的产生[27]. 该物理本构模型主要探讨孪生和滑移机制的相互作用, 及其对材料力学行为的影响. 因此, 当TWIP钢材料在变形过程中发生大量马氏体相变时, 本模型并不适用.

为了说明模型中孪生对滑移的影响, 根据计算结果作了不同分解下塑性应变率 $\begin{array}{c}{{\displaystyle \dot{\epsilon }}}_p\end{array}$, 孪生体积分数变化率、滑移引起的应变速率及位错密度随真应变的变化曲线, 结果如图9和10所示. 由图可见, 孪生速率和滑移速率负相关, 孪生速率增大时, 滑移速率减小, 而孪生趋于饱和时, 孪生速率逐渐减小为零, 滑移速率迅速增加, 该模型很好地反映了孪晶与滑移机制的相互作用. 而3种分解方式下, 孪生速率和孪生体积分数是一样的, 因为根据式(12), 其值主要取决于整体塑性变形的大小, 而与整体塑性变形的分解无关. 由于发生塑性变形后, 弹性应变速率远小于塑性应变速率, 故塑性变形速率只有在弹性区域和塑性变形前期间发生了变化, 而之后较稳定, 只略小于总应变速率. 因为前期未发生孪生, 3种分解完全等效, 故前期滑移引起的位错密度及塑性应变速率曲线完全重合, 塑性变形只受滑移一种机制控制, 当应变达到einit=0.05时, 开始发生孪生, 孪生率开始增加, 滑移引起的应变速率des/dt开始下降, 在e=0.3左右时, 孪生率达到最大值, 而des/dt达到最小值, 此时变形机制以孪生为主, 孪生本身提供较大的变形率, 同时滑移较难发生, 位错密度增殖较缓慢, 孪生对滑移产生明显的阻碍作用. 当材料进一步变形后, 孪生体积分数接近饱和值, 孪生速率迅速下降, 晶粒中已很难产生新的形变孪晶阻碍滑移, 故此时塑性变形以滑移为主, 滑移速率趋近饱和值.

图10   应变速率为0.1 s-1时Fe-22Mn-0.6C钢的不同变形分解下的位错密度及孪晶体积分数

Fig.10   Evolution of dislocation density and volume fraction of twins for Fe-22Mn-0.6C steel subjected to three decompositions with strain rate of 0.1 s-1

4 结论

(1) 建立了关于应变速率的TWIP钢物理本构模型, 模型拟合的结果与实验吻合较好, 平均误差只有0.84%.

(2) 模型考虑了孪晶内部的滑移对整体塑性变形的贡献, 同时考虑了孪晶区域Taylor因子不同, 对塑性应变作了新的分解. 3种不同分解在未发生孪生及孪生初期时完全相同, 只有在塑性变形较大时, 孪晶基本饱和后才有所区别.

(3) 通过建立应变速率与材料屈服应力的关系将应变速率的影响纳入到模型中. 计算的屈服应力结果与实验数据吻合较好, 同时还发现应变速率只影响材料的屈服应力, 而对后续硬化率的影响不明显.

(4) 模型很好地反映了孪晶与滑移机制的相互作用, 孪生速率与滑移速率之间负相关, 塑性变形前期, 孪生速率增大而滑移速率减小, 孪生趋于饱和时, 孪生速率降低而滑移速率迅速增加.