镁合金筒形件拉深成形的临界压边力

杨柳，官英平，姚丹，段永川，高安娜

(1.燕山大学先进锻压成形技术与科学教育部重点实验室，河北秦皇岛，066004；2.燕山大学机械工程学院，河北秦皇岛，066004)

筒形件拉深成形是最基本的轴对称类板料成形工艺，而在成形过程中出现的破裂和起皱失稳现象也是最常见并且是最主要的缺陷。在拉深成形过程中，合理的设置压边力可有效避免破裂和起皱失稳缺陷的出现，因此，压边力的研究受到广泛关注。MOSTAFAPUR 等[1]研究了一种可用于筒形件拉深的新型脉动冲裁压边系统，通过采用适当的频率和间隙，可以提高拉深深度，改善厚度分布。QIN等[2]利用一定的假设条件，通过数值计算方法预测轴对拉深成形临界压边力，研究了不同皱纹形状对临界压边力的影响，提出可通过径向分块压边方法控制皱纹形状，进而提高成形性能。张志远等[3]根据Mise-Hill屈服函数和Tresca准则对筒形件拉深成形过程中各变形区域应力分布进行了推导，获得了临界压边力解析式，并进一步对影响临界压边力的因素进行了定性分析和研究，通过试验验证了理论计算的准确性。赵升吨等[4]运用高阶厚向异性屈服准则及能量法给出了筒形件拉深成形的有效压边力区间，在有效压边力区间内进行拉深试验，可获得合格的成形件。镁合金具有密度低、比强度高、散热和消震性好等特点，广泛应用于航空、航天、汽车和电子产品等行业部门[5-6]。与传统金属材质(钢铁，铝合金)的筒形件相比，镁合金筒形件所具备的性能优势明显[7-10]。在镁合金拉深成形研究方面：XIE等[11]通过筒形件拉深试验，研究了电脉冲对AZ31B 镁合金薄板延性的影响，发现随着温度的升高，拉深成形极限逐渐增大，当温度大于373 K时，拉深成形极限增加得更快。CAO等[12]将固体颗粒介质成形技术与超声振动塑性成形技术相结合应用到AZ31B 镁合金薄板成形中，提出了超声振动颗粒介质成形技术，发现超声振动能提高颗粒介质的压力传递性能和板料的成形性。官英平等[13]针对AZ31B镁合金方盒形件的成形过程进行了工艺实验，确定了最佳工艺参数组合。杨柳等[14]依据差温成形特点，通过分析塑性变形后的AZ31B 镁合金方盒形件组织性能，得出差温所导致的孪晶可有效提高成形深度。唐伟琴等[15]针对AZ31 镁合金板材深冲过程中的制耳现象，提出一种基于晶体塑性理论的镁合金板材冲压制耳行为的预测方法。李彩霞等[16]利用有限元技术分析了工艺参数对AZ31镁合金盒形件拉深性能的影响，发现温度和速度对盒形件的拉深成形均有很大的影响。通过以上研究成果可知，镁合金筒形件拉深成形研究主要集中在数值模拟、成形试验的层面，对于镁合金筒形件拉深成形过程的法兰变形区应力分布、变形特点的研究较少，难以给出理论层面的解析，没有涉及成形过程的临界压边力控制。本文作者选用AZ31B 变形镁合金薄板作为研究对象，运用Fields-Backofen 流变应力本构模型，对AZ31B 镁合金筒形件拉深成形过程中的法兰区变形规律进行研究，给出应力的理论解析表达式，依据法兰变形区的理论解析公式，推导出拉深失稳起皱和破裂时临界压边力的解析模型，并利用数值模拟和试验相结合的方法对理论解析进行验证。

1 理论解析

1.1 本构关系模型

材料的本构关系模型是通过数学的形式来描述材料力学性质的方程，建立正确的本构关系模型对于分析材料的塑性变形过程，提高成形性能有很大的作用。AZ31B 镁合金流变应力属于动态再结晶类型，会出现加工软化现象，针对这一特征，镁合金本构方程的建立需要考虑应变速率和温度的影响。可采用Fields-Backofen本构关系模型[17]：

式中：σ为应力；ε为应变；˙为应变速率；K为材料的强度系数；n为应变硬化指数；m为应变速率敏感系数。K，n和m均与温度和应变速率相关。

