轧制问题的有限单元法解

D0I:10.13374/j.issnl001-053x.1983.02.033 北京钢铁学院学报 1983年第2 轧制问题的有限单元法解 力学教研室钱仁根周宝昆，乔端 摘 要 本文在平面应变条件下，假设轧件为应变硬化的刚塑性体，轧辊为不变形的刚 体，轧辊与轧件之间的接触摩擦条件为粘着，即轧辊与轧件之间无相对滑动。用刚 塑性有限单元法计算了平板轧制过程的单位压力，金属流动速废和应变、应力分布 等，并对接触弧长、刚塑性交界面、前滑系数和中性角等的确定提出新的看法。 有限单元法计算程序是以刚塑性广义变分原理为基础，采用八节点曲边四边形 等参单元。根据在四辊轧机上轧制铝板的实测数据，对计算结果进行验证。 摘 要 平板轧制问题力能参数的计算方法有切片法(Karman),一般解析法(HiI1),近 似能量法(Sieble),下、上限解法(Green),滑移线场解法(Alexander,.Fribank), 以及近几年发展起来的有限单元法（玉野敏隆，Ra0,森等）。多数的解法由于采用了一些 简化的假设，虽然在一定条件下与实验结果比较符合，但是对轧制过程金属流动的规律，塑 性变形区中应变、应力分布规律，以及刚塑性交界面等问题都难于得到令人满意的结果。 有限单元法是一种数值计算方法【！。自山田嘉昭等人【]提出了弹塑性应力一应变矩阵 的显式表达式后，弹塑性问题的有限单元法得到较大的进展。但最初，它只是求解简单的弹 塑性问题【9·」，后来才用来分析金属塑性加工过程，如锻压，挤压，拉拔，轧制等。 1973年小林史郎等人提出刚塑性问题有限单元法【s],进一步推进了用有限单元法研究塑性 加工过程。由于用有限单元法可以得到其他分析方法不易得到的一些结论，所以，虽说其 历史较短，但已显示出它的优越性，并日益受到人们的重视。关于这方面的研究工作，国外 已发表了不少论文[-]，读者可以在北美金属加工会议(NAMRC)和日本塑性加工连讲 会等专业会议论文集中看到。在国内的研究则还是刚开始【一1】。因轧制过程边界条件复 杂，所以用有限单元法研究轧制问题比较困难。玉野敏隆1！，Rao等【1s),Zienkiew- cz1]和森谦一郎等【17引在这方面都作了有益的工作。 本文假设轧件处于平面应变条件下，并为应变硬化的刚塑性体，轧辊与轧件之间为粘 着，即轧辊与轧件之间无相对滑动。用刚塑性问题的有限单元法计算了平板轧制过程的单位 压力，金属流动速度和应变，应力分布等，并对接触弧长，刚塑性交界面，前滑系数和中性 角等的确定提出新的看法。根据实测数据，对计算结果进行了验证。 本文的计算程序是以刚塑性变分原理为基础，采用八节点曲边四边形等参单元。这种单 元具有刘分单元方便，输入数据少等优点。 137

