面向多工序金属板料成形的混合硬化模型*

陈江 陈文亮

(南京航空航天大学 机电学院， 江苏 南京 210016)

金属板料冲压成形广泛应用于各种薄壁零部件的生产加工中.一些汽车结构件曲面多、形状结构复杂，往往需要通过拉延、切边、整形、翻边等多道工序才能得到最终形状，其模具制造的难度大、成本高、周期长.利用计算机进行数值分析，可以预测冲压成形过程中的拉裂、起皱、成形不足等缺陷.采用能够准确描述成形过程中的弹塑性行为、硬化行为以及各向异性的精确的材料模型是得到准确的数值分析结果的关键所在.随动硬化模型和混合硬化模型都考虑了成形过程中后继屈服面相对于初始屈服面的位置变化，也就是背应力张量对回弹结果的影响.很多研究者对这个问题进行了研究，并提出了多种背应力张量的表达公式.Frederick等[1]提出非线性随动硬化模型——Armstrong-Frederick模型，该模型可以表达材料成形过程中的鲍辛格效应中的永久软化，但是在循环拉伸压缩实验中显示，对于鲍辛格效应中的瞬时硬化行为的预测精度与实验数据存在差异.余海燕等[2]基于Chaboche理论的混合硬化模型考虑了不同的预应变水平，同样可以得到与实验结果非常吻合的本构模型，但是其需要标定的参数很多，获取准确的材料参数比较困难.庄京彪等[3]对DC06和DP600在不同预应变水平下做了拉伸实验，利用混合硬化系数来表达鲍辛格效应，但没有说明如何确定混合强化系数.Yoshida等[4- 5]考虑了材料在多次循环拉伸压缩过程中的鲍辛格效应，提出了Yoshida-Uemori(Y-U)混合硬化模型，通过一个U形件的回弹预测实验证明，对于包括高强钢在内的很多材料都能得到较高的回弹预测精度.王梦寒等[6]利用Y-U模型，成功预测了高强钢DP760在U形件成形后的回弹.张志强等[7]利用自制实验装置对TRIP800进行拉伸压缩实验，验证了Y-U混合硬化模型相比等效硬化模型的优势.He等[8]将Y-U用于FLD曲线的预测，取得了良好的效果.Fu等[9]在Y-U混合硬化模型的基础上作了改进，该模型通过对屈服面突变的表达，准确分析了Weldox960的冷弯曲和回弹过程.在金属板料多工序成形过程中，板料的应变路径非常复杂，首先在每一工序都可能经历循环加载、卸载过程，其次在不同工序的载荷加载方向可能会不同，另外先前工序的初始应变应力状态也会影响当前工序.王文平等[10]指出在复杂应力状态下应力张量与应变张量之间的本构关系有别于一般材料本构模型，但并未提出适用于复杂应力状态下的本构模型.

Y-U混合硬化模型在很多金属板料成形的数值分析中得到应用，但是在复杂的多轴向、多路径的多工序金属板料成形的数值分析中并不能得到令人满意的结果.一个原因是其仅允许屈服面的随动硬化，没有考虑屈服面的各向异性对称轴在不同工序的旋转变化；另一个原因是Y-U模型利用材料参数线性地表达其边界面的各向同性硬化，在多工序金属板料成形过程中，材料屈服应力会随着塑性应变的增加而增加，但是先前工序中的应力的变化往往会引起后续工序屈服应力的强化或者弱化.

为了准备表达不同工序的载荷加载方向下的应力应变响应，同时考虑先前工序导致的初始应变的影响，文中在Y-U硬化模型的基础上提出了多工序板料成形混合硬化模型.在经过改进的实验平台上，进行了高强钢DP600的单向拉伸压缩实验以及多轴向的拉伸压缩实验，得到应力应变曲线，然后通过两组实验进行有限元数值模拟.通过比较实验数据，验证多工序混合硬化模型用于多工序板料成形的数值模拟的有效性.此外，在某车型A柱的算例中，对比了 Y-U模型和多工序混合硬化模型的回弹预测精度.

1 初始的Y-U混合硬化模型

