耦合孪生晶体塑性模型硬化参数的灵敏度分析

工程科学学报，第37卷，第8期：1076-1083,2015年8月 Chinese Journal of Engineering,Vol.37,No.8:1076-1083,August 2015 D0l:10.13374/j.issn2095-9389.2015.08.016:http://journals..ustb.edu.cn 耦合孪生晶体塑性模型硬化参数的灵敏度分析 杨 竞”，黄杰12》，孙朝阳)四，郭宁”，王斌” 1)北京科技大学机械工程学院，北京1000832)宝山钢铁股份有限公司，上海201900 ☒通信作者，Email:suney@usth.cdu.cn 摘要基于经典晶体塑性理论，建立了耦合孪生的晶体塑性本构模型并进行了全隐式积分的数值实现.该本构模型采用 饱和硬化法则，并采用孪生阻力与滑移硬化之间的正比关系来描述孪生对滑移硬化影响及孪生硬化行为，针对该本构模型 的13个参数，结合各参数物理意义提出了参数的分类确定方法.以孪生诱导塑性(TWIP)钢Fe-22M0.6C为例，着重对硬 化参数的局部灵敏度进行了分析，研究了各硬化参数对宏观力学响应、孪生激活和演化的影响，根据变形机制的不同宏观变 形过程可区分为孪生硬化阶段和孪生硬化失效阶段，进而给出了硬化参数确定的步骤及其建议取值范围.结果表明：初始滑 移阻力与屈服极限线性相关，取值范围在80~160MP之间：孪生硬化指数增大使得孪生硬化阶段减弱，其取值范围应在0~ 3之间：李生阻力与滑移阻力比值增大，则孪生增长率降低，硬化率拐点后移，直至拐点消失，其取值范围在1~1.3之间 关键词孪生：晶体塑性：本构模型：灵敏度分析：硬化参数 分类号TG302 Sensitivity analysis of hardening parameters in the crystal plasticity model exhibiting deformation twinning YANG Jing,HUANG Jie,SUN Chao-yang,GUO Ning,WANG Bin 1)School of Mechanical Engineering,University of Science and Technology Beijing,Beijing 100083,China 2)Baoshan Iron Steel Co.Ltd.Shanghai 201900,China Corresponding author,E-mail:suncy@ustb.edu.cn ABSTRACT A crystal plasticity model exhibiting deformation twinning is introduced based on the classical crystal plasticity theory, and its numerical implementation is conducted in which a fully implicit integration procedure is employed.In this constitutive model, the saturation-type hardening law is adopted and the direct proportion relationship between twin resistance and slip hardening is used to describe the effect of twinning on slip hardening and twin hardening.By considering the physical meaning of all the 13 parameters in this model,the classification methodology is presented for these parameters.Taking Fe-22Mn-0.6C twinning induced plasticity (TWIP)steel as an example,a local sensitivity analysis of hardening parameters is investigated emphatically.The effects of hardening parameters on the macro-mechanical response,the activation and evolution of twinning and the strategy for determining the hardening parameters are discussed.According to the difference of deformation mechanisms,the macroscopic deformation are divided into a twin hardening stage and a twin hardening failure stage.And then the determination method and the suggestive ranges for the hardening parameters are given.It is found that the initial slip resistance is linearly associated with the yield limit,and the value of initial resist- ance stress ranges from 80 MPa to 160 MPa.The twin hardening stage weakens when the twin hardening index increases,the range of values allowed for the twin hardening index is 0 to 3.When the ratio between twinning resistance and slip resistance increases,the 收稿日期：2015-01-27 基金项目：国家自然科学基金资助项目(51105029,51575039)：国家自然科学基金委员会-中国工程物理研究院联合基金资助项目 (U1330121):非线性力学国家重点实验室开放基金资助项目(LNM201512):北京市自然科学基金资助项目(3112019)：北京市教委重点基金 资助项目(KZ201010005002)

工程科学学报，第 37 卷，第 8 期: 1076--1083，2015 年 8 月 Chinese Journal of Engineering，Vol． 37，No． 8: 1076--1083，August 2015 DOI: 10． 13374 /j． issn2095--9389． 2015． 08． 016; http: / /journals． ustb． edu． cn 耦合孪生晶体塑性模型硬化参数的灵敏度分析 杨 竞1) ，黄 杰1，2) ，孙朝阳1) ，郭 宁1) ，王 斌1) 1) 北京科技大学机械工程学院，北京 100083 2) 宝山钢铁股份有限公司，上海 201900  通信作者，E-mail: suncy@ ustb． edu． cn 摘 要 基于经典晶体塑性理论，建立了耦合孪生的晶体塑性本构模型并进行了全隐式积分的数值实现． 该本构模型采用 饱和硬化法则，并采用孪生阻力与滑移硬化之间的正比关系来描述孪生对滑移硬化影响及孪生硬化行为． 针对该本构模型 的 13 个参数，结合各参数物理意义提出了参数的分类确定方法． 以孪生诱导塑性( TWIP) 钢 Fe--22Mn--0. 6C 为例，着重对硬 化参数的局部灵敏度进行了分析，研究了各硬化参数对宏观力学响应、孪生激活和演化的影响，根据变形机制的不同宏观变 形过程可区分为孪生硬化阶段和孪生硬化失效阶段，进而给出了硬化参数确定的步骤及其建议取值范围． 结果表明: 初始滑 移阻力与屈服极限线性相关，取值范围在 80 ～ 160 MPa 之间; 孪生硬化指数增大使得孪生硬化阶段减弱，其取值范围应在 0 ～ 3 之间; 孪生阻力与滑移阻力比值增大，则孪生增长率降低，硬化率拐点后移，直至拐点消失，其取值范围在 1 ～ 1. 3 之间． 关键词 孪生; 晶体塑性; 本构模型; 灵敏度分析; 硬化参数 分类号 TG302 收稿日期: 2015--01--27 基金项目: 国家自然科学基金资助项目( 51105029，51575039 ) ; 国家自然科学基金委员会--中国工程物理研究院联合基金资助项目 ( U1330121) ; 非线性力学国家重点实验室开放基金资助项目( LNM201512) ; 北京市自然科学基金资助项目( 3112019) ; 北京市教委重点基金 资助项目( KZ201010005002) Sensitivity analysis of hardening parameters in the crystal plasticity model exhibiting deformation twinning YANG Jing1) ，HUANG Jie1，2) ，SUN Chao-yang1)  ，GUO Ning1) ，WANG Bin1) 1) School of Mechanical Engineering，University of Science and Technology Beijing，Beijing 100083，China 2) Baoshan Iron ＆ Steel Co． Ltd． ，Shanghai 201900，China  Corresponding author，E-mail: suncy@ ustb． edu． cn ABSTＲACT A crystal plasticity model exhibiting deformation twinning is introduced based on the classical crystal plasticity theory， and its numerical implementation is conducted in which a fully implicit integration procedure is employed． In this constitutive model， the saturation-type hardening law is adopted and the direct proportion relationship between twin resistance and slip hardening is used to describe the effect of twinning on slip hardening and twin hardening． By considering the physical meaning of all the 13 parameters in this model，the classification methodology is presented for these parameters． Taking Fe--22Mn--0. 6C twinning induced plasticity ( TWIP) steel as an example，a local sensitivity analysis of hardening parameters is investigated emphatically． The effects of hardening parameters on the macro-mechanical response，the activation and evolution of twinning and the strategy for determining the hardening parameters are discussed． According to the difference of deformation mechanisms，the macroscopic deformation are divided into a twin hardening stage and a twin hardening failure stage． And then the determination method and the suggestive ranges for the hardening parameters are given． It is found that the initial slip resistance is linearly associated with the yield limit，and the value of initial resist￾ance stress ranges from 80 MPa to 160 MPa． The twin hardening stage weakens when the twin hardening index increases，the range of values allowed for the twin hardening index is 0 to 3． When the ratio between twinning resistance and slip resistance increases，the

