D0I:10.13374/j.issn1001-053x.1982.02.037 北京钢铁学院学报 1982年第2期 金属塑性压缩过程三维流动的分析 力学教研室乔端林拥 摘 要 本文根据刚塑性体的广义变分原理采用20节点60自由度曲六面体等参单元，计 算了长、宽、高不等的六面体塑压过程的流动速度场、应变场、外力等。比较全面 准确地反映了金属塑压过程三维流动的情况，并将计算结果与实验进行了比较。 一、前 金属塑性加工过程所接触到的问题绝大多数是三维变形问题。严格地说，属于完金平面 应变的问题是极少的。如果说对于平面应变条件下理想塑性体的流动问题还可以用滑移线法 进行分析的话。那么，对于一般空间的三维流动问题的精确计算则只有在有限单元法发展的 今天才有可能。· 用有限单元法求解空间问题与解平面问题相比，由于单元和节点数目的急剧增加，而且 形成单元刚度矩阵所需的时间比平面问题也长得多，特别是对于变形量较大的金属塑性加工 过程这类非线性问题需要选代求解，其计算量的增加往往关系到是否可以用来实际求解。所 以，除要求计算机贮存和运算速度增加以外，·在选择计算模型、单元类型和计算技术上也必 需有更进一步的考虑。 用有限单元法对金属塑性加工过程三维流动进行分析，迄今有长松昭男〔1)〔2〕在1973年 采用弹塑性有限元法计算了正方形断面六面体的压缩，得到塑性区扩展的情况，以及开始屈 服时的变形、外力等。由于计算时间的限制，只算到压下率为万分之几便停止了。关于锻压 过程宽展的研究，1977年北原義之等〔3)将高度方向的变形视为均匀，仅考虑x、y方向的宽 展，用有限单元法计算得出若干结果。 实际上，对锻压技术而言，除希望知道被加工件中部的“宽展”外，还希望能知道侧面 变形的形状，也就是说希望获知变形区内各点的三维流动情况。 在作者等的早先一篇文章〔4)中已经指出：用刚塑性有限单元法计算金属塑性加工过 程，每次压下量可以比较大，因而是一种有力的计算工具。 对金属塑性加工过程，由于材料的塑性变形比较大，略去弹性变形部分，而视为钢塑性 体，其准确度是足够的。本文根据刚塑性体广义变分原理，采用20节点曲六面体等参元、计 ◆本文食在中国机械工程学会领压半会第二届塑性加工理论学术会议上宣读(1981年10 月) 132
北 京 钥 铁 学 院 学 报 年第 期 金属塑性压缩过程三维流动的分析 ’ 力学教研 室 养端 林拥 摘 要 本文根据 刚塑性体的广义 变分原理采用 节点 自由度曲六面体等参单元 ， 计 算了长 、 宽 、 高木等的六面体塑压过程 的 流动速度场 、 应变场 、 外力等 。 比较全面 准确地反映 了金属塑压过程三维 流动的情 况 ， 并将计算结果 与实验进行了比较 。 月呀 污 金属 塑性加工 过程所接触到的问题绝大多数是三维 变形 问题 。 严 格地说 ， 属于 完全平面 应 变的 问题是极少 的 。 如果 说对于 平面应 变条件下理想塑性体的 流动 间题还可 以用 滑移 线 法 进行分析的话 。 那 么 ， 对于一般 空 间的兰雏 流动问题的精确计算则只 有在有限 单元 法发展 的 今天才有可能 。 · 、 六 用有限单元 法求解空间问题与解平面问题相 比 ， 由于单元和 节点数 目的急剧 增加 ， 而且 形成单元 刚度矩 阵所需的时间比 平面 问题也长得多 ， 特别是 对于 变形量 较大的金 属 塑性加工 过程这类非线性问题需要选代求解 ， 其计算量 的增加 往往关系到是否可 以用 来实际求解 。 所 以 ， 除要求计算机贮存和 运算逮度增加以外 ， · 在选择计算模 型 、 单元 类型和 计算技术 上也必 需有更进一步 的考虑 。 用 有限单元 法对金属塑性加工 过程三维 流动进行分析 ， 迄今有长松 昭 男 〔 〕 〔幻 在 年 采用弹塑性有 限元法计算宁芷方形断面六面体的压缩 ， 得到塑性 区扩展 的情况 ， 以 及开 始屈 服时的变形 、 外力等 。 由子计算时间的限制 ， 只算到压下率为万分之 几便停止 了 。 关于锻压 过程宽 展 的研究 ， 年北原羲之等 〔 〕将高度方向的变形视为均匀 ， 仅考虑 、 方 向的宽 展 ， 用 有限单元 法计算得出若干结果 。 实际 上 ， 对 锻压技术而言 ， 除希望知道被加工件中部的 “ 宽展” 外 ， 还希望能知道侧面 变形的形状 ， 也就是说希望获知变形区 内各点 的三维 流动情况 。 在作者等的早先一篇文章 〔 〕 中已经指 出 用 刚 塑性有限 单元 法计算金属 ’ 塑性加工过 程 ， 每次压下量可以比较大 ， 因而是一种有力的计算工具 。 对金属 塑性加工过程 ， 由宇材料的塑性变形比较大 ， 略去弹性变形部分 ， 而视为钢 塑性 体 ， 其准确度是 足够的 。 本文根据刚塑性体广义 变分原理 ， 采用 节点 曲六面体等参元 、 计 今 水文 金在 中国 机械工程 学会镶 压 学会第二 届 塑性加 工 理论 学术 会议 上 宣 读 年 月 DOI ：10．13374／j ．issn1001－053x．1982．02．037
壶￥ 算了长、宽、高不等的六面体的压缩，压缩率达12.35%。继续的压缩实际上已无任何困 难。 采用二十节点曲六面体等参单元具有输入数据少，分单元容易，而且能很好地跟踪变形 过程单元畸变的形状等优点，更主要的在于具有较高的插值次数，可以使划分的单元和节点 数比较少而保证有较高的精度，这对于减少空间问题的计算量是具有关键意义的。 本文计算了具有应变硬化工业纯铝的矩形块体在压缩过程各个时刻变形区内各点的三维 流动速度场，并由此得出变形前矩形网格在变形过程各个时刻的畸变图。计算了单元各个积 分点的等效应变，可以看出变形区各部分变形不均匀程度的图象。计算了随着变形的增加， 工具作用在变形物体上的力，并与实验作了比较，取得较好的结果。 二、张本原理和计算方法 计算原理在参考文献〔4)〔5)中已有详细说明，这里简要叙述如下： 根据刚塑性体的广义变分原理，在忽略体积力的情况下速度场的正确解使得泛函 Φ=∫，gedv-∫s,TrV ds+-∫iTcdv (1) 取得极小。其中第一项积分表示物体的变形能，第二项积分是代表外力作的功，第三项表示 体积不可压缩条件。 将连续体离散后，（1)可以写成节点速度U-1和关于△U和入()的线性方程组 [a-81[]-[8w]-61 (2) 利用（2)式重复迭代计算，直到收敛。这时的速 度场即认为是速度场的正确解。 ⑧ 采用二十点曲六面体等三单元，单元为60个自由 图 度。在局部坐标(，n,)下，考察一中心在原点， 边长为2的立方体单元，将8个顶点及二十条棱边的中 点都取为节点，共有二十个节点，如图1所示。 取插值函数为如下形式 U=a,+a2ξ+aan+a45+agξ2+agn +a,52+aξn+an5+a1055+a1:52n ③ +a12525+a1sn2ξ+a14n25+a1552号 ① 迈 +a1.52n+a17ξn5+a18ξ2n5 +a10n25G+a2o℃25n (3) 图1 20节点曲六面体单元 可以知道，此插值函数可以保证在局部坐标下满足相容性条件。 记节点的函数值为U:,可将插值函数写成 20 U=∑Ni(5,n,U, (4) 1=1 :其中Ni(5,n,)(i=1,2,…20)为相应的形状函数，其丧达式为 138
安釜， 算了长 、 ’ 宽 、 高不等的六面体的压缩 ， 压缩率达 。 继续的压缩实际上 已无任何困 难 。 采 用二十节点 曲六面体等参单元 具有输入数据 少 ， 分单元 容易 ， 而且能很好地跟踪 变形 过程单元 畸变的形状等优点 ， 更主要 的在于具有较高的插值 次数 ， 可以使划分 的单元 和节点 数比较少而保证有较高的精度 ， 这对于减少 空 间问题 的计算 是具有关键意义 的 。 本文计算了具有应 变硬 化工 业 纯铝的矩形块体在压缩过程 各个时刻变形 区 内各点 的三维 流动速度场 ， 并由此得出变形前矩 形网 格在 变形过程 各个时刻的畸变图 。 计算了单元 各个积 分点 的等效应 变 ， 可 以 看出变形 区 各部 分变 形不 均匀程度的图 象 。 计算了随着变形的增加 ， 工 具作用 在变形 物体上的力 ， 并 与实验作 了比较 ， 取得较好的结果 。 毖本 原 理 和 计算方法 计算原 理在参考文 献 〔 〕 〔 〕中 已有详细说 明 ， 这里简要叙 述 如下 根 据 刚塑性体的广义 变分原 理 ， 在 忽略体积力 的情况下速度场 的正 确解使得泛 函 ， 二 孟 一 ’ ‘ 丁 ‘ 云 ’ ‘ · 取得极小 。 其 中第一项积 分表示 物体的变形能 ， 第二项 积 分是代表外 力作 的功 ， 第三项表示 体积不可 压缩条件 。 一 ‘ 将连续体 离散后 ， 可 以 写成 节点逮度 卜 和关于 △ 和 入《 。 的线性方程组 「 卜 〕 一 △卫卫 〕 「 卜 〕 厂 〕 ‘ 一 入 卜 利用 式 重 复迭代计算 ， 直 到收敛 。 这时 的速 度场 即认为是速度场 的正 确解 。 采用二十点 曲六面体等三单元 ， 单元为 个 自由 度 。 在局 部坐标 毛 ， ，， 切 下 ， 考察一 中心 在原点 ， 边长为 的立 方体单元 ， 将 个顶 点 及二十条棱边 的 中 点都取为节点 ， 共有二十个节点 ， 如图 所示 。 取插 值 函数为如下形式 色 月 ‘ 之 。 友 “ 。 月犷 屯 “ 七” 邑 。 息七 七 “ ” 息 “ 七 艺 邑 ‘ “ 屯 ， 乙 “ 邑 乙 名 ， 息 套 。 邑 “ 乙 ， 。 “ 息乙 。 “ 七” 岁性 图 节点 曲六 面体单元 可 以 知道 ， 此插值 函数可 以 保证 在局 部坐标下 满足 相 容性条件 。 一 记节点 的 函数值为 ‘ ， 可将插 值 函数 写成 乏 ‘ ， 。 ， 。 ‘ 二 ， … … 为相应的形状 函数 ， 其表 达式为 书 歇” 中 “ “ ， ” ， “ ’ “ 冬笋乒
音(1+ξ：)(1+n:n)(1+5:t)(E,5+n:n+℃：C-2) (i=1,3,5,7,13,15,17,19) +(1-*)(1+nrn)(1+C:C) N,(5,n,t)=.(i=2,6,14,18) +(1=n2)(1+5:)(1+t:) (i=4,8,16,20) +(1-t2)(1+5:)(1+n:n) (i=9,10,11,12) (5) 二十节点的空间等参单元在整体坐标下，其棱边是一个二次曲线，而其侧面则是由两族 这样的二次曲线交织而成的曲面。正由于此，对物体弯曲的边界，二十节点曲六面体单元在 边界附近能有较好的近似。 单元的应变速率为 0 0 0x 8 Ey 0 ay 0 0 8 Ez 0 8z t= 1 18 18 V2Y √28y√28x 0 (7) 1 18 V2 Yr: 0 18 √28z 8y 1 18 √20z 0 √28x 将插值公式(4)代入(7)式中得 =(B). (8) 其中 8 8x 0 0 0 6 0 0 0 6 N:00N200…N2000 8z B= 18 18 0 0N,00N20…0N2.0 √28yV2x t00N:00N2…00N2o 0 18 1 √ 8z √2 ay 10 8 V2 8z 0 1 √/2 0x (9) 每个单元的Q,F,P,H,采用十四点求积公式进行计算，这个积分公式为 rc,n cdtdnde 14
‘ ， 。 ， 气 舌 七 ‘ 七 介‘ ， 仁‘ 仁 毛 ， 七 ” ‘ 仁 ‘ 仁一 ， ， ， ， 一 ， ， 十 一 七 ， ， 仁 ‘ ‘ ， ， ， 去 ， 七 ‘ 七 七 ‘ 仁 ， ， ， 十 一 ‘ 七 ‘ 七 ‘ ， ， ， 二十节点 的空 间等参单元在整体坐标下 ， 其棱边是一个二次 曲线 ， 而 其侧面 则是 由两 族 这样的二次曲线交 织而成 的 曲面 。 正 由于此 ， 对物体弯曲的边界 ， 二十节点 曲六面体单元 在 边界附近能有较好 的近似 。 单元 的应变速率为 日 口 亿一 口 日 、 、 ‘ 、 、 ‘ … 卫 一立旦 砂一， 卜 ︸ 亿 一口目一抢在 亿了丫 ” 亿了丫” 亿 畜 识了 将插值公式 代入 式 中得 了 少、 、沪 ‘」卫 一 其中 一 ， 一叙。 。 侧 口 侧 一 召 百 侧 识 … … … …一 一 一 一 记 每个单元 的 ， ， ， ， 采用 十 四点求积公式进 行计算 ， 这个积分公式为 丈 ， 丈丈 “ “ ， ” ， “ ，‘ “ ” ‘ 冬笋洛
6 121 361 ∑fg,, E)+ 320 361 fE,, E.) (10) i=1 i=1 其中积分点的坐标分别是 （-,-a,-a) （， -a,-a) (-u, 2,-a) ,m,E)= (a, 以， (11) (-a, -4, a) (以，-， a) (-", ) （, a) 玫 (0,0,-B) (0,0, ) ,, (0,-B,0) )= (0, B,0) (12) (-B,0,0) (B,0,0) 61 图2单元上积分点的分布 而 -√ (13) 十四个积分点的分布如图2所示。 三、计算实例及与实验结果的比较 本文计算了长×宽×高=80×48×40毫米工业纯铝六面体的压缩，假设工具和被加工件 表面为粘着的摩擦条件，即相对无滑移。由实验〔8)回归出这种材料单向拉伸应力应变曲线 为： Y=3.23+11.60e0.58 (14) 考虑到对称性，取块体的1/8为计算对象（图3(）)。这样，将该部分分成图3(b) 所示的8个单元。 .135
带 艺 乙 礼 乙 纂 矛 曹 ‘ ， 育 ‘ ， 臼 其 中积 分点 的坐标分别是 一 仪 ， 〔又， 一 尤， ， 一 ， 茗尤 ， 一 口 ， 工， 一 以 ， 一 仪 ， ，， 了、矛 内下阳了 一 ， 一 足 一 一 一 一 以 ︸氛 ︸爪 ︸氨 一 誉 ‘ 育 ‘ ， 乙 ， ， 一 ， ， 日 ， 一 日 ， ， 日 ， 一 日 ， ， 日 ， ， “ …龙弓犷 ， 习 广 少 匕 … 口 图 单元上 积分 点的分布 而 ” 仅 二二 － 日蕊 所示 。 ‘ 」旦 十四个积 分点 的分布如图 三 、 计算实例及 与实验结果 的 比较 本文计算了长 宽 高 毫米工业 纯铝六面体的压缩 ， 假设工具和被加工件 表面 为粘着 的摩擦 条件 ， 即相对无 滑移 。 由实验〔 〕 回归 出这种 材料单向拉伸应 力应变曲线 甘吸七，‘ 、 口舀几 为 考 虑到对 称性 ， 取 块体的 所示的 个单元 。 · ‘ 为计算对 象 图 。 这 样 ， 将该部 分 分成图
⑧ ⑧ 2 ④ ① 80 (a) (b) 图3单元的划分 初始速度场可以为均匀场，也可以由弹性解给出。本文的作法是先将图3所示的整个变 形区分成2个单元，以均匀场作为初始速度场进行计算，求得收敛的速度场后，再将变形区 细分成8个单元，中间节点的速度值由初次计算的收敛速度场内插得出。这样做可以方便地 找到接近于真实速度场的初始速度场值，节省计算时间。图4分别用2个单元和8个单元计 算长×宽×高为10×10×20毫米块体压缩时在z=0截面上，边界点速度场值。由图中可以看 到，仅用两个单元计算时，得到的收敛速度场已与 Y 真实的速度场相当接近。因此用它通过内插得到更 细化单元的初速度场是合适的。同时，由两种单元 分划法计算出来结果的相近程度也说明，在本问题 中仅将物体分为8个单元计算的精度是足够的。 计算是在M150计算机上进行的，形成单元刚 。2单元计算结采 度矩阵约需30秒，每次选代时间约器4分40秒。开 义8单元计算结果 始时每次压下率取2/1000，几次后每次压下率给 (u值扩大100倍) 1/100左右。每增加一次压下，选代1一3次即可以 收敛。计算出收敛速度的场后，将各点的位移值迭 加在节点坐标上，即可得出坐标网格的畸变图。图 0 10 X 5为当高度方向压缩率为12.35%时变形区各个截 图41/2高裁面上用不同单元 面坐标网格的畸变图。 的计算结果 9 24 ⑨ 112 20 (a)z=0截面 136
双习 一⑥丫二⑧ 口 丫刁 丰 面 兰口夕「感 口 口 ② 冈口厂闻 凡 一口 诊 刁 图 单元 的划分 初始速度场可 以为均 匀场 ， 也可 以 由弹性解给 出 。 本文的作法是先将图 所示 的整个变 形 区分成 个单元 ， 以 均匀场作为初始速度场进行计算 ， 求得 收敛的速度场后 ， 再将变形 区 细分成 个单元 ， 中间节点 的速度值 由初次计算的收敛速度场 内插得 出 。 这样做可 以方便地 找到接近于真实速度场 的初 始速度场值 ， 节省计算时间 。 图 分别用 个单元和 个单元计 算长 宽 高为 。 毫米 合块体压缩时 在 。截面 上 ， 边界点速度场值 。 由图 中可 以 看 到 ， 仅用 两个单元计算时 ， 得到的收敛速度场 巳与 真实的速度场相 当接近 。 因此用 它通 过 内插得到更 细化单元 的初速度场 是合适的 。 同时 ， 由两种单元 分划法计算出来结果 的相近程度也说 明 ， 在本问题 中仅将物体分为 个单元计算的精度是足够 的 。 计算是在 计算机上进行的 ， 形成单元 刚 度矩 阵约需 秒 ， 每次迭代时 间 约需 分 秒 。 开 始时每次压下率 取 ， 几 次后每 次压下率给 左右 。 每增加一次压下 ， 迭代 一 次即可 以 收 敛 。 计算出收敛速度的场后 ， 将各点的位移值迭 加在节点坐标上 ， 即可得出坐标网 格的畸变图 。 图 为当高度方向压缩串为 时变形区各个截 面 坐标网 格的畸变图 。 卜 一 一 一 一 一 一 〕 单元计算结果 单元计算结 果 。 值扩大 倍 图 高截面上 用不 同单元 的计算结果 ’ 寻 ， 夕犷 一 一 ‘ 少 ‘ 下 臼 口 月卜 ， 一 乙一一一 户 醉 一 二二卫 ‘ 曰 口 一洲亡月尸 一 一 日成产山 乙一一一一 “ 巴 一 一 ‘ 灿 了 ‘ 卜洲尸 截 面 争尽