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

工程科学学报，第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 (<0或f>) 第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 ·