杨竞等：耦合孪生晶体塑性模型硬化参数的灵敏度分析 *1077· growth rate of twinning decreases and the turning point of the hardening stage moves backward until it is disappeared.The range of val- ues for the ratio between twinning resistance and slip resistance will be 1 to 1.3. KEY WORDS twinning:crystal plasticity:constitutive models;sensitivity analysis:hardening parameters 将非线性有限元与晶体塑性理论相结合，可以很 感系数对应变硬化的影响.以上研究说明了部分模型 好地再现材料的宏观力学响应”、各向异性、织构的形 参数对宏观力学响应及应变硬化等有显著影响，但仍 成演化-以及制耳特征，因而被广泛地用于金属塑 然没有建立各参数与宏观力学响应之间的联系，不能 性成形的模拟.经典的晶体塑性理论是以位错滑移理 为参数范围的确定及参数获取提供指导 论为基础的，其主要是考虑滑移对塑性变形的贡献. 本工作以FCC结构Fe-22Mn0.6C型孪生诱导 但是，滑移并不是多晶材料塑性变形的唯一方式，特别 塑性钢为例，建立其耦合孪生的晶体塑性本构模型和 是对一些低层错能金属，形变孪晶在塑性变形过程中 数值实现方法；基于参数物理意义，对其进行分类，并 发挥着重要的作用. 在此基础上着重对硬化参数的局部灵敏度进行分析， Kalidindi考虑形变孪晶对晶粒转向和应变硬化 研究各参数对宏观力学响应、李生的激活和演化的影 行为的影响，把孪生机制引入到晶体塑性本构模型中， 响，为各参数取值范围的确定和参数获取提供依据. 更好地描述了高锰孪生诱导塑性(twinning induced plasticity,TWIP)钢5、Mg合金m、a-T合金等 1耦合孪生的晶体塑性本构模型 金属的织构演化及制耳特征.但由于孪生与滑移耦合 1.1晶体塑性本构关系 的硬化模型比较复杂，目前对于孪生与滑移的相互作 对于诸如孪生诱导塑性钢、镁合金等多晶材料，晶 用尚无定论，需要根据材料滑移和李生变形的特点来 体的变形可以分解为两部分：由晶格畸变引起的弹性 进行选择和简化.如Kalidindi0采用的是一种饱和类 变形和刚体转动，以及由滑移和孪生引起的塑性变形 型的硬化法则，其着重考虑了孪生对滑移的影响，而忽 根据晶体塑性理论和变形梯度的乘法分解，可把晶体 略了其他方面对硬化的影响，孪生抗力在变形过程中 的变形梯度F分解为两部分： 与滑移抗力成正比，其硬化参数和饱和值都取决于晶 F=F·FP (1) 粒内的孪晶体积分数.该饱和硬化模型反映滑移和孪 式中，F为由塑性变形引起的变形梯度，F为弹性变 生相互作用机制是基于四个重要的假设：(1)一个晶 形和刚性转动合成的变形梯度 粒内形变孪晶的发生和发展会对所有的滑移系和孪生 考虑孪生对塑性变形的影响，塑性速度梯度L可 系产生硬化：(2)李晶边界与晶粒边界类似，会减小未 表示为回： 发生孪生区域中的滑移长度，这样通过Hall-Petch机 制就会增加应变硬化率@：(3)李晶要比基体组织更 L=mF-1=∑yS+∑jyS. (2) 难变形：(4)当变形程度较大时，李生会达到饱和，李 式中：L为塑性速度梯度；F为塑性变形梯度变化率： 晶体积分数的值不再增加.通过该模型，首次捕获滑 FP-为塑性变形梯度的逆：a和B分别为第a个滑移 移和李生之间复杂的相互作用关系圆 系和第B个孪生系；y为第α个滑移系的剪切速率： 然而，耦合李生的晶体塑性本构模型参数繁杂，如 为第B个孪生系的李晶体积分数变化率：y为孪生 Kalidindi的模型共l6个待定参数，且多数参数都无 引起的剪切量，对于FCC型晶体结构，y为常数，取 法直接与宏观力学响应相对应，难以直接确定.给出 y=0.707@:S“和S分别为第a个滑移系和第B个 的参数常常难以令分析结果具有唯一性或可靠性，目 孪生系的Schmid张量，即 前研究如Salem等和Wu等回主要是通过拟合不同 S=m⑧n (3) 变形方式，如单轴拉伸、单轴压缩、平面应变压缩及简 S9=m⑧n. (4) 单剪切条件下的宏观力学曲线来分步确定各个参数的 其中m°和n“分别为第α个滑移系的滑移方向和滑 值.但是，参数的选取范围、参数确定的先后顺序以及 移面法向单位向量，m和分别为第B个孪生系的 参数之间是否相互影响，没有定论 孪生方向和李生面法向单位向量 尽管国内外学者在晶体塑性本构模型参数对宏观 单晶体的弹性本构关系可表示为网 力学响应的影响方面做了许多研究."-网，试图解决 T=C:E (5) 上述问题.如i等如分析了率敏感系数对宏观力学 式中，T为第二Piola-Kirchhoff应力，E为与T能量共 性能的影响，研究表明率敏感系数越大，相同应变下的 轭的Green应变张量，而C为四阶弹性张量，是一个材 流变应力越大，其与应变速率对流变应力的影响一致. 料常数矩阵 郑华雷等四研究了α-金属中孪晶体积分数和率敏 第二Piola-Kirchhoff应力T可表示为

杨 竞等: 耦合孪生晶体塑性模型硬化参数的灵敏度分析 growth rate of twinning decreases and the turning point of the hardening stage moves backward until it is disappeared． The range of val￾ues for the ratio between twinning resistance and slip resistance will be 1 to 1. 3． KEY WOＲDS twinning; crystal plasticity; constitutive models; sensitivity analysis; hardening parameters 将非线性有限元与晶体塑性理论相结合，可以很 好地再现材料的宏观力学响应［1］、各向异性、织构的形 成演化［2--3］以及制耳特征，因而被广泛地用于金属塑 性成形的模拟． 经典的晶体塑性理论是以位错滑移理 论为基础的，其主要是考虑滑移对塑性变形的贡献． 但是，滑移并不是多晶材料塑性变形的唯一方式，特别 是对一些低层错能金属，形变孪晶在塑性变形过程中 发挥着重要的作用． Kalidindi［4］考虑形变孪晶对晶粒转向和应变硬化 行为的影响，把孪生机制引入到晶体塑性本构模型中， 更好地 描 述 了 高 锰 孪 生 诱 导 塑 性 ( twinning induced plasticity，TWIP) 钢［5--6］、Mg 合金［7］、α--Ti 合金［8--9］ 等 金属的织构演化及制耳特征． 但由于孪生与滑移耦合 的硬化模型比较复杂，目前对于孪生与滑移的相互作 用尚无定论，需要根据材料滑移和孪生变形的特点来 进行选择和简化． 如 Kalidindi［4］采用的是一种饱和类 型的硬化法则，其着重考虑了孪生对滑移的影响，而忽 略了其他方面对硬化的影响，孪生抗力在变形过程中 与滑移抗力成正比，其硬化参数和饱和值都取决于晶 粒内的孪晶体积分数． 该饱和硬化模型反映滑移和孪 生相互作用机制是基于四个重要的假设: ( 1) 一个晶 粒内形变孪晶的发生和发展会对所有的滑移系和孪生 系产生硬化; ( 2) 孪晶边界与晶粒边界类似，会减小未 发生孪生区域中的滑移长度，这样通过 Hall-Petch 机 制就会增加应变硬化率［10］; ( 3) 孪晶要比基体组织更 难变形; ( 4) 当变形程度较大时，孪生会达到饱和，孪 晶体积分数的值不再增加． 通过该模型，首次捕获滑 移和孪生之间复杂的相互作用关系［8］． 然而，耦合孪生的晶体塑性本构模型参数繁杂，如 Kalidindi 的模型共 16 个待定参数［4］，且多数参数都无 法直接与宏观力学响应相对应，难以直接确定． 给出 的参数常常难以令分析结果具有唯一性或可靠性，目 前研究如 Salem 等［8］和 Wu 等［9］主要是通过拟合不同 变形方式，如单轴拉伸、单轴压缩、平面应变压缩及简 单剪切条件下的宏观力学曲线来分步确定各个参数的 值． 但是，参数的选取范围、参数确定的先后顺序以及 参数之间是否相互影响，没有定论． 尽管国内外学者在晶体塑性本构模型参数对宏观 力学响应的影响方面做了许多研究［4，11--12］，试图解决 上述问题． 如 Li 等［11］分析了率敏感系数对宏观力学 性能的影响，研究表明率敏感系数越大，相同应变下的 流变应力越大，其与应变速率对流变应力的影响一致． 郑华雷等［12］研究了 α--Ti 金属中孪晶体积分数和率敏 感系数对应变硬化的影响． 以上研究说明了部分模型 参数对宏观力学响应及应变硬化等有显著影响，但仍 然没有建立各参数与宏观力学响应之间的联系，不能 为参数范围的确定及参数获取提供指导． 本工作以 FCC 结构 Fe--22Mn--0. 6C 型孪生诱导 塑性钢为例，建立其耦合孪生的晶体塑性本构模型和 数值实现方法; 基于参数物理意义，对其进行分类，并 在此基础上着重对硬化参数的局部灵敏度进行分析， 研究各参数对宏观力学响应、孪生的激活和演化的影 响，为各参数取值范围的确定和参数获取提供依据． 1 耦合孪生的晶体塑性本构模型 1. 1 晶体塑性本构关系 对于诸如孪生诱导塑性钢、镁合金等多晶材料，晶 体的变形可以分解为两部分: 由晶格畸变引起的弹性 变形和刚体转动，以及由滑移和孪生引起的塑性变形． 根据晶体塑性理论和变形梯度的乘法分解，可把晶体 的变形梯度 F 分解为两部分: F = Fe ·Fp ． ( 1) 式中，Fp 为由塑性变形引起的变形梯度，Fe 为弹性变 形和刚性转动合成的变形梯度． 考虑孪生对塑性变形的影响，塑性速度梯度 Lp 可 表示为［6］: Lp = F · p ·Fp － 1 = ∑α γ ·α Sα + ∑ β f ·β γtw Sβ ． ( 2) 式中: Lp 为塑性速度梯度; F · p 为塑性变形梯度变化率; Fp － 1为塑性变形梯度的逆; α 和 β 分别为第 α 个滑移 系和第 β 个孪生系; γ ·α 为第 α 个滑移系的剪切速率; f ·β 为第 β 个孪生系的孪晶体积分数变化率; γtw为孪生 引起的剪切量，对于 FCC 型晶体结构，γtw 为常数，取 γtw = 0. 707［10］; Sα 和 Sβ 分别为第 α 个滑移系和第 β 个 孪生系的 Schmid 张量，即 Sα = mα nα ． ( 3) Sβ = mβ nβ ． ( 4) 其中 mα 和 nα 分别为第 α 个滑移系的滑移方向和滑 移面法向单位向量，mβ 和 nβ 分别为第 β 个孪生系的 孪生方向和孪生面法向单位向量． 单晶体的弹性本构关系可表示为［9］ T = C∶ Ee ． ( 5) 式中，T 为第二 Piola-Kirchhoff 应力，Ee 为与 T 能量共 轭的 Green 应变张量，而 C 为四阶弹性张量，是一个材 料常数矩阵． 第二 Piola-Kirchhoff 应力 T 可表示为 · 7701 ·