北 京 钢 铁 学 院 学 报 年第 轧制问题的有限单元法解 力 学教研 室 锐仁根 周 宝眼 乔 端 摘 要 本文在平 面 应 变条件下 ， 假设 轧件为应 变硬化 的刚塑性体 ， 轧辊为不 变形 的刚 体 ， 轧辊与轧件之 间 的接触摩擦条件为粘着 ， 即 轧辊与轧件之 间无 相对 滑动 。 用 刚 塑 性有 限单元 法 计算了平板 轧制过程 的单位压 力 ， 金 属流动速度和应 变 、 应力分布 等 ， 并对接触弧 长 、 刚塑性交界面 、 前 滑系数和 中性角等的确定 提 出新 的看 法 。 有 限单元 法 计算程序是 以 刚塑性广义 变分原理 为墓础 ， 采用 八 节点曲边 四 边 形 等参单元 。 根 据在 四 辊 轧机上 轧制铝板 的实测数 据 ， 对 计算结果进行验证 。 摘 共 平板轧制问 题 力能参数 的计 算方法 有切片法 ， 一般解析法 ， 近 似能量法 ， 下 、 上限解法 ， 滑 移线场解法 ， ， 以 及近几 年发展 起来的有 限单元法 玉野 敏隆 ， 。 ， 森 等 。 多数 的解法 由于采 用 了一些 简化的 假设 ， 虽然 在一定 条件下与实验结果 比较符合 ， 但是 对轧制过程 金属 流动的规律 ， 塑 性变形区 中应 变 、 应 力分布规律 ， 以 及刚塑性 交 界面 等问题都难 于得 到令人 满意的结果 。 有 限单元法 是 一种数值计 算方法 〔 ‘ 】 。 自山 田 嘉 昭 等人 】 提 出 了弹塑性应 力一应 变矩 阵 的显式 表达 式后 ， 弹塑性问 题 的有 限单元法 得 到较大的进 展 。 但最初 ， 它只是求解简单的弹 塑性问 题 ， ’ ‘ ， 后来 才用 来 分 析 金 属 塑性 加工 过 程 ， 如 锻压 ， 挤 压 ， 拉 拔 ， 轧 制等 。 年小林史 郎等人 提 出刚 塑 性 问题 有 限单元 法 〔 ” ， 进一步推进 了用 有 限单元法 研究塑性 加工 过程 。 由于用 有 限单元法 可 以得 到 其他分析方法 不 易 得 到 的 一些 结 论 ， 所 以 ， 虽 说其 历史较短 ， 但 巳显示 出它的优越性 ， 并 日益受 到人们 的重 视 。 关于这方面 的 研究工 作 ， 国外 已发表了不 少论文 〔卜 ， 读者可 以在北美金属加工 会议 和 日本塑性加工 连讲 会等专业会议 论文 集中看到 。 在 国 内 的研究则还是 刚开 始 卜 ’ 】 。 因轧制过 程边界条件复 杂 ， 所 以 用有 限 单元法 研究轧制问 题 比较 困难 。 玉野 敏隆 【 ’ 毛 ， 。 等 ‘ “ ， 一 ’ 和森谦一 郎等 ‘ 在这方面都作 了有益 的工 作 。 本文假设轧件 处于 平面应 变 条件下 ， 并为 应 变硬 化的刚塑性体 ， 轧辊与轧件之 间为粘 着 ， 即轧辊与 轧件之 间无 相对滑 动 。 用 刚 塑性问题 的有 限单元法计算了平板 轧制过 程 的 单位 压 力 ， 金属 流 动速度 和应 变 ， 应 力分布 等 ， 并对 接触弧长 ， 刚塑性 交 界面 ， 前滑系数 和 中性 角等的 确定提 出新 的看法 。 根据实测数据 ， 对计 算结果 进 行 了验证 。 本文的计 算程 序是 以 刚塑性变分原理 为基础 ， 采 用八节点 曲边四 边形 等参 单元 。 这种 单 元县有划分 单元 方便 ， 输 入数据少 等优点 DOI ：10．13374／j ．issn1001－053x．1983．02．033

