板材成形回弹数值分析的静力隐式方法

中国机械工 程第 13卷第 17期 2002年 9月上半月 文章编号 ：1004—132 X(2002)1 7—1458—04 板材成形回弹数值分析的静力隐式方法 蔡中义 李明哲 李湘吉 摘要 ：给 出了板材成形加 载过程、卸载过程数值模拟的静 力隐式有限元 方法，通过平衡 迭代 求解有限元方程 ，并采用隐式算法进行应 力积分及接触 压力计算。这种 方法应力计算准确 ，因而回弹计算结果具有较 高的精度。数 值计算表明，纯弹性板大变形卸载后能够完全恢复到原始形状，说明该方法 对几何非线性问题的分析是非常准确的。最后对板材柱面、球面成形的回弹 进 行 了分析 。 关键词 ：板材成形 ；回弹 ；静力隐式 ；弹塑性有限元 中图分类号 ：TG386；TB115 文献标识码 ：A 回弹是板材成形 中不 可避免 的现象 ，是决 定 零件最终形状的重要因素 。如果回弹量不能有效 地控制 ，将严重影响零件的几何精度 ，成为成形缺 陷 。回弹计算是板材成形数值模拟 的一个十分重 要课题[】 ]。传统的回弹数值模拟通常采用动力显 式有限元方法L3]，由于显式时间积分格式不形成 总刚度矩阵，不进行力的平衡迭代，应力计算不准 确，因而无法给出准确的回弹计算结果。 回弹计算是在加载过程应力分析的基础上进 行的。为精确计算加载应力必须采用静力隐式有 限元格式L4]。这种分析格式采用迭代方法求解非 线性有限元方程 ，迭代过程中不断更新切线刚度 ， 直至节点不平衡力及位移增量小 于给定值为止 。 由于每一增量 步都进行平衡迭代 ，因而加 载过程 的应力计算精度高，以此为基础进行的卸载回弹 计算也 比较准确。 收稿 日期 ；2001一O1— 16 基金项目；“十五”国家科技攻关课题(2001BA203B11)；教育部跨 世纪优秀人才专项基金项 目 蔡 中义 副教授 1 静力隐式有限元理论 1．1 有限元方程 板 材成形是一种拟静态 的变形过程 ，这一变 形过程由各物质点在空 间的运 动构成 ，各物质点 的运动可 由位移来描述 。在任意时刻 ，物质点都处 于平衡状态。采用有限元格式空间离散 ，并考虑到 工具与板材间的接触作用后 ，平衡方程可简单表 示为 R + R“一 R“‘+ R (1) 式 中，R 为节点的 内力 向量 ；R“‘为作用在 节点 的外力 向 量； 及R 分别为由接触作用产生的节点内力向量及外 力向 量 。这些 向量 都 是 由各单 元 对节 点 的作 用 累加后 得 到的 。 在板材成形分析 中，最令人感兴趣的是变形 。 因此 ，需要从上述方程 中求解 出位移 。然而 ，由 于 R 及 R“是位移的复杂非线性 函数 ，方程不能 直接求解 ，必须将外载荷划分 成一 系列的载荷增 量 ，将非线性方程组线性化 。 从 t时刻给定 的平衡状 态出发 ，施加载 荷增 量后产生变形，t+ At时刻新构形下的平衡方程 F4J 彭志辉 ，余旭凡．不锈钢覆铝板成形极限的理论分 析 和实验 验证．中国有 色金属学 报 ，1999，9(2)：305 ～ 312 F5J 周磊．薄壁管液压胀形成形理论及试验研究：[硕士 学位 论文 政研论文下载论文大学下载论文大学下载关于长拳的论文浙大论文封面下载 ]．秦皇岛；燕山大学．2002． [6] 王祖唐，关延栋．金属塑性成形理论．北京：机械工 业 出版社 ，1989：47～140 [7] 杨玉英．大型薄板成形技术．北京：国防工业 出版 社 ，1996：146～163 ·1458· (编辑 周佑启 ) 作者简介：张 庆，男，1 949年生。燕山大学(河北省秦皇岛市 066004)机械工程学院教授。研究方向为材料加工工艺与装备、特 种成形技术。获省部级科技进步二等奖 3项、三等奖 3项。周 磊 ．男 ，1 977年生。燕 山大学机械工 程学院硕士研究 生。赵 长财 ， 男 ·1 965年生。燕 山大学机械工程学 院副教授 。仇 平。男 ，1969 年生 。燕 山大学 机械工程学院硕士 研究生 。 维普资讯 http://www.cqvip.com 板材成形回弹数值分析的静力隐式方法——蔡中义 李明哲 李湘吉 为 ( + AU)+ ’( + AU)一 + △ + + △ 。 (2) 式(2)左边展开成Taylor级数，略去高阶项， 得到关 于位移增量 △ 的线性方程 ( T+ Kc)AU 一 岸 + 苒 一 一 ；。 (3) 式 中 ，K 一 ； 一 一 Rim( 一 Rci( 一 + △ ； 一 + △ 。 采用 Newton— Raphson方法求解式(3)，通 过迭代使不平衡力 ” 一 + 一 ( + △ )一 ( + AU) (4) 满足给定的收敛条件R ≤Tol ，Tol 为检验迭 代 过程收敛的不平衡力残余值允许值，通常取 7 ol 一 10一 ～ 1O一。。 1．2 应力积分的隐式算法 每一增量步进行平衡迭代时，都必须首先确 定 变形状态 ，从而计算 出应力 。首先 假定 弹性状 态 ，计算出 t+ 时刻 的试算应力 r + im 一 + D △ (5) 式 中， 为 t时 刻的应力 列 向量 ；D 为 弹性模 量 ，△ 为应 变增量 。 当材料进入塑性变形状态后 ，按 弹性状态计 算的试算应力点将落到屈服面之外。这时，应力点 应修正到屈服面上，真实应力应为嘲 + 一 +』：fD (d —d ) (6) 式 中， 为 t时 刻的塑性 应变 。 直 接计算 式(6)是 比较困难 的，这 里给出隐 式计算方法。式(6)可写成嘲 + 一 G(△ ) (7) 式 中，G( )一 (，+ Z~CH)一。； 为服 从 von— Mises屈 服准则的等向强化材料的系数矩阵； 为塑性因子增量， △ 一△ ／ ； 为等效应力。 另外，t+ZXt时刻的真实应力应满足屈服函 数，即 F(△ )一 ( +rim ) G G t+iA 一 ( )一 0 (8) 式中， 为材料的单轴屈服应力；己为等效塑性应变。 式 (8)中的△ 可由Newton— Raphson迭代 方法求解 ，计算过程为 △ +I一△天一 F(△ t)(aF／a△ ) ] OF 一 一 ( ) G ( D G + } (9) HGD ) G 一 ， j 式中，H 为材料的强化模量。 由上式求得塑性 因子增量 后 ，由式 (7)即 可计算 t+ ZXt时刻的应力 + 。 1．3 处理接触边界的 Radial— Return方法 接触边界的处理是板材成形数值模拟的另一 关键技术问题。接触约束条件由罚函数法f6 引 入，并假定接触界面满足 Coulomb非经典摩擦条 件[7]。在罚 函数法 中允许 接触界面有法 向穿透量 g 及微小的切向微观滑移 gr。法向接触压力与穿 透量成正比，在粘着接触状态下切向接触压力与 微观滑移成正比。 在 每一增量步 中，首先按粘着接触状态计算 切向接触压力的试算值，当接触进入宏观滑动状 态后，则采用 Radial—Return方法修正到滑动屈 服面上 。按粘着接触 的罚函数法 ，t+ 时刻的接 触压力为 f 一f + ⅣAg^_1 一 f + △gr 7 。’ 式中， 女t分别为法向、切向罚参数。 式(1O)的第 一式为法 向接触力 ，第二式为切 向接触力 (摩擦力 )的试算值 。 根据 Coulomb定律 ，定义滑动函数 (fr，t ， )：一 llt丁ll一 It 『 (11) 式 中 ， 为摩擦 系数 。 如果 ( ， ，f )<0，接触点处于粘着状 态 ，则 f 一 i ；如果 ( ， ，f )≥ 0，则接 触点进入宏观滑动状态 ，需要对接触力进行修正 ， 使 (_『+ ， ，f )一0，f f 分别为t+zXt时 间的法向、切向接触力。修正过程见图 1。 J ： H： 0 J L、 ～ ＼ t 1、t 2—— 接触力 沿 2个切方向的分力 图 1 计算接触压力的 Radial— Return方法 此时修正的切向接触压力为 ir fT+ 一 ； I (12) II ‘v ll 2 回弹计算方法与步骤 为准确反映卸载过程的真实情况，采用逐步 卸载 的方法 ，考虑几何非线性效应 ，并且在卸载计 算中采用弹塑性本构关系，这样可以计算可能出 现 的卸载屈服 问题 。卸载 回弹分析仍采 用平衡迭 代方法求解有限元方程，每一增量步的迭代过程 如下 ： (1)边界条件初始化，解除约束，变量赋初 值 。1一 i，前一增量步的结果赋予 叮。及 ￡f等 ；计算 ·1459· 维普资讯 http://www.cqvip.com 中国机械工程第 13卷第 1 7期 2002年 9月上半月 应变矩阵 8 及切线刚度矩阵 K，等。 (2)计算卸载位移增量 AU 、应变增量 △毛 AU．一 一 K7 △R “ △ — B。AU． (3)确定变形状态 ，计算弹塑性应力 ffi+1。 (4)计算节点不平衡力 R ’ r I — I口 +ldV J V R 1一 R + AR； l — l — ’ i 1 (5)收敛检验 。如果 R ≤ Tol ，迭代收敛 ， 转向(1)，进入下一增 量步 ；否则 + 1一 ，转 向 (2)，继续平衡迭代过程 。 3 计算实例 根据上述静力隐式算法开发了有限元专用软 件吲，进行了考题计算及实 际的模具成形过 程与 多点成形过程的数值分析。 3．1 纯弹性板大变形的回弹 采用弹性本构关系可分析板材弹性非线性大 变 形 问题。无论受 载后变形多大 ，变形情况如何 ， 根据弹性 变形 的特 点 ，卸载后都应恢复到原始 的 未受力状态 。 (1)球形模 具作 用下板大变形的回弹 半径 为 250mm的球形模具作用于四角的简支方板 ，尺 寸为 140mm × 140mm × 2．0mm；材料弹性模量 为 E一 0．2TPa，泊松 比 一 0．3。利用对称性 ，计 算时只取 1／4，划分成 5× 5个四结点等参壳单 元。采用 50个增量步加载 ，并采用 同样增量步数 卸载。图 2给出了加载过程结束时板的变形及卸 载过程结束后(回弹后)板的形状。 (a)加载后的变形 (b)回弹后恢复初始形状 图 2 球面成 形加载 及回弹后 的形 状 表 1给出了加载结束时及卸载结束后中面z 方向的最大位 移、上表 面的最大等效应力及最大 等效应变的计算结果 。 表 1 球面成 形 回弹前 、后的计 算结果 方 向的最大 最大等效 最 大等效 位移 (mm) 应力 (MPa) 应变 加载结束 一 19．14 457，26 1．29× 10 回弹后 4．94× 10— 1 6．59 3．40× 10 (2)多点柱模作 用下板大 变形的回弹 在板 ·146O · 材的多点模具成形 中，接触问题 比整体模具成形 时更为复杂 ，接触边界处理对 回弹结 果的影响也 更大 。 图 3为一圆柱形多点模具作用的方板，多点 图 3 圆柱形 多点 模具作 用的方板 模 具 由 10× 10个基本体组成 ，基本体群的圆柱 形包络面曲率半径为200mm；方板尺寸为140mm × 140mm × 1．5mm，材 料 的 弹性 模 量 为 E 一 1．0TPa，泊松比 v一 0．3。 计 算 时 取 1／4，划 分成 10× 10个 四结 点等参 壳 单 元 ，采 用 100 个 增量 步加 载 ，并 采用 同样增量步数 卸载 。图 4为加 载 结束时及卸载结束 ＼ — ／ (a)加 载后的变 形 (b)卸载后恢复初始形状 图 4 柱面 多点成 形加载及 回 弹 后的形状 后(回弹后 )板的形状 。 表 2为加载结束时及卸载结束后中面z方向 的最大位移、上表面的最大等效应力及最大等效 应变。 表 2 柱面多点成形回弹前、后的计算结果 方 向的最 大 最大 等效 最 大等效 位移 (ram) 应 力 (MPa) 应变 加载结束 一 l1．51 985．85 6．75× 10一 回弹后 2．90× 10一 5O．O2 3．3O× 10 从表 1、表 2的计算结果 可以看出 ，回弹后 z 方 向的最大位移、最大等效 应力及最大等效应变 都 非常小。2个计算实例 回弹后 的最大位 移仅为 加载结束时的 0．025％ 左右，最大等效应力分别 为加载结束时的3．6 oA及 5 ，最大等效应变分别 为加载结束时的2．6 及4．8％。这 2个计算结果 表明，对于纯弹性的大变形问题 ，卸载后几何形状 能够完全恢复到未受力的初始状态 ，应力、应变等 也能基本恢复初始状态 ，这说明对于纯弹性问题 的回弹计算是 比较准确的。 3．2 板材成形的回弹 在板材成形过程中材料产生 弹塑性 大变形 ， 卸载后有塑性变形存在 ，不能恢 复到板材未受力 的初始状态 ，但卸载后应力、应变等分布都将发生 变化 。 维普资讯 http://www.cqvip.com 板材成形回弹数值分析的静力隐式方法——蔡中义 李明哲 李湘吉 采用基于静力隐式算法开发的有限元专用软 件分别计算板材的圆柱形模具成形及多点成形的 回弹 问 题。成 形 的板 材 为 140ram × 140ram × 2．0ram的方板，材料的弹性模量为 E：0．2TPa， 屈服应力 一IOOMPa，泊松 比 =0．3，材料为等 向线性强化，强化模量 H 一 1OGPa。 (1)柱形模具板材成 形的回弹 圆柱形模具 如图 5a所示，曲率半径为 250mm。根据对称性计 算只取1／4，划分成 5×5个等参壳单元，见图5b。 (a)柱 面成形模具 (b)有 限兀网格 图 5 整体模具 的板材成 形 加载过程采用 100个增量步 ，并 采用 同样增 量步数卸载。 图 6为加载过程结束时及卸载 回弹后圆柱形 件 的几何形状 。表 3列 出了卸载 回弹前后沿 轴 线上图 5b的 6个位置点 z方 向位移的计算结果 。 (a)加 载结束时 的变形 (b)卸载 回弹后 的形状 图 6 板材柱 面成形加载 及 回弹后 的形 状 表 3 柱面成形 回弹前后 方向位移 单位 ：mm l 位置点 1 2 3 4 5 6 加载结束 一 9．73 — 9．38 — 8．18 — 6．24 — 3．51 O I 回弹后 一 7．65 — 7．3O 一 6．19 — 4．89 — 2．82 O (2)板材 多点成形的回弹 多点成形的实际 零件见图 7，划分成 12× 12个等参壳单元。采用 200个增量步加载 ，并采用同样增量步数卸载 。图 7为加载过程结束时及卸载回弹后零件几何构形 的计算结果。 尽 (a)加载结 束时的变形 (b)卸载 回弹后的形状 图 7 多点成形 零件加载 及 回弹 后的形状 从图 6、图 7及表 3可以看 出板材成形卸载后 的回弹情况。 4 结束语 (1)回弹计算是板材成形数值模拟的重要内 容之一。本文给出了板材成形数值模拟的静力隐 式弹塑性有限元格式，采用平衡迭代方法求解有 限元方程，采用隐式算法计算应力积分及接触压 力 。 (2)由于静力隐式方法应力的计算精度高，对 纯弹性板材大变形问题得到了比较准确的回弹计 算结果 。 (3)板材柱面成形与多点成形回弹的数值分 析表明本文给 出的方法完全可以实际应用。 参 考文献 ： [1] Tang S C．Analysis of Springback in Sheet Form— ing Operation．In：Lange K，ed．Advanced Technolo— gy of Plasticity．Berlin：Springer— Verlag，1 987： 193～ 196 [2] Pourboghrat F，Chung K，Richmond O．A Hybrid Membrane／Shell Method for Rapid Estimation of Springhack in Anisotropic Sheet Metals．Journal of Applied Mechanics，1998，65：671～ 684 [3] Lee S W，Yang D Y．An Assessment of Numerical Parameters Influencing Springback in Explicit Fi— nite Element Analysis of Sheet Metal Forming Pro． cess．Journal of Materials Processing Technology， 1998．80— 81：6O～ 67 [4] Keck P，Wilhelm M，Lange K．Application of the Finite Element Method tO the Simulation of Sheet Forming Processes：Computation of Calculations and Experiments．International Journal Numerical M ethod in Engineering，1990，30：1415～ 1430 [5] 蔡中义．板材多点成形过程的数值模拟与最优成形 路 径 研究：[博 士学 位论 文]．长春 ：吉 林 大学， 2000 [6] Hunek I．On a Penalty Formulation for Contact— Impact Problems．Computers& Structures，1 993， 48(2)：193～ 203 [7] 蔡中义，李明哲，李广权．金属成形数值模拟 中的 接触单元法．中国机械工程，2000，l1(8)：933～ 935 (编辑 苏卫国 ) 作者简介：蔡中义，男，1963年生。吉林大学(长春市 130025)辊 锻工 艺研究所 副教授 、工学博 士。研究方 向为塑性加工 过程数 值 模拟、辊锻工艺及有限元法工程应用。获部科技进步三等奖 1项． 发表论文 3O余篇。李明哲．男．1 951年生。吉林大学辊锻工艺研 究所教授、博士研究生导师。李湘吉．男．1 975年生。吉林大学辊 锻工 艺研 究所助 教。 ·1461 · 维普资讯 http://www.cqvip.com