·1078· 工程科学学报，第37卷，第8期 T=F-1.{det (F)a).F-T. (6) 式中，F-l为弹性变形梯度的逆，det(F)为弹性变形 ”=(-安)∑ (14) 梯度的行列式的值，o为Cauchy应力，Fe-T为弹性变 式中，“和“分别为第α个滑移系的滑移阻力及其变 形梯度转置的逆. 化率，h:和S:分别为第α个滑移系的硬化率和滑移 Green应变张量E可表示为 阻力饱和值. E-i(FF-D. 滑移系硬化率h:可表示为 (7) hg=h.(1+w) (15) 式中，FT为弹性变形梯度的转置，I为单位张量 式中，h,为无李晶时的硬化率，地和b均为李生硬化参 四阶弹性张量C可表示为 数，其中b为李晶体积分数的指数 [Cu CI2 C1 0 0 01 滑移阻力饱和值S可表示为 C2C1C200 0 S.=Sn+Sfs (16) 0 0 0 式中，S。为无李晶时滑移阻力饱和值，S表征了Hall- C= (8) 0 0 0 Petch效应对硬化的贡献，指数o.5”通过推导Hall- Ca 0 0 0 0 0 0 Petch关系gxla5(σ表示屈服极限，l表示晶粒度) Cu 0 及孪晶体积分数与孪晶平均间距的反比关系∫∝得 0 0 00 0 C) 到，其中1为滑移长度，为李晶平均间距 式中，C、C2和Cu为材料常数. 变形过程中的李生变形阻力则设为与滑移阻力成 根据Schmid定律，第a个滑移系的分解剪应力 正比的关系： “为 s=6s“ (17) “=rS (9) 式中，s为第B个孪生系的李生阻力，δ为孪生与滑移 同理，第B个孪生系的分解剪应力为 阻力比值 =r:S8. (10) 1.2流动法则 2本构模型的数值实现与参数选取 对于率相关的晶体塑性本构模型，滑移系的剪切 2.1晶体塑性本构模型的数值实现 速率采用黏塑性理论中的幂指数关系： 对于率相关晶体塑性本构关系，有半隐式积分方 7=|F“e). (11) 案和全隐式积分方案.积分方案的区别主要在于建立 非线性方程组时选择的独立变量的不同.Kalidindi 式中，y°为参考塑性剪切率，s“为第α《个滑移系的滑 提出的全隐式积分方案，独立变量为塑性变形梯度和 移变形阻力，m为材料的率敏感系数，sign(r)为r的 滑移系变形阻力：Sarma和Zacharia提出的全隐式积 符号函数 分方案中独立变量为弹性变形梯度；Marin等提出的全 由孪生引起的塑性剪切量表示为 隐式积分方案的独立变量是刚性转动张量的.本工 y=yre (12) 作针对晶体塑性率相关本构模型，以第二Piola一Kirch- 式中，y√为第B个孪生系的孪生剪切量，为第B个李 of应力T和滑移阻力s“为独立变量，推导其全隐式 生系的李晶体积分数.设定∫为所有李生系中李晶体 积分过程，并基于ABAOUS/UMAT平台开发其子程 积分数的累积和，即∫=∑户当∫达到某一阀值6 序.其数值更新过程如下 时，李晶体积分数达到饱和，即不再发生李生四 第1步：输入已知量，ln时刻的变形梯度F1n 时刻F。T.s和F,T。和s:分别为前一增量步结 孪生具有方向性，其孪晶体积分数变化率可以表 示为阳 束时的收敛值，Schmid张量S和S; 第2步：判断总孪晶体积分数∫是否达到阈值f。, 1/m (>0), 如果超过阈值，则不再发生孪生且∫=，否则直接进 (13) 入第3步： P=0 () 第3步：根据式(1)、式(5)和式(7)，求解t.1时 1.3滑移与孪生耦合的硬化模式表征 刻的状态变量试探值，即 目前应用比较广泛的硬化模型有线性硬化模型、 Fe=FF, 幂指数模型和饱和硬化模型8.本文采用饱和 硬化模型来描述孪生对滑移的硬化.滑移阻力的演化 E=2(FF-D, 规律如下： T=C:E

工程科学学报，第 37 卷，第 8 期 T = Fe － 1·{ det( Fe ) σ}·Fe － T ． ( 6) 式中，Fe － 1为弹性变形梯度的逆，det( Fe ) 为弹性变形 梯度的行列式的值，σ 为 Cauchy 应力，Fe － T 为弹性变 形梯度转置的逆． Green 应变张量 Ee 可表示为 Ee = 1 2 ( FeTFe － I) ． ( 7) 式中，FeT为弹性变形梯度的转置，I 为单位张量． 四阶弹性张量 C 可表示为 C = C11 C12 C12 0 0 0 C12 C11 C12 0 0 0 C12 C12 C11 0 0 0 0 0 0 C44 0 0 0 0 0 0 C44 0 0 0 0 0 0 C                  44  ． ( 8) 式中，C11、C12和 C44为材料常数． 根据 Schmid 定律，第 α 个滑移系的分解剪应力 τα 为 τα = σ∶ Sα ． ( 9) 同理，第 β 个孪生系的分解剪应力 τβ 为 τβ = σ∶ Sβ ． ( 10) 1. 2 流动法则 对于率相关的晶体塑性本构模型，滑移系的剪切 速率采用黏塑性理论中的幂指数关系: γ ·α = γ ·0 τα s α 1 /m sign( τα ) ． ( 11) 式中，γ ·0 为参考塑性剪切率，s α 为第 α 个滑移系的滑 移变形阻力，m 为材料的率敏感系数，sign( τα ) 为 τα 的 符号函数． 由孪生引起的塑性剪切量表示为 γβ = γtw f β ． ( 12) 式中，γβ 为第 β 个孪生系的孪生剪切量，f β 为第 β 个孪 生系的孪晶体积分数． 设定 f 为所有孪生系中孪晶体 积分数的累积和，即 f = ∑ β f β ． 当 f 达到某一阈值 f0 时，孪晶体积分数达到饱和，即不再发生孪生［13］． 孪生具有方向性，其孪晶体积分数变化率可以表 示为［4］ f ·β = γ · 0 γ ( tw τβ s β ) 1 /m ( τβ ＞ 0) ， f ·β = 0 ( τβ ＜ 0 或 f ＞ f0 ) { ． ( 13) 1. 3 滑移与孪生耦合的硬化模式表征 目前应用比较广泛的硬化模型有线性硬化模型、 幂指数模型［14］和饱和硬化模型［4，8--9］． 本文采用饱和 硬化模型来描述孪生对滑移的硬化． 滑移阻力的演化 规律如下: s ·α = hα s ( 1 － s α Sα ) s ∑α γ ·α ． ( 14) 式中，s α 和 s ·α 分别为第 α 个滑移系的滑移阻力及其变 化率，hα s 和 Sα s 分别为第 α 个滑移系的硬化率和滑移 阻力饱和值． 滑移系硬化率 hα s 可表示为 hα s = hs( 1 + wfb ) ． ( 15) 式中，hs 为无孪晶时的硬化率，w 和 b 均为孪生硬化参 数，其中 b 为孪晶体积分数的指数． 滑移阻力饱和值 Sα s 可表示为 Sα s = Ss0 + Spr f 0. 5 ． ( 16) 式中，Ss0为无孪晶时滑移阻力饱和值，Spr表征了 Hall￾Petch 效应对硬化的贡献［8］，指数“0. 5”通过推导 Hall￾Petch 关系 σ∝l － 0. 5 ( σ 表示屈服极限，l 表示晶粒度) 及孪晶体积分数与孪晶平均间距的反比关系 f∝l － 1 tw 得 到［4］，其中 l 为滑移长度，ltw为孪晶平均间距． 变形过程中的孪生变形阻力则设为与滑移阻力成 正比的关系: s β = δs α ． ( 17) 式中，s β 为第 β 个孪生系的孪生阻力，δ 为孪生与滑移 阻力比值． 2 本构模型的数值实现与参数选取 2. 1 晶体塑性本构模型的数值实现 对于率相关晶体塑性本构关系，有半隐式积分方 案和全隐式积分方案． 积分方案的区别主要在于建立 非线性方程组时选择的独立变量的不同． Kalidindi［4］ 提出的全隐式积分方案，独立变量为塑性变形梯度和 滑移系变形阻力; Sarma 和 Zacharia［14］提出的全隐式积 分方案中独立变量为弹性变形梯度; Marin 等提出的全 隐式积分方案的独立变量是刚性转动张量［15］． 本工 作针对晶体塑性率相关本构模型，以第二 Piola--Kirch￾hoff 应力 T 和滑移阻力 s α 为独立变量，推导其全隐式 积分过 程，并 基 于 ABAQUS /UMAT 平 台 开 发 其 子 程 序． 其数值更新过程如下． 第 1 步: 输入已知量，tn + 1时刻的变形梯度 Fn + 1，tn 时刻 Fn、Tn、s α n 和 Fp － 1 n ，Tn 和 s α n 分别为前一增量步结 束时的收敛值，Schmid 张量 Sα 和 Sβ ; 第 2 步: 判断总孪晶体积分数 f 是否达到阈值 f0， 如果超过阈值，则不再发生孪生且 f = f0，否则直接进 入第 3 步; 第 3 步: 根据式( 1) 、式( 5) 和式( 7) ，求解 tn + 1 时 刻的状态变量试探值，即 Fe tr n + 1 = Fn + 1·Fp － 1 n ， Ee tr n + 1 = 1 2 ( Fe trT n + 1·Fe tr n + 1 － I) ， Ttr n + 1 = C∶ Ee tr n + 1 ; · 8701 ·

杨竞等：耦合孪生晶体塑性模型硬化参数的灵敏度分析 ·1079· 第4步：设定初值为To,=T1,进入Newton跌代 值.本文选取C,=198GPa,C2=125GPa,C4= 法两层迭代阿求解关于T。和s非线性方程组： 122 GPa 第5步：如果综合判断T:”和s4”已满足收敛 (2)晶体结构参数.Fe-22Mn-0.6C的李生诱导 条件，跳到第6步，否则k=k+1并跳到第2步继续迭 塑性钢为面心立方结构，其有12个滑移系和12个李 代计算： 生系.滑移系统为{11}〈110)，孪生系统为{111} 第6步：计算其他状态变量，根据式(6)的逆运算 112)刃. 计算tn,时刻的Cauchy应力oa (3)李晶体积分数阈值参数.孪晶体积分数阈值 2.2晶体塑性本构模型参数的分类及选取 。,可通过检测不同变形量下孪晶体积分数获得，而李 为了使用该模型计算滑移和李生导致的宏观塑性 晶体积分数测试方法参考Renard和Jacques图的工 力学行为及其微观结构的演化，还需要确定该模型中 作.当应变达到一定值时，李晶体积分数达到饱和基 的相关材料参数.该模型中参数较多，共13个待定参 本不变，可认此时的李晶体积分数即为阈值。，本文取 数.为此，本工作则以Fe-22Mn0.6C钢为例，考虑各 为0.29-20 种参数的物理意义，提出了一种晶体塑性本构模型参 (4)率相关参数.对于准静态变形时的力学行 数的分类和确定方法 为，参考塑性剪切率。1设为0.001s1,率敏感系 (1)弹性参数.四阶弹性张量C中的三个材料常 数m设为0.024. 数C、C2和C4,该材料常数的值一般可在工具书中 (5)滑移和李生硬化参数{s,h.,So,b,,Se, 查到.对于FCC结构Fe-22Mn0.6C的李生诱导塑性 δ}参考王伟华Fe-22Mn0.6C钢的研究工作，其具 钢，可通过文献6]中给出的方法确定三个参数的 体取值见表1. 表1Fe-22Mn0.6C的孪生诱导塑性钢品体塑性本构模型材料参数 Table 1 Material parameters in the crystal plasticity constitutive model for Fe-22Mn-0.6C TWIP steel Cu/GPa Cp/GPa Ca/GPa Yo/s-1 h,/MPa So/MPa b w S/MPa s8/MPa fo 198 125 122 0.0010.024300 300 220 100 90 0.21.16 800T 3本构模型硬化参数的灵敏度分析 700 在FCC结构晶体塑性本构模型的基础上，以 Fe-22Mn0.6C钢为例，着重对硬化参数的局部灵敏 0和。w、泽 600 50 度进行了分析，研究了初始滑移阻力s:、孪生硬化指数 44444444 --=50 MPa b、饱和硬化参数S和孪生阻力与滑移阻力比值δ对 300 ★s。=l00MPa 宏观力学响应、孪生激活和演化的影响，并给出其建议 200 8=150 MPa 7-=200NPa 取值范围. 100 3.1初始滑移阻力s的影响 0.1 02 0.3 0.4 0.50.6 等效应变 式(14)只给出了滑移阻力的变化率，因此还需要 图1不同初始滑移阻力下等效应力一应变曲线 给定一个初始的滑移阻力s。当分解剪应力“≥$。 Fig.I Equivalent stress-strain curves at different initial slip resist- 时，滑移系开始启动，材料才能产生塑性变形.由 ances 式(9)可知，在Schimd因子变化很小时，分解剪应力与 Cauchy应力呈线性关系.当“=s6时，滑移系统刚开 120MPa,Schmid因子约为sg/o,=0.42,与文献21]所 始滑移，材料屈服产生塑性变形，此时的Cauchy应力 测结果基本一致.同时，初始滑移阻力影响硬化率的 则为屈服极限，因此屈服极限与初始滑移阻力也呈线 幅度大小，但不改变其变化趋势 性关系.图1给出了s:在50~200MPa之间变化时等3.2孪生硬化指数b的影响 效应力一应变曲线. 式(15)中，b为李晶体积分数的指数，直接影响滑 图2为不同初始滑移阻力下应变硬化率(dσ/移系硬化率h:值的大小，如图4所示，当b=0时，相 dε)，图3为屈服极限与初始滑移阻力之间的关系.随 当于完全忽略了李生硬化的作用，h:为常数，因此b 着初始滑移阻力的增加，计算的等效应力曲线的屈服 应取大于0的数值 极限呈线性增加，当s:=50MPa时，屈服极限σ，约为 图5为等效应变硬化率曲线，可将曲线分成两个

