随着汽车产业的快速发展和能源的不断消耗，汽车轻量化的问题越来越重要^[1]，而实现汽车轻量化的有效方式之一是采用轻量化材料，如高强钢^[2-3]。在汽车产业中，相当一部分部件是通过冲压加工而成的。在金属板料冲压成形过程中，材料成形效果的评价非常重要，工业应用中通常采用成形极限曲线(forming limit curve，FLC)作为材料是否发生破裂的判断标准。成形极限曲线的概念最早由Keeler和Goodwin提出^[4-5]，因为应用效果好而沿用至今^[6-7]。成形极限曲线以第2主应变为横坐标轴，以第1主应变为纵坐标轴，曲线呈“V”字形。当产品的所有单元的应变对都位于曲线以下时认为产品安全，当产品中存在某处单元的应变对位于曲线以上时则认为产品会出现破裂。

自成形极限曲线被提出以来，其获取方法得到了很大的发展。最初，人们主要通过大量实验得到不同材料的成形极限曲线^[8-9]，但是实验过程繁琐、费时费力。经过大量实验经验的积累，研究者们总结出经验性的Keeler公式^[4]，即已知低碳钢的硬化系数n和厚度t，就可以得到材料的近似成形极限曲线，但此曲线与真实成形极限曲线还有一定的出入。常用的获取成形极限曲线的方法还有理论方法^[10-11]和有限元方法^[12-14]，而现在的大部分理论方法也需要借助有限元方法来实现，因此有限元仿真在成形极限曲线预测上的应用越来越广泛。

虽然对高强钢的成形极限曲线的研究比较多，但是其中大部分采用的是各向同性材料或者仅仅考虑厚向异性^{[10, 12, 15-17]}，少数涉及到面内各向异性的研究工作也很少采用有限元虚拟试验进行预测^{[8, 11, 18]}。本文考虑材料的各向异性，对高强钢Q490C的成形极限曲线进行了仿真预测，并和试验结果、Keeler公式计算结果进行了对比，为板料冲压成形等后续工作奠定基础。

1 仿真模型

根据国标GB/T 15825.8—2008《金属薄板成形性能与试验方法第8部分：成形极限图(FLD)测定指南》的要求，参照实际试验中模具和板料的几何尺寸，在软件UG中分别建立凸模、凹模、压边圈和板料的几何模型，然后将几何模型导入有限元分析软件。选取的试件如图 1所示，其宽度b取40、60、80、100、120、140、160、180 mm，以分别代表实际试验中各种规格的试件的成形情况。所有试件的厚度均为3 mm。仿真所用的模具如图 2所示，有限元模型如图 3—7所示。板料轧制方向均沿试件的长边方向。

仿真采用板料冲压仿真专用有限元软件Dynaform。凸模、凹模和压边圈都采用刚体单元，拉延筋采用等效拉延筋。板料采用完全积分壳单元，厚度方向采用7个积分点。网格大小为2 mm，具体的单元数目因试件大小而有所不同。接触处理采用罚函数法，摩擦方法采用Coulomb摩擦。冲压过程中，压边力为210 kN，凸模速度为2 m/s。摩擦系数为0.125。仿真过程分为两个工步：

1) 压边圈闭合。此时凹模与凸模都保持静止，板料在压边圈作用下向凹模运动，直到压边圈和凹模距离为3 mm时，此工步结束。

2) 拉延工步。此时压边圈闭合，并施加恒定压边力210 kN。与此同时，凸模向凹模运动，带动板料完成成型过程。给定一个较大的凸模行程，当达到这一行程时仿真结束。

2 材料模型和屈服准则

为了采用有限元程序进行有效仿真，需要测试所用材料的力学特性。本文根据国家标准，分别测试与轧制方向成0°、45°、90°的3个方向的应力应变曲线，然后将其转化为真实应力应变曲线。每个方向测试3个试件。在弹性阶段，求出各试件的弹性模量。因为本文采用的是Dynaform中的36号材料，假定材料在弹性阶段是各向同性的，但是在塑性阶段是各向异性的^[19-20]，所以对每个方向的弹性参数求平均值，得到整体平均弹性模量为217.2 GPa，平均Poisson比为0.26。板料的密度为7 850 kg/m³。

