板料快速成形预测方法与流程

本发明涉及一种板料快速成形预测方法，属于板料成形技术领域。

背景技术：

板料成形是用板料等薄壁型材等作为原材料进行冲压加工的方法。板料成形是汽车覆盖件制造的一个重要部分，它是一个涉及几何形状、材料物理以及接触摩擦三重非线性的复杂问题。早期依靠工程师直觉的经验法和试验法，虽然能在一定程度上满足制造的需求，但是具有周期长、费用高等缺点，因此人们试图通过数学和力学的工具来对这个问题进行研究。在众多的数学方法中，由于具有坚实的理论基础和良好的实践检验，有限元法成为模拟板料成形的一种有效数值方法。

根据计算原理的不同，有限元法可以分为动力显式有限元法和静力隐式有限元法。动力显式有限元法，虽然能够较为精确地模拟板料厚度等物理量分布，但是考虑到数值计算的稳定性因素，必须将计算过程中的时间步长设置地足够小，导致它需要大量的运行时间，满足不了模具设计阶段的快速响应需求；静力隐式有限元法尽管可以将时间步长设置地足够大，但是存在平衡方程数值求解的收敛性困难的问题，造成该算法在实际应用中受到很大限制。因此，研究实现一种具有良好收敛性的静力隐式有限元法，来快速准确地预测板料成形具有十分重要的理论意义和应用价值。

技术实现要素：

本发明要解决技术问题是：克服上述技术的缺点，提供一种能够满足设计阶段的及时响应需求、可以在较短的时间内实现板料成形的快速准确模拟的方法。

为了解决上述技术问题，本发明提出的技术方案是：一种板料快速成形预测方法，包括以下步骤：

第一步、对凸模、凹模、压边圈和板料建模后进行网格划分；

第二步、将板料的四周边界通过压边圈固定，然后将凸模从板料处下压并运动一预设行程步使板料作弯曲变形；通过接触搜索算法寻找板料与凸模发生穿透的网格节点以及板料与凹模接触的网格节点作为约束节点，板料的其余网格节点作为自由节点；根据约束节点采用网格变形方法来求解当前行程步的板料预示形状作为滑移约束面；

第三步、将板料节点约束在所述滑移约束面上，使板料可在滑移约束面上流动，即将板料整体作拉伸变形；

根据虚功原理，内外力所做的虚功相等建立平衡方程，对该平衡方程使用迭代法进行求解，迭代收敛后得到板料当前行程步的所有节点的位移增量，从而得到当前行程步的拉伸变形形状；

第四步、根据当前行程步的拉伸变形形状计算当前行程步的板料每个网格单元的应力应变，并根据体积不变原理，即板料的每个网格单元变形前后的体积不变，计算得到变形后每个网格单元的厚度；

第五步、重复执行第二步到第四步，直到凸模和凹模两者合模时为止，最终得到板料成形的过程及变形情况。

本发明带来的有益效果是：本发明将板料的成形过程分解为弯曲变形和拉伸变形两个独立的变形过程，弯曲变形时板料的边界(即压边部分)固定不动(即边界形状不发生变化)，与凸模接触部分随着凸模前进某一行程步而发生弯曲变形形成滑移约束面，拉伸变形时在滑移约束面上进行滑动，拉伸变形前后板料的边界形状由于滑动发生了改变。

本发明将板料的成形过程分解为弯曲变形和拉伸变形两个独立的变形过程后，简化了计算过程，大大提高预测效率。传统的计算方法通常每个节点需要计算六个自由度，而本发明在弯曲变形求解时每个节点只需要求解三个自由度，拉伸变形求解时只需要求解两个自由度，从而使得本发明节约了大量的计算时间，提高了预测速度，能够在较短的时间内进行板料成形的模拟分析，可以满足设计阶段的及时响应需求。

另外，本发明可采用大步长隐式求解方法，即前面所述的行程步可采用较大的步长，而传统的显式求解步长非常小，效率较低，本发明的大步长也节约了较多的计算时间。

上述技术方案的进一步改进是：第三步中，所述平衡方程为：

其中，[k]是单元刚度矩阵，{δu}是位移增量向量，是第e个网格单元的节点内力向量，是第e个网格单元的等效节点外力向量，1≤e≤n，n为板料的网格单元数量；对所述平衡方程使用newton-rapshon迭代法进行求解。

上述技术方案的再进一步改进是：第四步中，计算当前行程步的板料每个网格单元的应力应变时，在比例加载的条件下得到弹塑性本构方程：

其中，

式中，r为板料平均面内的厚向各项异性系数，m为板料的强化系数，ε0为网格单元的初始应变，n为板料的硬化指数，为网格单元的应力，为网格单元的应变。板料平均面内的厚向各项异性系数r为板料试件单向拉伸试验中宽度应变与厚度应变之比，是板料固有的参数之一。

附图说明

下面结合附图对本发明作进一步说明。

图1是本发明实施例的流程示意图。

图2是本发明实施例的初始状态示意图。

图3是板料弯曲变形的示意图。

图4是板料拉伸变形的示意图。

图5是板料的成形过程分解原理示意图。

具体实施方式

实施例

