基于GTN细观损伤模型的激光拼焊板成形极限预测

马向东　官英平

燕山大学先进锻压成形技术与科学教育部重点实验室,秦皇岛,066004

基于GTN细观损伤模型的激光拼焊板成形极限预测

马向东官英平

燕山大学先进锻压成形技术与科学教育部重点实验室,秦皇岛,066004

建立了从细观损伤角度预测拼焊板成形极限的Gurson-Tvergaard-Needleman(GTN)损伤模型,用有限元逆向法确定了损伤模型中的各损伤参数。采用有限元软件ABAQUS耦合基于Mises屈服准则的弹塑性GTN损伤模型，对拼焊板半球凸模胀形过程进行了数值模拟。设计了拼焊板半球凸模胀形物理试验，试验过程中通过改变试件的宽度得到了不同应变状态下完整的拼焊板成形极限图，并与GTN细观损伤模型预测到的拼焊板成形极限图进行对比分析，验证了GTN细观损伤模型预测拼焊板成形极限图的准确性。

拼焊板；GTN模型；成形极限图；半球凸模胀形

0　引言

拼焊板(tailor-welded blanks,TWB)成形技术是采用不同厚度、不同材质或不同表面涂层的板材，经下料后焊接成整体预制板坯，再进行冲压成形的工艺方法[1]。拼焊板成形技术在减轻车身质量、增强车身结构安全性等方面具有重要作用[2]。近年来，随着社会对环境保护、节能减排等方面的要求越来越高,拼焊板成形技术在汽车制造业中的应用越来越广泛。但是,拼焊板母材或厚度不同或强度不同，且存在焊缝和热影响区，导致拼焊板在成形技术上与常规单一板材相比有很大不同,有很多新的技术问题亟待解决，如焊缝移动和成形极限降低等问题[3-4]。

成形极限图用来表征弹塑性材料发生塑性变形时所能达到最大变形的能力，是衡量板材成形性能最直接、最有效的方法之一。Mohebbi等[5]基于修正的M-K理论分别对拼焊板薄厚两侧板料的成形极限进行了预测。Korouyeh等[6]通过物理实验与理论计算相结合的方法研究了板料厚度比对拼焊板成形极限的影响,建立了拼焊板极限板料厚度比的理论计算模型。Lee等[7]基于Yld2000和Cazacu各向异性屈服准则，并结合半球凸模胀形实验对五种不同类型拼焊板的成形性能进行了研究。以上研究均从宏观角度出发，以应变成形极限图来衡量拼焊板的成形性能，而应变成形极限图与加载路径有很大的相关性[8-9],实验得到的近线性路径下拼焊板的应变成形极限图用于复杂加载路径时会存在较大误差。Arrieux[10-11]通过实验得到了单向拉伸和双向拉伸预加载条件下板材的应力成形极限图，实验结果表明,应力成形极限图与加载路径无关,不同应变路径下的应力极限曲线是唯一的。这一结论后来也得到了Matin等[12]和Panich等[13]的验证。Bandyopadhyay等[14]通过研究发现，对于复杂应变状态下的拼焊板成形，用应力成形极限图来描述拼焊板的成形性能要比用应变成形极限图来描述拼焊板的成形性能更为准确。但是,应力成形极限图仅能通过实验极限应变和塑性理论结合屈服准则间接求得，无法通过物理实验方法直接获取，其结果的准确程度与极限应变的测量精度有很大关系。

Gurson-Tvergaard-Needleman(GTN)细观损伤模型[15-17]是一种从细观损伤力学角度出发，描述板材成形性能最常用、最有效的方法之一,该模型不仅能较好地预测金属板材的成形极限,而且与应变加载路径无关,此外,对于塑性较差或者断裂时没有明显颈缩现象的金属板材,GTN细观损伤模型具有显著的优势。Liu等[18]基于GTN损伤模型,采用数值模拟与半球凸模胀形实验相结合的方法，对AA5052/聚乙烯/AA5052复合层板的成形极限图进行了预测。Chen等[19]建立了基于塑性各向异性的GTN细观损伤模型，并将该模型应用到铝合金板材的成形中。上述研究均以塑性较差或者成形极限较低的复合板材或者合金板材为研究对象。拼焊板由于其母材厚度或者强度不同，且存在焊缝和热影响区，从而导致拼焊板的成形性能和成形极限较差，而目前从细观损伤角度研究拼焊板成形极限的研究还很少。