在强化阶段，采用幂函数σ=K(ε+b)ⁿ进行拟合(其中：K、b、n分别为拟合曲线的参数，σ为单向拉伸试验的真实应力，ε为单向拉伸试验的真实应变)，并针对每个方向求3个试件的平均值，得到3个方向的拟合参数。

为了考虑厚度方向的影响，分别测试了与轧制方向成0°、45°、90°的3个方向的各向异性系数r(具体为r₀、r₄₅、r₉₀)。r为单向拉伸试验时，试样宽度方向和厚度方向真实的塑性应变的比值。同样地，每个方向测试3个试件，并求平均值。

通过以上的两组试验，得到仿真所需要的参数如表 1所示。

通常认为轧制钢板的力学特性为正交各向异性，此种情况常用的屈服准则是Hill1948屈服准则和Barlat1989屈服准则，对于各向异性系数小于1的板料，后者更精确，因此本文采取了使用Barlat1989屈服准则的Barlat三参数塑性模型(即Dynaform中的36号材料)，其屈服准则为

$ \begin{array}{l} f = a{\left| {{k_1} + {k_2}} \right|^M} + a{\left| {{k_1}-{k_2}} \right|^M} + \\ \;\;\;\;\;\;\;\;\;\;\;\;c{\left| {2{k_2}} \right|^M} = 2\sigma _{\rm{e}}^M. \end{array} $

(1)

其中σ_e为屈服应力。M为整数，此处M=6。系数k₁和k₂可由式(2)计算，