本实施例以一方盒零件为例进行说明，凸模总行程为51.89mm，所采用的板料crdq钢板的参数如下表所示：

表1-板料材料参数

本实施例的板料快速成形预测方法，如图1所示，包括以下步骤：

第一步、如图2所示，对凸模1、凹模2、压边圈和板料3建模后进行网格划分。

第二步、如图3所示，将板料3的四周边界通过压边圈固定，然后将凸模1从板料3处下压并运动一预设行程步使板料3作弯曲变形；通过接触搜索算法寻找板料与凸模发生穿透的网格节点以及板料与凹模接触的网格节点作为约束节点，板料的其余网格节点作为自由节点；根据约束节点采用网格变形方法来求解当前行程步的板料预示形状作为滑移约束面。

接触搜索算法为现有技术，可参考相应技术文献，其原理是：运用划分空间网格的方法来搜索与板料节点距离最近的模具(根据情况选择凸模或凹模)主表面节点，然后通过模具主表面上的节点与网格单元的拓扑关系来确定投影单元的范围，进而确定投影点。在确定了投影点后，利用几何计算来判断节点是否穿透凸模主表面，以及是否与凹模接触。

网格变形方法也为现有技术，可参考海量的现有文献。比如采用以下方法：通过建立一个基于二阶离散lbo算子的方程组：

这里ik是k×k维单位矩阵，l²是(n-k)×(n-k)维矩阵，d＝(d1,…,dn)^t是待求解的n×3维的未知矩阵，c1,…,ck是指定的约束节点位置，l＝m^-1ls，而m是一个对角矩阵，且ls是一个对称矩阵，

其中，是节点vi的邻接节点集合，而

这里，αij和βij是边eij所在三角形中的两个对角，ai为围绕顶点vi的voronoi三角形的面积。

可以证明，上述的方程组是对称线性方程组，这里使用lu分解法来进行数值求解，依次求解接触后的线性方程组，即可得到滑移约束面的形状。

第三步、将板料节点约束在所述滑移约束面上，使板料可在滑移约束面上流动，即将板料整体作拉伸变形；

根据虚功原理，内外力所做的虚功相等建立平衡方程，对该平衡方程使用迭代法进行求解，迭代收敛后得到板料当前行程步的所有节点的位移增量，从而得到当前行程步的拉伸变形形状。

所述平衡方程为：

其中，[k]是单元刚度矩阵，{δu}是位移增量向量，是第e个网格单元的节点内力向量，是第e个网格单元的等效节点外力向量，1≤e≤n，n为板料的网格单元数量；对所述平衡方程使用newton-rapshon迭代法进行求解。

根据材料流动理论法则，求解板料拉伸变形的形状结果，如图4所示，由于板料产生流动运动，因此它的外轮廓线(边界)形状发生了改变。

第四步、根据当前行程步的拉伸变形形状计算当前行程步的板料每个网格单元的应力应变，并根据体积不变原理，即板料的每个网格单元变形前后的体积不变，计算得到变形后每个网格单元的厚度。

计算当前行程步的板料每个网格单元的应力应变时，在比例加载的条件下得到弹塑性本构方程：

其中，

式中，r为板料平均面内的厚向各项异性系数，m为板料的强化系数，ε0为网格单元的初始应变，n为板料的硬化指数，为网格单元的应力，为网格单元的应变

第五步、重复执行第二步到第四步，直到凸模和凹模两者合模时为止，最终得到板料成形的过程及变形情况。。

本实施例将板料的成形过程分解为弯曲变形和拉伸变形两个独立的变形过程，其原理如下，如图5所示，将坐标系固定在节点切平面所在的局部坐标系中，图5中oxyz为整体坐标系，oxyz网格为单元局部坐标系，为网格节点随体坐标系，则每个行程步中需要求解下述方程组：

其中，δu1和δu2是待求解的节点位移增量向量，u0为指定节点在节点随体坐标系中的位移向量，f2是由于外力在切平面上引起的摩擦力，r1和r2分别是指定节点位移在节点所在的法线方向和切平面上产生的力，kmij和kgij(i,j＝1,2)是相应的刚度矩阵。

假设网格节点和网格单元的的法线方向之间的夹角可以忽略，这时，只需独立地求解下列两组有限元列式：

km11(u0)·δu1＝-r1(u0)(弯曲变形)

km22(u0)·δu2＝f2-r2(u0)(拉伸变形)

因此，板料成形过程就被分解为弯曲变形和拉伸变形两个独立的过程，从而简化了计算，提高了板料成形的预测效率。

采用本实施例方法和传统的ls-dyna软件的模拟结果相比较，两者的厚度分布趋势和主应变变化趋势相一致，说明了本方法的正确性。但是，两者的运行时间差距很大，在相同的计算机硬件配置下进行验证，本实施例方法的运行时间为30.56s，ls-dyna软件的运行时间为180.42s，可见本实施例方法大幅提高了计算效率，对复杂形状的板料成形预测有着积极的意义。

本发明不局限于上述实施例所述的具体技术方案，除上述实施例外，本发明还可以有其他实施方式。对于本领域的技术人员来说，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等形成的技术方案，均应包含在本发明的保护范围之内。