本文在上述研究的基础上,首先建立了GTN细观损伤模型,通过有限元逆向法确定GTN模型中的各损伤参数，使其能很好地描述拼焊板的损伤演化过程；其次，基于有限元数值模拟软件ABAQUS，耦合GTN损伤模型对拼焊板半球凸模胀形过程进行数值模拟；最后，设计了获取拼焊板成形极限图的半球凸模胀形物理试验，并将该试验得到的拼焊板成形极限图与GTN细观损伤模型预测到的拼焊板成形极限图进行了对比分析。

1　GTN细观损伤模型

材料的宏观断裂与细观损伤之间有着本质联系。在金属塑性加工领域，材料多以塑性破坏为主，且内部本身存在一定的缺陷。如图1所示，金属材料损伤演化过程大致包括以下四个步骤:①材料本身内部缺陷引起孔洞形核(图1a)；②随着塑性变形孔洞长大(图1b)；③相邻孔洞聚合形成微裂纹(图1c);④微裂纹扩展导致宏观断裂(图1d)。

(a)步骤①　(b)步骤②(c)步骤③(d)步骤④图1　金属材料损伤演化过程示意图

GTN损伤模型可用来描述基体材料在塑性变形过程中由于孔洞的形核、长大、聚合而导致的宏观破裂现象。与传统塑性变形过程中体积不变这一假设相比,GTN损伤模型考虑了宏观静水应力对屈服面的影响，其具体表达式如下[15-17]：

[1+q3(f*)2]=0

(1)

当q1=q2=q3=1时,GTN损伤模型退化为最初的Gurson损伤模型。当f*=0(初始孔洞体积分数f0=0)时,材料没有损伤,GTN损伤模型退化为标准的Mises屈服函数。损伤变量f*可表示为

(2)

式中，fc为孔洞开始聚合时的临界孔洞体积分数;κ为孔洞长大加速因子;ff为材料最终破裂时的孔洞体积分数。

(3)

基体材料总的损伤演过包括两部分:原有孔洞长大和新孔洞形核引起的孔洞体积分数的变化,即

df=dfgrowth+dfnucleation

(4)

假设基体材料是不可压缩的，根据质量守恒定律,原有孔洞长大的变化量dfgrowth依赖于宏观的塑性体积变形，即与静水压力分量有关：

dfgrowth=(1-f)dεp∶I

(5)

式中，I为二阶单位张量。

孔洞形核的变化量dfnucleation既可由应力控制，也可由应变控制。本文采用由应变控制的形核准则，其表达式如下：

(6)

(7)

式中,A为应变控制的孔洞形核系数;fN为形核粒子的孔洞体积分数;εN为孔洞形核时的平均塑性应变;SN为形核应变的标准偏差。

GTN损伤模型中,材料力学性能参数和损伤参数的准确确定是合理描述材料行为的前提。其中，材料的基本力学性能参数，如弹性模量E、应变硬化指数n、各向异性系数r、屈服强度σs和抗拉强度σb等，可通过单向拉伸试验获得。而损伤参数的确定目前还没有一套完整的普遍接受的方法。

2　拼焊板制备及母材力学性能和损伤参数的确定

本文选取厚度h分别为0.8mm、1.2mm的ST12钢板作为拼焊板母材进行研究。ST12钢板的化学成分如表1所示。拼焊板制备在型号为HC-AW300的脉冲激光焊接机上进行(图2),焊接电流为130A,激光频率为20Hz,脉宽为4.5ms,焊接速度为120mm/min。拼焊板焊接方向均与母材轧制方向平行。