$ \left\{ \begin{array}{l} {k_1} = \frac{{{\sigma _{11}} + h{\sigma _{22}}}}{2}, \\ {k_2} = {\left[{{{\left( {\frac{{{\sigma _{11}}-h{\sigma _{22}}}}{2}} \right)}^2} + {p^2}\sigma _{12}^2} \right]^{\frac{1}{2}}}. \end{array} \right. $

(2)

σ₁₁和σ₂₂分别是平面内两个方向的主应力, σ₁₂为剪应力。a、c、h为材料相关的参数，可以由不同方向的各向异性系数r₀、r₄₅、r₉₀得到，

$ \left\{ \begin{array}{l} a = 2-c = 2-2\sqrt {\frac{{{r_0}}}{{1 + {r_0}}}\frac{{{r_{90}}}}{{1 + {r_{90}}}}}, \\ h = \sqrt {\frac{{{r_0}}}{{1 + {r_0}}}\frac{{1 + {r_{90}}}}{{{r_{90}}}}} . \end{array} \right. $

(3)

将表 1中的3个方向的各向异性系数r(具体为r₀、r₄₅、r₉₀)代入式(3)可以得到a=1.194 9, c=0.805 1, h=0.969 8。

任意方向的各向异性系数r_ϕ可以由式(4)得到，

$ {r_\phi } = \frac{{2M\sigma _{\rm{e}}^M}}{{\left( {\frac{{\partial f}}{{\partial {\sigma _{11}}}} + \frac{{\partial f}}{{\partial {\sigma _{22}}}}} \right){\sigma _\phi }}}-1. $

(4)

其中σ_ϕ是对应方向的单向拉伸真实应力。参数p通过拟合45°方向的单轴拉伸测试数据，按式(4)由软件自行进行迭代计算得到^[21]。

3 仿真失效准则和仿真数据处理

试件采用壳单元进行模拟计算，需要引入专门的失效准则^[22]。常用来判断薄板冲压集中失稳的失效准则有最大凸模力准则、应变梯度变化准则和应变速率变化准则^[23]。其中，最大凸模力准则通过有限元仿真得到凸模与板料的冲击力相对于时间的曲线(如图 8所示)，当材料失稳时凸模力会突然下降，因而此时板料上的最大应变就可以作为极限应变。本文试验中用到的失稳判断准则也是此准则，因此仿真分析也以此准则作为失效准则以获取不同规格试件的极限应变。板料失稳时截面线上的主应变分布如图 9所示。

通常，对于成形极限曲线，横坐标轴负的一侧的试验数据呈明显的线性特征，横坐标轴正的一侧呈明显的二次曲线特征，因此仿真的成形极限曲线的左侧采用线性拟合，右侧采用二次函数拟合，两个拟合公式分别如式(5)所示^[24]。采用最小二乘法使得总误差取最小值，具体通过数值计算软件MATLAB实现。

$ \left\{ \begin{array}{l} {\varepsilon _1} = {\rm{FL}}{{\rm{C}}_0} + {\varphi _1}{\varepsilon _2}, \;\;\;\;\;{\varepsilon _2} \le 0;\\ {\varepsilon _1} = {\rm{FL}}{{\rm{C}}_0} + {\varphi _2}\varepsilon _2^2, \;\;\;\;\;{\varepsilon _2} > 0. \end{array} \right. $

(5)

其中：ε₁为第1主应变，ε₂为第2主应变，FLC₀为成形极限曲线与纵坐标的截距，φ₁和φ₂为拟合曲线的系数。

4 仿真结果讨论和结果验证

为了验证仿真结果的正确性，采用厚度为3 mm的Q490C板材，按照国标GB/T 15825.8—2008《金属薄板成形性能与试验方法第8部分：成形极限图(FLD)测定指南》开展试验。试件的形状如图 1所示。试验机和试验后的试样如图 10所示。采用GMASystem网格应变自动测试系统测试得到颈缩区网格的第1主应变和第2主应变，绘制在以第2主应变为横坐标轴、以第1主应变为纵坐标轴的坐标系中。因为实际测量的是颈缩区或者临近颈缩区的网格的应变，所以会有某些数据点的安全余量很大，应去除安全余量过大的数据点，用剩下的数据点根据式(5)拟合得到试验的成形极限曲线，结果如图 11所示。

针对不同规格的试件，分别进行了有限元仿真，得到的典型的Mises等效应力云图、等效塑性应变云图、减薄率云图分别如图 12—14所示。在Mises等效应力云图 12中，在试件的球形凸起上，等效应力最大且呈环形分布，为971.731 MPa；同时，球形凸起与压边圈压实的板料之间的过渡圆角处，应力也较大。在等效塑性应变云图 13中，同样在试件的球形凸起上，等效塑性应变最大但呈两个“瓣”状分布，为0.538。因为两个“瓣”状处的塑性应变较大，所以试件多在“瓣”状对应位置处破裂而非球形凸起的中线。对于减薄率云图 14，在试件的球形凸起上，减薄率最大的位置也呈两个“瓣”状分布，为37.065%，这种分布也与等效塑性应变分布一道解释了试件的破裂位置。

采用第3节提到的仿真方法分别得到极限应变，并进行分段拟合，得到的结果如图 15所示。可以看出，在单向拉伸区域(图 15中曲线的左侧)和平面应变区域(图 15中横坐标为0附近的曲线)，采用最大凸模力准则预测得到的极限应变与仿真拟合的成形极限曲线吻合得较好；在双向拉伸区域(图 15中曲线的右侧)采用最大凸模力准则预测得到的极限应变与仿真拟合的成形极限曲线的偏差稍大。

每一次仿真相当于一次虚拟试验，对8个不同规格的试件进行仿真得到8组极限应变，然后采用分段拟合的方法可以得到仿真的成形极限曲线，相当于对不同的试验结果进行统计意义上的平均，可以得到更为真实可靠的结果。如图 15所示，仿真结果拟合得到的成形极限曲线与试验得到的成形极限曲线在左侧和中间位置很接近，其中试验拟合得到的成形极限曲线在纵坐标轴的截距为0.434 7，而仿真结果拟合得到的成形极限曲线在纵坐标轴的截距为0.452 3；仿真结果拟合得到的成形极限曲线在右侧略高于试验得到的成形极限曲线，在第2主应变为0.25时，仿真得到的成形极限曲线比试验得到的成形极限曲线高0.043 4。

然后，对仿真结果、试验结果和经验公式进行比较。在冲压仿真软件Dynaform中，经验性的成形极限曲线与第1主应变轴的截距(%)由Keeler公式^[4]计算得到：

$ \left\{ \begin{array}{l} {\rm{FL}}{{\rm{C}}_0}/\% = \frac{{\left( {23.3 + 14.13t} \right)n}}{{0.21}}, \;\;\;0 < t < 2.54;\\ {\rm{FL}}{{\rm{C}}_0}/\% = \frac{{\left( {20 + 20.67t-1.94{t^2}} \right)n}}{{0.21}}, 2.54 \le t < 5.33;\\ {\rm{FL}}{{\rm{C}}_0}/\% = 75.125n/0.21, t \ge 5.33. \end{array} \right. $

(6)

其中：t为板料厚度；n为材料的硬化系数，当n>0.22时，取n=0.22。本材料3个方向的硬化系数(n₀、n₄₅、n₉₀)由3个方向的拉伸曲线得到(如表 1所示)，并由式(7)求得平均值。

$ \bar n = \frac{{{n_0} + 2{n_{45}} + {n_{90}}}}{4} = 0.1518. $

(7)

可以得到当板料厚度为3 mm时的经验性的成形极限曲线与第1主应变轴的截距为

$ {\rm{FL}}{{\rm{C}}_0} = 46.66\% . $

(8)

其截距换算成小数为0.466 6。对于板料厚度为3 mm的试件，根据试验数据拟合得到的截距为0.434 7，经验公式比试验拟合曲线的截距高7.34%。根据仿真结果得到的截距是0.452 3，经验公式比其高3.16%。可见，仿真结果介于经验公式与试验结果之间。

5 与各向同性仿真的对比

本文在充分考虑材料各向异性的基础上对成形极限曲线进行了仿真预测。为了对比其与各向同性仿真的差别，分别采用了表 1中0°、45°和90°的材料参数作为各向同性的参数，采取本文所述的方法进行成形极限曲线预测，得到的结果如图 16所示。

可以看到，考虑了各向异性的仿真结果与试验结果吻合得最好。以0°和45°的参数作为各向同性参数时得到的成形极限曲线比试验值偏高很多；以90°的参数作为各向同性参数时得到的成形极限曲线在大部分区域比试验值偏低较多；而各向异性本构综合了3个方向的材料特性，更能反映材料真实的特性，可以明显提高预测精度。相比其他各向同性仿真得到的曲线，考虑各向异性的仿真结果总体误差最小。

举例来讲，在第2主应变为0.25时，试验得到的成形极限曲线的第1主应变为0.492 5。采用0°的材料参数作为各向同性参数进行仿真预测得到成形极限曲线的第1主应变为0.610 9，比试验值高24.04%；采用45°的材料参数作为各向同性参数进行仿真预测得到成形极限曲线的第1主应变为0.659 3，比试验值高33.87%；采用90°的材料参数作为各向同性参数进行仿真预测得到成形极限曲线的第1主应变为0.448 4，比试验值低8.94%；而采用各向异性的模型进行仿真预测得到成形极限曲线的第1主应变为0.535 9，比试验值高8.81%。

6 结论

本文针对各向异性高强钢Q490C，采用最大凸模力准则对板料的极限应变进行了有限元仿真，得到以下结果：

1) 得到了高强钢Q490C有限元仿真预测的成形极限曲线；

2) 进行试验验证，有限元仿真预测得到的结果与试验数据得到的曲线吻合较好；

3) 将有限元仿真预测得到的曲线的截距与经验性的Keeler公式给出的截距进行比较，二者一致性较好；

4) 将各向同性的成形极限曲线的仿真结果与考虑各向异性的仿真结果进行了比较，发现考虑各向异性的成形极限曲线仿真可以明显提高预测精度。

本文对各向异性材料Q490C的成形极限曲线的研究，有助于进一步开展此高强钢的冲压成形开裂预测工作。