李军超, 谢锋, 赵泽, 龚鹏程
重庆大学 材料科学与工程学院, 重庆 400044
收稿日期：2018-01-09
基金项目：国家自然科学基金资助项目(51575066); “控制科学与工程”浙江省重中之重学科开放基金资助项目。
作者简介：李军超(1979-)，男，河南平顶山人，重庆大学副教授。

摘要：针对板材渐进成形破裂预测对应变路径的过分依赖且无法实时预测的问题, 首先, 运用FORTRAN语言, 通过ABAQUS有限元软件的材料子程序接口, 将于忠奇破裂准则引入DC56D+Z钢板的VUMAT材料子程序中；其次, 分别将于忠奇模型与ABAQUS自带Von Mises模型进行单向拉伸和渐进成形的模拟；最后, 结合渐进成形实验与有限元模拟分析, 以实验结果为标准, 根据有限元模拟的结果, 逆向寻找出于忠奇准则下渐进成形模型的临界破裂积分值I.研究结果表明:两种有限元模型的应力应变大小数量级相同, 分布一致, 从而说明了该子程序的有效性; 渐进成形模拟过程中最大破裂积分值I出现的位置在零件的侧壁, 与实验结果一致, 积分值I=17可以作为预测渐进成形过程中板材是否破裂的有效条件.
关键词：渐进成形有限元分析于忠奇破裂准则临界破裂积分值破裂预测
Numerical Simulation and Fracture Prediction of Incremental Sheet Forming of Metals
LI Jun-chao, XIE Feng, ZHAO Ze, GONG Peng-cheng
College of Materials Science and Engineering, Chongqing University, Chongqing 400044, China
Corresponding author: LI Jun-chao, E-mail: ljchust@163.com
Abstract: In order to solve the problems in the fracture prediction of incremental sheet forming of metals, in which the prediction strongly depends on the strain path and a real-time prediction cannot be achieved, a numerical model has been built. To perform that, the Yu Zhongqi ductile fracture criterion was firstly added into ABAQUS software through the interface subroutine of the VUMAT for DC56D+Z steel. Then, the uniaxial tensile and incremental forming processes were simulated by the new user-defined model and ABAQUS built-in Von Mises model, respectively. Finally, according to the results of finite element simulation, the critical fracture integral value I of the incremental forming model under the Yu Zhongqi ductile fracture criterion can be reduced reversely based on the actual experimental results. The results show that the calculated stress and strain from the two models have the same order of magnitude and similar distribution, which prove the validity of the subroutine. The maximum fracture integral value I in the simulated sheet forming process appears on the side wall of the workpiece, which is consistent with the experimental ones. The critical fracture integral value I=17 can be used as an effective way to predict the fracture in the incremental sheet metal forming process.
Key words: incremental formingfinite element analysisYu Zhongqi ductile fracture criterioncritical fracture integral valuefracture prediction
金属板料数控渐进成形工艺是一种基于计算机、数控和塑性成形技术基础之上的先进柔性制造工艺.它的特点是引入快速原型制造技术“分层制造”的思想, 将复杂的三维形状沿z轴方向离散化, 实现了设计与制造一体化的柔性快速制造.这种对板料进行分层渐进成形的方法, 无需模具, 零件的形状和结构也不受约束^[1].在小批量生产、产品定制、缩短研发周期和生产成本等方面具有优越性^[2-3].目前, 对于汽车和航空航天等制造业中要求产品多品种、小批量、快速更新以及其他领域中薄壳类新产品的开发具有巨大的应用价值^[4].
自渐进成形技术问世以来, 板材减薄及预测一直是研究的重难点.周六如^[5]通过实验, 验证了渐进成形材料板厚变化符合正弦定律, Hussain等^[6]基于正弦定律提出了一种测试板材减薄极限的新方法, 并指出成形极限与零件特征和工艺参数有密切联系.Ambrogio等^[7]用板材成形过程中的成形力的大小来预测板材是否断裂.Petek等^[8]同样运用成形力的大小作为预测板材破裂发生的时间和位置的判断依据.但是在不同应力状态下, 相同材料的极限应力也会不同, 且只能通过破裂位置的极限应变间接得到极限应力, 因此用该方法预测破裂十分困难.Hussain等^[9]建立了渐近成形过程的经验成形极限图(FLD), 并验证了当层间距和工具头直径在一定范围内时, 所建FLD图的有效性.但由于渐近成形过程中弯曲和剪切变形的存在, FLD图只适用于特定零件, 具有一定局限性.李磊等^[10]运用Oyane破裂准则证明渐进成形的成形极限远高于传统板材成形工艺.于忠奇等^[11]成功将韧性破裂准则引入到数值模拟中, 并准确地预测了铝合金板材的成形极限.但该方法需要先输出易破裂区单元的历史应力应变信息, 而后进行计算判断, 且要求易破裂区位置容易估计, 因而对板材的破裂预测具有一定的局限性.目前, 用于渐进成形方面的破裂准则相关研究不够深入, 可用于指导生产的相关破裂准则少之又少.因此, 找到较为准确的预测渐进成形破裂情况的破裂准则并建立相关数值模型至关重要.本文基于ABAQUS平台, 开发子程序并引入于忠奇韧性断裂准则, 对板料渐进成形过程的破裂行为进行数值模拟分析; 将渐进成形实验与有限元模拟相结合, 运用逆向寻找的方法, 确定临界韧性断裂积分值.
1 材料本构模型实验所用的材料为DC56D+Z钢板, 厚度为0.8 mm.由于该型钢板对应变速率及温度不敏感, 根据文献[12]将其本构关系描述为