拼焊板母材的材料力学性能参数可通过单向拉伸试验获取。沿板材轧制方向、与板材轧制方向成45°方向和垂直于轧制方向三个方向用线切割切取单向拉伸试样,其几何尺寸参照GB/T228-2002的规定。通过单向拉伸试验得到的拼焊板母材的力学性能参数如表2所示。

表1　拼焊板母材化学成分(质量分数)　%

图2　拼焊板激光焊接过程

表2　拼焊板母材力学性能参数

确定GTN损伤参数的常用方法大致可归纳为金相法、代表性体积单元法和有限元逆向法[20-22]。金相法得到的实际材料的孔洞大小并不能完全等同于GTN损伤模型中的孔洞体积分数参数。代表性体积单元法确定GTN损伤模型中的相关形核参数时需给定受力条件，过程比较复杂。有限元逆向法是结合数值模拟和实验数据，采用最优化的参数识别法来确定各参数的。本文采用有限元软件ABAQUS分别对拼焊板薄厚两侧母材的单向拉伸过程进行数值模拟，通过调整各损伤参数，使数值模拟得到的单向拉伸名义应力-名义应变曲线(或者力-行程曲线)与实验测得的曲线相吻合，从而确定各损伤参数。

图3　单向拉伸有限元模型

图4　单向拉伸试验与模拟的名义应力-名义应变曲线(h=1.2 mm)

根据单向拉伸试样的对称性，取试样四分之一建立有限元模型，如图3所示。单向拉伸数值模拟过程用ABAQUS/Explicit动力显式算法求解器计算，试样的基本力学性能参数见表2，单元类型采用考虑沙漏控制的8节点线性减缩积分体单元C3D8R,厚度方向设置2层单元。多次调整损伤参数得到拼焊板母材ST12(h=1.2 mm)的名义应力-名义应变曲线如图4所示。与试验结果对比发现，数值模拟模型及损伤参数可以很好地描述试验过程。通过有限元逆向法确定的拼焊板薄厚两侧母材的GTN各损伤参数如表3所示。由表2和表3可以看出，两种厚度的ST12拼焊板母材虽然均产自同一厂家，但是其力学性能参数和GTN损伤参数仍有一定的差异。为了体现损伤模型应用于预测拼焊板成形极限分析时所特有的性质，下文在进行拼焊板半球凸模胀形过程数值模拟时，将拼焊板薄厚两侧母材的力学性能参数和GTN损伤参数分别赋予有限元模型。

表3　拼焊板母材的GTN损伤参数

3　获取成形极限图的物理试验和数值模拟

3.1半球凸模胀形

本文通过半球凸模胀形物理试验来获取拼焊板及其母材的成形极限图。半球凸模胀形试验在型号为BCS-50A的板材成形性能试验机上进行，如图5所示。半球凸模胀形的模具结构示意图见图6。考虑到拼焊板母材有厚度差且胀形过程中伴随着焊缝移动，压边圈采用阶梯式设计，且在焊缝处留有一定的缝隙，防止厚板向薄侧移动时，与薄侧压边圈侧壁挤压导致焊缝处非正常破裂。试验过程中， 凸模速度为35 mm/min， 压边力为110 kN,板料与模具之间采用聚乙烯薄膜包裹凡士林的润滑方式。为获取不同应变状态下拼焊板及其母材的极限应变，试验板材的长度保持180 mm不变,改变其宽度,具体尺寸如图7所示，拼焊板焊缝方向与面内主应力方向垂直。试验前，在被测板材的表面用电腐蚀打标机印制2 mm×2 mm的正方形网格，板材胀形完成后，用网格应变仪测量临近破裂区域网格的面内主次应变即可获取板材在该应变状态下的极限应变。连接各应变状态下的极限应变点即可得到完整的成形极限图。

图5　BCS-50A板材成形性能试验机

3.2有限元数值模拟