Y-U硬化模型为双曲面模型[4]，分别由尺寸大小保持不变而空间位置会改变的屈服面和尺寸大小和空间位置都会改变的边界面组成.屈服面方程f如式(1)所示：

f(σ-α)-Y2=0

(1)

式中，α表示背应力张量，Y表示屈服应力，σ为流动应力.边界面方程F由边界面半径以及背应力张量描述,如式(2)所示：

F(σ)=f(σ-β)-(B+R)=0

(2)

式中：β为偏应力空间中边界面的中心点位置，其初始值为0；B表示边界面的初始尺寸；R表示各向同性硬化分量.

α*表示屈服面中心与边界面中心的相对位移,如式(3)所示，其演化方程如式(4)所示：

α*=α-β

(3)

(4)

a=B+R-Y

(5)

(6)

式中，Rsal为R的饱和值的材料参数，m为控制各向同性硬化速率的材料参数.边界面的背应力的演化方程如式(7)所示：

(7)

式中，b同样为材料参数.

2 多工序混合硬化材料模型

原始的Y-U模型并没有指定如何选择初始屈服面，为了准确表达材料的各向异性特性，本模型选用Barlat-Lian[11]屈服准则定义屈服面，同时为不同工序设置不同的材料参考坐标系来模拟屈服面的旋转.另外，针对不同工序的成形条件，同时考虑先前工序对于当前边界面的背应力的影响，采用新的定义方式可以为不同工序定义不同的硬化系数.

2.1 屈服准则

在本模型中，假定初始的材料参考坐标系C0与初始的各向异性主轴一致，即X轴的方向为材料轧制方向，Y轴为材料所在平面内与X轴垂直方向，Z轴为材料厚度所在方向，那么在i工序的材料参考坐标系Ci如式(8)所示：

Ci=Ci-1Ti

(8)

其中，Ti表示参考坐标系的变换矩阵.假定在当前工序的参考坐标系下，Barlat-Lian[11]屈服准则的定义如式(9)所示：

(9)

其中:α初始值为0;d、c、K1、K2均为各向异性相关的材料参数;n为材料参数，面心立方晶体结构材料的为8，体心立方晶体结构材料的为6;σY为初始屈服应力.

2.2 硬化模型

初始的Y-U模型中采用了等效应变相关的线性方程表达边界面的等效各向同性应力分量,这种表达方式多用于单工序成形分析，在多工序成形分析中往往不够精确.文中采用了一种新的边界面各向同性应力分量的表达公式，如式(10)所示：

Ri=ξiK(ε0+εi)N

(10)

(11)

ai=B+Ri-Y=B-Y+ξiK(ε0+εi)N

(12)

当前工序Ri的演化方程则可以定义为

(13)

利用Ls-Dyna提供的用户自定义材料模型的接口，可以实现多工序混合硬化模型子程序，在完成编译链接后，便能得到新的Ls-Dyna求解器.在多工序金属板料成形过程中，由于先前工序中的应变路径的影响，不同工序的材料硬化曲线并不完全一致，本模型的应变路径影响因子用于表达先前公共应变路径对当前工序的材料模型的强化或者弱化，高强钢DP600材料参数如表1所示，其中r00和r90分别为0°和90°方向的厚向异性系数，在不同应变路径影响因子下的平面拉伸压缩的应力应变响应有所区别，如图1所示.

表 1 DP600材料参数Table 1 Material property of DP600

3 模型验证

3.1 平面拉伸压缩实验及其材料模型参数标定

Yoshida等[5]利用图2所示的标准试样通过单向拉伸压缩实验完成Y-U硬化模型参数的标定.文中为了验证在多工序成形过程中经常会发生的多轴向拉伸的情形，选用了可以进行两个轴向进行拉伸压缩实验的试样，如图3所示.由于薄板在平面压缩过程中容易发生翘曲而失稳，文中在参考Y-U实验装置的基础上，结合实验室已有设备，利用一套槽形板作为夹持装置.该装置结构简单，易于操作，可以安装在试样两侧，在试样拉伸压缩过程中，槽形板发生相对运动，避免试样在压缩过程中失稳，槽形板结构如图4所示.另外，实验时试样与槽形板中间涂上了凡士林以减少摩擦.