以上述镁合金本构方程为基准，利用AZ31B 镁合金板材在不同温度、不同应变速率下所进行的单向拉伸试验[18]，确定本构方程中K，n和m参数与温度的关系式：

1.2 法兰区应力解析

为简化AZ31B镁合金筒形件法兰区的应力分析，本文进行如下假设：板料为厚向异性材料；法兰变形区厚向应力σz=0；等效应力与等效应变关系符合

根据以上假设，在拉深的某一瞬间，取筒形件法兰区一部分基元体作受力分析，如图1所示。

图1 筒形件法兰微元体受力分析Fig.1 Mechanics analysis on flange microelement of deep drawing cylinder parts

法兰区任一点处的瞬时应力平衡微分方程和塑性方程分别为

式中：σρ为变形质点在半径为ρ处的径向应力；σθ为变形质点在半径为ρ处的周向应力；ρ为变形质点的半径。

法兰外边缘可近似为单向压缩应力状态，其等效应变等于周向应变的绝对值，有

式中：R0为毛坯的初始半径；Rw为变形质点法兰外缘瞬时半径。

联立式(4)和(5)，得

联立式(3)和(6)可得

利用边界条件，当ρ=Rw时，σρ=0，则有法兰区任意半径ρ处变形质点的径向应力为

联立式(6)和式(8)，有

式中：C为积分常数，为作用在法兰外缘的压边力Q所引起的附加径向拉应力，μ为法兰区板料上下表面的总平均摩擦因数；t为板料的厚度。

此时，法兰区域的应力分布情况为：

1.3 破裂临界压边力模型

AZ31B 镁合金筒形件拉深时，在筒形件侧壁与凸模圆角相切处危险断面上的瞬时拉应力σmax为

式中：σwd为毛坯绕过凹模圆角区时所受的弯曲及反弯曲应力，σs为屈服强度；rd为凹模圆角半径；σwp为拉深初期毛坯在凸模圆角处的弯曲应力，rp为凸模圆角半径；eμα为毛坯沿凹模圆角区滑动所受摩擦阻力对筒壁拉应力的影响可用系数；包角

对于AZ31B 镁合金筒形件拉深，其极限承载能力σlim[19]为

式中：σb为抗拉强度；r为各向异性系数。

为使拉深试件不因破裂失稳而失效，应满足危险断面处的最大拉应力小于承载极限，即σmax≤σlim。依据式(12)和(13)可得出不发生破裂的临界压边力Qrup满足下式：

1.4 起皱临界压边力模型

依据法兰区起皱能量变化以及不发生起皱临界条件，有

式中：Uw为波纹隆起需要的弯曲功；Uθ为应力释放能量；UQ为压边力消耗功。

推得起皱数学模型为[20]

式中：y为质点挠度；y0为单波最大挠度；φ为单波任一弧段所对圆心角；φ0为单波所对圆心角；Rd为法兰的内缘半径(凹模内半径)。

假设任意ρ处，因周向应力σθ作用，周长缩短ΔS，释放的能量为Uθ，有

将式(11)和(18)代入式(17)中，可得单波释放的能量：

起皱失稳时，波纹隆起需要的弯曲功Uw可依据梁弹性弯曲的能量方法进行计算，用塑性切线模数D代替弹性模量E，有

式中：dI为半径ρ处厚度t、宽dρ的惯性矩，dI=dx为单波微分段在x轴的投影，dx=ρdφ；l为单波总长度在x轴投影，l=ρdφ0。

因为x=ρφ，dx=ρdφ，由式(16)可知：

定义塑性切线模数D为真实应力应变曲线的切线斜率：

令D=f2(ρ)，将式(23)和(24)代入式(20)，得

设总压边力为Q，共有N个波纹，则有N=2π/φ0。因在法兰外边缘处波纹挠度最大，压边力作用在法兰外缘部分，故每个单波上压边力做功为

将式(19)，(25)和(26)代入式(15)，可以求得压边力Q为

将式(27)对φ0进行微分，并令求得最小压边力下的φ0为

将式(28)代入式(27)，可得到起皱临界压边力为

2 数值模拟及试验

2.1 有限元模型建立

为验证上述理论分析推导的临界压边力公式可行性，利用有限元分析软件eta/DYNAFORM 对AZ31B镁合金筒形件拉深成形过程进行数值模拟，建立有限元模型如图2所示。采用Belytschko-Tsay(BT)单元类型和3参数Barlat材料模型，运用动力显式计算算法，凸模、凹模和压边圈均为不变形刚体。