杨 竞等: 耦合孪生晶体塑性模型硬化参数的灵敏度分析 第4 步: 设定初值为 T( 0) n + 1 = Ttr n + 1，进入 Newton 跌代 法两层迭代［15］求解关于 Tn + 1和 s α n + 1非线性方程组; 第 5 步: 如果综合判断 T( k + 1) n + 1 和 s α( k + 1) n + 1 已满足收敛 条件，跳到第 6 步，否则 k = k + 1 并跳到第 2 步继续迭 代计算; 第 6 步: 计算其他状态变量，根据式( 6) 的逆运算 计算 tn + 1时刻的 Cauchy 应力 σn + 1 ． 2. 2 晶体塑性本构模型参数的分类及选取 为了使用该模型计算滑移和孪生导致的宏观塑性 力学行为及其微观结构的演化，还需要确定该模型中 的相关材料参数． 该模型中参数较多，共 13 个待定参 数． 为此，本工作则以 Fe--22Mn--0. 6C 钢为例，考虑各 种参数的物理意义，提出了一种晶体塑性本构模型参 数的分类和确定方法． ( 1) 弹性参数． 四阶弹性张量 C 中的三个材料常 数 C11、C12和 C44，该材料常数的值一般可在工具书中 查到． 对于 FCC 结构 Fe--22Mn--0. 6C 的孪生诱导塑性 钢，可通过文献［16］中给出的方法确定三个参数的 值． 本 文 选 取 C11 = 198 GPa，C12 = 125 GPa，C44 = 122 GPa［6］． ( 2) 晶体结构参数． Fe--22Mn--0. 6C 的孪生诱导 塑性钢为面心立方结构，其有 12 个滑移系和 12 个孪 生系． 滑 移 系 统 为{ 111} 〈110〉，孪 生 系 统 为{ 111 } 〈112－ 〉［17］． ( 3) 孪晶体积分数阈值参数． 孪晶体积分数阈值 f0，可通过检测不同变形量下孪晶体积分数获得，而孪 晶体积分数测试方法参考 Ｒenard 和 Jacques［18］ 的工 作． 当应变达到一定值时，孪晶体积分数达到饱和基 本不变，可认此时的孪晶体积分数即为阈值 f0，本文取 为 0. 2［19--20］． ( 4) 率相关参数． 对于准静态变形时的力学行 为，参考塑性剪切率 γ ·［3--4，11］ 0 设为 0. 001 s － 1，率敏感系 数 m 设为 0. 024． ( 5) 滑移和孪生硬化参数{ s α 0，hs，Ss0，b，w，Spr， δ} 参考王伟华［15］Fe--22Mn--0. 6C 钢的研究工作，其具 体取值见表 1． 表 1 Fe--22Mn--0. 6C 的孪生诱导塑性钢晶体塑性本构模型材料参数 Table 1 Material parameters in the crystal plasticity constitutive model for Fe--22Mn--0. 6C TWIP steel C11 /GPa C12 /GPa C44 /GPa γ · 0 /s － 1 m hs /MPa Ss0 /MPa b w Spr /MPa s α 0 /MPa f0 δ 198 125 122 0. 001 0. 024 300 300 2 20 100 90 0. 2 1. 16 3 本构模型硬化参数的灵敏度分析 在 FCC 结 构 晶 体 塑 性 本 构 模 型 的 基 础 上，以 Fe--22Mn--0. 6C钢为例，着重对硬化参数的局部灵敏 度进行了分析，研究了初始滑移阻力 s α 0、孪生硬化指数 b、饱和硬化参数 Spr和孪生阻力与滑移阻力比值 δ 对 宏观力学响应、孪生激活和演化的影响，并给出其建议 取值范围． 3. 1 初始滑移阻力 s α 0 的影响 式( 14) 只给出了滑移阻力的变化率，因此还需要 给定一个初始的滑移阻力 s α 0 ． 当分解剪应力 τα ≥s α 0 时，滑移 系 开 始 启 动，材 料 才 能 产 生 塑 性 变 形． 由 式( 9) 可知，在 Schimd 因子变化很小时，分解剪应力与 Cauchy 应力呈线性关系． 当 τα = s α 0 时，滑移系统刚开 始滑移，材料屈服产生塑性变形，此时的 Cauchy 应力 则为屈服极限，因此屈服极限与初始滑移阻力也呈线 性关系． 图 1 给出了 s α 0 在 50 ～ 200 MPa 之间变化时等 效应力--应变曲线． 图 2 为 不 同 初 始 滑 移 阻 力 下 应 变 硬 化 率( dσ/ dε) ，图 3 为屈服极限与初始滑移阻力之间的关系． 随 着初始滑移阻力的增加，计算的等效应力曲线的屈服 极限呈线性增加，当 s α 0 = 50 MPa 时，屈服极限 σy 约为 图 1 不同初始滑移阻力下等效应力--应变曲线 Fig． 1 Equivalent stress--strain curves at different initial slip resist￾ances 120 MPa，Schmid 因子约为 s α 0 /σy = 0. 42，与文献［21］所 测结果基本一致． 同时，初始滑移阻力影响硬化率的 幅度大小，但不改变其变化趋势． 3. 2 孪生硬化指数 b 的影响 式( 15) 中，b 为孪晶体积分数的指数，直接影响滑 移系硬化率 hα s 值的大小，如图 4 所示，当 b = 0 时，相 当于完全忽略了孪生硬化的作用，hα s 为常数，因此 b 应取大于 0 的数值． 图 5 为等效应变硬化率曲线，可将曲线分成两个 · 9701 ·