(1)

其中:σ为等效塑性应力; ε为等效塑性应变; a, b, c为待求材料常数.
本文通过单向拉伸实验来确定式(1)中的常数以及DC56D+Z钢板的力学性能.实验过程中, 为确保测量材料性能参数的准确性, 将板料按与轧制方向成0, 45°, 90°三个方向取样, 每组3个试样, 分别进行室温单向拉伸实验.试样尺寸如图 1所示.
图 1(Fig. 1)

图 1 单向拉伸实验试样尺寸(mm)Fig.1 Test sample size for a uniaxial tensile test

依据实验得到的各项结果, 并通过得到的真实应力-应变曲线拟合得到其本构模型中的相关参数.得到的DC56D+Z的本构模型如式(2)所示:

(2)

根据计算, DC56D+Z材料模型的其他性能参数如表 1所示.
表 1(Table 1)

表 1 材料DC56D+Z相关力学性能参数Table 1 Mechanical parameters of DC56D+Z steel sheet 弹性模量 屈服强度 密度 硬化指数 泊松比 GPa MPa kg·m-3 207 136.19 7 850 0.24 0.28

表 1 材料DC56D+Z相关力学性能参数 Table 1 Mechanical parameters of DC56D+Z steel sheet

2 于忠奇韧性断裂准则于忠奇韧性断裂准则^[13]如式(3)所示:

(3)

式中:σ_m, σ分别为平均应力和等效应力; ε, ε_f分别为等效应变和断裂产生处等效应变; p, C为材料常数.在应用韧性断裂准则预测金属板材在成形过程中的破裂时, 必须先确定材料常数p, C，对于不同的材料、不同的加载条件, 准则中的材料常数也不尽相同.为减小误差, 本文采用单向拉伸和平面应变拉伸实验所得的破裂时的极限应变, 再结合适合各向异性的Hill屈服准则求解出于忠奇韧性断裂准则中的材料常数p, C.
平面应变拉伸实验选用厚度为0.8 mm的DC56D+Z钢板为研究对象, 实验设备及条件均与上文单向拉伸实验一致, 且同单向拉伸实验一样将板材取样方向按与轧制方向0, 45°, 90°三个方向取样, 每组3个试样.平面应变拉伸试样尺寸如图 2所示.
图 2(Fig. 2)

图 2 平面应变拉伸试样(mm)Fig.2 Plane strain tensile specimen

根据实验的结果, 取真实应力最大时刻的应变为破裂时的极限应变, 并对3组试样结果取平均值作为最终结果.实验得到单向拉伸时的极限应变为0.217, 平面应变拉伸时的极限应变为0.428.
获得破裂时的极限应变后, 继而解决式(3)中的项的求解.根据Hill屈服准则, 在不考虑平面各向异性差别的情况下, 应力三轴度和等效应变与单向拉伸第一主应变的比值可分别简化为式(4)和式(5).
单向拉伸状态时:

(4)

平面拉伸状态时:

(5)

式(4)与式(5)中的r为塑性应变比, 其值可由式(6)求得.

(6)