图2 有限元模型Fig.2 Finite element model

2.2 模拟参数设定

工艺及模具参数为：虚拟拉深速度ν=2 000 mm/s；最大拉深高度h=30 mm；凸模半径Rp=25 mm，凹模内半径Rd=26.10 mm；凸、凹模圆角半径rp=rd=5 mm；毛坯半径R0=50 mm；毛坯厚度t=1.2 mm；摩擦因数μ=0.125；温度为250 ℃。

2.3 拉深成形试验

为验证理论解析和数值模拟结果的准确性，对AZ31B 镁合金板料进行了筒形件拉深试验。依据试验需要，以可实现压边力控制的BCS30D通用板材成形性试验机为平台，设计了模具及温度控制系统。试验拉深速度ν=30 mm/min，润滑剂采用耐高温镍基润滑油脂(0.10 ≤μ≤0.15)。

3 结果与讨论

以筒形件的侧壁最大减薄率和法兰外缘最大增厚率作为破裂和起皱的衡量指标，当筒形件的侧壁厚度低于板料初始厚度的30%认为破裂，当筒形件的法兰外缘厚度超过板料初始厚度的10%认为起皱。以恒定压边力加载，通过数值模拟可知临界起皱压边力为3 kN，破裂临界压边力为12 kN，厚度变化率分布分别如图3和图4所示。

依据理论解析与有限元数值模拟分别得出AZ31B 镁合金筒形件在拉深过程中破裂和起皱临界压边力的变化情况，如图5所示。

图3 压边力3 kN下的厚度变化率分布Fig.3 Distribution of thickness variation under blank holder force of 3 kN

图4 压边力12 kN下的厚度变化率分布Fig.4 Distribution of thickness variation under blank holder force of 12 kN

从图5可以看出：有限元数值模拟与理论解析所得临界压边力曲线相比，理论解析临界破裂最低点值高于有限元数值模拟值，相对误差为9.62%，其主要原因可能是拉应力分布不均导致危险区域应力过于集中而发生拉裂。图6所示为筒形件侧壁拉应力分布不均的现象，可见理论解析临界起皱最高点值低于有限元数值模拟值，相对误差为6.76%，主要原因可能是拉深过程中实际法兰上的等效应变不与瞬时径向坐标成简单的反比例关系。图7所示为筒形件等效应变分布不均的现象。

从图5还可以看出：在给定安全区范围内的恒定压边力作用下，由于有限元数值模拟过程中设置的参数与试验中的拉深工艺参数存在一定的差异(应力状态、变形状态、摩擦等)，实际试验中上限压边力附近有拉深破裂趋势，拉深成功试件所需的上限压边力整体降低。

图5 理论解析与数值模拟临界压边力变化情况Fig.5 Theoretical analysis and numerical simulation of critical blank holder force change

图6 筒形件侧壁拉应力分布Fig.6 Tensile stress distribution on side wall of cylindrical parts

图7 筒形件等效应变分布Fig.7 Equivalent strain distribution of cylindrical parts

图8所示为部分成形试验及有限元模拟结果的对比。从图8可以看出：在压边力F为1 kN时，由于压边力过小导致法兰部分板料周向压应力增大，板料径向流动受阻，出现起皱现象，拉深高度h小于12 mm；当压边力F为13 kN 时，因压边力过大使得法兰部分板料径向应力增大，向凹模流动受阻，凸模圆角位置发生破裂，拉深高度h不足14 mm；当压边力F为6 kN时，试件成形效果较好，拉深高度h接近30 mm。

从总体上看，理论分析计算结果与模拟和试验之间的差别不是很大，说明塑性理论分析计算具备一定的可行性和可靠性。

图8 试验及有限元模拟成形件Fig.8 Experiment and finite element simulation of forming parts

4 结论

1)基于平面应力假设条件，依据Fields-Backofen流变应力本构模型，对AZ31B 镁合金筒形件拉深成形过程中的法兰区应力分布进行解析研究。

2)运用应力解析公式及能量法，推导出AZ31B镁合金筒形件拉深成形的破裂和起皱临界压边力表达式。

3) 平面应力假设条件下所推得的理论解析结果与有限元数值模拟和试验验证的结果具有良好一致性，证明了理论解析计算的正确性。