半球形零件拉伸成型的数值分析.doc

北京航空航天大学毕业设计(论文) 第 78 页 1 绪 论1.1 背景及目的拉深是一种用平板毛料制成容器状零件的冲压成形方法，在飞行器、汽车、家用电器制造业得到了广泛的应用，并且可以作为评价板料成形性能的一种重要试验方法。早在三十年代人们就开始了对筒形件拉深成形的理论和实验研究，但系统和全面的研究，是在筒形件拉深被作为一种重要试验方法提出来之后才大量开展的。Swift在1939年提出以拉深试验作为评价材料成形性能好坏的基本试验方法之一。由于此试验方法的独特优点，已经被冲压界定为评定材料成形性能的必不可少的试验方法。筒形件拉深的研究是从理论和实验两方面进行的。实验研究的目的是为探讨拉深试验和其他较为简单的工艺试验的相关性，能否用简单工艺试验的试验参数来评定材料的拉深成形性能，以及试验上确定板料拉深性能的方法。Swift最早对拉深成形的力、功、应变以及能够成功地拉深成形的条件进行了系统的试验研究，为筒形件拉深成形的理论分析提供了可靠的基础；同时还研究了凹模圆角上的弯曲所引起的拉伸和变载问题，给出了摩擦系数的确定方法。肖夫曼在实验上得出摩擦和几何条件对材料的厚度变化分布有显著影响。在探讨单向拉伸试验测得的基本应力应变性质和拉深性能的关系方面，M.H.Sommer做了最早的努力，其后C.Arbel继续了这方面的工作，发现了应变强化模数和拉深性能的相关性。R.L.Whiteley对22种钢板、A.N.Bramley和P.B.Mellar对钛、锌板进行了试验，以研究拌料的各向异性对拉深性能的影响，认为代表各向异性的单向拉伸试验之宽向和厚向应变之比R越大，其拉深性能越好。在理论研究方面，人们一直在努力对此成形过程进行理论分析并计算成形力。Pomp和Sachs对此进行了最早的工作，为以后的理论研究奠定了基础。Hill接着分析了径向拉延问题1，他考虑的是平面应力、应变的极限情况。最早对拉深成形问题进行全面理论分析的是Chung和Swift，他们求解了变形区域的应力应变分布以及最大凸模载荷、总凸模功率和载荷行程曲线，讨论了凹模圆角上的弯曲和摩擦的影响，改进了弯曲计算模型；以后，Siebel和Panknin又继续寻找简单的方法以计算成形力。在利用单向拉伸试验的结果来推算各种材料的拉深性能方面，Norris、G.sachs、C.Arbel等人进行了深入的研究。G.G.Moore和J.F.Wallace，K.M.Frommann和W.F.Barday等人从理论上分析讨论了材料的各向异性指数R以及应变强化指数n对拉深成形的影响，其结论是R值对拉深成形性能有显著的影响，而n的影响较小；R值大，材料的性能好。六十年代电子计算机的普遍使用，为数值计算分析带来了极大的方便，也为板料成形的模拟和分析展示了广阔的前景。1963年D.M.Woo2提出了板料轴对称成形过程的一般分析方法，随后对液压胀形、球底和平底的杯形件拉深成形进行了完整细致的分析。E.I.Odell和W.E.Clausen运用增量理论与弯曲理论相结合的方法，考虑了刚塑性材料的拉深成形问题。D.C.Chiang和S.Kobayashi运用从边界条件出发的迭代法，着重分析了各向异性和应变强化特性对杯形件拉深的应力应变分布的影响，得到的结论是应变强化特性和各向异性对环向应变影响小，对径向应变和厚向应变影响大；径向应力和成形力随厚向异性指数R的减小而增加，而随n值的增加其变化程度减小。山田又提出了级数求解方法，并对比了数值解和级数解的结果，认为级数解在估算摩擦系数方面很好，且对某些问题也是足够精确的。尽管人们对拉深问题进行了大量的工作，理论和实验上都比较完备，但时之今日，仍然有很多问题需要解决，人们也没有停止过在这方面的研究。半球形零件在板料拉深成形领域中占有相当的比例，研究其成形机理具有重要的意义。在汽车、钢铁、航空航天、船舶、兵器等国民经济部门，随处可见半球形零件以及各种更复杂曲面零件的身影，而且，对零件外形质量的要求也越来越高。进行半球形零件拉深成形的数值分析，建立精确计算成形过程的应力、应变分布以及零件外形、厚度分布的力学模型，是研究起皱和破裂的基础，为预测成形过程可能出现的起皱、破裂等缺陷，并确定其中的一些重要工艺参数提供了理论依据。另外，进行半球形零件拉深过程的分析研究，具有典型性，是认识复杂曲面零件成形的基础。因为，半球形拉深件具有法兰区域、与凸凹模圆角相接触的区域以及不与模具相接触的悬空区域，这些典型区域可以组合出很多复杂曲面零件，研究复杂曲面零件的成形，一样需要分析这些区域。本课题是对半球形拉深的理论分析。半球形件拉深与一般平底筒形件拉深不同，文献中针对其特点，采用数值算法进行详细的理论分析的论述尚不多见。因此本文通过对半球形零件拉深成形过程分析研究，建立了一个方便快捷的数值算法对其成形过程进行模拟计算，得到零件的应力应变和外形厚度分布。这个算法可以用来研究筒形件、锥形件等成形的规律以及工艺参数对零件成形的影响，有助于进一步研究起皱和破裂，可以应用于教学，展示各种具有上述特征区域的零件的成形应力应变。另外，也希望能应用于钣金件智能化生产，对半球形零件拉深过程中的压边力等进行实时控制。1.2 国内外研究现状对于半球形零件拉深成形这样的工程塑性力学问题，计算方法有有限元法、解析方程-解析解、解析方程-数值解。有限元法数值计算是目前分析塑性成形问题最有力的工具，但对于板料成形因涉及到几何物理及边界多重非线性并存问题，计算量大，进行系统的分析需要大量计算时间。目前国内这方面的研究其中有哈尔滨工业大学金属精密热加工国防科技重点实验室的刁法玺、张凯锋进行的动力显式有限元模拟3，他们建立了适合于三维板料成形分析的显式算法的有限元数学模型,采取集中质量矩阵,用动力显式算法来分析, 基于时间中心差分格式, 使位移计算显式化, 避免了由材料、几何、边界条件等高度非线性因素引起的计算收敛问题。并根据该模型开发了动力显式算法的板料成形过程模拟的有限元分析程序。解析方程-解析解需要对成形模型进行很多假设，早期北京航空航天大学704教研室的曹增强、胡世光进行的半球形拉深件成形研究4，推导出了计算每一瞬时的应力公式，解析方法直观，便于分析影响应力的因素，容易总结出经验公式。但这种算法进行很多假设，比如板料厚度变化均匀、简单加载、理想塑性材料等条件，有些部分与实际情况不符合。解析方程-数值解虽然和解析方程-解析解一样，需要写出较符合模型的力学方程、本构方程和几何方程，但不推导最后的计算表达式，而是用计算机进行数值计算。这种算法在进行力学分析的时候，只需做很少的假设或者不做假设。通常力学方程是微分形式的，本构、几何方程是代数形式的，可以把力学方程离散化，全部变成代数方程来求解，比如D.M.Woo的算法。这种算法对函数的要求低，但它通常写为隐式的，在分析像板料成形这样高度非线性过程时，初值难以确定和计算耗时的问题尤为突出。可以把本构和几何方程微分化，全部微分形式来计算，这种算法对函数要求高，如果函数的一阶导数不连续可导，就可能得到奇异解，判断计算结果的真伪和修正十分麻烦。这种算法可以写成显式的，也可以写成隐式的，比如俄罗斯国立科技大学的研究人员，就提出了的显式算法，不需要求解大型联立方程组，也不需要迭代，适合于大规模的强非线性问题的分析，但要求研究对象的一阶导数连续。本课题对Woo D M的算法感兴趣，它容易理解，基本能满足课题的应用要求。这个算法在板料与凹模圆角相接触的区域的方程比较复杂，还可以进一步化简，更方便掌握和应用。目前，具有悬空区的轴对称零件研究中仍存在一些问题。在零件的悬空段，通常假设为锥面。如果这样假设，这部分应该是单向应力状态，但一些研究仍然算有切向应力，例如文献4 、5。悬空区实际上是具有拐点的S型曲面，计算其型面也是本课题的目标之一。在零件与凸模接触的头部区域，由于存在摩擦，而且比较大，所以应力状态不是双向等拉。通常研究要假设为双向等拉或者部分材料双向等拉，这样才能计算出结果。此区域的实际应力状态究竟如何，尚存在疑问。不做任何假设，探究头部区域的应力状态，也是本课题研究的重点和难点。1.3 研究方法鉴于半球形件的轴对称性，本文基于增量理论的逐步逼近数值法，在分析其拉深成形过程的基础上，综合考虑板材厚度变化、厚向异性、加工硬化以及摩擦力、压边力等多种因素的影响，建立成形过程不同阶段各变形区的力学模型，推导出用以计算的基本方程，编写程序，对拉深过程中应力应变状态和零件外形厚度分布进行完整的数值解析。对计算模型进行理论推导时，首先，要进行分区，建立各个区域的混合方程组；其次，对凹模圆角接触区、悬空区、球头区引入参数，尽可能使平衡和几何方程的表述更方便；最后，在基本方程的推导过程和建立数值算法的时候，将符合增量理论的本构方程和微分平衡方程离散化处理，在编程求解过程中，根据实际情况只将拉深过程分成了一个阶段计算，因此实际上采用了全量理论进行数值模拟。在整个数值模拟建立以后，通过将计算结果与其它方法（解析法、有限元模拟以及试验等）的对比来验证本算法的准确性与有效性，之后就可以用来模拟在各工艺条件下的拉深成形结果，以便我们研究工艺参数对拉深成形的影响。基于Woo D M的解析方程数值算法，尽可能不做应力应变状态的近似假设，求解半球形件拉深成形过程中的悬空区的外形以及头部球形区的应力应变分布，是本课题研究的重点，也是难点。1.4 论文构成和研究内容在数值计算分析的文章中，以D.M.Woo的分析最为全面，但在他的文章中没有考虑材料的各向异性，特别是凸缘部分和平底拉深部分分开考虑，使应力、应变分布曲线间断，误差较大。本文以半球形零件的拉深过程为研究对象，通过对其进行数值模拟计算，全面且深入的分析研究D.M.Woo的算法，改进了其中的一些精度问题，并提出了改进逐步逼近数值解法的方法。提出了采用弦截法自动求解有压边力的凸缘区域的厚度以及确定分界点。提出了确定材料在两个圆角上包角的方法。论证了采用增量理论进行解析方程数值解时产生的困难和不实用性，论证了采用全量理论来处理轴对称零件是足够精确的。比较精确地求解出了半球形拉深件的外形，确定了悬空区的S 形曲面。确定了凸模球头区的应力状况。之后与各种计算方法进行了对比分析，以确定算法的有效性，研究了各种工艺参数对拉深成形的影响，得到了一些对生产实际和理论分析有意义的结论。2 半球形零件拉深成形的数值模拟2.1 计算模型的建立2.1.1 板料拉深过程的模拟对于半球形零件拉深成形过程，首先将原始板料虚拟地划分为一定数量的环形单元(i=1，2，M)，再将连续的成形过程离散化，使整个拉深过程成为一系列的变形瞬时阶段(j=1，2，N)。图2.1是半球形件拉深过程中的一个瞬时模型，它展示了板料变形后的外形和上述的划分方法。根据这一模型，可以将整个零件划分为5个不同的区域：压边力作用下的法兰边缘区域AB，不与压边装置接触的法兰内缘区域BC，与凹模圆角处相接触的区域CD，不与模具接触的悬空区域DE，与冲模接触的球形区域EF。这5个区域的受力状态不尽相同，算法也有所差异，将分别单独分析。1-冲模；2-压边圈；3-板料；4-凹模图2.1 半球形零件拉深成形模型2.1.2 基本计算方程当进行一个数值分析的时候，必须注意的最重要的因素之一就是所用的本构方程能够精确的记录材料的行为。目前的各种材料，如铝合金等很明显展示出各向异性，那么所选的材料模型能够在成形期间恰当的记录其行为就是很重要的。因此，本次研究选用了Hill6的各向异性屈服准则。对于板料变形的情况，一般假设板料在板面内各向同性，只有厚向异性，本研究也假设板料具有厚向异性指数R，而其他方向各向同性。对于半球形拉深这种轴对称成形问题，其应力和应变增量之间的关系为 （2.1）其中，为应力强度 （2.2）应变强度增量为 （2.3）式中 、分别为板料的径向、环向和厚向应力，、分别为板料的径向、环向和厚向应变增量。R为板料的厚向异性指数。根据式（2.1），可以进一步推导出如下的应力应变关系： （2.4） （2.5）应力强度和应变强度之间的关系由实验确定，理论近似表达式可以采用Hollomon6的指数型应变硬化式 （2.6）式中 A为强化系数，n为应变硬化指数，e为应变强度；。2.2 微分平衡方程与体积不变条件2.2.1 法兰区域在拉深成形中，为防止凸缘材料起皱，需要加一防皱压边装置，在此只考虑压边力恒定的情况。由于法兰区域的材料环向受压，径向受拉，厚度增加，因此产生了与压边装置接触的等厚区域AB和不受压边力的内缘区域BC，下边将分开讨论。（1）区域AB该区域某一微体的受力状态如图2.2所示，其中、为主应力，对应方向为主方向。设此区域的摩擦系数为，并忽略摩擦力对应力主方向的影响，则也取为主应力，对应的方向为主方向。由径向平衡条件可以得到因很小，并忽略二阶及其以上小项，上式可化简为 （2.7）图2.2 有压边力的区域AB之微体受力状态图设区域AB的变形后的厚度为，将上式写成积分形式为 （2.8）设材料塑性变形时体积保持不变，由此条件有 （2.9）（2）区域BC本区域是上区域结束点到凹模圆角圆心半径之间的区域，其材料变形仍为环向受压缩，径向受拉伸，所以材料厚度由外向内逐渐减小。其区域材料不与压边装置接触，材料处于平面应力状态。此区域某一微体的受力状态如图2.3所示。由径向平衡条件可以得到因很小，并忽略二阶及其以上小项，上式可化简为 （2.10）图2.3 不与压边装置接触的法兰内缘区域之微体受力状态图将上式写成积分形式： （2.11）该区域的几何方程为 （2.12）式（2.8）、式（2.9）、式（2.11）、式（2.12）中的、以及的定义见图2.1。2.2.2 凹模圆角区域材料经上两个区域的拉伸成形后，进入凹模圆角区域CD。此区域微体的受力状态如图2.4所示，是此区域的摩擦系数。有些文献认为，例如文献7，此区域的材料经历了由直变弯和由弯变直两次弯曲变形，因此应该对径向应力进行修正，其修正量就是肖夫曼所设的分量是凹模的圆角半径（见图2.1）。如果时，这一项的影响不大，在正常的拉深条件下进行应力分析计算的时候可以忽略不计5。然而，在此区域的微体受力状态全面无误的前提下，通过微体的受力平衡方程从一个点的应力应变状态求解出的下一点的应力应变状态就也应该是正确无误的，对点的应力强加一个修正是没有道理的。这个修正量应该是宏观计算板料经过圆角弯曲后的应力变化时才考虑的6，是从变形功的角度得到的结果，一般应用于解析解。因此，本文认为采用D.M.Woo的方法进行数值计算时，再考虑肖夫曼分量是不对的，在此区域应该考虑的是如何处理主应力在厚度方向分布的不均匀性。比如，板料上表面不受力，下表面受凹模的挤压，因此厚向应力是不均匀的；径向应力和周向应力由于板料的弯曲在厚向分布上也是不均匀的6。因为板料厚度相对来说很小，像大部分文献那样，本文也假设三个主应力在厚度方向是均匀分布的，但厚向应力没有忽略，而是采用了板料上下表面外应力的平均值进行计算。忽略摩擦因素对主应力、应变方向的影响，由圆角切向和法向力的平衡条件可以得到平衡方程：其中，是凹模圆角圆心C到板料中面的距离 （2.13）事实上，我们就是采用板料中面上的点的应力应变分布以及厚度高度等来描述整个板料的。因为、很小，所以有、，又， 图2.4 凹模圆角区域微体的受力状态图并略去二阶及其以上高阶小量，消去N，可化简整理为： （2.14）另外，将法向平衡方程两端除以微体贴模面积（这里用板料中面的面积近似微体贴模面积，可以证明，它们差个高阶小量），可以进一步得到下面的公式： （2.15）上式中的就是板料下表面的表面应力。将式（2.14）进行积分，得到积分式 （2.16）式中，、变形瞬时板料的环带单元在凹模圆角处的包角。板料一个环带单元在圆角区域的体积计算可以用单元中面的面积与平均厚度乘积。由积分公式已知，任意曲线的旋转表面积，等于该曲线的长度乘以曲线重心的旋转周长8，再根据体积不变条件，就可以得到下面的几何方程 （2.17）假设材料在凹模圆角上的包角是，显然，利用脱模点D的应力状态以及几何参数 ，可以计算拉深力F为 （2.18）2.2.3 悬空区域区域DE处于凸凹模之间，材料不与模具接触。该区域的实际外形是具有拐点的比较复杂的S形曲面，为了求解这个曲面形状，必须采用足够精确的微体模型进行受力分析。我们知道，某时刻拉深力F其实可以用悬空区域中任一点的应力状态和几何信息计算，形式同式（2.18）。显然F的方向是在拉深方向（竖直方向）上的，微体上截面的径向应力必须有向上的分量，那么，该区域的微体模型只有图2.5（a）和（b）这两种形式。图2.5 悬空区域的微体形态图2.6是该区域微体的受力状态图，图2.6（a）、图2.6（b）分别对应于图2.5的（a）、（b）。该区域的受力与区域CD相似，但没有厚向应力和摩擦力。根据微体（a）的受力状态，由力的平衡条件可以得到：切向方程法向方程根据微体（b）的受力状态，由力的平衡条件可得切向方程法向方程因为、很小，所以有、，又，（对应图2.6（a）、（对应图2.6（b），并忽略二阶及其以上小项，将上述的切向方程化简，可以得到一个统一的切向方程，它适用于整个悬空区域DE： 图2.6 悬空区域微体受力状态图。（a） 反向双曲微体的受力状态（对应图2.5（a）；（b） 同向双曲微体的受力状态（对应图2.5（b）。可以看到，两种微体的中面都是双曲面，其差异就在于，图2.5（a）微体中面的弯曲中心在曲面的两侧，而图2.5（b）微体中面的弯曲中心在曲面的同一侧。假设图2.5（a）微体的一个曲率半径为负的，图2.5（b）微体的为正的，那么它们的法向平衡方程也可以化简为一个统一的形式 （2.19）将切向平衡方程进行积分，写成积分形式 （2.20）本区域环带单元的体积近似用锥台的体积计算，建立体积不变条件表达式为 （2.21）以上曾提到拉深力F可以用悬空区任意一点的径向应力和几何状态求解，事实上悬空区不受任何外力，所以都与公式（2.18）得到的拉深力F相等，这就得到一个整体平衡方程，可以求解任一点的法线与竖直方向的夹角（图2.5中的），其表达式为 （2.22）2.2.4 凸模球头区域凸模球头区域EF同凹模圆角区域CD的情况很相似，其某一微体的受力状态如图2.7所示。设此区域摩擦系数为，忽略摩擦对应力应变主方向的影响，则由力的平衡条件可得其中是冲模球头球心到板料中面的距离 （2.23）因为、很小，、，且，略去二阶及其以上高阶小项，消去N可以得到 （2.24）对于法向方程，照区域CD做同样的处理，也可以得到 （2.25）上式中的是板料上表面的表面应力。图2.7 凸模球头区域微体的受力状态图将式（2.24）进行积分，得到积分公式 （2.26）同区域CD的考虑方式，此区域的体积不变条件为 （2.27）2.3 数值计算方法的建立2.3.1 应变计算计算中，把应变取为有限应变增量来考虑，对于j阶段任一单元i边界的应变为： （2.28）应变增量表示为： （2.29）2.3.2 计算方法和编程技术关键成形过程中，由于各区域所建的平衡方程、几何方程和边界条件不同，因此分析计算的方法和技术问题的处理也不同，应分区考虑计算。算法的基本思想是，在已知环带单元上边界点的应力应变以及几何关系的条件下，采用了逐步逼近数值解法，求解下一个点的应力应变。其计算流程如下所述。（1）区域AB由于此区域的厚度不变，先假设一值。在当前板料上取环带单元步长，设此步长为（该值应比较小，比如取0.01mm），则下一点的半径就可以由得到。利用式（2.9）可以得到该点在原始板料上对应的点半径，然后点的应变就可以由式（2.28）和式（2.29）给出。利用应变可以从式（2.3）、式（2.6）得到该点的应变强度和应力强度。将上边得到的这些量代入到式（2.4）、式（2.5）中，就可以求出点的两组主应力之差：、。为了进一步求出每一个主应力，用梯形公式积分后，整理得 (2.30)经过以上的计算，式（2.30）右边的量都是已知的，就可以求出点的径向应力，进而知道周向应力和厚向应力。需要说明的是，板料的起始点A点的径向应力为零，所以不需要用式（2.30）就可以求出A点所有的应力，进而循环求解其他点的应力应变。当时应停止循环，因为区域AB有压边力作用，厚向应力小于等于零。像求解点应力应变一样，我们也可以用点的应力应变求解这两点中间任一点的应力应变，所以可以将看作是r（）的函数，而点的应力应变可以看作是函数的参数。在处，其函数值，在处，函数值，从这两点开始，采用弦截法（本人对于本文的问题而言，二分法不够方便）就可以循环求解出函数的零点，这个点就是区域AB和区域BC的分界点B点。另外，如果（的定义见图2.1），表示此点落到了凹模圆角区域CD，而B点不可能位于此区域，最多和C点重合，所以这种情况不用步长求，直接让等于并结束循环。在计算出点和的应力状态后，可以求解出压边装置施加在这个环带单元上的压边力，用梯形公式表达为 (2.31)程序走到B点后，可由下式 (2.32)计算出总压边力，并与预先给定的压边力Q比较，如果相等，则表示假定的是对的，否则就要改变的值，重新计算整个区域AB。求解正确的时，也可以让程序用弦截法自动求解。上述的整个过程，可以看作是由自变量经过一个带很多参数的函数求解的过程。记 (2.33)对应于每一个都有一个，问题就变成求使等于零的值。当然，用弦截法时要先给函数值一正一负的两个点。可以采用的一个点是，因为法兰区域的材料周向收缩，厚度一定是变厚的，只有压边力非常大的时候才可能使厚度不变，在实际拉深中一般没有那么大的压边力，所以在处的一定大于零；另一个点可以这样得到：假设没有压边力，则区域AB退变成一个点A，也就是点A处的厚度，将代到式（2.5）中就可以解出一个，显然这个值使小于零。图2.8显示的是该区域的计算框图。（2）区域BC经过上一区域的计算，分界点B的位置以及应力应变状态都已经知道，以此为边界条件，可以继续求解区域BC其他点的信息。由于每个区域的变形状态不同，所以计算精度的要求也应该不同，比如此区域厚度变化，所以步长应该比区域AB小。以后为了根据情况方便调整精度，每个区域都有各自的步长。记该区域步长为，则下一点半径，初设，由几何条设定式（2.9）：式（2.28）、式（2.29）：、式（2.3）、式（2.6）：、式（2.4）、式（2.5）：、式（2.30）：、或定出、求是否是否下一区域图2.8 区域AB的计算框图件(2.12)可以得到该点在原始板料上对应的点半径，由式（2.28）和式（2.29）给出其应变，再根据式（2.3）和式（2.6）求得应变强度和应力强度。将上边得到的这些量代入到式（2.4）、式（2.5）中，就可以求出点的两组主应力之差：、，由于这个区域的厚向应力等于零，所以可以进一步得到和，记。另外，用梯形公式改写积分，得 (2.34)将已知条件代入到式（2.34）后可以得到一个值。如果和满足一定的误差条件（视所取精度而定），则设置合理，可顺序求解下一节点。否则由式（2.5）当时，将代入得： (2.35)改变的初设值，重复计算直到达到精度要求。当时不用步长求这个点，直接让计算出边界点C的应力应变后终止该区域的计算。此区域的计算框图见图2.9。（3）区域CD凹模圆角区域CD的计算思想同区域BC相似。为了方便计算，这里考察的自变量不用点的半径，而用材料中面上的点与凹模圆角圆心的连线和竖直方向的夹角，由于板料厚度变化相对于凹模圆角半径来说很小，所以也可看作是该点处法线方向与竖直方向的夹角。显然边界等于零。设积分步长为，则下一点。初设该点厚度，由几何条件（2.17）可以得到该点在原始板料上对应的点半径，而该点现在的半径可由初设式（2.12）：式（2.28）、式（2.29）：、式（2.3）、式（2.6）：、式（2.4）、式（2.5）：、式（2.34）： 式（2.35）：是否是否下一区域图2.9 区域BC的计算框图 (2.36)得到，将代入到式（2.13）可以得到上式中的。然后由式（2.28）和式（2.29）给出其应变，再根据式（2.3）和式（2.6）求得应变强度和应力强度。将上述求得的条件代入到式（2.4）中，得到主应力之差。根据梯形公式将式（2.16）化为数值积分的形式，并加以变形，使方程的右边只包含目前已知的应力差，则式（2.16）成为： (2.37)其中 (2.38) (2.39) (2.40) (2.41)以及 (2.42) (2.43)经过以上的计算后可以得到一个值。由于在这个区域考虑了厚向应力，所以要采用与区域BC不一样的处理方法，但基本思想是一样的，这一点会在本节的最后详细讨论。这里直接用得到点的径向应力：，然后代到上边求得的中解出：，将、以及上边求得的一些几何参数代入法向方程（2.15）中去，得到该点对应的板料下表面外应力。根据假设，该点的厚向应力。现在记，注意到利用点的应变可以直接从本构方程式（2.5）得到一个值，如果和满足一定的误差条件，则设置合理，可顺序求解下一节点。否则将代入到本构方程（2.5）得： (2.44)改变的初设值，重复计算直到达到精度要求。关于该区域的计算终止条件，可先假定板料在凹模圆角上包角的大小，当时不用步长求这个点，直接让计算出边界点D的应力应变后终止该区域的计算，最后，利用式（2.18）求出拉深载荷F。假定的是否合理，需要在后续区域EF进行判断。此区域的计算框图见图2.10。4）区域DE悬空区域DE虽然具有复杂的形状，但其受力状态和区域BC十分相似，从式（2.11）和式（2.20）可以看到，其积分方程相同。设该区域的积分步长为，则下一点的半径，初设，另外注意到几何方程（2.21）中用到环带单元两节点值的平均值，所以还要设，然后就可以从式（2.21）中得到该点在原始板料上对应的点半径。此后的计算同（2）中所述相同，另外，在用修正的同时，还可以由式（2.22）用来修正： (2.45)当重复计算达到精度要求之后，可以利用式（2.19）得到该点的曲率，期望这个曲率能更加精确的确定悬空区域的外形。此部分计算的停止条件是初设式（2.17）：，式（2.13）、式（2.36）： 式（2.28）、式（2.29）：、公式（2.3）、式（2.6）：、式（2.4）、式（2.37）（2.43）、式（2.15）：、式（2.5）： 式（2.44）：是否是否下一区域图2.10 区域CD的计算框图边界点E的计算同区域AB的B点计算方法相同，构造的函数是。记E点处的值为，这是材料在凸模上的包角，也是下一区域的边界条件。悬空区域DE的计算框图如图2.11所示。（5）区域EF凸模球头区域EF是半球形拉深件的最后一个区域。通常研究要假设该区域应力应变状态为双向等拉或者一部分双向等拉，这样才能计算出结果。本文的主要目的就是采用最为符合实际的模型，较精确的计算拉深件的应力应变分布状态，因此这里为求解区域EF的实际应力应变状态，不采用过多的假设。该区域的受力状态和区域CD十分相似，先设积分步长为，则下一点。初设该点厚度，同（3）中的考虑方式， 采用与区域EF对应的方程，确定顺序求解或改变值重新计算，求解出点的应力应变状态以及厚度。其具体采用的对应的方程见计算框图2.12。用梯形公式改写本区域的积分式（2.26）如下： (2.46)其中 (2.47) (2.48) (2.49) (2.50)下一区域初设、式（2.21）：式（2.28）、式（2.29）：、式（2.3）、式（2.6）：、式（2.4）、式（2.5）：、式（2.34）： 式（2.35）（2.45）：、是否是否定出E点图2.11 悬空区域DE的计算框图此区域的终止条件理论上应该是等于零，但零点没办法计算应变，所以相对于半径70mm左右的板料来说，可以用作为终止条件，终止点F的半径可以近似为0.1mm，另外根据球头半径可以换算出终止角。区域EF的计算框图如图2.12所示。在步骤（3）讨论凹模圆角区域的计算时，曾假设了一个包角大小，这个值是否正确，可以从两方面判断。第一，在整个拉深件的数值计算过程中，每个环带单元的体积都要求在原始板料上有一块体积与之对应。板料的总体积是固定的，如果值设置的不恰当，就会导致材料不够用或者过多，这样的话模拟计算在半径到达0.1mm之前就无法继续了，所以应该调整值，使程序能计算到终止点F而不出现数据溢出等错误；第二，在球头的顶点处，径向应力和周向应力相互垂直，并且都垂直于冲模的中心轴，根据零件的轴对称性，这两个应力一定相等。作为替代顶点的0.1mm处的终止点F的径向应力和周向应力也应该近似相等，如果计算到最后，不满足条件，就要重新设置的大小，并返回步骤（3）重新计算。2.3.3 几点讨论本文的主要任务之一就是深入研究D.M.Woo的数值算法。以上对半球形拉深件的数值模拟，可以体现出逐步逼近数值解法的基本思想，可以说，D.M.Woo的算法也只是这种思想的体现， 而不是全部，深刻体会这种思想，可以寻找到更多的、在各方面都更加好的解析方程数值解法。首先，我们知道，有一种方程求根的计算方法叫不动点迭代法9：对于求解方程，可以改为求解，函数被称为迭代函数。初设一个x值，代入迭代函数后得到一个函数值，将这个函数值作为新的x，并一直迭代下去，所有函数值组成了一个数列。如果所给初值以及迭代函数满足一定的条件，那么这个数列是可以收敛的，记该数列收敛于，则有，这个点就是，也是方程 初设式（2.27）：，： 式（2.28）、（2.29）：、公式（2.3）、（2.6）：、公式（2.4）、（2.46）（2.50）、（2.25）：、公式（2.5）： 式（2.44）：是否是否结束图2.12 球头区域EF的计算框图的解。从区域BC开始，主要的计算过程可以概括为，初设一个厚度值t，经过一系列的运算后得到一个新的厚度值，然后以代替t迭代反复运算直到达到精度要求。可以将这一系列的运算看作是迭代函数，逐步逼近数值解法实质上就是迭代法。比如，在区域BC的计算中，D.M.Woo以及本文都是将代入到本构方程（2.5）中得到应变增量，然后再得到一个，这是一种迭代函数；如果直接用除以从式（2.5）得到的，也可以得到一个，这也是一种迭代函数，本研究通过编程计算已经证明这种方法也是收敛的。再比如，计算区域CD时考虑了厚向应力，因此采用了与D.M.Woo算法不同的计算过程，这也算一种新的迭代函数。能否迭代求解出正确的x值，以及收敛速度的快慢，是取决于的好坏的，如果要从计算速度和适用性等角度来改进逐步逼近数值解法，可以考虑寻找更好的迭代过程。进行半球形拉深件等轴对称零件以及更复杂的零件的数值分析，需要用到很多的方程，包括本构方程、几何方程、微分平衡方程，如果在众多的方程中组合出最好的迭代过程，是值得研究和探讨的。其次，D.M.Woo算法中经常用一个乘积进行循环控制，认为计算的厚度值如果是真实值，从不同的方程中求得的应该相等。本文认为，如果用于数值解析的解析方程组不是只有唯一解，这是不合理的。比如除了实际的和t的乘积满足与相等外，可能还有另一组和t也满足，用这个乘积作为循环控制条件，显然可能得不到真实解。通常认为轴对称等拉深件的解析方程组只有唯一解，本文也倾向于这种提法，在此前提下，可以有很多种循环控制条件。比如，1）区域CD的计算中，控制条件是应力差；2）前已叙述，此种算法可以看作是不动点迭代法，当然就可以直接拿新算的厚度值和之前的厚度值比较，如果收敛，计算前后的厚度值应该相等，所以厚度t也可以做循环控制条件；3）从当前厚度t可以计算应力应变并由基本方程可以得到一个硬化模数，另外，用微分平衡方程并从本构方程中可以解出一个硬化模数，比较这两个值也可以进行循环控制。不同的控制条件得到的计算精度不同，相同的误差条件下，可能用硬化模数进行循环控制比用厚度进行控制计算的结果更为精确，而且可以应用于非轴对称的零件拉深计算中。值得一提的是，误差条件都应该采用相对误差。最后，本算法仍遗留了一个问题，让程序自动求解正确的应力应变状态分布存在困难。在计算区域CD的时候，假设了一个包角大小，计算完区域EF后，也给出了判断该值正确与否的条件。第一个条件中提到，如果的值偏离真值过多，程序就会出错，比如循环不收敛而导致数据溢出或者数据范围错误等，这就给让程序自动求解的值带来了困难，因为数据溢出了程序就停止了，从而不能自动改变初始值重新计算。另一个困难是，这个包角的初值很难给出，只要与真值偏差0.1弧度甚至更小，算法都很难全部运算完。本文为了求解更接近实际的应力应变状态，一切从实际受力状态出发，不能有过多的假设，比如在球头区域EF没有假设有双向等拉区域，这种做法必然导致对包角的初值要求很高。本算法仍是根据计算结果手动改变的大小，直到满足判定条件。如何解决上述的两个困难，或者采用其他方法使算法能自动求解出正确的结果，也是值得人们继续研究下去的。2.4 半球形零件的型面计算经过以上的数值计算，半球形拉深件的应力应变分布和厚度分布已经知道，也知道板料中面上的点的法线方向和曲率，用这些信息可以较精确地？算出零件成形后的型面。基于零件的轴对称性，零件的纵截面就可以代表整个零件，所以只需要找到一种表示零件纵截面的方法。首先，将法兰的下表面定为高度的零点，对于任一半径处，用高度y定出零件纵截面中线上的点，然后沿着该点的法线方向向上和向下各偏移，就可以得到零件的上表面和下表面上的点，这样就表示出了整个零件。前边已经计算出了中线上点的法线方向与竖直方向的夹角以及厚度t，所以这一部分只需要给出高度y的计算方法。在法兰区域，即区域AB和区域BC，厚度的一半就是中面上的点的高度。在凹模圆角区域，可以用凹模圆角的圆心坐标来定： (2.51)上式中的是凹模圆角圆心的坐标，显然。在悬空区域DE，已知的几何条件有法线方向、厚度以及曲率，但这些条件并不能确定点的高度，点的曲率中心无法确定，这里只能用差分计算。由角可以知道该点的斜率，所以用于计算该区域中线高度的公式如下： (2.52)原本期望能用来精确计算外形的曲率对这种方法没有帮助，不过能证明半球形拉深件的外形曲线是二阶可导的，所以零件外形上不应该有尖点或突变。经过悬空区域的计算，边界点E的高度是已知的，记为，同区域CD相似，我们利用球头的球心来计算高度： (2.53)上式中的是E点与球心的距离。 以上是半球形拉深件数值模拟的过程和一些技术问题的处理，按照上述的计算步骤和计算框图可以编写C语言程序进行计算。该计算程序由五部分组成，能够快捷的计算出零件的应力应变分布和外形，可对各工艺参数进行研究。3 计算结果与对比分析3.1 增量理论与全量理论的实际应用问题板料拉深成形时，材料基本处于受两向以上的应力作用的复杂应力状态，其塑性变形区的应力应变关系是相当复杂的，目前常用的有全应变理论和增量理论两种应力应变关系。增量理论表示在塑性变形的某一个瞬间应变增量与主应力之间的关系，经过积分便可以把变形过程的特点反映出来，它比较接近于实际的情况。全量理论仅仅表示塑性变形终了时的主应力和主应变之间的关系，它不能反映出变形过程中应力与应变的变化过程所产生的影响。假如塑性变形过程中的主应力方向不变，而且各应力间的比例也保持不变，全量理论和增量理论的计算结果是一致的，所以在这种情况下完全可以应用全量理论。此外，在单调的塑性变形过程中也可以应用全量理论。增量理论在计算上引起的难度很大，尤其在冷变形硬化时，计算就更复杂了。本文在数值模拟的理论推导部分也采用了增量理论，但在编程计算的时候，发现有很多困难，如果严格按照D.M.Woo的算法进行，又发现很不实用。如果采用增量理论，首先，不能在当前变形的板料上取步长，因为板料还在拉深过程中，所有的点的状态都是不断变化的，所以只能在原始板料上取积分步长，这样就丧失了很多采用当前步长为计算带来的巨大方便。其次，所有的分界点十分难以确定，比如在j时刻确定的分界点D，到了时刻就不再是分界点，就完全没用了，只好在原始板料上重新插入一个分界点（且不说这个新分界点如何找到），并从第一阶段重新计算，十分浪费计算时间。如果不这么关心分界点，每一个j阶段都近似用原始板料上的某个节点代替分界点，那么，划分的阶段越多，累计的误差越大，与采用增量理论的目的相背，所以分界点的问题是无法回避的。再次，从上一章的分析讨论中知道，材料在凹模圆角上的包角十分重要，在j时刻，计算想要顺利进行，包角大小的初值要十分接近真值，但这个值在下一时刻就不能保证程序进行到底了。将变形阶段划分的很细小，保证在时刻也能满足程序进行到底，则可以利用顶点双向等拉的条件自动修正逐步逼近正确的，这是让计算机自动求解材料在凹模圆角上包角大小的一种方法，但实际上D.M.Woo的算法对这个值要求十分严格，阶段划分的要非常多，而阶段划分的越多，寻找边界点的工作就越复杂。文献7采用增量理论进行了筒形件的拉深计算，其处理技术问题的方法与上述方法相似，为求解某材料的极限拉深系数，将拉深过程划分为八个阶段，仍然需要进行十几个小时的计算，所以增量理论的求解方法在实际生产中推广运用有一定困难。筒形件和圆锥形件具有筒底的双向等拉区域，但半球形零件没有这样的特殊区域，半球形零件的球头区域的数值模拟计算要继续在变化中靠近零半径点，所以如果不做额外的假设，用增量理论进行半球形零件的数值模拟计算更难以在实际生产中推广运用。文献7进行了筒形件的拉深试验和板料的单向拉伸试验，并用增量理论和全量理论进行了数值模拟。无论是在筒形件拉深的条件下还是在单向拉深的条件下，两种理论数值模拟的理论结果都和试验结果符合的很好。两种理论对应力的计算误差小于2%，对应变的计算小于6.6%，所以相对于全量理论，采用增量理论在结果上并没有实质性的改善，以全量理论代替增量理论，进行拉深成形计算是可行的。虽然相对于筒形件，半球形零件拉深过程的数值模拟更难以运用增量理论，但其变形复杂程度以及应力状态的复杂程度却更简单，更接近于单调的塑性变形过程，所以用全量理论代替增量理论进行半球形零件的拉深成形计算，从结果上说也一定是可行的。而且，采用全量理论大大减少了计算的复杂性。本研究虽然利用增量理论推导了方程，但编写程序的时候采用的是全量理论，上一章节中的一些技术问题处理方式，也是只适用于全量理论的，该程序运行一遍只需要几秒的时间，调节出正确的值并得到最终正确结果也只花费最多五分钟的时间。另外，采用全量理论后可以更自如的控制各个变量，求解边