式中:b₀, b分别代表试样变形前后的宽度; t₀, t分别代表试样变形前后的厚度.
结合式(3)~式(5), 计算得到于忠奇韧性断裂准则厚度为0.8 mm的DC56D+Z板材的材料常数C=0.2, p=0.006.
3 于忠奇韧性断裂准则子程序模型的建立与验证3.1 于忠奇韧性断裂准则子程序模型的建立本文利用ABAQUS软件进行渐近成形过程模拟.为了使用于忠奇韧性断裂准则实时准确地预测渐近成形过程的破裂情况, 并实现板材破裂的可视化, 本文通过ABAQUS的VUMAT接口编写于忠奇韧性断裂准则材料模型, 进而实现该软件的二次开发.
将式(3)的于忠奇韧性断裂准则变为式(7)的形式:

(7)

式中, I为准则的积分值, 当该值大于某一临界值时, 材料将发生破裂.
由于有限元软件可以得到变形单元每步对应的应力应变值, 故考虑采用梯形法计算式(8)的积分值.梯形法的计算原理如下:

(8)

对应于忠奇准则, 式(8)中f(x)即为将式(7)的于忠奇准则写成梯形法的形式, 如式(9)所示:

(9)

接下来便可用Fortran语言将式(9)编入子程序中完成于忠奇准则的编写.设置积分值I为子程序中判定材料是否破裂的标准.为了实现材料破裂的可视化, 在子程序中指定控制单元删除标志的状态变量, 当该变量值为1时, 材料点为激活状态; 若某单元的破裂积分值大于临界值时, 赋值该状态变量为0, 材料点便被删除, 从而实现板料破裂的可视化.
3.2 子程序模型的验证为了验证开发的VUMAT子程序的正确性, 利用ABAQUS软件对板料进行单向拉伸实验模拟.模拟分为两组进行.第一组模拟调用含有于忠奇准则的VUMAT用户子程序.为了与实际情况相符, 将程序中判断板料破裂的临界积分值I设为1^[10].第二组采用软件自带的Mises模型.通过对比两组模型拉伸模拟发生缩颈时的应力应变, 验证于忠奇模型的有效性.图 3为模拟后的应力分布云图.可见, Mises模型模拟后的应力区间为138.7~485.6 MPa, 于忠奇准则模拟后的应力区间为138.1~493.8 MPa.同样的, 在图 4中, Mises模型的应变区间为0 ~0.525 4, 于忠奇准则的应变区间为0.006 478~0.520 6.可见两种模型的应变大小数量级相同, 分布一致, 这说明该子程序编写的VUMAT子程序能够有效地分析DC56D+Z钢板的变形规律.
图 3(Fig. 3)

图 3 单向拉伸模拟应力分布云图Fig.3 Simulated stress contour of the sample in uniaxial tensile test (a)—Mises模型结果; (b)—于忠奇准则模型结果.

图 4(Fig. 4)

图 4 单向拉伸模拟应变分布云图Fig.4 Simulated strain contour of the sample in uniaxial tensile test (a)—Mises模型结果; (b)—于忠奇准则模型结果.

观察图 3与图 4, 两组模拟的应力应变分布规律相同, 最大值均出现在板料的中间发生颈缩处, 与实际相吻合.总体来说, 于忠奇准则模型子程序对单向拉伸实验模拟结果与ABAQUS自带Mises准则模拟结果一致, 同时二者与实际实验结果相差微小, 由此可证明编写的于忠奇模型子程序是有效的.
3.3 子程序模型在渐进成形过程中的验证验证子程序有效后, 应进一步在渐近成形过程中去验证其有效性, 本文采用尺寸为160 mm×140 mm×0.8 mm的DC56D+Z钢板作为成形坯料, 成形如图 5所示的圆锥制件.
图 5(Fig. 5)

图 5 需成形制件尺寸(mm)Fig.5 Size of the target workpiece

设置的工艺参数如下:工具头半径2.5 mm, 层间距1 mm, 上下压边圈的压边力为2 000 N, 工具头加工速度为2 500 mm/s, 板料与工具头之间的摩擦系数为0.05, 板料与其他部件间的摩擦系数为0.25.此外, 设置模拟的放大系数为1 500.渐进成形模拟的模型装配图如图 6所示.
图 6(Fig. 6)

图 6 板材数控渐进成形CAE模型Fig.6 CAE model for NC incremental sheet forming for metals

Mises模型与于忠奇准则及程序模型模拟结果如图 7, 图 8所示.由图 7可见, Mises模型计算出的应力大小区间分布为13.81~661.9 MPa; 于忠奇准则模型计算出的应力分布区间为11.27~691.2 MPa.两组模拟的应力结果数量级相同, 大小相近.同时, 两组模拟结果的应力分布结果大致相同.
图 7(Fig. 7)