图1 多工序混合硬化模型中DP600在不同的应变路径影响因子下的应力应变响应

Fig.1 Strain stress response of DP600 produced by improved Y-U model with various strain path coefficients

图2 Y-U模型中材料参数标定时所用试样(单位：mm)

Fig.2 Specimen used for parameters calibration of original Yosheda-Uemori model(Unit：mm)

图3 验证多工序混合硬化模型采用的试样

Fig.3 Specimen for verification of multi-stagecombined hardening model

图4 槽形板夹具示意图Fig.4 Schematic diagram of channel shape clamp

所取试样为单片厚度为1.0 mm的DP600高强钢，如图4中所示，材料参数如表1所示，沿轧制方向以及90°方向获取.为了得到带有残余应力的试样进行其他轴向的拉伸压缩实验，文中在试样完成轧制方向的拉伸压缩实验后，利用同一片试样再次进行90°方向(横向)的拉伸压缩实验.实验步骤如下：

步骤1 沿轧制方向进行最大塑性应变为0.06的单次循环拉伸压缩实验，即首先在轧制方向拉伸试样直到发生0.06的塑性应变，然后开始卸载，完全卸载后反向加载到-0.06的塑性应变，接着再次加载到应变为0的状态.弹性阶段拉伸速率为0.5 mm/min，塑性阶段以及反向拉伸阶段的拉伸速率为1.5 mm/min.

步骤2 取下试样，再次沿横向进行材料的最大塑性应变0.06的单次循环拉伸压缩实验.

参考Standerm的迭代优化算法[12]，基于实验得到的塑性应力-应变曲线，利用LSTC公司开发的商业有限元分析软件Ls-Dyna以及优化工具包Ls-Opt，可以标定Y-U材料模型参数(如表 2所示)，以及多工序混合硬化材料模型参数(如表3所示).

表2 Y-U模型参数Table 2 Material parameters of Y-U model

表3 多工序混合硬化模型参数Table 3 Material parameters of combined hardening model

3.2 平面拉伸压缩实验的有限元分析

为了简化数值分析，文中采用只包含一个四边形壳单元的有限元模型，如图5所示.壳单元尺寸为100 mm×100 mm的正方形，定义在XY平面内，X方向为轧制方向，厚度为1 mm，厚度方向的积分点个数为7，单元类型选用Belytschko-Tsay.N1、N2、N3、N4为壳单元的4个节点，其中N1为固定节点，N2可以在X方向平移，N4可以在Y方向运动，N3可以在X和Y两个方向平移.首先对N2和N3设置节点运动边界条件，使单元模型在轧制方向进行拉伸压缩，同样对N3和N4设置节点运动边界条件，可以使单元模型在90°方向拉伸压缩.

图5 有限元实验示意图Fig.5 Schematic diagram of FEM experiment

经过二次开发的Ls-Dyna可以支持多工序混合硬化模型，同时Ls-Dyna已经支持包括Y-U模型和Swift模型在内的多种材料本构模型.利用Ls-Dyna求解器进行显示计算，在轧制方向分别选用Y-U模型、多工序混合硬化模型和Swift模型，图6显示了这3种材料模型在试样拉伸压缩过程中的数值分析结果与实验结果的对比.由图6可以看出，试样在轧制方向的拉伸压缩实验中，Y-U模型、多工序混合硬化模型和实验所得的应力应变响应值基本一致，而Swift非线性硬化模型因为不能描述鲍辛格效应，所以在试样压缩阶段与实验结果存在一定偏差.

将轧制方向应该拉伸压缩的试样再次在90°方向拉伸压缩，Y-U模型、多工序混合硬化模型的数值分析的结果以及实验结果如图7所示.由于Swift模型在轧制方向上已经与实验结果存在一定偏差，所以没有进一步被用于90°方向的分析.由图7可以看出，在经过一次拉伸压缩实验后，对于不同轴向的拉伸压缩的应力应变响应值的预测，多工序混合硬化模型的精度更高.从图6和图7的对比可以看出，在轧制方向的拉伸压缩强化了横向的应力应变响应，这一点在Y-U模型和多工序混合硬化模型都能得到体现，但是Y-U模型中预测的应力值偏大.

