D0I:10.13374/j.issn1001-053x.1999.05.038 第21卷第5期 北京科技大学学报 Vol.21 No.5 1999年10月 Journal of University of Science and Technology Beijing 0ct.1999 非圆截面管辊弯成形过程的变形 和应力场的FEM模拟 徐树成”王先进) 刘才) 1)北京科技大学材料科学与工程学院,北京1000832)燕山大学,秦皇岛066004 摘要采用三维大变形弹塑性有限元法,考虑材料和几何双重非线性,基于Prandtl-Reuss流 动规律和Mss屈服准则,对圆管的辊弯二次成形过程中金属的流动规律及其应力分布进行了 模拟分析,分别得到非圆截面管辊弯成形过程中金属的流动规律及应力的信息, 关键词辊弯:弹塑性有限元:金属流动:应力分布:变形 分类号TG335.7 近年来,有限元法开始应用于金属压力加 元的刚度迭加起来,可得到整个变形体中各节 工中,早期只是日本学者木内学做了关于这种 点的位移,继而得到了应力应变等信息. 成形方式的实验研究,如成形力的测定等.近 本文采用U.L大变形弹塑性有限元法来模 些年,有限元法也开始应用到辊弯成形过程.早 拟非圆管二次成形从非稳态到稳态的全过程: 在1987年以小野田义富等为代表的日本学者 在非圆管的有限元模拟过程中,采用了最小二 就已经开始这方面的有限元法研究,但他们仅 乘法将孔型多段圆弧进行回归,使工件上与孔 采用了刚塑性有限元法,不计变形过程中的弹 型圆弧接触处尖点的接触点的归属问题得到了 性部分,当弹性部分较大时会出现大的误差.虽 解决. 然这种处理使问题得到了简化,然而对接触区 按Prandtl-.Reuss弹塑性增量理论和Mises 以外的变形,以及出口外的弹性恢复等一系列 屈服准则得到弹塑性本构关系方程: 问题不能给出正确的解, △a=DR△E克 (1) 本文采用三维大变形弹塑性有限元法,考 其中:D为弹性本构矩阵,其表达式为: 虑了材料和几何的双重非线性,因此允许用较 30a'0u (2 大的应变增量进行计算,从而大大地缩短了弹 D呢F 261+2(1+DH 3 塑性有限元法求解的全过程,并用来描绘了圆 式中:D品为与小变形情况相同的弹塑性本构矩 管二次成形时质点的流动规律和应力分布规 律,这对合理选择和设计设备、控制产品的质量 阵,成为Konkr系:数C2为刚性模量: 都有着重要的意义. ō为等效应力:H”为等效应力对于塑性应变曲 线的斜率. 1理论基础 按虚功原理每个单元的外力功等于内力 功,即δW-δW=0可得: 有限元法是通过计算变形体内有限个节点 Knm=K码+Kn+K (3) 的位移而求得整个变形体的解的一种方法.首 其中:K=BDuBundv是与小变形弹塑性有 先将变形体划分成若干个有限单元,并由节点 限元法刚度矩阵相同的变形矩阵,反映了 连接在一起.每一个单元中节点上的力与位移 由线变形引起的刚度;=(NiOON- 之间的关系可以通过应力和应变之间的关系和 2B,doBn)dy是对物体大变形时刚性转动所引 应变与位移的关系及虚功原理求得,把所有单 起的应力变化部分进行校正的应力校正矩阵, 1999-03-14收稿徐树成男,35岁,讲师,博士后
第 21 卷 第 5期 1 9 9 9 年 10 月 北 京 科 技 大 学 学 报 J o u r n a l o f U n vi e r s iyt o f S e i e n e e a n d Te e h n o l o yg B e ij i n g V O I . 2 1 N 0 . 5 o e t . 1 9 9 9 非 圆截面 管辊 弯成 形 过 程 的 变形 和 应 力场 的 F E M 模拟 徐树成 ” 王 先进 ` , 刘 才 2 , 1) 北京 科技 大学 材料科 学 与工程 学院 , 北京 10 0 0 8 3 2 ) 燕 山 大学 , 秦 皇岛 0 66 0 0 4 摘 要 采用 三维大 变形 弹 塑性有 限元 法 , 考 虑材 料和 几何 双重 非线 性 , 基 于 rP an dlt 一 R e us s 流 动规 律和 M is es 屈服准 则 , 对 圆管 的辊弯 二次成 形过 程 中金 属 的流动 规 律及其 应力分 布进行 了 模拟分 析 , 分 别得 到非 圆截面 管辊 弯成 形过程 中金 属 的流动 规律 及 应力 的信息 . 关 键词 辊 弯 : 弹 塑性 有限元 ; 金属 流动 ; 应 力分 布 : 变 形 分类 号 T G 3 35 .7 近 年来 , 有 限元法 开始应用 于 金 属压 力加 工中 . 早 期只 是 日本 学者木 内学做 了关于 这种 成形方式 的实验研究 , 如成形 力 的测 定等 〔1 . 近 些年 , 有限 元法也 开 始应用 到辊弯成 形 过程 . 早 在 19 8 7 年 以小野 田 义 富等 LZI 为代表 的 日本 学者 就 已经开 始这 方面 的有 限元 法研 究 , 但他 们仅 采 用 了刚塑性有 限元法 , 不 计变形过程 中的弹 性部分 , 当弹性部分较大时会出现 大的误差 . 虽 然这 种处 理 使 问题得到 了简化 , 然而 对接 触区 以外 的变 形 , 以及 出 口 外 的弹性恢复等 一系列 问题 不 能 给 出 正确 的 解 . 本文 采用三 维 大 变形弹塑 性 有 限元法 , 考 虑 了 材料和 几何 的 双重非线性 , 因此 允许用较 大的 应变增量进 行计算 , 从而大大地 缩短 了弹 塑性 有限元法求 解的全过程 , 并用 来描绘 了圆 管二 次成 形 时质 点 的流动 规 律 和 应 力分 布 规 律 , 这对合理选择和 设计设备 、 控制产品 的质量 都有着重 要的意义 . 元 的刚度迭 加起来 , 可 得到整 个变 形体 中各节 点 的位移 , 继 而得到 了 应力应变等信 息 . 本文采用 .U L 大变形弹 塑性有 限元法来模 拟非 圆 管二 次成形 从 非 稳态 到稳态 的 全过 程 ; 在非 圆 管的有 限 元模拟过程 中 , 采用 了最 小二 乘法将 孔型 多段 圆弧进行 回归 , 使工 件上 与孔 型圆弧接触 处尖 点的接触 点的归属 问题得到 了 解 决 . 按 P arn id l 一 R e us s 弹塑性增量 理论和 M i se s 屈 服准则得 到弹塑 性本构关 系方程 : △端刁私战 ( 1) 其 中 : 刀界 才为弹 性本构矩 阵 , 其 表达式 为 : 刀粉佘… 。 ` 。 。、 (击) 万二了二粤奇溉画{ 、2 ) L ` U L ` , 一二厂一 ) J 1 理论 基础 有限元法是通过计算变形体 内有 限个节 点 的位移而求得整个变 形体 的解的 一种方法 . 首 先将变形 体划 分成 若干个有限单元 , 并 由节 点 连接在一起 . 每一 个单元 中节 点上的力与位移 之间的关系可 以通过应力和应变之间的关系和 应变与位移的关系及虚功 原理求得 . 把所有 单 19 9 一 03 一 14 收稿 徐树 成 男 , 35 岁 , 讲 师 , 博 士后 式 中 : D 孔为与小变形情况相 同的弹 塑性本构矩 阵 3[] : 汉为 。 nek e r 系数 ; G 二不具不 为刚性模量 ; 阶 ` u 少/ , ` ~ “ “ 一 “ ` 小 琳 ’ ~ 2 ( 1 + v ) / 廿 ’ “ 切 ’ 一 ’ 八 ~ ’ 云为等 效应 力 ; H ` 为等 效应 力对 于 塑性应 变 曲 线的斜率 . 按虚 功 原 理每 个 单 元 的 外 力 功等 于 内力 功 , 即 6琳 一 6 磷= 0可 得 : 无时 , = 川盆汁川毒+ 川舔 (3 ) 其 中 : 聪 一 工B涵 。迢kun dy 是 与小变形弹 塑性有 限 元 法 刚 度 矩 阵 相 同 的变形 矩 阵 , 反映 了 由 线 变 形 引 起 的 刚 度 ;心 , 一 工呱底办.kN * - ZB * 药两刀、 d) 、 是对物体大变形时刚性转动所 引 起 的 应力变 化部分 进行校 正 的应力校 正 矩 阵 . DOI: 10. 13374 /j . issn1001 -053x. 1999. 05. 038
·456· 北京科技大学学报 1999年第5期 K=K[Bd∫Bd-∫BBdy是对 沿2个轧辊连心线方向且正向指向上方,}轴指 体积不变约束引起的单元刚度进行了校正的等 向根据其他两坐标轴形成右手系, 体积变形矩阵,K30-2列 E 将圆管的1/4分成1216个节点、810个单 元.如图2所示. 狰水应力的求解方法:静水应力可以按体 积应变增量△E与体积弹性模量K进行求解,即: △=K△8b (4) 而在实际的有限元法计算中,用这种方法 求得的静水应力会产生很大的误差,因此本文 采用一种间接求解静水应力的方法,实践证明, 这种方法求得的静水应力有较好的精度 这种方法是把应力张量分为偏应力张量和 静水应力分量两部分,根据平衡方程可知: 0(a'+0)00x00=0 图2工件的初始网格 dx dy dz (5) g+o++g=0 Ox dy Oz (6) 3计算结果与讨论 景+号+ao-0 (7) 图3表示第1道次进入稳态后的速度场, 由于本文采用了间接求静水应力的方法, 图中清楚地反映出辊弯成形时金属的流动规 因此能可靠的显示工件的应力分布信息, 律.在刚开始咬入到孔型最低点处,金属向下方 2计算条件 向和扩展的方向流动,在最低点到出口处,即工 件处于弹复阶段时,金属向上方向和扩展的反 21孔型及其几何尺寸 方向流动,这与实际成形过程是相吻合的.从图 该道次只有上下2个水平辊没有立辊,2个 4和图5中的2个分速度图中可以看到,在顶部 水平辊安装好后形成的孔型如图1所示.其上 的Z方向速度比Y方向的速度大,向侧边的扩 下辊的最小半径为134mm,曲率半径为 展过程中Z方向的速度逐渐减小,Y方向的速 158.37mm. 度逐渐增大,而到边部Z方向的速度比Y方向 的速度小, 下辑 图3稳态后的速度场 图1第1道次孔型图 2.2试料尺寸及特性 母管半径R=44.683mm,壁厚t=7.407mm, 摩擦因数为μ=0.15,弹性模量E=196GPa,泊淞比 v=0.3,摩擦层厚度FIL=1m,压下量为28.09%, 材料的屈服函数=512.2(c+0.6)°26MPa,轧辊的 位移增幅为△a-0.03490658rad, 根据圆管的对称性和其变形特点,只取管 材截面的1/4部分为研究对象.坐标原点取在管 图4Z方向的速度场 件末端截面的中心处,X轴正向为轧制方向轴
. 4 5 6 . 北 京 科 技 、 域令如 ilm dy )( 执剧 一 扭、 汕 ,et] 对 体积 不变 约 束引 起的单元刚 度进行 了校正 的等 体积 变形 矩 阵 , =K 市髦刁 . 静水应力 的求解方法 : 静 水应力可 以按体 积应 变增量 △命与体积弹性 模量尤进行求解 , 即 : △ * = 尤△£* ( 4 ) 而在 实际 的有 限 元法计 算 中 , 用这种 方法 求得 的静 水应 力会 产生 很 大的误差 , 因此 本文 采用一 种 间接求解静水应力 的方法 , 实践证 明 , 这 种 方法 求得 的静 水应力 有较好 的精 度 `4] . 这种方法是把应 力张量分为 偏应 力张 量和 静水应 力 分量两 部 分 , 根据 平衡 方程可 知 : 大 学 学 报 1 9 , 9 年 第 5 期 沿 2 个轧辊连心线 方 向且正 向指 向上方 , ) 尹 轴指 向根据其他 两 坐标轴 形成右 手系 . 将 圆管 的 1 ,4/ 分 成 1 21 6 个节 点 、 81 0 个单 元 . 如图 2 所示 . 图 2 工 件 的初 始 网格 、户. 、少`2. 、 砚 、I, J了 O, z 厂 ` , z .J、了r . 由于 本 文采用 了 间接求 静水应 力 的方法 , 因此 能可靠 的显 示工 件 的应力 分布信息 . 2 计算 条件 2 . 1 孔型及 其 几 何尺 寸 该道次只 有上 下 2 个 水平辊没有立辊 , 2 个 水平 辊安装好 后形 成 的孔型如 图 l 所 示 . 其 上 下辊 的 最 小 半 径 为 13 4 m m , 曲 率 半 径 为 15 8 . 3 7 m m . 3 计算结果 与讨论 图 3 表 示第 1 道 次进入稳态 后的速度场 , 图 中清 楚地 反映 出辊 弯 成形 时金 属 的流 动规 律 . 在刚开始咬入 到孔型最低点处 , 金属向下 方 向和 扩展 的方 向流动 , 在最 低点到 出 口 处 , 即工 件 处 于 弹复阶段 时 , 金属 向上 方 向和扩 展的反 方 向流动 , 这 与实际成形过程是相吻合的 . 从 图 4 和 图 5 中的 2 个分速度 图中可 以看到 , 在顶部 的 Z 方 向速度 比 Y 方 向的速度大 , 向侧边 的扩 展过程 中 Z 方 向的速度逐渐减 小 , Y 方 向的速 度逐渐增大 , 而到边部 Z 方 向的速度 比 Y 方 向 的速度小 . 图 3 稳 态后 的速 度场 图 1 第 1 道次 孔型 图 2 . 2 试料 尺 寸及 特性 母 管半径尺 = 4 4 . 6 8 3 m , 壁 厚 r = 7 . 4 0 7 m m , 摩擦 因数为产= 0 . 巧 ,弹性模量=E 1% G P a , 泊 淞 比 v 二 .0 3 , 摩擦层厚度 IF L = 1 林m , 压 下 量为 2 & 09 % , 材料 的屈 服 函 数=0 5 12 . 2 (扮 0 . 6 ) 。 · 。 , ` M P a , 轧 辊的 位移增幅为△a = 0 . 0 3 4 9 0 6 5 8 r a d · 根据 圆 管的对称性和 其变 形特 点 , 只 取 管 材截面的 14/ 部分为研究对象 . 坐 标原点 取在管 件末端截面的中心 处 ” 丫轴 正 向为轧制方 向 2轴 图 4 2 方 向的速 度场
Vol.21 No.5 徐树成等:非圆截面管辊弯成形过程的变形和应力场的FEM模拟 ·457· (a) 300 300 200 200 100 100 0 0 5-100 -100 edW/'o -200 -200 图5Y方向的速度场 4080 图6为稳态成形后计算管壁厚的增厚率和 实测值的情况,本文模拟计算的材料与文献[2] 轧制方向位移mm 40 120160 宽度方向位移mm 0 0 实验用的材料相同.从图中可以看出,从顶部开 (b) 始壁厚逐渐增大,当到达工件刚刚脱离辊子时, 300 壁厚增厚率至最大,然后增厚率开始减小,但变 300 200 200 化不大,从理论值和实验值可以看出,在顶端的 100 100 增厚率相差较大,因为在计算时没有考虑焊缝 0 0 处的材质与其他位置的差异而导致的误差.其 -100 -100 增厚率的最大值为1%左右,而最小值仅为0.2 -200 -200 -300 -300 %左右,与实测值接近 -400 -400 16 80120 60 轧制方向位移 宽度方向位移mm 40 0 0 m 0 图7轧制方向的单元矩心G的三维分布, (a)内层单元:(b)外层单元 -8 一理论值 ·实测值 (a) -161LL1L上 0 0 010203040506070 沿外层单元表面的长度mm -100 -100 图6管壁增厚率T理论与实测值比较 -200 -200 EdW/ 图7(a),b)分别表示轧制方向上.的内外层单 -300 -300 元矩心应力σ的三维分布,从计算结果可以看 出,工件与轧辊接触区域显示为压应力状态,而 4080 非接触区域为拉应力状态,也就是说,从轧辊与 4080120160 轧制方向位移/mm 宽度方向位移mm 020 工件刚开始接触点到工件的顶部皆为压应力状 (b) 态,从接触点附近开始到工件的侧表面呈现拉 应力,从图中亦可看出,在成形过程中,当工件 0 进入稳态成形后,未进入孔型的工件开始有点 -100 -100 凸起,显示出内压外拉的应力状态 美 -200 图8(a),(b)分别表示成形过程中宽展方向的 -300 6 内外层单元矩心应力σ,的三维分布情况.从计 -300 算结果可以看出,在整个变形区内呈现出压应 4080120 60 力状态,但其最大压应力值出现在工件与轧辊 轧制方向位移 6 宽度方向位移mm 20 开始接触点附近,从图中可以看出,上流侧的压 m 0 应力值大于下流侧呈现出的压应力值. 图8宽展方向的单元矩心0,的三维分布. (a)内层单元:(b)外层单元
6 V 1 . N 2 1 0 . 5 徐 树成 等 : 非 圆截 面 管辊弯 成 形过程 的变 形和 应 力场 的 模拟 M F E 一 7 5 4 - 芝、招` d 芝、戈` d 芝、月` d 0 0 3 2 00 00 1 0 一 00 1 一 2 00 一 00 3 一 00 4 . 6 芝、 d 图 方 向的速 度场 S Y 6 图 为稳态 成形 后计算 管壁 厚 的增 厚率和 实测 值的 情况 , 本文 模拟计算 的 材料 与文 献【2] 实验 用 的材料相 同 . 从图 中可 以看 出 , 从顶部 开 始壁 厚逐渐增大 , 当到达工 件刚 刚脱离辊子 时 , 壁厚增 厚率至最 大 , 然后 增厚率开 始减 小 , 但变 化不大 . 从理论值和 实验值可 以看 出 , 在顶端 的 增厚 率相 差较 大 , 因 为在计 算时没有考 虑焊缝 处 的材质 与其 他位置 的差 异而 导致 的误差 . 其 增 厚率 的最大值为 1 % 左 右 , 而最小值仅 为 .0 2 % 左 右 , 与实测值接近 . 日ùà 攀、卜 一 理论 值 . 实测 值 图 7 轧制 方 向的单 元矩心 饥的三维 分布 . a( ) 内层 单元 ; (b) 外层 单元 . 芝、d目泞 。 ùō俨 一 1 6 ` ~ ~ 山- ~ 盛~ se ` . es 盛 . e 目 . . . 曰 . . . 曰 .」 0 1 0 20 30 4 0 50 6 0 7 0 沿外层 单元 表面 的长度 /m m 图 6 管壁增 厚率 T 理论 与实 测值 比较 一 10 0 一 20 0 . 芝、d ó 曰 图 7 ( a) , (b) 分别 表示 轧制方 向上 的 内外层 单 元矩心 应 力氏 的三 维分 布 . 从计算结果 可 以看 出 , 工件与轧辊接触区域显 示为压 应 力状态 , 而 非 接触区域 为拉 应力状 态 . 也 就是 说 , 从轧辊与 工件 刚开 始接触点到工 件的顶部 皆为压应力状 态 , 从 接触点 附近 开 始到 工件 的侧表面 呈现拉 应 力 . 从 图 中亦 可看 出 , 在成 形 过程 中 , 当工 件 进入 稳态 成形 后 , 未进入孔 型 的工 件 开 始有 点 凸起 , 显 示 出 内压 外 拉 的 应力 状态 . 图 8 (a) , b( )分别表示 成形过程 中宽展方 向的 内外层 单元矩心 应 力坏 的三维 分布情况 . 从 计 算结果可 以看 出 , 在整 个 变形 区 内呈现 出压 应 力状态 , 但其最大压 应力 值出 现在 工件 与轧辊 开 始接触点附近 . 从 图 中可 以看 出 , 上 流侧 的压 应力 值大于 下 流侧呈 现 出 的 压 应力值 . 一 3 0 0 鸯“ 碱 d芝. 、时 图 8 宽展方 向的单 元矩心 巧的 三维 分布 . a( ) 内 层单元 ; (b) 外 层单 元
·458· 北京科技大学学报 1999年第5期 4结论 寸方实验的研究一第7报一构形口一北仁女石角管 0成形(4).生产研究,1978,30(12:444 (1)从速度场图可以知道,在形成棱角越远 2小野田义富,大柴茂.刚塑性有限要素法仁为角 处金属的流动越大,由此可知角部的成形比边 钢管”口一儿成形了口七天”之子之3y.见: 部困难,因此根据模拟结果来修改孔型可降低 第38回塑性加工连合讲演会.富山,1987.285 成本,提高生产率,增加效益. 3徐树成.辊弯成形的三维大变形有限元模拟:[学位 论文].秦皇岛市:燕山大学,1997 (2)从应力场图可知,在棱角附近应力变化 4 Liu C.Simulation of Strip and Slab Rolling Using an Elas- 较大,容易产生裂纹,这样可根据模拟结果修改 tic-Plastic Finite-Element Method:[dissertation].Bir- 孔型,不断提高产品质量,生产出合格产品. mingham:Department of Mechanical Engineering,Uni- versity of Birmingham,1985 参考文献 】木内学,新谷贤,户澤正孝.异形管”口一北成形关 Deformation and Stress Fields of Non-circular Tubes in Roll Forming Process Xu Shucheng,Wang Xianjin,Liu Car 1)Material Science and Engineering School,UST Beijing,Beijing 100083,China 2)Yanshan University,Qing Huangdao 066004,China ABSTRACT A simulation on metal flow rule and stress distribution of non-circular tubes in roll forming process by using 3D large deformation elastic plastic finite element method with considering of material and geometrical non-linearity is presented.This method is based upon Prandtl-Reuss flow ruler and mises yield cri- terion.The simulation results of grids deformation flow pattern and stress distribution are obtained. KEY WORDS roll-forming;elastic-plastic finite element;metal flow;stress distribution;deformation
4 5 8 . 北 京 科 技 4 结论 ( 1) 从速度场 图 可 以 知 道 , 在 形 成棱角越远 处 金属 的流 动越大 , 由此可知 角部 的成形 比边 部 困难 , 因此 根据模拟 结果来修 改孔型可 降低 成 本 , 提高 生产 率 , 增加效 益 . (2 ) 从应力场 图可 知 , 在棱角 附近应 力变化 较大 , 容易产生 裂纹 , 这样可根据模拟结果修改 孔 型 , 不断提 高产 品 质量 , 生 产 出合格产 品 . 参 考 文 献 1 木 内学 , 新谷 贤 , 户泽 正孝 . 异形管 。 口 一 沙 成形 忆 关 大 学 学 报 19 9 , 年 第 5 期 寸 朽实 验 的研究一第 7 报一构形 口 一 沙 忆 上 石 角管 。 成 形 ( 4 ) . 生 产研 究 , 1 97 8 , 30 ( 1 2) : 4 4 4 2 小 野 田 义 富 , 大 柴茂 . 刚塑 性有 限要 素法 忆 巧 石 角 钢 管 。 口 一 羚 成形 弓 口 七 久 。 夕 之 二 卜 少 习 夕 . 见 : 第 38 回塑 性加 工连 合讲 演会 . 富 山 , 1 98 7 . 2 85 3 徐树 成 . 辊 弯成 形 的三 维大变 形 有 限元模 拟 : [学位 论 文 ] . 秦 皇 岛市 : 燕 山大 学 , 199 7 4 L i u C . Si m u lat i o n o f Sitr P an d S lab R o l li n g U s i n g an E l a s - ti e 一 Pl a s t i e F i n iet 一 E l e m e n t M e ht o d : [ d i s s e rt a ti o n ] . B i r - m i n gh am : D e P别rt m e in o f M e e han i e a l 助g in e e inr g , nU i - V e r s iyt o f B i n n l n gh am , 19 85 D e fo mr a t i o n a n d S tr e s s F i e ld s o f N o n 一 e i r e u lar uT b e s i n R o ll F o mr i n g P r o e e s s Xu hS u c h e 馆 , ), 肠 gn iX a nj in , ), Li u aC 产, l ) M at e r i a l S e i e n e e an d E n g i n e e r i n g S e h o o l , U S T B e ij i n g , B e ij in g 10 0 0 8 3 , C h i n a Z) 丫恤n sh an U n i v e sr iyt, Qi n g H u an g da o 0 6 6 0 04 , C h i n a A B S T R A C T A s im u l a t i o n o n m e t a l fl o w ur l e an d str e s s d i s itr b u t i o n o f n o n 一 e i r e u l ar t u b e s i n r o ll fo mr i n g P r oc e s s by u s i n g 3 D l a gr e d e fo mr at i o n e l a s t i e Pl a s t i e 加it e e l em e nt m e ht o d w iht e o n s id e ir n g o f m at e r i a l a n d g e o m e itr e a l n o n 一 li n e a r iyt 1 5 Pr e s e nt e d . T h i s m e ht o d 1 5 b a s e d uP o n P r an id l 一 R e u s s fl o w ur l e r an d m i s e s y i e ld e ir - t e ir o n . T h e s im u l a t i o n r e s u lt s o f gr id s d e fo mr a t i o n fl o w Pat em an d str e s s d i s itr b ut i o n ar e o b t a i n e d . K E Y W O R D S r o ll 一 fo mr i n g : e l a st i e 一 P l a st i e if in et e l e m e nt ; m e at l if o w ; s etr s s d i s itr b ut i o n : d e fo mr at i o n