二、基本理论 在金属塑性加工过程中，被加工物体弹性变形比起塑性变形要小得多，忽略弹性变形也 不会引起很大的误差，因而可以假设轧件为刚塑性材料。解决刚塑性问题必须满足基本方 程I1)。 根据最小位能原理，在所选择的运动许可的速度场ⅴ：中，能使泛函 ⑨=∫edv-∫sr下v:ds (1) 取最小值的Y:是问题的正确解。其中g为广义应力， 0=（11,) (2) 若采用Mises屈服条件，并设os,k分别为拉伸屈服点和剪切屈服点，则 o=0.=v 3 k, (3) e为广义应变速率， (4) 亚，为应力边界Sr上给定的表面力，v为体积。 为了便于选择速度场，可以把体积不可压缩条件，利用Lagrange乘子法，引入泛函式 (1)中，并考患到(3)式和(4)式，则泛函变为 ∮-∫7(n)声+∫.A8dv-∫ 亚：v;ds (5) ST 其中，入为Lagrange乘子，81为Kro necker记号。所有满足v:和et1之间的几何关系和 速度边界条件的速度场中，能使泛函式（5)取最小值的速度ⅴ：是正确解，这就是刚塑性 不完全广义变分原理（以后简称为刚塑性广义变分原理）。 为了获得有限单元法的最终矩阵方程，首先将刚塑性广义变分原理的泛函用矩阵表示， ∮=∫oedv+∫，A&cdv-∫，Tvds (6) 其中，、V分别为应变速率列阵和速度列阵；T为物体表面S上给定的外力列阵；C为矩 阵记号，它满足体积不可压缩条件eC=0,对平面应变何题CT=（110),矩阵右下角 的T表示矩阵转置。 关于有限单元法离散化和非线性矩阵方程的线性化步骤，在文献【：！中已有详细推导。 最终可以得到矩阵方程 [s].[8]-[R-:] (7) 其中，〔S。-1)和〔Ra-1)分别为整体刚度矩阵和等效载荷列阵，它们是由第n-1次迭代的速 度场解un-1计算到到的。(7)式则是关于△un和入n的线性方程组。给定运动许可速度场 138

二 、 基 本理 论 在金属塑性加工过 程 中 ， 被加工 物 体弹性 变形 比起塑性变形要小得多 ， 忽略弹性变形 也 不 会引起 很大 的误 差 ， 因 而可 以 假设轧件为刚 塑性材料 。 解决刚塑性问 题 必须满足 基本方 程 ’ 。 根据最小位 能原 理 ， 在所选 择 的运 动许可的速 度场 ‘ 中 ， 能使泛 函 歹 ‘ 了万百‘ 一 了厂 ‘ 一 取最小值的 ‘ 是问 题 的 正确解 。 其 中云 为广义应 力 ， 。 ‘ ，·‘ ， 专 一 百、、 、 一 一 若采 用 屈 服 条件 ， 并设 。 ， 分别为拉伸屈服 点和剪切屈服点 ， 则 二 亿 飞厂 ， 一 一 。 为广义 应 变速率 ， 二 。 二 ， 、 士 ‘ 、 巴 － ， 一 ‘ 一 、 一 ‘ 一 ‘ 犷，为应 力边界 上给定的 表面 力 ， 为体积 。 为 了便 于选择速 度场 ， 可 以把 体积 不可压缩条件 ， 利 用 乘 子法 ， 引入泛 函 式 中 ， 并考虑 到 式和 式 ， 则 泛函变为 户 ， 二 、 女 ， 少 亿 一 气 ‘ ， “ ‘ ， 入 ‘ ，各‘ ， 一 ，‘ ‘ 一 一 、 ， 一 一 其中 ， 入为 乘子 ， 各门 为 。 。 记号 。 所有满足 ‘ 和 后 ， ，之 间的 几何关系和 速度边 界条件的速度场 中 ， 能使泛 函 式 取 最 小 值的速 度 ‘ 是正确解 ， 这就是刚塑性 不 完全 广义 变分原理 以后 简称为刚 塑性广义 变分原理 。 为 了获得有 限单元法的 最终矩 阵方程 ， 首先将刚塑性广义 变分原理的泛 函用 矩 阵表示 ， 歹了万百‘ 二 了 ‘ “ ‘ 一 了 · ‘ 其 中 ， 、 分别为应 变速 率列阵和 速度列阵， 为物 体表面 上给定的外力列 阵， 为矩 阵记号 ， 它满足体积 不可压缩 条件 。 二 。 ， 对平面应 变问题 ， ， 矩 阵 右 下角 的 表示矩 阵转置 。 关于有 限 单元法 离散化和非线性矩 阵方程 的线性化步骤 ， 在文献 ’ “ 】 中已有详细推导 。 最终可 以得到矩 阵方程 「一 〕 〔 △ 〕 〔 一 〕 其 中 ， 〔 。 一 〕 和 。 一 〕分 别 为整体 刚度 矩 阵和 等效载荷列 阵 ， 它们 是 由第 一 次迭 代 的 速 度场解 卜 计算到到的 。 式 则 是 关于△ 和 入 。 的线 性 方程组 。 给定运 动许可速 度场