图 7 渐进成形模拟应力分布云图Fig.7 Calculated stress contour in incremental sheet metal forming simulation (a)—Mises模型结果; (b)—于忠奇准则模型结果.

图 8(Fig. 8)

图 8 渐进成形模拟应变分布云图Fig.8 Calculated strain contour in incremental sheet metal forming simulation (a)—Mises模型结果; (b)—于忠奇准则模型结果.

两组模拟的应变结果如图 8所示.Mises模型所得最大应变为2.155, 于忠奇准则模型所得最大应变为2.436.二者的差距在可接受的范围之内.此外, 与应力分布相同的, 二者的应变分布结果也一致.成形后的板料顶部与法兰部分未参与变形, 应变均为0, 与实际相符.
4 于忠奇韧性断裂准则对渐进成形过程破裂的预测本文采用渐进成形实验测定出破裂时的成形高度, 根据有限元模拟出破裂高度所在位置处的应变，应变值逆向求解得出破裂时的临界积分值I.
4.1 实验参数为了寻找破裂积分值I, 本文设计了3组模拟, 并有3组工艺参数相同的实验与之对应.每组的工艺参数如表 2所示.
表 2(Table 2)

表 2 工艺参数Table 2 Parameters in processing 方案 层间距p/mm 工具头半径r/mm 成形角α/(°) 加工速度v/(mm·s-1) 1 80 2 0.5 4 85 2 000 3 90

表 2 工艺参数 Table 2 Parameters in processing

实验目标制件形状及尺寸, 除成形角外, 与上节模拟目标圆锥制件一致.本文实验设备为重庆大学NH3532三轴数控渐进成形机.
4.2 实验和模拟结果分析与讨论实际所得制件和相同参数下的3组模拟如图 9所示.可见, 每个模拟的结果有一个共同点, 破裂积分值I从上至下大体呈现先增大后减小的趋势, 这与渐进成形过程中板料厚度减薄率的变化情况相同^[14].
图 9(Fig. 9)

图 9 实验和有限元模拟结果对比Fig.9 Comparison between experimental and finite element simulation results (a)—成形角=80°; (b)—成形角α=85°; (c)—成形角α=90°.

各个模拟的最大破裂积分值I所在单元对应的Z向坐标值及同参数下加工制件破裂位置对应加工高度如表 3所示.由表中数据模拟所得最大破裂积分值I对应坐标与实际发生破裂的位置基本相符.说明模拟结果是符合实际的, 模拟所得的最大破裂积分值即为临界破裂积分值.分析各组模拟所得最大I, 判断临界破裂积分值在17左右.考虑到实际破裂位置为制件破裂后一段时间所测得, 并且数值模拟所得值为单元瞬时值, 故模拟值应略高于实际值, 从而确定临界破裂积分值I=17能作为有效地预测渐进成形过程中板材破裂是否发生的判据.
表 3(Table 3)

表 3 模拟结果与实验制件破裂位置相关数据比较Table 3 Comparison of the correlation data between the simulation results and the fracture position from the experiment 变量 80° 85° 90° 实际加工高度/mm 8 7.5 6 数值模拟Imax出现位置/mm 7.94 7.67 6.23 最大破裂积分值 16.30 16.92 18.44

表 3 模拟结果与实验制件破裂位置相关数据比较 Table 3 Comparison of the correlation data between the simulation results and the fracture position from the experiment

5 结论1) 结合单向拉伸实验和平面应变拉伸实验确定了于忠奇韧性破裂准则中的材料常数值为:C=0.2, p=0.006.
2) 基于ABAQUS平台建立了于忠奇模型子程序, 并通过采用该子程序模型和软件自带的Von Mises模型模拟单向拉伸的变形过程, 对比模拟结果与实验所得数据, 验证了子程序的有效性.
3) 通过模拟板材渐进成形过程, 发现两种模型计算出的应力应变一致, 且发现其他参数一定的情况下, 随着成形深度的增加, 破裂积分值I先增加后减小, 与实际相符, 进而证明了于忠奇模型子程序在模拟渐进成形过程中的有效性.
4) 通过多组实验和模拟分析发现:渐进成形过程中的最大破裂积分值I出现在零件侧壁, 且破裂位置对应的最大破裂积分值I均与17接近, 故可认为I=17为判断板材渐进成形过程中是否破裂的临界条件.
参考文献

