0 引言
在冲压成形过程中，由于板料与模具之间的摩擦或自身塑性变形，使得金属工件上的温度会有一定程度的升高，与周围环境相比有明显的温差[1]。常用的低碳钢塑性较好，变形过程中内应力较小，导致同等应变下产生的温升很小，对板料的影响几乎可忽略不计。但是随着先进高强钢的广泛应用，由于高强钢的屈服强度高，在变形过程中应力变化大，金属内部晶格会产生较大塑性应变能[1]，进而导致金属板料的温度有大幅度升高，且增幅为传统低碳钢的几倍到十几倍，故有必要考虑冲压过程中的自发热问题。
冲压过程中板料局部温度升高，局部区域的材料塑性增强，会直接导致同等条件下材料内应力变小。而随着材料属性的变化，局部区域的发热量也随之相对动态变化。传统冷冲压数值模拟无法考虑变形过程中的温升问题，复杂的自发热问题的介入必然导致其模拟结果与实际冲压结果存在较大差异，而考虑自发热效应的完全动态耦合问题研究能够较为准确地计算每一个时间历程下的温度变化情况，必然会使数值模拟结果更为准确。
目前，更多的研究集中在依靠外部加热的温热成形研究，且主要集中于先进超高强钢热成形和镁铝合金温成形，而对于双相钢的温成形和成形过程中的自发热问题研究较少。CAI等[2]利用温拉伸试验得到了DP590在不同温度下的多向应力变化关系，建立了与温度相关的屈服准则；SUNG等[3]基于Hollomon-Voce模型对比研究了DP590、DP780和DP980钢的温变属性。对于变形热和摩擦热，虽然早在1978年AUTILA等[4]就通过红外测温仪器发现高强钢金属试样拉伸断裂时的温度变化超过38℃，但后续针对变形热和摩擦热的研究并不完整。文献[5-8]研究表明，变形做功和摩擦做功产生的能量并不是完全转化为热量，因此在计算时往往需要增加热效应转换系数进行调节。JUNGHWAN等[9]利用高速摄影机和红外监测仪测量汽车侧围件冲压时的温升，并建立了热力耦合模型进行对比，但是没有考虑与温度相关的真应力模型，其结果存在可探讨性。
本文选取双相钢DP780钢板为研究对象，利用马弗炉和热拉伸设备进行了温拉伸试验，针对自发热效应问题，建立了温度区间为20～140 ℃的DP780动态塑性本构方程，并且基于此本构方程，采用完全热力耦合数值模拟方法，提取壳单元建立冲压模具与板料的接触模型，应用显式向前差分时间整合算法，建立了考虑变形热和摩擦热的新式冲压成形研究方法，将热效应计算模型应用于U形件的回弹数值模拟中，并将该新式研究方法模拟结果、传统冷冲压模拟结果以及试验试制结果进行对比，验证了考虑自发热效应的新式冲压成形研究方法的准确性。
1 变形热及摩擦热理论
1.1 变形热机理及计算
1.1.1 应变能产生机理
当板料发生冲压变形时，原子平面相互滑动而产生微观摩擦热，原子摩擦热逐渐积累，进而产生不可忽略的温升。在这个运动过程中，原子晶格以及一些缺陷相互滑动，改变它们的相对位置，正是这种错位导致了塑性变形，而且在塑性变形过程中常伴随着金属的温升。
金属均匀变形过程中，按照增量理论，总应变分为弹性应变和塑性应变：
(1)
式中，dεij为某一节点应变总增量；
为单节点塑性应变增量；
为单节点弹性应变增量。
按照外力加载的做功条件和定积分原理，应变能分为弹性应变能和塑性应变能，可得到塑性应变能的表示方式为
(2)
式中： Wp为塑性应变能；W为总应变能；We为弹性应变能；σij为节点的应力值。
通常通过拉伸试验可得到金属某一状态下的真应力-应变曲线。在几何学上，总应变能即真应力曲线下包络的面积如图1所示。其中，弹性应变能在外载荷卸载之后会由于弹性回复而抵消，塑性应变能则会转化为金相组织转化的潜热以及金属局部温升，由于本文考虑环境温度不超过150 ℃，板料几乎无相变，因此潜热带来的影响可忽略不计。
图1 单向均匀拉伸应变能分布
Fig.1 Uniform tensile strain energy distribution
1.1.2 变形热转化机理
变形过程中伴随着热传导、热辐射和热弹性效应，变形热即塑性应变能转变的热量。与弹性应变功相反，塑性应变功在卸载后会保留在工件中。塑性应变功的85%～95%会转变为变形热，这使得工件的温度在塑性变形过程中升高。
对于变形热带来的温升，有
(3)
式中，J为机械热当量；ρ为金属密度；Cp为常量压力下的质量热容；T为温度；
为有效塑性应变；σ为应力值；E为弹性模量；λ为热耗散系数；α为热膨胀系数；
为与变形热相关的Taylor-Quinney热效应系数；
为有效弹性应变率。
实际上，塑性成形过程在绝热、温度多变还是等温方式下发生，材料的温升有着显著的不同，因此式(3)引用了变形热效应系数
作为调节热量转换的重要参数。
根据与位错密度相关的Zehnder模型，基于金属流变应力曲线满足硬化准则的前提，可得到
(4)
式中，n为硬化指数；ε0为初始塑性应变值。
考虑到工程应用和该系数变化区间极小，参照文献[10]，根据经验定义后续仿真中变形热效应系数
在塑性变形过程中，由于向工具和环境导出的热量会使工件的温升减小(在多变关系中)，变形速率越大，导出的这部分热量越少，因此大多数情况下，如果
则可忽略式(3)左侧的导热部分，即近似忽略热导率带来的影响。
通常，热弹性效应带来的温度影响|ΔTe|<0.5 ℃对结果影响极小，故对式(3)右侧潜热部分取零简化处理。
在近似绝热条件下，式(3)可简化为
(5)
其中，式(5)左侧代表热效应带来的温升变化，右侧代表塑性应变功做功效率。
采用宏观的温升速率无法计算微小步长下的温升变化，定义等效塑性应变差值Δε=εi+1-εi极小时，推导出瞬间局部温升为
(6)
1.1.3 变形热应用机理
某一节点下一时刻的温度可表示为
(7)
由式(7)可以看出，某一节点下一时刻的温度不仅与当前时刻的温度相关，而且与当前温度下的金属塑性本构关系以及应力、应变有关。
考虑到温度和应变速率的影响，传统硬化准则不能够满足拟合要求，因此本文根据试验曲线和误差估计，采用Norton-Hoff模型对双相钢温变塑性本构模型进行描述：
(8)
n=n0exp(-Cn(T-T0))
m=m0exp(-Cm(T-T0))
式中，
为应变率；C为温度相关系数；m为应变速率敏感系数；K为硬化指数模量；Cm、Cn分别为与应变速率敏感系数相关的系数和与硬化指数相关的系数；T0为室温。
将式(8)代入式(7)中，可得到单节点某时刻宏观温升与应变的计算方程：
(9)
整理得
(10)
定义
是只与材料属性和应变属性相关的系数)，则
(11)
由式(11)可见，在只考虑变形热的情况下，某节点的瞬间温升与对应步长下的变形增量和温度直接相关。
1.2 摩擦热机理及计算
1.2.1 摩擦热产生机理
金属拉延时，不但有变形热产生，还存在大量的摩擦热。通常板料边由压板压住，部分零件需要增加拉延筋以控制板料的流入速度，压边力和拉延筋的存在都会导致板料与模具摩擦力增大，经拉延产生摩擦热，导致摩擦后的板料存在一定的温升。
传统的冷冲压仿真并未考虑板料与模具摩擦发热的问题，一方面是因为存在技术难度，需要针对不同材料发热和散热机理进行描述，另一方面，低碳钢硬度和屈服强度都较低，摩擦热带来的影响几乎可忽略不计，但高强钢屈服强度较大，不易变形，成形过程中压机压力较大，局部产生极大的接触应力，所以生产过程中，摩擦热效应影响较大。在连续生产中，变形热和摩擦热效应带来的温升可达180 ℃[11]，温度与材料本身相互影响，必将带来计算结果的差异。本文将摩擦做功作为工具与接触区中的热源来考虑。
模具与板料接触面以相对速度v滑动时，如果不考虑磨损的影响，滑动干摩擦接触面上产生的总摩擦功为
(12)
式中,μi为当前时刻的摩擦因数；r为微元体的位置；t为摩擦时间；p为局部接触压力；v为相对移动速度；A为接触面积。
在数值模拟计算中，需要依据瞬态能量转化率及瞬态温升进行热量相关的计算，故当Δt=ti+1-ti极小时，考虑摩擦热引起的热流密度简化计算公式为[12]
(13)
式中，η为依赖于边界条件的来自摩擦能量的热转换率；τ为瞬态摩擦应力；Δs为时间步长的增量滑移。
在本文的研究中，研究对象均为双相钢冲压钢板，研究目标一致，故取
与变形热效应系数对应，以保证数值模拟的简便与准确性。
传统仿真过程的冲压摩擦因数常设置为固定值，大小一般为0.15～0.3，但研究表明：随着压力的变化，摩擦因数也存在着不同程度的变化[13-14]，另外，摩擦因数也随温度的升高而增大[15]。实际工作条件下，金属局部区域的摩擦力生成极其复杂，如U形件的冲压，在进料的圆角区域，板料与模具上模圆角充分接触，成形力的存在导致对应位置的接触压力较大，摩擦因数与温升和应力变化相互影响。
本文结合与温度和压力相关的摩擦因数研究，将其应用于考虑摩擦热的冲压过程。
1.2.2 变摩擦力模型
目前，用于CAE仿真模拟的传统摩擦力模型包括库仑准则和恒定摩擦因数模型。当接触压力较小时，模拟结果准确性可满足要求，但高强钢成形过程中接触和摩擦行为复杂，传统的恒定摩擦因数模型不再适用，有必要建立更加符合实际的变摩擦因数模型。
赵玉璋等[13]考虑板料边界润滑条件下接触界面载荷对摩擦因数的影响，采用如下变摩擦因数模型进行表述：
μ(p)=μ0(p/p0)np-1
(14)
式中，p0为参考载荷；μ0为p0下所对应的摩擦因数；np为模型指数(p0>0，0.5≤np≤1)。
采用MMW-1型立式万能销-盘式摩擦磨损试验机测量摩擦因数，所用试样为DP780裸钢板，摩擦副采用冷模具钢DC53，试验结果如图2所示，图中θ=T-273 ℃。采用最小二乘法进行拟合，得到适用于变摩擦因数模型的参数np=0.793,p0=2.05,μ0=0.125。
图2 不同温度下的变摩擦因数模型
Fig.2 Variable friction coefficient model at
different temperature
另外，笔者基于不同温度下的热摩擦试验，建立了DP780在20～160 ℃温度下的温变摩擦模型[16]：
μ(θ)=aμ(T0)(T/T0)nT+b
(15)
式中， μ(T0)为室温下测得摩擦因数；a、b分别为拟合参数；nT为模型指数。
通过MG-2000型销-盘式摩擦试验机得到不同压力下的摩擦因数，进行最小二乘拟合后，得到a=0.042 89，nT=1.697，b=0.123 1，拟合曲线见图3。
图3 不同压力下的变摩擦因数模型
Fig.3 Variable friction coefficient model at
different pressures
为了更加准确地描述双相钢的摩擦因数，综合式(14)与式(15)的变摩擦因数模型，建立双因素变摩擦因数简易模型如下：
μ(p,θ)=aμ0(p/p0)np-1(θ/θ0)nT+b
(16)
冲压过程中，只有局部(拉延筋、成形圆角等)位置存在摩擦力的畸变，其他区域的摩擦力与压边力直接相关。本文考虑在得到室温下的原始摩擦因数后，参考式(16)建立的温变摩擦模型进行冲压过程的摩擦描述，计算每一个步长下的摩擦因数随温升和压力变化的动态数值。
由于摩擦因数与压力和温度相关，对式(13)进行修正，摩擦热效应热流密度计算公式为
(17)
1.3 热耗散与热传递机理
在笛卡儿坐标系描述下，对于连续均匀的介质，三维非稳态导热微分方程可由下式描述：
(18)
式中,k为材料的热导率；
为内热源率。
非稳态热传导问题中温度与时间t有关，导热微分方程需要一定的条件才能进行求解，主要包括初始条件和边界条件。式(18)要求有一个定解则需要一定的初始条件，即
T(x,y,z,t=0)=T0(x,y,z)
(19)
式中，T0(x,y,z)为时间为零时所规定的的温度分布。
目前主要应用的边界条件有三类：
(1)第一类边界条件是指某给定表面上质点的温度值在传热过程中保持不变，即温度值为给定的边界条件，设该表面为S1，则表示为
T(x,y,z,t)=T0 (t>0,S∈S1)
(20)
(2)若物体表面S2给定热流密度q，则称之为第二类边界条件，即
(21)
式中，χ表示S2表面上任意点的外法线方向。
当模具与板料接触并且有相对位移时，就会导致热流密度增大，如果材料相同，则生成的摩擦热会有一半进入板料中，但往往模具材料与板料材料的属性存在一定的区别。
通常热流密度分配计算公式为
(22)
式中，c1、c2，ρ1、ρ2，k1、k2分别为板料和模具的质量热容、密度和热导率。
通过材料属性值可计算得到：q1/q2=1.76。
(3)第三类边界条件。板料与模具温度升高后，高于室温的边界面与空气接触，即开始与环境对流换热，即
(23)
式中，h为导热物体边界面与周围环境流体件的表面传热系数；T∞为环境温度；As为板料与空气接触面积。
冲压成形过程中自由表面上的边界条件为第三类边界条件，与模具接触的表面同时存在第二类和第三类边界条件。
1.4 数值模拟应用方法
对于传热问题来说，根据变分原理，传热方程如下：
(24)
Q=Qp+Qc+Qd+Qf
(25)
式中，Kconv为该节点的传热系数矩阵；T和
分别为节点的温度矢量和节点的温度变化率矢量，为总体热传导矩阵；C为总体热容矩阵；Q为总体热流矢量；Qp、Qc、Qd、Qf分别为变形热源热流矢量、对流和辐射损失热流矢量、模具接触传导热损失热流矢量和摩擦生热热流矢量。
求解热流问题的差分形式，t时刻温度Tt与t+Δt时刻的温度Tt+Δt有如下关系：
(26)
将式(24)左右两边同时乘以1-φ和φ并相加，结合式(26)消去对时间的导数
得到
(27)
式中，Tt、TΔt+t分别为t时刻和t+Δt时刻热流矢量。
可以由t=0时刻的已知T0，用式(27)求出TΔt，然后由TΔt计算出T2Δt，依此类推，解决温度场耦合问题。
由此推出考虑自发热效应的热力耦合方法：①给定初始温度场T0；②计算对应的速度场
③更新节点坐标和有关变形场量；④由步骤②，利用式(27)求解本步的温度场Tt+Δt；⑤在Tt+Δt温度状态下求解本步的速度解
⑥重复③～⑤步，直到变形结束。
依据热力耦合算法，应用显式向前差分时间整合算法，把传热方程集成于显式中心差分整合算法的力学求解中，可得到考虑变形热和摩擦热的求解思路如图4所示。
图4 考虑变形热和摩擦热的求解思路
Fig.4 Solutions considering deformation heat and
friction heat
2 双相钢热拉伸试验
2.1 热拉伸试验过程
试验材料选用双相钢DP780，其成分如表1所示。利用江苏大学MT5105微机控制电子万能试验机(图5)进行恒温条件下的热拉伸试验。
表1 DP780化学成分(质量分数)
Tab.1 Chemical composition of DP780
(mass fraction) %
质量分数w(C)w(Si)w(Mn)w(P)w(S)w(Al)数值0.090.192.140.0120.0050.041
图5 MT5105微机控制电子万能试验机
Fig.5 MT5105 microcomputer controlled electronic
universal testing machine
试验过程如下：①将图6所示试样安装在试验机夹具上，微调上下端距离进行预紧，合上马弗炉；②分别设置温度θ为20 ℃、70 ℃、140 ℃，进行加热，根据马弗炉中上、中、下3个温敏传感器监测内部温升情况；③设置0.1 s-1、0.01 s-1、0.001 s-13个不同的应变率进行拉伸，直至样件拉断。
图6 拉伸试样尺寸
Fig.6 Size of the tensile specimen
在恒温马弗炉中，炉壁、夹具、内部空气、试样均加热到相应温度，且本次试验拉伸率并非高应变速率，试样有充分的时间将变形热传递到炉内的夹具和空气中，故近似认为拉伸过程为恒温过程。
2.2 热拉伸试验结果
在获取试验的原始力-位移曲线后，按照国标要求进行数据处理，通过热拉伸试验得到20～140 ℃温度范围内不同应变率下的真应力-应变曲线(图7～图9)。由图7可以看出，随着应变率的增大，真应力曲线也呈现上升趋势，材料塑性降低，断裂伸长量也对应减小；而随着温度的升高，材料塑性明显增强，应力曲线呈现下降趋势，但是在该区间内的延伸率却呈现减小的趋势，该现象不同于文献[12，15，17]中铝合金中的流变应力曲线趋势，考虑可能是双相钢含碳量高，低温状态下，温度升高，碳元素导致金属脆性增强。
图7 θ=20 ℃条件下真应力-真应变曲线
Fig.7 True stress-strain curve at 20 ℃
图8 θ=70 ℃条件下真应力-真应变曲线
Fig.8 True stress-strain curve at 70 ℃
图9 θ=140 ℃条件下真应力-真应变曲线
Fig.9 True stress-strain curve at 140 ℃
2.3 考虑温度的塑性本构方程的求解
采用Norton-Hoff模型对材料曲线进行拟合，对式(8)左右两边同时取对数，得到
(28)
利用MATLAB对
等曲线进行数据拟合，求解各个参数，对DP780试样在20～140 ℃范围内、应变率在0.001～0.1 s-1范围内的应力-应变曲线进行描述，最终拟合的本构方程为
n=0.215 4exp(-0.000 9(T-T0))
m=60.021 3exp(-0.007 9(T-T0))
3 数值模拟研究
3.1 变形热效应数值模拟研究
3.1.1 变形热数值模拟过程
为研究变形热的产生对试样性能的影响程度，采用完全热力耦合方式，同时对拉伸过程中的力和热进行耦合求解。将试样模型导入ABAQUS软件中，左端固定，右端施加强制位移，拉伸直至试样断裂，其中考虑变形发热及传热、辐射等能量传导。网格属性定义为C3D8RT，能够较为准确地计算三维模型的热力耦合结果。材料属性定义如表2所示，考虑到此后的冲压需要用到断裂界定准则，选择使用成形极限图(forming limit diagrams，FLD)破裂准则作为仿真破裂准则，数据参考文献[18]中的DP780不同温度下的FLD模型。
表2 材料参数
Tab.2 Material parameters
密度(kg/m3)弹性模量(GPa)热导率(W/(m·K))质量热容(J/(kg·K))膨胀系数(l/K)7 900210524801.9×10-5
3.1.2 变形热效应结果讨论
考虑变形热效应的拉伸最终结果如图10和图11所示，最大应力点位于试样中心位置，这与试验结果相符；最大温升也处于中间区域，最大温度可达48 ℃，这与文献[4]所监测到的温度在量级上是一致的。
图10 断裂前一刻的应力分布
Fig.10 Stress distribution before fracture
图11 断裂前一刻的温度场分布
Fig.11 Temperature distribution before fracture
考虑到变形热效应系数大小的影响，设置不同变形热效应系数对结果进行对比，如图12所示。图中列出了不同热效应系数β下的温升情况以及考虑热效应影响的应力应变曲线对比。由图12发现，随着热效应系数β的减小，温升也存在一定程度的减小，符合实际生产情况，但是当热效应系数差值
时，断裂时刻的温升变化ΔT<3 ℃，对试样属性影响极小，应力曲线基本重合，此处不再加以讨论，并规定后续
图12 试样中心点的真应力-真应变曲线
与本构曲线对比
Fig.12 Comparison of true stress-strain curve and
constitutive curve at the center of the specimen
另外，由图12可以看出，不考虑变形热拉伸结果与室温(20 ℃)下试验结果近似相同(图中重合)，在考虑变形热之后，弹性变形阶段因温升极小，拉伸结果与室温下相同；到达塑性阶段后，随着温度升高，应力曲线逐渐低于室温下的应力曲线，且越来越低，这是由于受变形发热量影响，塑性变形区域温度不断升高，金属塑性变形不断增强，曲线中的应力值随着应变值的增大越来越小，最终表现为图12中的应力状态曲线。
3.2 热效应数值模拟研究
3.2.1 摩擦热数值模拟过程
为模拟研究钢板从压边面和下模之间抽出时的温度变化情况，首先进行简单的摩擦热力耦合模拟确定相关参数，并且研究摩擦热变化情况。
在上下模具压住钢板之后，施加一定的压边力，将钢板从模具(模具设为刚体，无形变)之间抽出，抽出的钢板均存在不同大小的温升，如图13所示。板料材料参数按照表2输入，模具材料参数见表3，均采用C3D8RT网格形式。为避免变形热的交叉影响，将板料设置成线弹性体。
(a)t=0.025 s (b)t=0.05 s
(c)t=0.075 s (d) t=0.1 s
图13 试验过程中温度的变化
Fig.13 Temperature change during the test
表3 模具参数设置
Tab.3 Mold parameter settings
密度(kg/m3)弹性模量(GPa)热导率(W/(m·K))质量热容(kJ/(kg·K))膨胀系数(K-1)7 4002100.0134601.2×10-5
3.2.2 摩擦热效应结果讨论
选取板料中间线作为选取线，以板料右端作为基准，板料右端被模具压住，从右至左逐步发生位移，以距离右端的距离L作为参考依据，L越小的区域，摩擦时间越长，产热越多。
图14 板料选取点温度随时间的变化
Fig.14 Temperature change with time at the
point of sheet
由图14可以看出：①选取点的温度随着时间延长而上升，说明摩擦热逐渐产生，且升温速度近似相同；②距离右端大于或等于10 mm的选取点，温度变化会出现一个平台，不再上升，这是由于超过10 mm的选取点，在拉伸一段时间后，会被拉出模具的压边范围，不再有摩擦热产生，又因为设置时间极短，温度散失较小，所以会出现一个趋于等温的平台。
将冷摩擦模型下的钢板拉出模拟实验结果与温变摩擦模型下的实验结果进行对比，模拟实验采用强制位移的方式拉动钢板，过程中位移变化均匀，选择拉出钢板所需要的自反力作为衡量标准，对比两种摩擦模型下的自反力结果，如图15所示。
图15 不同摩擦模型下的自反力
Fig.15 Reflexive forces under different
friction models
由图15可以看出：在0.005～0.05 s期间，自反力增大到一定值，两种模型下的自反力均呈现波动式变化，且温变摩擦模型下的自反力增长更快；在0.05～0.1 s期间，常温摩擦模型下的自反力曲线趋于平稳，而温变摩擦模型下的自反力曲线却仍在上升，且波动较大，这是由于温度升高导致摩擦力增大，相同应变下，产生的应力更大，另外塑性的不同使得各个区域的工况变得更加复杂，导致出现选点的应力变化随之出现较大跳动变化。
3.3 考虑热效应的冲压成形仿真
3.3.1 新式冲压成形方法
在研究变形热和摩擦热各自对金属变形过程的影响之后，现同时考虑两种热效应影响，采用完全热力耦合研究方法，将其应用于简单U形件的冲压数值模拟分析中，如图16所示。由于U形件左右对称，为加快计算速度，选用一半的板料以及模具进行计算，得到冲压结果后，根据成形结果计算其回弹量。
图16 U形件冲压模具示意图
Fig.16 U-shaped stamping die pictures
设置板料和模具的初始温度为20 ℃，将凸模、凹模以及压边圈均设置为刚体，不产生变形热，但与板料之间存在摩擦热，而板料存在变形热以及摩擦热。材料塑性本构模型按照得到的Norton-Hoff模型输入，摩擦力模型按照式(16)输入。
另加，增加一组冷冲压模型下的冲压成形数值模拟计算，该模型不考虑变形热和摩擦热带来的影响，将其与两种模型下的结果进行对比研究。
3.3.2 新式冲压成形结果讨论
3.3.2.1 冲压过程温度变化
设置分析步为与温度相关的显式动态分析过程，即Dynamic-Temp-Disp-Explicit分析步，得到热力耦合模型下的U形件冲压过程的温度变化如图17所示，图中箭头方向为板料流动方向。
(a)t=0.25 s (b)t=0.5 s
(c)t=0.75 s (d) t=1.0 s
图17 完全热力耦合冲压过程温度变化
Fig.17 Temperature change of the stamping process
with fully thermal coupling
由图17可看出，在冲压过程中，侧壁靠近上圆角处(图17d中②处)的发热量最大，一方面该区域变形较大，产生变形热最多，另一方面，圆角区域板料与模具贴合力较大，接触面之间产生较大的摩擦力，产生较多摩擦热，最终导致该区域板料温升最快；其次，板料与压边面接触区域(图17d中③处)产生的温升也相对较大，该区域产生的摩擦热较多，对冲压过程造成的影响显著。
3.3.2.2 凸模力结果对比
在U形件的冲压过程中，随着凸模不断下降，板料的变形增大，产生的变形抗力不断增大，采集两种冲压模型的凸模力进行对比研究，结果如图18所示。
由图18可知，在两种冲压模型下，凸模力都随着时间延长而增大，但是由于存在较大的波动，不易观察规律，故建立两种模型下的凸模力的趋势线，以此进行对比研究。
图18 凸模力变化趋势对比
Fig.18 Comparison of punch force changes
根据趋势线的变化程度发现，在初始成形阶段，热效应给凸模力带来的影响微乎其微，此时尚未有较大温升；在冲压成形中后期，存在较多的塑性变形热和摩擦热积累，热效应带来的温升明显提高了材料的塑性，相同应变下，凸模产生的力较小。
3.3.2.3 回弹结果对比
在得到两种模型下的冲压成形残余应力分布结果后，建立新的静力隐式模型，将最后一个分析步导入作为初始状态(initial state)，约束顶部表面，计算回弹。
冲压成形应力最终结果如图19所示，由图19b可明显看出，在不考虑热效应的情况下，成形最终阶段的最大内应力相对较大，最大内应力差值为17.2 MPa。而图19a中，由于热效应温度升高促使材料塑性增强，在相同的变形条件下，板料内部产生的应力较小，进一步验证了温变塑性模型的准确性。
按照残余应力分布计算得到U形件回弹分布(图20)，选取剖面图中的16个选点(图20a)进行回弹对比，探究热效应对回弹的影响。
(a)考虑热效应
(b)不考虑热效应
图19 最终应力分布对比
Fig.19 Comparison of final stress distribution
(a)考虑热效应
(b)不考虑热效应
图20 回弹分布侧视图对比
Fig.20 Side-view comparison of springback distribution
图21 回弹计算结果对比
Fig.21 Comparison of springback calculation results
由图21可知，在两种计算模型下，回弹结果出现明显的差异，在靠近顶部圆角处(图20中1～5选点)，回弹量的差值较小；在侧壁靠近法兰处以及法兰边上，回弹量都随着选点位置增大逐渐增大，且两种模型下的差值逐渐增大。显然，考虑热效应的模型与传统冷冲压分析模型存在明显区别，且考虑热效应的新式冲压模型的回弹量普遍较小。
为验证新式冲压模型的准确性，采用试制的方法得到该U形件的实际冲压结果，如图22所示。按照图20中的选点对比较难操作与测量，容易产生极大误差，因此采用法兰边和侧壁的回弹角度进行对比，结果如表4所示。对比表4可知，考虑变形热显然更加符合实际的回弹情况。
图22 冲压U形件及最大回弹测量
Fig.22 Stamping U-shaped part and maximum
springback measurement
表4 回弹角度大小对比
Tab.4 Comparison of springback angle size
侧壁回弹角法兰回弹角冷冲压模型(°)5.9013.32热效应模型(°)5.3810.91试制U形件(°)5.1610.36冷冲压模型误差(%)14.3428.57热效应模型误差(%)4.265.31
4 结论
(1)考虑变形热条件下的拉伸结果明显与恒温下的拉伸结果不同，热效应导致流变应力曲线低于恒温下的应力曲线，变形热导致金属塑性增强，且不可忽略不计。
(2)金属钣金冲压过程中摩擦区域的摩擦热近似与摩擦面的相对位移成正比增长，一旦脱离摩擦面，温升不再产生，产生的摩擦热对金属塑性存在明显的影响。
(3)将变形热和摩擦热同时考虑进模型中，通过钣金冲压热力耦合分析发现，同时考虑变形热和摩擦热的情况下，U形件的厚度变化更小，最终步长下的残余应力相对较小，且冲压所需的凸模力也更小。
(4)对计算模型下的U形件进行试制，对比两种模型下的回弹结果与实际回弹量贴合性,发现传统的冷冲压模型和考虑自发热效应的新式模型下分析的回弹结果差异明显，考虑自发热效应的模型下，侧壁回弹的准确率提高10.08%，法兰处回弹准确率提高23.26%，明显更加贴合实际生产的测量结果。
[1] NADERI M. Constitutive Relationships for 22MnB5 Boron Steel Deformed Isothermally at High Temperatures[J]. Materials Science and Engineering：A， 2008， 478: 130-139.
[2] CAI Zhengyang, WAN Min, LIU Zhigang, et al. Thermal-mechanical Behaviors of Dual-phase Steel Sheet under Warm-forming Conditions[J]. International Journal of Mechanical Sciences, 2017, 126:79-94.
[3] SUNG J H , KIM I H, WAGONER R H. A Plastic Constitutive Equation Incorporating Strain, Strain-rate, and Temperature[J]. International Journal of Plasticity， 2010， 26(12):1746-1771.
[4] AUTILA M , LANTTO E , ARKKIO A . Determination of Forces and Linearized Parameters of Radial Active Magnetic Bearings by Finite Element Technique[J]. IEEE Transactions on Magnetics, 1998, 34(3):684-694.
[5] RUSINEK A， KLEPACZKO J R.Experiments on Heat Generated during Plastic Deformation and Stored Energy for TRIP Steels[J]. Materials and Design， 2009， 30:35-48.
[6] HODOWANY J, RAVICHANDRAN G, ROSAKIS A J， et al. Partition of Plastic Work into Heat and Stored Energy in Metals[J]. Experimental Mechanics， 2000， 40:113-123.
[7] KNYSH P, KORKOLIS Y P. Determination of the Fraction of Plastic Work Converted into Heat in Metals[J]. Mechanics of Materials 2015， 86:71-80.
[8] PREZ-CASTELLANOS J L. Temperature Increase Associated with Plastic Deformation under Dynamic Compression:Application to Aluminium Alloy Al 6082[J]. Journal of Theoretical and Applied Mechanics， 2012， 50: 377-398.
[9] JUNGHWAN C, BEOMJOON C, SOOJEONG H. Numerical Modeling of the Thermal Deformation during Stamping Process of an Automotive Body Part[J]. Applied Thermal Engineering， 2018，128:159-172.
[10] MAJIDI O, BARLAT F, LEE M G. Effect of Slide Motion on Springback in 2-D Draw Bending for AHSS[J]. International Journal of Material Forming, 2016, 9: 313-326.
[11] PEREIRA M P, ROLFE B F. Temperature Conditions during ‘Cold’ Sheet Metal Stamping[J]. Journal of Materials Processing Technology， 2014，214:1749-1758.
[12] LI Daoming, AMIT G. Biaxial Warm Forming Behavior of Aluminum Sheet Alloys[J]. Journal of Materials Processing Technology， 2004， 145：281-293.
[13] 赵玉璋，王凯，王武荣，等. 高强度双相DP780钢板冲压成形变摩擦系数模型及其使用[J]. 上海交通大学学报, 2015, 49(10):1446-1451.
ZHAO Yuzhang, WANG Kai, WANG Wurong, et al.Application of Variable Friction Coefficient Model in Forming of Advanced High-Strength Steel[J]. Journal of Shanghai Jiaotong University, 2015, 49(10):1446-1451.
[14] WANG Wurong, ZHAO Yuzhang. A Study on Variable Friction Model in Sheet Metal Forming with Advanced High Strength Steels[J]. Tribology International, 2016, 93:17-28.
[15] LI Daoming, AMIT G. Effects of Temperature and Blank Holding Force on Biaxial Forming Behavior of Aluminum Sheet Alloys[J]. JMEPEG， 2004， 13:348-360.
[16] 聂昕，谭广，乔晓勇. 基于热-力耦合和变摩擦因数的高强钢冷冲压成形性[J].中国机械工程, 2018, 29(16):1996-2002.
NIE Xin, TAN Guang, QIAO Xiaoyong. Formability of HSS in Cold Stamping Based on Variable Friction Coefficient and Thermal-mechanics Coupling[J]. China Mechanical Engineering, 2018, 29(16):1996-2002.
[17] LI Daoming, AMIT G. Tensile Deformation Behavior of Aluminum Alloys at Warm Forming Temperatures[J]. Materials Science and Engineering：A， 2003， 352：279-286.
[18] 刘大海，许钢碧，常春. DP780 高强钢板温热成形极限图及其计算模型[J]. 塑性工程学报, 2017, 24 (2):192-197.
LIU Dahai, XU Gangbi, CHANG Chun. Forming Limit Diagram and Its Calculation Model for DP780 High Strength Steel Sheet at Elevated Temperatures[J]. Journal of Plasticity Engineering, 2017, 24 (2):192-197.