的初始值u,利用(7)式迭代让算，直到获到稳定的速度场。解的收歙尔件可按下式判定 I△uI/IuI<8 (8) 其中， :I△uI=V∑△u?,IuI=V∑u},8为-个 小的正数。 如文【1所述，取双二次插值函数，对八节点四边形等参单元可以在划分单元较大的情 况下，仍有足够的精度。 三、单位压力计算与实验结果比较 实验是在中110/中320×300mm四辊冷轧机上进行的1」。在工作辊上安装四组(8根) 倾斜测压针，测定轧制铝板时沿接触弧分布的单位压力。由压上螺丝处的压头测定总轧制 力。本文计算结果是与位于板宽中央附近一对测压针的实测数据进行比较。 工作辊是直径D=110mm,材料为GCr15,未经淬火的软面辊。轧辊转速是35r·p·m。 轧件材料为工业纯铝，通过单向拉伸试验可求得实验曲线及其回归方程【2]。 根据单一曲线假设，轧件的广义应力一广义应变关系为 0=3.2310+0.9048e06530 (9) 轧件宽为B=50mm,轧件原始厚度H=2.90mm,本文计算了压下率分别为13.07%、 20.77%、29.79%三种情形。 单元的划分如图1所示。由于下上对称，只取下半部进行计算。考虑到在靠近轧辊表面 的速度有较为剧烈的变化，因此，沿Y方向上、下分为高度比为1·2的二层单元，即靠近 轧辊的一层较为扁平一些。在入口和出口截面两外侧的单元，其X方向长度为接触弧长L的 一半。由计算实践可知，在被划分单元的区域以外的范围，金属接近刚性运动。 H/2 图1。单元划分 图2()、(b)、(c)为轧件原始厚度不变，改变压下率时，用有限单元法计算的单位压 力分布与实验【2！结果的比较。表1则是计算得到的平均单位压力”和总轧制力P与实验值 的比较。 表1 有限单元法计算结果与实测值比较 板 轧前 轧后 压下 接触弧长 平均单位压力 径 厚度 厚度 率e L(单m) p(kg/mm 2) 总轧制力P(T) B 6 (mm)(mm)(mm)(mm)(% 计算实测计算 实测 计算 实测 误差 110 50 2.87 2.49513.07 4.5625.3012.95 9.77 2.954 2.59 12% 110 50 2.84 2.25 20.775.7226.94 15.629.80 4.468 3.40 24% 110 50 2.922.05 29.796.9358.46 19.6315.60 6.806 6.60 3% 139

的 初始值 。 ， 利用 式迭 代计 算 ， 朋△ 直 到获到稳定 的速 度场 。 乙 训 名 △ 了 解的 收款条件可 按下 式判 定 其 中 ， 皿△ 小的 正数 。 如文 ‘ ” 所述 ， 取双 二 次插值 函数 ， 况 下 ， 仍有 足够 的精度 召 公 矛 ， 乙 为一个 “朋 声 对八 节点四 边形 等参单元可 以 在划 分 单元 较大 的情 三 、 单位 压 力计算 与实验 结果 比较 实验是 在 中 小 幻。 四 辊冷轧机 上进 行的 ‘ ’ 。 在工 作辊上安装 四组 根 倾斜测压针 ， 测 定 轧 制 铝 板 时沿接触弧分布的 单位压力 。 由压 上螺丝处的压 头测定总轧制 力 。 本 文计 算结果 是 与 ， 位于板宽中央 附近 一对测压针的实测数据 进 行 比较 。 工 作辊是直径 ， 材料为 ， 未经 淬火 的软面 辊 。 轧辊转 速是 · ， 。 轧 件材料为工 业纯 铝 ， 通过 单向拉伸试验可求得实验 曲线 及其 回归方程 〔 ” 。 根 据 单一 曲线假设 ， 轧件的广义应 力一广义 应 变关系 为 于 百 · “ “ “ 轧 件宽为 ， 轧件原 始厚度 二 ， 本文计 算了压 下 率 分 别 为 、 、 三种 情形 。 单元 的 划 分如 图 所示 。 由于下 上对称 ， 只取 下半部进 行计 算 。 考虑 到在靠近 轧 辊表面 的速度有较为剧 烈 的变 化 ， 因此 ， 沿 方向上 、 下分为高度比为 的二层 单元 ， 即靠近 轧辊的 一层 较为扁平 一些 。 在 入 口 和 出 日截面 两外 侧的 单元 ， 其 方 向长 度为接 触 弧 长 的 一半 。 由计 算实践可 知 ， 在被划分单元 的 区域 以外 的范 围 ， 金属接近 刚性运 动 。 ‘ ’ 卜 协 图 、 力分布 与实验 的比较 。 厅 表 图 单元 划分 、 为轧件原始厚度不 变 ， 改 变压下 率 时 ， 用有限 单元 法计 算的单 位压 “ ‘ 结果 的 比较 。 丧 则是计 算得到 的 平 均单 位压 力， 和 总轧 制 力 与实验值 有 限单元法计 算结果 与实测值 比较 径 轧后 厚度 压 下 率 。 总轧制力 ‘ 吸一‘ ﹄备 … 板宽 辑凰前度 接触弧长 平均单位压 力 生粤粤 ， 巫且边塑虹 计算 实测 计算 …实 测 实，。 误 差 月城︸了一 行八内才 勺﹄ 计 算 户卜 ︸﹄ 、

p(kg/mm) p(kg/mm) H=2.87mm Hs2.84m血 h=2.495mm 20 h=2.25mm 20 c=13.07% g=20.77% 0 10 计算曲线 计算曲线 一●一实测曲线 一。一一实测曲线 87 5 98 4 L(mm) L(mm) (a) (b) 单位压力分布曲线，除了入口截 p(kg/mm+) 30 面附近以外，其他地方一般与实验结 H=2.92mm h=2.05mm 果符合程度尚好。至于入口截面附近 t=29.79% 的压力分布现象，在其他文献【22·2！ 和森谦一郎最近的工作2」中，也有 同样的结果。 计算得到的单位压力高于实验 值，其原因是两方面的。一是实验本 20 身的误差，二是计算时按照粘着摩擦 条件，导致单位压力偏高。由于实测 的接触弧长与计算值之差较大，因而 平均单位压力p较之总轧制力P有更 大的误差。 在轧制问题的分析方法中，通常 10 总是人为地规定塑性变形区。但是用 有限单元法进行分析计算时可以避免 这种不真实的假设。用有限单元法计 计算曲线， 算时，把轧件入口截面后和出口截面 实测曲线： 前取一定的长度。由计算得知，考虑了 入口和出口截面两外侧区域，增加了 单元的数目，因而增加了计算时间， 109 654321 但能较真实反映轧制过程的应力、应 1(mm) (c) 140 图2单位压力分布

盆 之 ‘ 》 二 五二 七 二 留 也 口 卜 ‘ 二 盯 ‘ 一 计算曲线 － 。 －实测曲线 一 ‘一 计算曲线 一 一 实侧曲跳 口 单位 压 力分布 曲线 ， 除 了入 口 截 面 附近 以外 ， 其 他地 方一般 与实验结 果 符 合程 度尚好 。 至 于 入 口 截面 附近 的压力分布现象 ， 在 其他文献 【 ‘ “ 和 森谦一 郎最近 的 工 作 ‘ 中 ， 也有 同样的结 果 。 计 算得到 的单 位 压 力 高于实验 值 ， 其原 因是 两 方面 的 。 一是实验 本 身的误 差 ， 二 是计 算时按照粘 着摩擦 条件 ， 导 致单位压 力偏高 。 由于实测 的接触 弧长 与计 算值之 差较大 ， 因而 平均单位压 力乡较 之 总 轧制 力 有更 大的误 差 。 在轧制问 题 的分析方法 中 ， 通 常 总是人为地 规定塑性变形 区 。 但是用 有 限单元法 进 行分析计 算时可 以避免 这种不真实的假设 。 用有限单元法计 算时 ， 把轧件 入 口 截面 后 和出 口 截面 前取一定 的长度 。 由计算得知 ， 考虑 了 入 口 和 出 口 截面 两外侧区域 ， 增加 了 单元 的 数 目 ， 因而增加 了计 算时间 ， 但 能较 真 实反 映轧 制过程 的应 力 、 应 二 。 二 石 ‘ 二 您 万 ， 计算曲舰 实橄曲线 图 单位压 力分布

变状态。如果只考虑入口和出口截面之问的区域，则与真实状态有较大的误差。本文对这两 种区域划分都作过计算，计算结果如炎2所示。 表2 两种区域计算轧制力比较 辊 板 轧前 轧后 总轧制力P(T) 厚度 厚度 D 不考虑 考 虑 mm)(mm) mm) (mm) 外侧区域 外侧区域 误差 110 50 2.87 2.495 2.298 2.954 22.2% 110 50 2.84 2.250 3.580 4.468 19.9% 根据本文计算结果， 取两外侧区域长度为L/2时，两端面已经十分接近为刚体运动，因 ·此计算结果较好。 四、其他计算结果和讨论 用有限单元法计算轧制问题，除了能求得单位压力外，还能求得诸如轧件各点的流动速 度，应变、应力分布等结果。根据速度场和应变等值线可以确定刚塑性交界面，前滑系数等 以及对接触弧长问题提出新的看法。· 1.速度分布 通过计算可以得到轧件内各点的速度分布。图3(a)为轧件原始厚度H=2.84mm,轧 后厚度h=2.25mm时速度分布图。图中表示速度矢量的箭头在Y方向被放大了，图3(h) 为与X轴垂直的11个截面上速度分量ⅴx沿轧件厚度上的分布。图3(b)中横座标轴表示速 度vx(m/sec),纵座标轴表示轧件厚度方向。 Y轴 0.2m/Sec 2 3-34-45.-56-67-78-89-91-1011=11 X轴 0.25 色 0.20 0.15 010 0.05 X轴 图3(a)轧件各点速度 141