利用有限元分析软件ABAQUS,基于物理试验装置建立拼焊板半球凸模胀形的有限元模型。考虑到拼焊板的对称性，取二分之一模型建模，如图8所示。胀形数值模拟过程用ABAQUS/Explicit动力显式算法求解器计算。拼焊板母材的基本力学性能参数和GTN损伤参数见表2和表3，单元类型采用考虑沙漏控制的8节点线性减缩积分体单元C3D8R,厚侧板料的厚度方向设置3层单元，薄侧板料的厚度方向设置2层单元，焊缝处做刚性连接处理，且焊缝方向垂直于主应力方向。板料与模具之间的接触方式为面面接触，摩擦服从库仑摩擦定律，摩擦因数取0.125，与物理试验的摩擦因数基本一致。当试样的孔洞体积分数达到破裂孔洞体积分数ff时，则认为被测试样破裂，测得此时临近破裂区域的面内主次应变即可得到成形极限图。

图8　拼焊板半球凸模胀形有限元模型

4　结果与讨论

图9　胀形破裂后的试件形状

图10　拼焊板及其母材成形极限图

通过半球凸模胀形物理试验得到胀形破裂后的试件形状如图9所示。用网格应变仪测量破裂临近区域的应变即可获得拼焊板及其母材的成形极限图，如图10所示。由图10可以看出,拼焊板薄侧母板的成形极限略高于厚侧母板的成形极限,差别并不大,但是二者的成形极限均高于拼焊板的成形极限,特别是在拉-拉变形区即成形极限图的右半部分,薄厚两侧母板的成形极限要明显大于拼焊板的成形极限。这是因为本文所研究的两种厚度的ST12拼焊板母材均产自同一厂家，且它们的基本力学性能参数差距很小，所以成形性能差别也不大。而拼焊板由于焊缝及热影响区的存在，且薄侧母材与厚侧母材相比，其承载能力较低，胀形过程中厚侧母材基本不发生塑性变形，或者仅有较小的塑性变形，变形主要集中在拼焊板薄侧母材，从而导致胀形提前破裂。这也进一步解释了拼焊板与母材相比，成形性能降低的原因。

以宽度为180 mm的拼焊板为例，当胀形高度为25 mm时，通过数值模拟得到的孔洞体积分数云图如图11所示。从图11中可以看出，空洞体积分数最大的位置在薄侧母材临近焊缝处，此处也是胀形时拼焊板破裂的起源。另外,除法兰区由于压边导致拼焊板的孔洞体积分数有所增大外,孔洞体积分数发生明显增大的区域集中在拼焊板与胀形凸模接触的顶端，且厚侧板料的孔洞体积分数明显小于薄侧的孔洞体积分数,这也说明了在胀形过程中,薄侧母材与厚侧母材相比，其承载能力较低，成形过程中发生了较大的塑性变形，容易率先发生破裂。

图11　孔洞体积分数分布云图

图12所示为胀形高度为25 mm时,宽度为180 mm拼焊板及其母材的孔洞体积分数随焊缝垂直距离(或穹顶中心距离)的变化曲线。从图12中可以看出，当胀形高度为25 mm时,孔洞体积分数在拼焊板薄侧板材距离焊缝2 mm处达到最大值,最大值为0.06,此处也是胀形时拼焊板的破裂位置。拼焊板厚侧板材的孔洞体积分数随着焊缝距离的增大而减小,但最大值不超过0.03,与破裂孔洞体积分数(0.05)相差较远,胀形过程中发生较小的塑性变形，不易破裂。两种厚度的拼焊板母材的孔洞体积分数随着穹顶中心距离的增大而减小，且孔洞体积分数均小于破裂孔洞体积分数(0.05和0.06),这也说明了当胀形高度为25 mm时,拼焊板已经发生破裂，而其母材仅发生均匀塑性变形，并没有明显局部厚度减薄或者破裂，这也进一步说明了拼焊板与其母材相比，成形性能明显降低。

图12　孔洞体积分数随焊缝距离的变化曲线