图6 试样DP600在轧制方向拉伸压缩实验的应力应变响应

Fig.6 Strain stress response of specimen DP600 of uniaxial tension-compression test in roller direction

图7 试样DP600在横向拉伸压缩时的应力应变响应

Fig.7 Strain stress response of specimen DP600 of tension compression test in transverse

4 应用实例

某车型A柱采用DP600高强钢，成形过程的数值模拟需要4工序，分别是OP05(重力分析)、OP20(拉延成形)、OP30(切边冲孔)、OP40(翻边整形)，最后进行回弹分析，如图8所示.采用商业有限元分析软件eta/DYNAFORM进行数值模拟.初始坯料尺寸为2 180 mm×445 mm×1 mm，采用单动拉延模，初始定位于压边圈，压边力为1.95 MN，压边速度为2 m/s，拉延速度为5 m/s，因为坯料表面覆盖了一层湿膜润滑剂，根据经验设定坯料与工具的摩擦系数为0.075.eta/DYNAFORM的自动设置功能结合一种高效灵活的自动定位方法[13- 15]，可以方便地完成多工序成形模拟过程.

图8 A柱多工序成形过程Fig.8 Multi-stage forming process of A-pillar

A柱零件实际成形结果如图9所示，由于是高强钢材料，并且在成形过程中多次被拉延或者弯曲，出现了较大回弹.另外，该零件需要4个工序才最终成形，每一工序都会引起回弹，所以最终回弹量的形成过程非常复杂.为了比较Y-U模型和多工序混合硬化模型在该零件多工序成形过程中回弹预测的精度，文中利用三坐标光学扫描仪对实际零件的外表面进行测量，将测量结果与数值分析结果利用最小二乘法完成最佳匹配后进行对比分析.如图10所示，采用Y-U模型的回弹预测结果的最大误差为1.52 mm，采用多工序混合硬化模型的回弹预测结果的最大误差为0.52 mm.图11为实验结果和两个模型的数值分析结果的一个截面形状，实验结果与两个模型的数值分析结果基本一致，在该截面线的左侧对齐的情况下，由右侧边缘的局部放大图可以看出，Y-U模型的回弹分析结果相比于实际测量结果偏大，多工序混合硬化模型的结果与实验结果更加接近.由此可见，采用Y-U模型和多工序混合硬化模型都能比较准确地预测高强钢DP600的回弹结果，对于多工序成形的回弹预测，多工序混合硬化模型具有更高的精度.

图9 A柱实际成形结果Fig.9 Final forming result of A-pillar

图10 A柱回弹预测偏差分析Fig.10 Springback deviation analysis of A-pillar

图11 A柱截面线上实验与数值分析结果比较

Fig.11 Comparison of measured and calculated section profiles of A-pillar

5 结语

针对多工序板料成形过程中板料会受到多轴向、多种预应变水平、循环加载卸载影响的特点，提出了一种改进的Y-U模型——多工序混合硬化模型.通过实验数据与多种有限元材料模型的比较，证明该模型不仅可以准确描述材料在单向拉伸压缩过程中的鲍辛格效应，而且能准确表达复杂情况下的应力应变响应.准确的材料模型是得到精确的回弹分析的基础，某车型A柱的成形过程的数值分析以及实验结果表明，该模型在多工序成形的数值分析中具有更高的回弹预测精度.

文中通过理论和实验相结合的方法对高强钢DP600的多工序成形提出了一种混合硬化模型，在今后的工作中，一方面可以将本模型应用到高强钢回弹补偿的实例中，另一方面需要将本模型推广应用到其他具有包辛格效应的新材料中，进一步地验证该模型在多工序成形过程中的有效性.

[1] FREDERICK C O,ARMSTRONG P J.A mathematical representation of the multiaxial Bauschinger effect [J].Materials at High Temperatures,2007,24(1)：1- 26.

[2] 余海燕,王友.一种基于Chaboche理论的混合硬化模型及其在回弹仿真中的应用 [J].机械工程学报,2015,51(16)：127- 134.