变状态 。 如 果 只 考虑 入 口 和 出 「截面 之 间的 区域 ， 则 与真实状态 有较大 的误 差 。 卞文 对这两 种 区域 划 分都 作过计 算 ， 计 算结果 如 丧 所 示 。 表 两 种 区域计 算轧 制 力比较 轧厚后度 总轧 制 力 虑 区 域 误 差 径辊 板宽 辐凰前度 ︸︹︸ 匀口产 ” 根据 本文计 算结果 ， 取 两外 侧 区域 长 度为 时 ， 两 端而 已经 十分 接近 为刚体运 动 ， 因 一 此计算结 果较好 。 四 、 其 他 计算结果 和 讨论 用有 限单元法计 算轧 制问 题 ， 除 了能求得 单位压 力外 ， 还能求得诸如 轧件各点的 流动速 度 ， 应 变 、 应 力分布 等结果 。 根据 速 度场和应 变 等值线可 以 确定 刚塑性 交 界面 ， 前滑系数等 以 及对接触弧长问 题 提 出新的 看法 。 速度分布 通过计 算可 以得到轧件内各点的速度分布 。 图 为轧件原始厚度 ， 轧 后 厚度 时速度分布图 。 图 中表示速度矢量 的 箭头 在 方 向被 放大 了 ， 图 为与 轴垂直 的 个截面 上 速 度 分 量 二 沿轧件厚度 上的 分布 。 图 中横座 标 轴 表示速 度 ， 纵座标 轴 表示轧件厚度方 向 。 “ 、 、 ‘ ，轴 一 一 二 曦一 一 一 · 了 一 一 一 ， 轴 弓 么组 址绝 ， 城的功拐、月目口 抽 图 乳件各 点速度