基于GTN细观损伤模型,通过有限元数值模拟预测得到拼焊板的成形极限图见图13。从图13中可以看出,预测得到的成形极限图与物理试验所获得的成形极限图吻合较好，特别是在拉-压变形区(成形极限图的左半部分)。在拉-拉变形区(成形极限图的右半部分),预测结果与试验结果有一定的差距，最大相对误差为7.7%,这是因为胀形过程中当拼焊板处于双拉应变状态时，破裂位置离焊缝较近，应变状态受到焊缝的影响而产生了漂移。

图13　拼焊板成形极限图数值模拟预测与物理试验结果对比

5　结论

(1)建立了预测拼焊板成形极限的GTN细观损伤模型，通过有限元逆向法确定了GTN模型中拼焊板薄厚两侧母材的损伤参数，该模型及损伤参数能很好地描述拼焊板的损伤演化过程。

(2)拼焊板薄侧母材的成形极限略高于厚侧母材的成形极限，但差别并不大;拼焊板整板的成形极限明显低于其母材的成形极限,这也是拼焊板与母材相比，成形性能降低的原因。

(3)基于GTN细观损伤模型,通过有限元数值模拟预测到拼焊板的成形极限图与半球凸模胀形物理试验得到的拼焊板成形极限图吻合很好，虽然在拉-拉变形区由于受到焊缝牵制，使应变路径产生了漂移，成形极限图的预测误差稍大，但最大相对误差也仅为7.7%。

[1]Merklein M,Johannes M,Lechner M,et al.A Review on Tailored Blanks-Production,Application and Evaluation[J].Journal of Materials Processing Technology,2014,214(2):151-164.

[2]Abbasi M,Ketabchi M,Ramazani A,et al.Investigation into the Effects of Weld Zone and Geometric Discontinuity on the Formability Reduction of Tailor Welded Blanks[J].Computational Materials Science,2012,59:158-164.

[3]朱元右,姜银方,杨继昌,等.拼焊板冲压成形中的焊缝移动分析[J].中国机械工程,2006,17(15):1618-1622.

Zhu Yuanyou,Jiang Yinfang,Yang Jichang,et al.Analysis of Weld-line Movement in Tailor-welded Blanks’ Formation[J].China Mechanical Engineering,2006,17(15):1618-1622.

[4]宋燕利,华林.车身覆盖件拼焊板冲压成形技术的研究现状及发展趋势[J].中国机械工程,2011,22(1):111-118.

Song Yanli,Hua Lin.Current Research Status and Trends of Tailor Welded Blanks Applied in Auto-body Stamping[J].China Mechanical Engineering,2011,22(1):111-118.

[5]Mohebbi M S,Akbarzadeh A.Prediction of Formability of Tailor Welded Blanks by Modification of MK Model[J]. International Journal of Mechanical Sciences,2012,64(1):44-51.

[6]Korouyeh R S,Naeini H M,Torkamany M J,et al.Experimental and Theoretical Investigation of Thickness Ratio Effect on the Formability of Tailor Welded Blank[J].Optics & Laser Technology,2013,51:24-31.

[7]Lee W,Chung K H,Kim D,et al.Experimental and Numerical Study on Formability of Friction Stir Welded TWB Sheets Based on Hemispherical Dome Stretch Tests[J].International Journal of Plasticity,2009,25(9):1626-1654.

[8]Chu C C.An Investigation of the Strain Path Dependence of the Forming Limit Curve[J]. International Journal of Solids Structure, 1982, 18(3):205-215.

[9]Zhalehfar F,Hashemi R,Hosseinipour S J.Experimental and Theoretical Investigation of Strain Path Change Effect on Forming Limit Diagram of AA5083[J].International Journal of Advanced Manufacturing Technology, 2015, 76(5/8):1343-1352.

[10]Arrieux R.Determination and Use of the Forming Limit Stress Diagrams in Sheet Metal Forming[J].Journal of Materials Processing Technology,1995,53(1/2):47-56.