·1080· 工程科学学报，第37卷，第8期 1500 3000 --3=50 MPa -★s=100MPa 1250 a。=150MPa 2500 o-b=0 7-=2001Pa ★b=1 1000 2000 △-=2 -=3 750 hy 4chind-ihy 1500 女年 50 10 aAA 0哈其。w、合含合合A 500 250 守0-00女* 01 0.2 0.5n 0.6 0.1 0.20.30.405 0.6 等效应变 等效应变 图2不同初始滑移阻力下应变硬化率 图5不同孪生硬化指数时等效应变硬化率曲线 Fig.2 Strain hardening rate at different initial slip resistances Fig.5 Equivalent strain hardening rate curves at different twin hard- ening indices 500 生硬化.当b=3时，孪生硬化的效应已失效，阶段A 400 与阶段B拐点基本消失.因此b应取小于3的值，故b 300 的大小应在范围0~3之间调整. 对于b的取值，并不影响孪生系的激活，如图6所 200 示，不同b取值时，均激活李生系11)12]和(1i1) I00 Di2],且两个孪生系的李晶体积分数演化规律保持 不变 50 100150 200 250 初始滑移阻力fMP 0.12 d0000D0-0g040 图3初始滑移阻力与屈服应力之间的关系 0.10 Fig.3 Relationship between initial slip resistance and yield stress 0.08 7000 0.06 00fr-0-0-0-0-0-0-0-0-0-0-0-0-0-0-0-0-0-0-0-0-0 4-(111112引 0.04 0-(1111121 5000 0-b=0 0.02 ★b=1 △-b=2 -3 0.1 020.30.40.50.6 等效应变 2000 图6不同孪生硬化指数时激活孪生系的孪品体积分数 1000 Fig.6 Twinning volume fraction of the activated twin system at dif- 0.1 0.2030.40.5 0.6 ferent twin hardening indices 等效应变 图4不同孪生硬化指数时h:值 3.3饱和硬化参数S的影响 Fig.4 h values at different twin hardening indices 在式(16)中，饱和硬化参数S是李晶体积分数的 系数，最终影响滑移阻力饱和值的大小.图7为滑移 阶段，即李生硬化阶段A和李生硬化失效阶段B.在 阻力饱和值在S不同取值时的变化曲线.随着S增 阶段A硬化率上升，孪生对滑移产生明显的硬化效应： 大，滑移阻力饱和值一起增大，在等效应变约0.21时 而在阶段B,硬化率持续下降，此时不再有新的孪晶产 达到最大值.图8为S不同取值时滑移阻力的变化曲 生，孪生硬化效果消失.应变硬化率(dσ/ds,σ一应 线.滑移阻力受S的影响较小，只有在等效应变较大 力，ε一应变)曲线体现为李生硬化阶段A的变化.当 时才有所区别.这是因为在式(14)中，虽然1-s“/S b≥1时，随着b逐渐增大，：逐渐减小，阶段A不断减 随着饱和值的增大而增大，但幅度较小，远没有h:对 弱，影响阶段A与阶段B拐点处的应变硬化率逐渐减 滑移阻力的贡献大.图9为应变硬化率与等效应变的 小，与Wu等回对a-Ti的研究结果相似.当b=0时， 关系.可以看出，S同样也影响阶段A与阶段B拐点 没有孪晶体积分数的影响，曲线一直下降，没有出现孪 处硬化率的大小，但没有参数b对阶段A的影响显著

工程科学学报，第 37 卷，第 8 期 图 2 不同初始滑移阻力下应变硬化率 Fig． 2 Strain hardening rate at different initial slip resistances 图 3 初始滑移阻力与屈服应力之间的关系 Fig． 3 Ｒelationship between initial slip resistance and yield stress 图 4 不同孪生硬化指数时 hα s 值 Fig． 4 hα s values at different twin hardening indices 阶段，即孪生硬化阶段 A 和孪生硬化失效阶段 B． 在 阶段 A 硬化率上升，孪生对滑移产生明显的硬化效应; 而在阶段 B，硬化率持续下降，此时不再有新的孪晶产 生，孪生硬化效果消失． 应变硬化率( dσ/dε，σ—应 力，ε—应变) 曲线体现为孪生硬化阶段 A 的变化． 当 b≥1 时，随着 b 逐渐增大，hα s 逐渐减小，阶段 A 不断减 弱，影响阶段 A 与阶段 B 拐点处的应变硬化率逐渐减 小，与 Wu 等［9］对 α--Ti 的研究结果相似． 当 b = 0 时， 没有孪晶体积分数的影响，曲线一直下降，没有出现孪 图 5 不同孪生硬化指数时等效应变硬化率曲线 Fig． 5 Equivalent strain hardening rate curves at different twin hard￾ening indices 生硬化． 当 b = 3 时，孪生硬化的效应已失效，阶段 A 与阶段 B 拐点基本消失． 因此 b 应取小于 3 的值，故 b 的大小应在范围 0 ～ 3 之间调整． 对于 b 的取值，并不影响孪生系的激活，如图 6 所 示，不同 b 取值时，均激活孪生系(1 － 11) ［1 － 12 － ］和( 11 － 1) ［11 － 2 － ］，且两个孪生系的孪晶体积分数演化规律保持 不变． 图 6 不同孪生硬化指数时激活孪生系的孪晶体积分数 Fig． 6 Twinning volume fraction of the activated twin system at dif￾ferent twin hardening indices 3. 3 饱和硬化参数 Spr的影响 在式( 16) 中，饱和硬化参数 Spr是孪晶体积分数的 系数，最终影响滑移阻力饱和值的大小． 图 7 为滑移 阻力饱和值在 Spr不同取值时的变化曲线． 随着 Spr增 大，滑移阻力饱和值一起增大，在等效应变约 0. 21 时 达到最大值． 图 8 为 Spr不同取值时滑移阻力的变化曲 线． 滑移阻力受 Spr的影响较小，只有在等效应变较大 时才有所区别． 这是因为在式( 14) 中，虽然 1 － s α / Sα s 随着饱和值的增大而增大，但幅度较小，远没有 hα s 对 滑移阻力的贡献大． 图 9 为应变硬化率与等效应变的 关系． 可以看出，Spr同样也影响阶段 A 与阶段 B 拐点 处硬化率的大小，但没有参数 b 对阶段 A 的影响显著， · 0801 ·

杨竞等：耦合孪生晶体塑性模型硬化参数的灵敏度分析 1081 而且S完全不影响拐点处对应等效应变值的大小，因 化率随着6增大而不断减小.在8较小时，孪晶体积分 此参数S取值范围较广，其取值应对硬化率曲线拐点 数变化率很大，孪晶体积分数能很快达到饱和，如 进行微调. 图11所示.同时，一般认为孪生阻力大于滑移阻力 800 因此，δ适合在范围1~1.3之间取值. =50 MPa 0.25 700 S=100 MPa 9=20011Pa 600 =500 MPa 0.20 5=1000MP 效 500 0.15 400 0-d-1.14 8-1.16 △-d=1.18 0.05 7-d=1.20 0.1 0.2 0.30.40.5 0.6 4-d=130 等效应变 KK444111444444KK+ 图7不同S时的阻力饱和值 0.1 02 03 0.4 0.5 0.6 等效应变 Fig.7 Resistance saturation values at different S values 图10不同8时李品体积分数 500 Fig.10 Twinning volume fraction at different 6 values --S=50 MPa 400 女S=100MPa 1.5 =200 MPa =500 MPa 1.2 300 0.9 200 0.6 0.3 0.1 0.20.30.4 0.5 0.6 等效应变 图8不同S时的滑移阻力 1.16 1.20 124 1.28 1.32 李生阻力和滑移阻力比值 Fig.8 Slip resistances at different S values 图116与李晶体积分数变化率之间的关系 250 Fig.11 Relationship between 6 values and the change rate of twin -o-S =50 MPa ★-S=100MP4 crystal fraction 20 B --=200 MPa S=500 MPa 在图12中，当8=1.2和8=1.3时，应变硬化率 =1000MP 44444-44-4444 都没有出现图5所说的阶段A和阶段B,而其余三个 值都有很明显的拐点.这是因为当8=1.2和δ=1.3 时，李晶体积分数均未达到饱和值，晶粒未发生转向. 当8<1.2时，随着8增大，拐点对应的应变值也不断 500 办8个各各各含者水8-各办-t 增大.这与图10孪晶体积分数达到饱和时应变不断 0.1 0.2 0.30.40.5 0.6 增大类似.而且，孪李晶体积分数达到饱和值时等效应 等效应变 变值与图12硬化率拐点处等效应变值基本相当.可 图9不同S时等效应变硬化率 见，δ不仅影响李生变化率，还影响应变硬化率拐点处 Fig.9 Equivalent strain hardening rate under different S values 对应等效应变值的大小.故可通过应变硬化率曲线拐 点处应变值确定δ，同时还可用孪晶体积分数随应变 3.4孪生阻力与滑移阻力比值δ的影响 的变化曲线来进行验证. 在图10中，随着8不断增大，孪晶体积分数达到 饱和时应变不断增大，直至不能达到饱和，如在δ= 4讨论 1.2时.在6=1.3时，孪晶体积分数甚至接近为零，因 由上面分析可知，各参数均能影响李生硬化行为， 此再增大8值已没有意义.孪晶体积分数饱和前的变 且参数可能同时对不同的阶段产生影响，如图2所示

杨 竞等: 耦合孪生晶体塑性模型硬化参数的灵敏度分析 而且 Spr完全不影响拐点处对应等效应变值的大小，因 此参数 Spr取值范围较广，其取值应对硬化率曲线拐点 进行微调． 图 7 不同 Spr时的阻力饱和值 Fig． 7 Ｒesistance saturation values at different Spr values 图 8 不同 Spr时的滑移阻力 Fig． 8 Slip resistances at different Spr values 图 9 不同 Spr时等效应变硬化率 Fig． 9 Equivalent strain hardening rate under different Spr values 3. 4 孪生阻力与滑移阻力比值 δ 的影响 在图 10 中，随着 δ 不断增大，孪晶体积分数达到 饱和时应变不断增大，直至不能达到饱和，如在 δ = 1. 2 时． 在 δ = 1. 3 时，孪晶体积分数甚至接近为零，因 此再增大 δ 值已没有意义． 孪晶体积分数饱和前的变 化率随着 δ 增大而不断减小． 在 δ 较小时，孪晶体积分 数变化 率 很 大，孪 晶 体 积 分 数 能 很 快 达 到 饱 和，如 图 11 所示． 同时，一般认为孪生阻力大于滑移阻力． 因此，δ 适合在范围 1 ～ 1. 3 之间取值． 图 10 不同 δ 时孪晶体积分数 Fig． 10 Twinning volume fraction at different δ values 图 11 δ 与孪晶体积分数变化率之间的关系 Fig． 11 Ｒelationship between δ values and the change rate of twin crystal fraction 在图 12 中，当 δ = 1. 2 和 δ = 1. 3 时，应变硬化率 都没有出现图 5 所说的阶段 A 和阶段 B，而其余三个 值都有很明显的拐点． 这是因为当 δ = 1. 2 和 δ = 1. 3 时，孪晶体积分数均未达到饱和值，晶粒未发生转向． 当 δ ＜ 1. 2 时，随着 δ 增大，拐点对应的应变值也不断 增大． 这与图 10 孪晶体积分数达到饱和时应变不断 增大类似． 而且，孪晶体积分数达到饱和值时等效应 变值与图 12 硬化率拐点处等效应变值基本相当． 可 见，δ 不仅影响孪生变化率，还影响应变硬化率拐点处 对应等效应变值的大小． 故可通过应变硬化率曲线拐 点处应变值确定 δ，同时还可用孪晶体积分数随应变 的变化曲线来进行验证． 4 讨论 由上面分析可知，各参数均能影响孪生硬化行为， 且参数可能同时对不同的阶段产生影响，如图 2 所示， · 1801 ·

·1082 工程科学学报，第37卷，第8期 1600 图5中阶段A,同时受参数s:、b、S和8的影响.其中 0=1.14 。对后面李生硬化阶段A的影响，可在对后续参数的 1200 ò=1.20 拟合中得到调整；参数b对阶段A的幅值影响最大：参 ■6-1.30 A■ 数δ对硬化阶段的长度影响最大，这是因为随着δ增 800 大，李生阻力不断增大，可在式(13)和式(17)中反映， 孪晶体积分数变化率则不断降低，如图11所示，因此 400 李晶体积分数达到饱和值越慢，孪生硬化阶段越长 S只对阶段A的幅值有较小影响.因此可先确定阶段 0.1 0.2 0.3 0.4 0.5 0.6 A的幅值，获得孪生硬化指数b,然后根据硬化阶段A 等效应变 的长度获得6，同时8对阶段A的幅值产生影响，则需 图12不同8时等效应变硬化率 要调整S来消除. Fig.12 Equivalent strain hardening rate at different 6 values 由前面的分析也可以进一步分析各硬化参数的取 值范围，也是通过这两个阶段来判定.如初始滑移阻 初始滑移阻力不仅与屈服极限线性相关，而且还影响 力s,其与屈服极限呈线性关系，其比值约为$/σ，= 其应变硬化率的大小.某一阶段还可能同时受不同参 0.42,孪生诱导塑性钢的屈服极限一般在200~ 数影响，如图5、图9和图12所示，参数s6、b、S和6对 400MPa之间2-，故孪生诱导塑性钢s6的取值范围 硬化阶段A都有影响.因此在参数确定时，需考虑参 在80~160MPa之间.参数b的取值范围则是通过孪 数确定的先后顺序.首先，确定弹性参数如CC2、 生硬化阶段是否显著来确定.如图5所示：当b=0 C4。和率参数如参考应变速率及率敏感系数，然后 时，没有孪生硬化阶段，其塑性变形完全由滑移引起： 根据宏观力学响应中应变由小到大的顺序，分两阶段 当b=3时，李生硬化阶段的拐点也不明显：当b在0~ 确定对应参数的值. 3之间时，其孪生硬化阶段明显，此时李生机制才起主 对于图5和图9的孪生硬化阶段A和李生硬化失 导作用，故b的取值范围确定为0~3之间. 效阶段B.对于孪生硬化阶段前的材料屈服主要受初 根据对模型硬化参数灵敏度的分析，重新调整了 始滑移阻力8的影响，此时孪生不是主要机制，滑移 孪生诱导塑性钢耦合孪生的晶体塑性本构模型硬化参 系的开启直接影响材料的屈服.孪生硬化阶段A,即 数的值，如表2所示 表2优化的Fe-22Mn0.6C李生诱导塑性钢本构模型材料参数 Table 2 Optimized model parameters for Fe-22Mn-0.6C TWIP steel Cu /GPa Cn2/GPa C/GPa Yo/s-1 m h,/MPa So/MPa b S /MPa so/MPa 198 125 122 0.0010.01 600 600 3 10 400 120 0.21.16 图13为真应力-真应变曲线及硬化参数优化前后 1200 文献23]实验结果 的计算结果.由图13可知，不同参数下真应力一真应 0-表1参数计算结果 -一表2参数计算结果 变曲线计算结果差距很大.参数未优化调整前，其计 1000 算结果与实验结果存在较大的偏差：调整硬化参数后， 800 计算结果与实验结果吻合较好.基于以上灵敏度分 600 析，调整硬化参数后的晶体塑性本构模型能较好地反 400 映李生诱导塑性钢的宏观力学性能，说明硬化参数灵 敏度分析对参数获取具有一定指导意义 200 5 结论 0.10.20.30.4050.60.7 真应变 建立了FCC结构李生诱导塑性钢耦合孪生的晶 图13参数优化前后的计算结果与实验结果对比 体塑性本构模型，通过全隐式积分对模型进行了数值 Fig.13 Comparison of results before and after parameter optimization 实现.在对本构模型特征参数分类的基础上，对硬化 with experimental results 参数进行了灵敏度分析，获得了各参数对宏观变形、李 实验结果吻合较好，更好地反映了孪生诱导塑性钢的 生激活和演化的影响大小，确定了各参数的取值范围 宏观力学行为 及其获取步骤.基于以上分析，优化调整后的参数与 (1)提出了本构模型参数的分类确定方法，尤其

工程科学学报，第 37 卷，第 8 期 图 12 不同 δ 时等效应变硬化率 Fig． 12 Equivalent strain hardening rate at different δ values 初始滑移阻力不仅与屈服极限线性相关，而且还影响 其应变硬化率的大小． 某一阶段还可能同时受不同参 数影响，如图 5、图 9 和图 12 所示，参数 s α 0、b、Spr和 δ 对 硬化阶段 A 都有影响． 因此在参数确定时，需考虑参 数确定的先后顺序． 首先，确定弹性参数如 C11、C12、 C44、f0 和率参数如参考应变速率及率敏感系数，然后 根据宏观力学响应中应变由小到大的顺序，分两阶段 确定对应参数的值． 对于图 5 和图 9 的孪生硬化阶段 A 和孪生硬化失 效阶段 B． 对于孪生硬化阶段前的材料屈服主要受初 始滑移阻力 s α 0 的影响，此时孪生不是主要机制，滑移 系的开启直接影响材料的屈服． 孪生硬化阶段 A，即 图 5 中阶段 A，同时受参数 s α 0、b、Spr和 δ 的影响． 其中 s α 0 对后面孪生硬化阶段 A 的影响，可在对后续参数的 拟合中得到调整; 参数 b 对阶段 A 的幅值影响最大; 参 数 δ 对硬化阶段的长度影响最大，这是因为随着 δ 增 大，孪生阻力不断增大，可在式( 13) 和式( 17) 中反映， 孪晶体积分数变化率则不断降低，如图 11 所示，因此 孪晶体积分数达到饱和值越慢，孪生硬化阶段越长． Spr只对阶段 A 的幅值有较小影响． 因此可先确定阶段 A 的幅值，获得孪生硬化指数 b，然后根据硬化阶段 A 的长度获得 δ，同时 δ 对阶段 A 的幅值产生影响，则需 要调整 Spr来消除． 由前面的分析也可以进一步分析各硬化参数的取 值范围，也是通过这两个阶段来判定． 如初始滑移阻 力 s α 0，其与屈服极限呈线性关系，其比值约为 s α 0 /σy = 0. 42，孪 生诱导塑性钢的屈服极限一般在 200 ～ 400 MPa之间［22--23］，故孪生诱导塑性钢 s α 0 的取值范围 在 80 ～ 160 MPa 之间． 参数 b 的取值范围则是通过孪 生硬化阶段是否显著来确定． 如图 5 所示: 当 b = 0 时，没有孪生硬化阶段，其塑性变形完全由滑移引起; 当 b = 3 时，孪生硬化阶段的拐点也不明显; 当 b 在 0 ～ 3 之间时，其孪生硬化阶段明显，此时孪生机制才起主 导作用，故 b 的取值范围确定为 0 ～ 3 之间． 根据对模型硬化参数灵敏度的分析，重新调整了 孪生诱导塑性钢耦合孪生的晶体塑性本构模型硬化参 数的值，如表 2 所示． 表 2 优化的 Fe--22Mn--0. 6C 孪生诱导塑性钢本构模型材料参数 Table 2 Optimized model parameters for Fe--22Mn--0. 6C TWIP steel C11 /GPa C12 /GPa C44 /GPa γ · 0 /s － 1 m hs /MPa Ss0 /MPa b w Spr /MPa s α 0 /MPa f0 δ 198 125 122 0. 001 0. 01 600 600 3 10 400 120 0. 2 1. 16 图 13 为真应力--真应变曲线及硬化参数优化前后 的计算结果． 由图 13 可知，不同参数下真应力--真应 变曲线计算结果差距很大． 参数未优化调整前，其计 算结果与实验结果存在较大的偏差; 调整硬化参数后， 计算结果与实验结果吻合较好． 基于以上灵敏度分 析，调整硬化参数后的晶体塑性本构模型能较好地反 映孪生诱导塑性钢的宏观力学性能，说明硬化参数灵 敏度分析对参数获取具有一定指导意义． 5 结论 建立了 FCC 结构孪生诱导塑性钢耦合孪生的晶 体塑性本构模型，通过全隐式积分对模型进行了数值 实现． 在对本构模型特征参数分类的基础上，对硬化 参数进行了灵敏度分析，获得了各参数对宏观变形、孪 生激活和演化的影响大小，确定了各参数的取值范围 及其获取步骤． 基于以上分析，优化调整后的参数与 图 13 参数优化前后的计算结果与实验结果对比 Fig． 13 Comparison of results before and after parameter optimization with experimental results 实验结果吻合较好，更好地反映了孪生诱导塑性钢的 宏观力学行为． ( 1) 提出了本构模型参数的分类确定方法，尤其 · 2801 ·

杨竞等：耦合孪生晶体塑性模型硬化参数的灵敏度分析 ·1083· 是详细给出了硬化参数确定步骤，即先根据屈服极限 hardening behavior in V-5Cr-5Ti alloy.Chin J Appl Mech, 确定初始滑移阻力s:,再根据应变硬化率拐点和阶段 2012,44(2):334 A分别确定8和孪生硬化指数b,最后对整个应变硬化 (余勇，潘晓霞，谢若泽，等.孪生诱发塑性对V5C5Ti合 率曲线微调确定饱和参数S 金应变硬化行为的影响.力学学报，2012,44(2)：334) 01] Li H W,Yang H,Sun Z C.A robust integration algorithm for (2)初始滑移阻力s8与屈服极限呈线性关系。孪 implementing rate dependent crystal plasticity into explicit finite 生硬化指数b增大使得孪生硬化阶段A减弱.孪生阻 element method.Int J Plast,2008,24(2)267 力与滑移阻力比值8增大，孪生增长率降低，且硬化率 [12]Zheng H L,Yang H,Li H W.Crystal plasticity finite element 拐点后移，直至8=1.2时拐点消失 modeling for uniaxial compression of commercially pure titanium. (3)获得了Fe-22Mn0.6C型孪生诱导塑性钢耦 JPlast Eng,2013,20(1):95 (郑华雷，杨合，李宏伟.纯钛压缩变形下的品体塑性有限 合孪生晶体塑性模型硬化参数的取值范围，其中s。为 元分析.塑性工程学报，2013,20(1)：95) 80~160MPa之间，b为0~3之间，8为1~1.3之间. 03] Bouaziz O,Allain S,Scott C P,et al.High manganese austenit- ic twinning induced plasticity steels:a review of the microstruc- 参考文献 ture properties relationships.Curr Opin Solid State Mater Sci, 2011,15(4):141 1]Zhang K S,Zhang G.Feng L Stress calculation of single crystal [14]Sarma G,Zacharia T.Integration algorithm for modeling the elastoviscoplastic response of polycrystalline materials.J Mech under finite slip plastic deformation.Acta Mech Sin,2002,34 (4):636 Phys Solid,1999,47(6):1219 [15]Wang W H.Study on Plastie Deformation Beharior Induced by (张克实，张光，冯露.单品体塑性滑移有限变形下的应力计 算.力学学报，2002,34(4)：636) Slipping and Ticinning for TWIP Steel [Dissertation].Beijing: Li D Y,Zhang S R.Peng Y H,et al.Finite element simulation University of Science and Technology Beijing,2012 (王伟华.TWP钢滑移和孪生耦合的塑性变形行为[学位论 of sheet metal stamping with polycrystalline plasticity.Chin 文].北京：北京科技大学，2012) Mech Eng,2008,44(1):190 06 (李大永，张少蓉，彭颖红，等.板材冲压成形的品体塑性有 Gebhardt T,Music D,Kossmann D,et al.Elastic properties of 限元模拟.机械工程学报，2008,44(1)：190) fce Fe-Mn-X (X=Al,Si)alloys studied by theory and experi- B3]Pi H C.Han JT,Zhang C G,et al.Simulation of the rolling tex- ment.Acta Mater,2011,59(8):3145 ture of pure Al using crystal plasticity finite element method.J [17]Sun C Y,Guo X R,Huang J,et al.Modelling of plastic deform- Univ Sci Technol Beijing,2007,29(9):920 ation on coupling twinning of single crystal TWIP steel.Acta Met- (皮华春，韩静涛，章传国，等.面心纯铝轧制织构的晶体塑 all Sin,2015,51(3):357 性有限元模拟.北京科技大学学报，2007,29(9)：920) (孙朝阳，郭样如，黄杰，等.耦合李生的TWP钢单品品体 4]Kalidindi S R.Modeling anisotropic strain hardening and deforma- 塑性变形行为模拟研究.金属学报，2015,51(3)：357) tion textures in low stacking fault energy fec metals.Int Plast, [18]Renard K,Jacques P J.On the relationship between work hard- 2001,17(6):837 ening and twinning rate in TWIP steels.Mater Sci Eng A,2012, [5]Barbier D.Favier V,Bolle B.Modeling the deformation textures 542:8 [19]Bouaziz O,Allain S,Scott C.Effect of grain and twin bounda- and microstructural evolutions of a Fe-Mn-C TWIP steel during ries on the hardening mechanisms of twinning-induced plasticity tensile and shear testing.Mater Sci Eng A,2012,540:212 [6]Dancette S,Delannay L,Renard K,et al.Crystal plasticity mod- steels.Scripta Mater,2008,58(6):484 [20] eling of texture development and hardening in TWIP steels.Acta Sun C Y,Huang J,Guo N,et al.A physical constitutive model Matr,2012,60(5):2135 for Fe-22Mn-0.6C TWIP steel based on dislocation density.Ac- [Wang N,Lei L P,Fang G,et al.Deformation analysis of magne- ta Metall Sin,2014,50(9)1115 sium alloy based on crystal plasticity theory.Chin Rare Met, (孙朝阳，黄杰，郭宁，等.基于位错密度的Fe一22Mn0.6C型 2008,32(6):766 TWP钢物理本构模型研究.金属学报，2014,50(9)：1115) (王娜，雷丽萍，方刚，等.镁合金变形的晶体塑性有限元分 221]Barbier D,Gey N,Allain S,et al.Analysis of the tensile behav- 析.稀有金属，2008,32(6)：766) ior of a TWIP steel based on the texture and microstructure evolu- [8]Salem AA,Kalidindi S R,Semiatin S L.Strain hardening due to tions.Mater Sci Eng A,2009,500(1-2)196 deformation twinning in a-titanium:constitutive relations and crys- [22] Grassel O,Kriiger L,Frommeyer G,et al.High strength Fe- tal-plasticity modeling.Acta Mater,2005,53(12):3495 Mn-(Al,Si)TRIP/TWIP steels development-properties-appli- 9]Wu X P,Kalidindi S R,Necker C,et al.Prediction of crystallo- cation.1 nt J Plast,2000,16(10-l1):1391 graphic texture evolution and anisotropic stress-strain curves during 23]Koyama M,SawaguchiT,Lee T,et al.Work hardening associ- large plastic strains in high purity aitanium using a Taylor-type ated with s-martensitic transformation,deformation twinning crystal plasticity model.Acta Mater,2007,55(2):423 and dynamic strain aging in Fe-17Mn-0.6C and Fe-17Mn- [10]Yu Y,Pan XX,Xie RZ,et al.Study on TWIP effect on strain 0.8C TWIP steels.Mater Sci Eng A,2011,528(24)7310