[1]	莫健华, 叶春生, 黄树槐. 金属板料数控无模成形及快速制模[J].电加工与模具, 2002(1): 15–18. ( Mo Jian-hua, Ye Chun-sheng, Huang Shu-huai. Numerical control die forming and rapid die making of metal sheet[J].Electromachining & Mould, 2002(1): 15–18.DOI:10.3969/j.issn.1009-279X.2002.01.005)
[2]	Martins P A F, Bay N, Skjoedt M, et al. Theory of single point incremental forming[J].CIRP Annals-Manufacturing Technology, 2008, 57(1): 247–252.DOI:10.1016/j.cirp.2008.03.047
[3]	Nimbalkar D H, Nandedkar V M. Review of incremental forming of sheet metal components[J].International Journal of Engineering Research & Applications, 2013, 3(5): 23–28.
[4]	李燕乐, 陈晓晓, 李方义, 等. 金属板材数控渐进成形工艺的研究进展[J].精密工程, 2017(1): 1–9. ( Li Yan-le, Chen Xiao-xiao, Li Fang-yi, et al. Research progress of numerical control incremental forming of sheet metal[J].Precision Engineering, 2017(1): 1–9.)
[5]	周六如. 板料数控渐进成形变形区厚度变化规律研究[J].机械工程学报, 2011, 47(18): 50–54. ( Zhou Liu-ru. Study on variation of thickness of asymmetrical deformation zone of sheet metal[J].Chinese Journal of Mechanical Engineering, 2011, 47(18): 50–54.)
[6]	Hussain G, Gao L. A novel method to test the thinning limits of sheet metals in negative incremental forming[J].International Journal of Machine Tools and Manufacture, 2009, 47(3/4): 419–435.
[7]	Ambrogio G, Filice L, Micari F. A force measuring based strategy for failure prevention in incremental forming[J].Journal of Materials Processing Technology, 2006, 177(1/2/3): 413–416.
[8]	Petek A, Kuzman K, Suhac B. Autonomous online system for fracture identification at incremental sheet forming[J].CIRP Annals-Manufacturing Technology, 2011, 60(1): 283–286.DOI:10.1016/j.cirp.2011.03.146
[9]	Hussain G, Gao L, Hayat N, et al. A new formability indicator in single point incremental forming[J].Journal of Materials Processing Technology, 2009, 209: 4237–4242.DOI:10.1016/j.jmatprotec.2008.11.024
[10]	李磊, 周晚林, HussainG. 金属板料单点渐进成形极限的数值模拟预测[J].机械工程学报, 2010, 46(18): 102–107. ( Li Lei, Zhou Wan-lin, Hussain G. Numerical simulation of single point incremental forming limit of sheet metal[J].Journal of Mechanical Engineering, 2010, 46(18): 102–107.)
[11]	于忠奇, 杨玉英, 王永志, 等. 基于韧性断裂准则的铝合金板材成形极限预测[J].中国有色金属学报, 2003, 13(5): 1223–1226. ( Yu Zhong-qi, Yang Yu-ying, Wang Yong-zhi, et al. Prediction of forming limit of aluminum alloy sheet based on ductile fracture criterion[J].The Chinese Journal of Nonferrous Metals, 2003, 13(5): 1223–1226.DOI:10.3321/j.issn:1004-0609.2003.05.034)
[12]	刘春达, 张旭, 刘健. 板材数控渐进成形过程数值分析与实验研究[J].锻压技术, 2012, 37(5): 35–39. ( Liu Chun-da, Zhang Xu, Liu Jian. Numerical analysis and experimental study on numerical control progressive forming process of sheet metal[J].Forging & Stamping Technology, 2012, 37(5): 35–39.DOI:10.3969/j.issn.1000-3940.2012.05.008)
[13]	于忠奇.基于Lemaitre损伤理论的韧性断裂准则建立及板料成形极限预测[D].哈尔滨: 哈尔滨工业大学, 2003: 21-28. ( Yu Zhong-qi.Establishment of ductile fracture criterion based on Lemaitre damage theory and prediction of sheet metal forming limit [D].Harbin: Harbin Institute of Technology, 2003: 21-28.http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=Y551986)
[14]	陶龙, 王进, 姜虎森. 单点渐进成形时工艺参数对成形能力的影响[J].锻压技术, 2012, 37(3): 19–22. ( Tao Long, Wang Jin, Jiang Hu-shen. Influence of process parameters on forming ability in single point progressive forming[J].Forging & Stamping Technology, 2012, 37(3): 19–22.DOI:10.3969/j.issn.1000-3940.2012.03.005)