[11]Arrieux R.Determination and Use of the Forming Limit Stress Surface of Orthotropic Sheets[J].Journal of Materials Processing Technology,1997,64(1/3):25-32.

[12]Matin P H,Smith L M,Petrusevski S.A Method for Stress Space Forming Limit Diagram Construction for Aluminum Alloys[J].Journal of Materials Processing Technology,2006,174(1/3):258-265.

[13]Panich S,Barlat F,Uthaisangsuk V,et al.Experimental and Theoretical Formability Analysis Using Strain and Stress Based Forming Limit Diagram for Advanced High Strength Steels[J].Materials and Design,2013,51:756-766.

[14]Bandyopadhyay K,Basak S,Panda S K,et al.Use of Stress Based Forming Limit Diagram to Predict Formability in Two-stage Forming of Tailor Welded Blanks[J].Materials and Design,2015,67:558-570.

[15]Gurson A L.Continuum Theory of Ductile Rupture by Void Nucleation and Growth:Part I-Yield Criteria and Flow Rules for Porous Ductile Media[J].Journal of Engineering Materials and Technology,1977,99(1):2-15.

[16]Tvergaard V,Needleman A.Analysis of the Cup-cone Fracture in a Round Tensile Bar[J].Acta Metallurgica,1984,32(1):157-169.

[17]Needleman A,Tvergaard V.An Analysis of Ductile Rupture in Notched Bars[J].Journal of the Mechanics and Physics of Solids,1984,32(6):461-490.

[18]Liu J G, Liu W,Xue W.Forming Limit Diagram Prediction of AA5052/Polyethylene/AA5052 Sandwich Sheets[J].Materials and Design,2013,46:112-120.

[19]Chen Z Y,Dong X H.The GTN Damage Model Based on Hill’48 Anisotropic Yield Criterion and Its Application in Sheet Metal Forming[J]. Computational Materials Science,2009, 44(3):1013-1021.

[20]Mahnken R.Aspects on the Finite-element Implementation of the Gurson Model Including Parameter Identification[J].International Journal of Plasticity,1999,15(11):1111-1137.

[21]Sun D Z,Siegele D,Voss B,et al.Application of Local Damage Models to the Numerical Analysis of Ductile Rupture[J].Fatigue & Fracture of Engineering Materials & Structures,1989,12(3):201-212.

[22]Steglich D,Siegmund T,Brocks W.Micromechanical Modeling of Damage Due to Particle Cracking in Reinforced Metals[J].Computational Materials Science,1999,16(1/4):404-413.

(编辑王艳丽)

Forming Limit Prediction of Tailor-welded Blanks Based on GTN Damage Model

Ma XiangdongGuan Yingping

Key Laboratory of Advanced Forging & Stamping Technology and Science(Yanshan University),Ministry of Education of China,Qinhuangdao,Hebei,066004

The GTN damage model was established to predict the forming limit of TWBs,and the damage parameters were determined by the use of finite element inverse method. Hemispherical punch bulging test process was simulated by the finite element software ABAQUS and coupled the Mises yield criterion with GTN damage model. Experiments of hemispherical punch bulging test were designed,which changed specimen width to obtain complete forming limit diagram of TWBs in different strain conditions.Meanwhile, the prediction and experimental results were compared to validate the accuracy of GTN damage model.

tailor-welded blanks(TWB);Gurson-Tvergaard-Needleman(GTN) model;forming limit diagram;hemispherical punch bulging test

2015-06-30

国家自然科学基金资助项目(51275444);河北省钢铁联合研究基金资助项目(E2014203271);高等学校博士学科点专项科研基金资助项目(20121333110003)

TG386DOI：10.3969/j.issn.1004-132X.2015.22.013

马向东,男,1986年生。燕山大学机械工程学院博士研究生。主要研究方向为板材成形。发表论文5篇。官英平(通信作者),男,1963年生。燕山大学机械工程学院教授、博士研究生导师。