杨 竞等: 耦合孪生晶体塑性模型硬化参数的灵敏度分析 是详细给出了硬化参数确定步骤，即先根据屈服极限 确定初始滑移阻力 s α 0，再根据应变硬化率拐点和阶段 A 分别确定 δ 和孪生硬化指数 b，最后对整个应变硬化 率曲线微调确定饱和参数 Spr ． ( 2) 初始滑移阻力 s α 0 与屈服极限呈线性关系． 孪 生硬化指数 b 增大使得孪生硬化阶段 A 减弱． 孪生阻 力与滑移阻力比值 δ 增大，孪生增长率降低，且硬化率 拐点后移，直至 δ = 1. 2 时拐点消失． ( 3) 获得了 Fe--22Mn--0. 6C 型孪生诱导塑性钢耦 合孪生晶体塑性模型硬化参数的取值范围，其中 s α 0 为 80 ～ 160 MPa 之间，b 为 0 ～ 3 之间，δ 为 1 ～ 1. 3 之间． 参 考 文 献 ［1］ Zhang K S，Zhang G，Feng L． Stress calculation of single crystal under finite slip plastic deformation． Acta Mech Sin，2002，34 ( 4) : 636 ( 张克实，张光，冯露． 单晶体塑性滑移有限变形下的应力计 算． 力学学报，2002，34( 4) : 636) ［2］ Li D Y，Zhang S Ｒ，Peng Y H，et al． Finite element simulation of sheet metal stamping with polycrystalline plasticity． Chin J Mech Eng，2008，44( 1) : 190 ( 李大永，张少睿，彭颖红，等． 板材冲压成形的晶体塑性有 限元模拟． 机械工程学报，2008，44( 1) : 190) ［3］ Pi H C，Han J T，Zhang C G，et al． Simulation of the rolling tex￾ture of pure Al using crystal plasticity finite element method． J Univ Sci Technol Beijing，2007，29( 9) : 920 ( 皮华春，韩静涛，章传国，等． 面心纯铝轧制织构的晶体塑 性有限元模拟． 北京科技大学学报，2007，29( 9) : 920) ［4］ Kalidindi S Ｒ． Modeling anisotropic strain hardening and deforma￾tion textures in low stacking fault energy fcc metals． Int J Plast， 2001，17( 6) : 837 ［5］ Barbier D，Favier V，Bolle B． Modeling the deformation textures and microstructural evolutions of a Fe--Mn--C TWIP steel during tensile and shear testing． Mater Sci Eng A，2012，540: 212 ［6］ Dancette S，Delannay L，Ｒenard K，et al． Crystal plasticity mod￾eling of texture development and hardening in TWIP steels． Acta Mater，2012，60( 5) : 2135 ［7］ Wang N，Lei L P，Fang G，et al． Deformation analysis of magne￾sium alloy based on crystal plasticity theory． Chin J Ｒare Met， 2008，32( 6) : 766 ( 王娜，雷丽萍，方刚，等． 镁合金变形的晶体塑性有限元分 析． 稀有金属，2008，32( 6) : 766) ［8］ Salem A A，Kalidindi S Ｒ，Semiatin S L． Strain hardening due to deformation twinning in α-titanium: constitutive relations and crys￾tal-plasticity modeling． Acta Mater，2005，53( 12) : 3495 ［9］ Wu X P，Kalidindi S Ｒ，Necker C，et al． Prediction of crystallo￾graphic texture evolution and anisotropic stress-strain curves during large plastic strains in high purity α-titanium using a Taylor-type crystal plasticity model． Acta Mater，2007，55( 2) : 423 ［10］ Yu Y，Pan X X，Xie Ｒ Z，et al． Study on TWIP effect on strain hardening behavior in V--5Cr--5Ti alloy． Chin J Appl Mech， 2012，44( 2) : 334 ( 余勇，潘晓霞，谢若泽，等． 孪生诱发塑性对 V--5Cr--5Ti 合 金应变硬化行为的影响． 力学学报，2012，44( 2) : 334) ［11］ Li H W，Yang H，Sun Z C． A robust integration algorithm for implementing rate dependent crystal plasticity into explicit finite element method． Int J Plast，2008，24( 2) : 267 ［12］ Zheng H L，Yang H，Li H W． Crystal plasticity finite element modeling for uniaxial compression of commercially pure titanium． J Plast Eng，2013，20( 1) : 95 ( 郑华雷，杨合，李宏伟． 纯钛压缩变形下的晶体塑性有限 元分析． 塑性工程学报，2013，20( 1) : 95) ［13］ Bouaziz O，Allain S，Scott C P，et al． High manganese austenit￾ic twinning induced plasticity steels: a review of the microstruc￾ture properties relationships． Curr Opin Solid State Mater Sci， 2011，15( 4) : 141 ［14］ Sarma G，Zacharia T． Integration algorithm for modeling the elasto-viscoplastic response of polycrystalline materials． J Mech Phys Solids，1999，47( 6) : 1219 ［15］ Wang W H． Study on Plastic Deformation Behavior Induced by Slipping and Twinning for TWIP Steel ［Dissertation］． Beijing: University of Science and Technology Beijing，2012 ( 王伟华． TWIP 钢滑移和孪生耦合的塑性变形行为［学位论 文］． 北京: 北京科技大学，2012) ［16］ Gebhardt T，Music D，Kossmann D，et al． Elastic properties of fcc Fe--Mn--X ( X = Al，Si) alloys studied by theory and experi￾ment． Acta Mater，2011，59( 8) : 3145 ［17］ Sun C Y，Guo X Ｒ，Huang J，et al． Modelling of plastic deform￾ation on coupling twinning of single crystal TWIP steel． Acta Met￾all Sin，2015，51( 3) : 357 ( 孙朝阳，郭祥如，黄杰，等． 耦合孪生的 TWIP 钢单晶晶体 塑性变形行为模拟研究． 金属学报，2015，51( 3) : 357) ［18］ Ｒenard K，Jacques P J． On the relationship between work hard￾ening and twinning rate in TWIP steels． Mater Sci Eng A，2012， 542: 8 ［19］ Bouaziz O，Allain S，Scott C． Effect of grain and twin bounda￾ries on the hardening mechanisms of twinning-induced plasticity steels． Scripta Mater，2008，58( 6) : 484 ［20］ Sun C Y，Huang J，Guo N，et al． A physical constitutive model for Fe--22Mn--0. 6C TWIP steel based on dislocation density． Ac￾ta Metall Sin，2014，50( 9) : 1115 ( 孙朝阳，黄杰，郭宁，等． 基于位错密度的 Fe--22Mn--0. 6C 型 TWIP 钢物理本构模型研究． 金属学报，2014，50( 9) : 1115) ［21］ Barbier D，Gey N，Allain S，et al． Analysis of the tensile behav￾ior of a TWIP steel based on the texture and microstructure evolu￾tions． Mater Sci Eng A，2009，500( 1 － 2) : 196 ［22］ Grssel O，Krüger L，Frommeyer G，et al． High strength Fe-- Mn--( Al，Si) TＲIP /TWIP steels development--properties--appli￾cation． Int J Plast，2000，16( 10--11) : 1391 ［23］ Koyama M，Sawaguchi T，Lee T，et al． Work hardening associ￾ated with ε － martensitic transformation，deformation twinning and dynamic strain aging in Fe--17Mn--0. 6C and Fe--17Mn-- 0. 8C TWIP steels． Mater Sci Eng A，2011，528( 24) : 7310 · 3801 ·