YU Hai-yan,WANG You.A combined hardening model based on Chaboche theory and its application in the springback simulation [J].Chinese Journal of Mechanical Engineering,2015,51(16)：127- 134.

[3] 庄京彪,刘迪辉,李光耀.基于包辛格效应的回弹仿真分析 [J].机械工程学报,2013,49(22):84- 90.

ZHUANG Jing-biao,LIU Di-hui,LI Guang-yao.Analysis of springback simulation based on Bauschinger effect [J].Chinese Journal of Mechanical Engineering,2013,49(22):84- 90.

[4] YOSHIDA F,UEMORI T.A model of large-strain cyclic plasticity describing the Bauschinger effect and work-hardening stagnation [J].International Journal of Plasti-city,2002,18(5/6):661- 686.

[5] YOSHIDA F,UEMORI T.A model of large-strain cyclic plasticity and its application to spring-back simulation [J].International Journal of Mechanical Sciences,2003,45(10):1687- 1702.

[6] 王梦寒，岳宗敏，王根田.考虑包辛格效应的高强钢U型件冲压回弹规律分析 [J].锻压技术,2016,41(2):58- 63.

WANG Meng-han,YUE Zong-min,WANG Gen-tian.Study on the springback law of high strength steel U-shaped part considering Bauschinger effect [J].Forging & Stamping Technology，2016,41(2):58- 63.

[7] 张志强,贾晓飞,袁秋菊.基于Yoshida-Uemori模型的TRIP800高强钢回弹分析 [J].吉林大学学报(工学版),2015,45(6):1852- 1856.

ZHANG Zhi-qiang,JIA Xiao-fei,YUAN Qiu-ju.Springback analysis of trip high strength steel based on Yoshida-Uemorimodel [J].Journal of Jilin University(Engineering and Technology Edition),2015,45(6):1852- 1856.

[8] HE J,XIA Z C，ZHU X H,et al.Sheet metal forming li-mits under stretch-bending with anisotropic hardening [J].International Journal of Mechanical Sciences,2013,75(10):244- 256.

[9] FU L，DENG C H,REN H L,et al.A modified Yoshida-Uemori constitutive model and its application to cold-bending in Weldox960 [J].Advanced Materials Research,2013,834/835/836:407- 415.

[10] 王文平,刁可山,吴向东,等.板料屈服行为及强化规律的研究进展 [J].机械工程学报，2013,49(24):7- 14.

WANG Wen-ping,DIAO Ke-shan,WU Xiang-dong,et al.Review on yielding and hardening behavior of sheet metal [J].Chinese Journal of Mechanical Engineering,2013,49(24):7- 14.

[11] BARLAT F,LIAN J.Plastic behavior and stretchability of sheet metals(Part I):a yield function for orthotropic sheet under plane stress condition [J].International Journal of Plasticity,1989,5(1):51- 66.

[12] CRAIG K J,STANDER N.On the robustness of a simple domain reduction scheme for simulation-based optimization [J].Engineering Computations,2002,19(3):431- 450.

[13] 陈江,陈文亮,鲍益东.基于约束投影的多工序板料成形自动定位 [J].华南理工大学学报(自然科学版),2016，44(2):53- 59.

CHEN Jiang,CHEN Wen-liang,BAO Yi-dong.An auto-position method of multi-stage sheet metal forming based on constrained projection [J].Journal of South China University of Technology(Natural Science Edition),2016，44(2):53- 59.

[14] 杜国康，陈文亮，鲍益东.基于近似板料构形的全工序工具定位方法 [J].机械科学与技术，2014，33(4)：522- 526.

DU Guo-kang,CHEN Wen-liang,BAO Yi-dong.Method of tools’ positioning in complete processes based on approximate blank configuration [J]. Mechanical Science and Technology,2014，33(4)：522- 526.

[15] CHEN Jiang,CHEN Wenliang,DU Guokang.An approach to predict blank shape for auto-position in multistage sheet metal forming [C]∥Proceedings of NUMISHEET 2014.Melbourne:AIP,2014:1028- 1031.