（材料加工工程专业论文）铝合金板料深冲成形的晶体塑性有限元模拟.pdf

铝合金板料深冲成形的晶体塑性有限元模拟 摘要 金属板材经过机械变形或热处理以后通常会显示出各向异性，这种各向异 性会对板材成形性能产生较大的影响。近年来，随着计算机技术和有限元技术 的飞速发展，板材冲压成形的数值模拟取得了很大的进展，金属薄板冲压成形 过程的计算机模拟技术已经逐渐从实验室走向实际应用。为了提高数值模拟的 准确性，已经有大量研究致力于如何将各向异性性能引入到有限元模拟过程之 中。 为了满足汽车轻量化的要求，各种轻量化材料得到了迅猛的发展。与汽车 用钢板相比，铝合金板具有密度小、强度高、耐锈蚀等优点，能满足汽车轻量 化的特殊要求，已成为汽车轻量化技术中替代钢板的主要材料。由于晶体结构 上的差别，传统钢板的成形理论和技术不能直接用到铝合金板上，需要探索出 适合汽车用铝合金板自身的冲压成形理论和技术。 本文对5 a 0 2 铝合金板料深冲过程进行了实验研究，并基于晶体塑性变形的 微观滑移机理和硬化机制，在率相关晶体塑性本构理论框架下，实现了5 a 0 2 铝合金板料深冲过程的晶体塑性有限元模拟。研究了深冲过程中的制耳和板厚 的分布情况，揭示了最大厚度和最小厚度的变化规律，分析了摩擦系数和压边 力对制耳轮廓的影响，并进一步预测了极限拉深系数值。结果表明：晶体塑性 有限元法模拟的制耳轮廓的高度和极限拉深系数值与实测值比较吻合，随着摩 擦系数和压边力增大，制耳轮廓高度也增加，板料的厚度变化与实际情况一致。 关键词：深冲；制耳；极限拉深系数；晶体塑性：铝合金 s i m u l a t i o no ft h ed e e pd r a w i n go fa l u m i n u ms h e e t u s i n gc r y s t a lp l a s t i c i t yf e m a bs t r a c t t h ea n i s o t r o p yo fs h e e tm e t a le x e r t sa ne v i d e n ti n f l u e n c eo nt h ef o r m a b i l i t y i nr e c e n ty e a r s ，g r e a tp r o g r e s sh a sb e e nm a d ei nh o wt oi n t r o d u c ea n i s o t r o p y p a r a m e t e r si n t ot h ep l a s t i c i t yf i n i t ee l e m e n ta n a l y s i so fs h e e tm e t a lf o r m i n g i no r d e rt os a r i s f yi n c r e a s i n gr e q u e s to ft h ea u t o m o b i l el i g h t w e i g h t ，m a n y k i n d so fl i g h t 。w e i g h tm a t e r i a lh a v eb e i n gs t u d i e da n dd e v e l o p e d c o m p a r i n gw i t h t h es t e e ls h e e t ，t h ea l u m i n u ma l l o ys h e e th a ss m a l ld e n s i t y ，h i g hs t r e n g t h ，c o r r o s i o n r e s i s t a n c ea n ds oo n ，a n dt h e r e f o r ec a ns a r i s f yt h es p e c i a lr e q u e s to ft h ea u t o m o b i l e l i g h t w e i g h t a l u m i n u ma l l o y sh a v eb e c o m et h ep r i n c i p a lm a t e r i a l st h a ts u b s t i t u t e s t e e ls h e e ti nt h ea u t o m o b i l el i g h t w e i g h tt e c h n o l o g y b e c a u s eo ft h ed i f f e r e n c eo f c r y s t a ls t r u c t u r e ，t h et r a d i t i o n a lf o r m i n gt h e o r ya n dt e c h n o l o g yo fs t e e ls h e e tc a n n o tb ed i r e c t l yu s e di nf o r m i n go ft h ea l u m i n u ma l l o ys h e e t ，a n dw en e e dt oe x p l o r e n e wf o r m i n gt h e o r ya n dt e c h n o l o g yt h a ta r es u i t e dt ot h ea l u m i n u ma l l o ys h e e tu s e d i na u t o m o b i l e t h ed e e pd r a w i n go f5 a 0 2a l u m i n u ms h e e tw a se x p e r i m e n t a l l yi n v e s t i g a t e d ， a n df u r t h e ri nt h ef r a m e w o r ko fr a t e d e p e n d e n tc r y s t a lp l a s t i c i t yt h e o r y ，t h ed e e p d r a w i n gp r o c e s sw a sn u m e r i c a l l ys i m u l a t e db yc r y s t a lp l a s t i c i t y f i n i t ee l e m e n t m e t h o d ( c p f e m ) t h ee a r i n gp h e n o m e n o na n dt h ev a r i a t i o no fs h e e tt h i c k n e s s w e r es t u d i e d ，a n dt h ee v o l u t i o no ft h em a x i m u ma n dm i n i m u mt h i c k n e s sd u r i n gt h e d r a w i n gp r o c e s sw a sd i s c o v e r e d t h ee f f e c t so ft h e f r i c t i o nc o e m c i e n ta n dt h e b l a n kh o l d i n gf o r c eo nt h ee a r i n gp r o f i l ew e r es t u d i e d ，a n df u r t h e rt h el i m i t d r a w i n gc o e f f i c i e n to f5 a 0 2a l u m i n u ms h e e tw a sa l s oo b t a i n e d t h ep r e d i c t i o n s f r o mt h ec p f e ms i m u l a t i o n sa r ef o u n di nr e a s o n a b l ea g r e e m e n tw i t h t h e e x p e r i m e n t a lr e s u l t s ，w h i c hi n c l u d et h eh e i g h to ft h ee a r i n gp r o f i l ea n d t h ev a l u eo f t h el i m i td r a w i n gc o e f f i c i e n t ，i n c r e a s i n gf r i c t i o nc o e f f i c i e n to rb l a n kh o l d i n gf o r c e w i l le n a b l ei n c r e a s i n gh e i g h to ft h ee a r i n gp r o f i l e k e y w o r d s ：d e e pd r a w i n g ；e a r i n g ；l i m i td r a w i n gc o e f f i c i e n t ；c r y s t a lp l a s t i c i t y ； a l u m i n u ma l l o y 插图清单 图2 - 1 晶体弹塑性变形的乘法分解9 图3 - 1 实验主要设备8 0 吨开式双柱可倾压力机2 0 图3 - 2 拉深模具及关键尺寸2 1 图3 - 35 a 0 2 拉深圆筒件2 2 图3 - 4 压边力对制耳轮廓高度影响2 3 图3 - 5 摩擦系数对制耳高度的影响2 3 图3 - 65 a 0 2 板料厚度与压边力的关系2 4 图3 - 85 a 0 2 拉深系数为0 5 9 6 时的筒形件2 6 图3 - 95 a 0 2 坯料直径过大造成的拉裂2 6 图4 15 a 0 2 铝合金拉伸应力位移响应模拟与实验结果对比2 8 图4 - 2 深冲模拟用的网格划分( a - 坯料b 一压边圈c 一凸模d 一凹模) 2 8 图4 3 深冲的有限元模型及关键部位尺寸2 9 图4 - 4 深冲过程圆筒杯应力云图3 0 图4 - 5 深冲过程圆筒杯应变云图3 1 图4 - 6 深冲过程圆筒杯应力云图3 2 图4 - 7 实验深冲的筒形件3 2 图4 9 深冲过程最大和最小厚度的变化3 3 图4 1 0 不同压边力对制耳轮廓的影响3 4 图4 - 1l 不同摩擦系数对制耳轮廓的影响3 5 图4 - 1 2 模拟预测筒壁厚度与实测值的比较3 5 图4 - 1 3 不同压边力对筒壁厚度的影响3 6 图4 - 1 4 不同摩擦系数对筒壁厚度的影响3 7 图4 15 板料应力云图( 坯料直径8 5 m m ) 3 7 图4 一1 6 板料应变云图( 坯料直径8 5 m m ) 3 8 图4 - 1 7 实验深冲的筒形件( 坯料直径8 5 m m ) 3 8 图4 1 8 板料应力云图( 坯料直径9 0 m m ) 3 8 图4 1 9 板料应变云图( 坯料直径9 0 m m ) 3 9 图4 - 2 0 实验深冲的筒形件( 坯料直径9 0 m m ) 3 9 图4 2 l 正面应力云图( 完全深冲，直径9 0 m m ) 4 0 图4 2 2 侧面应力云图( 完全深冲，坯料直径9 0 m m ) 4 0 图4 - 2 3 正面应变云图( 完全深冲，坯料直径9 0 m m ) 4 l 图4 - 2 4 侧面应变云图( 完全深冲，坯料直径9 0 m m ) 4 1 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作及取得的研究成果。据我所 知，除了文中特别加以标志和致谢的地方外，论文中不包含其他人已经发表或撰写过的研究成果， 也不包含为获得金目巴工业太堂或其他教育机构的学位或证书而使用过的材料。与我一同工作 的同志对本研究所做的任何贡献均已在论文中作了明确的说明并表示谢意。 p 加 学位论文版权使用授权书 本学位论文作者完全了解金起王些太堂有关保留、使用学位论文的规定，有权保留并向 国家有关部门或机构送交论文的复印件和磁盘，允许论文被查阅或借阅。本人授权金日墨王些太 兰l 可以将学位论文的全部或部分论文内容编入有关数据库进行检索，可以采用影印、缩印或扫 描等复制手段保存、汇编学位论文。 ( 保密的学位论文在解密后适用本授权书) 名：钎求岔新躲 签字日期：孑年咿月，占日 ” 触乌 签字日期：。q 7 年q 月7 日 电话： 邮编 致谢 在合肥工业大学浓厚的学术氛围中，我顺利地完成了硕士研究生阶段的学 业，并按时完成了我的硕士论文。在此，我谨向所有曾支持与帮助过我的老师、 同学和朋友致以深深的谢意! 首先，我最衷心地感谢我的导师陈忠家副教授，恩师严谨的治学态度、活 跃的学术思想、丰富的实践经验以及对知识不倦的追求精神不仅让我在学业上 跃上了一个新的台阶，更让我懂得了人生的许多道理，让我受益终生。从论文 的选题到论文的整理、从论文详细审稿到最后定稿无不倾注了陈老师的辛勤汗 水，借此机会向恩师表达我最诚挚的谢意和深深的祝福! 在进行论文的工作中，得到了刘全坤、陈文琳、张强和祖方道等老师的支 持和关心，借此机会表示深深的感谢。本届同门何艳生、师弟张铧、杜海龙、 师妹汪祥鹏、刘萍，他们对我的研究工作提出了很好的建议和意见； 深深的感谢我的父母家人，有了他们的鼓励和支持，我才能够安心的完成 学业，他们是我前进的精神动力。 最后，由衷地感谢各位评委在百忙之中抽出时间给我评审。 作者：韩永志 2 0 0 9 年4 月 第一章绪论 1 1 引言 先进塑性成形技术在促进国民经济发展，推动国防现代化和提高国家竞争 力方面发挥着越来越重要的作用。然而基于理论、经验和反复实验的研究方法， 现在已经难以满足先进塑性成形理论研究和成形技术的数字化、高技术化发展 的需求。计算机建模仿真与理论及实验研究有机结合，能够实现虚拟成形，将 大量反复实验在计算机上完成，可以进行一些理论和实验暂时还做不到的研究。 因此它已经成为研究和发展先进塑性成形技术，高质量、低成本、短周期地实 现塑性成形产品开发的强有力工具。 但是，这一强有力的分析工具的可靠性首先取决于材料本构模型的准确性。 半个世纪以来，众多学者运用统计学等理论对多种材料建立了多种宏观本构模 型。然而，航空、航天、汽车及一系列高新技术的发展向塑性加工技术提出了 精确控形和精确控性的双重要求。这就给计算机数值模拟提出了更高的要求： 不但要能够预测宏观成形规律，更要能揭示微观成形机制，预测产品的微观性 能。这就对宏观本构模型提出了挑战，而晶体塑性理论和细观力学的发展与成 熟，为建立材料细观结构、内部缺陷与宏观力学之间的定量关系提供了依据， 从而形成了一套在细观尺度上描述材料成形的新的本构理论框架。然而，晶体 塑性理论涉及力学、材料学、数学、计算科学等多学科知识的交叉，其模型复 杂，数值化难度大。所以，晶体塑性模型的数值化研究已成为先进塑性成形前 沿研究领域的重要课题，也日益成为国际上的研究热点之一。 1 2 国内外晶体塑性理论及应用的研究概况及评述 晶体的大应变塑性变形常常伴随着应变局部化、塑性损伤等现象，经历了 大应变塑性变形的物体，内部的应变场、应力场一般具有强烈的非均匀性，这 样的问题难以用通常的解析法研究，必须用有限元法进行数值模拟，至今国内 外尚未见任何基于晶体塑性理论内核的商用有限元软件发布，且现有商用软件 均是建立在宏观各向异性屈服函数的基础之上，不能描述上述应变局部化、塑 性损伤等现象，而晶体塑性有限元法( c r y s t a lp l a s t i c i t yf i n i t ee l e m e n tm e t h o d ， c p f e m ) 基于晶体塑性变形的微观滑移机理和硬化机制，采用了基于晶体塑性 理论的运动学描述和率相关的本构方程的晶体塑性有限元法，可以对这些现象 进行数值模拟。c p f e m 的基础理论最初由h i l l 和r i c e 提出 1 - 3 1 , 它主要依据的 是晶体塑性变形的s c h m i d 定律。 1 2 1 国内晶体塑性塑性理论及其应用的发展 上个世纪九十年代以来，我国的学者们也对晶体塑性有限元进行了广泛而 深入的研究。 1 9 9 2 年，梁乃刚等人【4 】提出了三维组集式本构模型。它是一种半物理半数 学的弹塑性本构模型，与数学塑性本构理论的不同之处在于由模型结构的宏观 力学效应来导出本构关系，避免了将变形分解成弹性变形和塑性变形两部分所 招致的理论上的麻烦。模型中的构元随宏观变形而伸缩和转动，使得几何非线 性效应对于材料性质的影响得到了反映。因此，这种模型不仅能模拟复杂加载 情况，而且也能解释非均匀局部变形现象。这个模型与经典物理塑性理论的不 同在于以简单的一维构元组成能量等价系统代替真实的细观结构，得到显式的 本构关系，保持了简洁的数学描述方式，可以较方便地应用于结构分析与工程 设计计算。基于此模型对晶体单轴拉伸中的宏观剪切带的分析，结果与实验吻 合。从而说明此本构模型能够模拟有限变形中的几何非线性效应和晶体材料塑 性变形中的宏观力学行为。本模型没有涉及到取向不同的纤维在变形中的相互 作用，因而无法反映材料内部组织的潜在强化效应。 1 9 9 3 年，仲政等人【5 】在等效夹杂的基础上，提出了均值自洽理论。并给出 两种求解速度梯度集中系数张量和应力率集中系数张量的方法：均值等效夹杂 法和自洽有限元法。此理论能有效预计由于晶粒形状非椭球和晶界滑错而导致 的不均匀变形对多晶本构的影响。这一理论将仅适用于固连椭球形状夹杂的 h i l l 自洽理论推广到适用于任意形状固连或界面可滑错夹杂的情况。 1 9 9 6 年，刘志宏等人【6 】在“三维组集式本构模型 研究的基础上，放弃了 “均匀应力”假设，发展了一种称之为“取向元”的概念，使这种新型的本构 理论更具有一般性。在一维“取向元 上考虑率敏感效应，也使得粘性的描述 变得非常简单，不必再引入“过应力 、“背应力”、“平衡应力”等容易混淆的 概念，方便工程应用。并由率相关的取向元推导出了多晶金属材料弹粘塑性的 本构方程，基本描述了粘塑性理论中的两大问题。与以往的粘塑性模型相比， 该模型不仅能同时考虑多晶金属材料的率相关性和路径相关性，而且由于引入 了物理机制，使模型具有较强的预测能力，数值算例也证实了该模型的合理性。 19 9 7 年，董湘怀等【7 】开发了基于晶体塑性理论的板料成形过程的动力显式 有限元分析程序，对单晶体和多晶体的塑性失稳以及织构对板材成形性能的影 响进行了深入地研究，提出了统一的单晶体硬化模型；在率相关晶体塑性模型 的计算中提出了忽略非活动滑移系的算法，大大缩短了计算时间。 2 0 0 0 年，扶名福等人【8 】讨论了变形梯度弹塑性乘积分解的唯一性。在不同 构形中给出了塑性应变在晶体塑性中的表述，建立了塑性滑移率与塑性应变及 应变率之间的关系，它们为有限塑性变形的宏微观桥梁关系及晶体塑性的进一 步研究提供了基础。 2 0 0 1 年，苏世忠等人【9 】开发了一套适用于f c c ( 面心立方) 金属的三维弹 塑性有限元程序，该程序将与应变速率无关的多晶体弹塑性模型引入到立方体 2 单元晶粒塑性应变增量的计算中。以a 1 板为例计算出常温下单向拉伸和压缩、 双向压缩、平面应变、面内剪切和厚向剪切等六种典型受力状态下的织构演变 和各向异性变化，并对a 1 圆管扭曲与拉伸复合加载和a 1 圆管冷径缩挤拔加工 进行了数值模拟，结果显示该模型能够较准确的从细观力学角度预测出多晶体 金属的织构演化，又能够从宏观角度预测出塑性变形行为。 2 0 0 2 年，冯露、张克实等人【i o 】对延性单晶在拉伸载荷作用下的应变局域化 和颈缩等非均匀变形过程进行了三维有限元数值模拟，将相关晶体塑性本构模 型及一种新的数值积分方法补充到a b a q u s 6 1 商用有限元软件中。该方法的 特点是，利用晶体塑性的动力学方程，获得一个关于晶体弹性变形梯度的演化 方程，采用半隐式积分方案进行求解，并推导出一种新的应力应变本构矩阵， 使计算速度和计算稳定性大大提高。加载方式、边界条件和变形程度等因素影 响着滑移系的启动状况，这是平面模型所不能预测的。他们还利用三维有限元 方法模拟了不同取向下滑移系的启动状况，全而地考虑了f c c 单晶材料1 2 个 可能滑移系在变形过程中的启动状况，合理地模拟了f c c 面心立方单晶沿不同 取向加载时晶轴旋转导致的应变局域化和颈缩等非均匀变形过程。2 0 0 3 年，冯 露、张克实等人【l l 】采用宏观h i l l 模型和晶体塑性模型对面心立方单晶( f c c ) 材料的非均匀变形进行了数值模拟，比较了两种不同尺度的模型对塑性各向异 性的描述能力的差异。为了使两种模型具有可比性，对于f c c 单晶材料提出一 种用晶体塑性模型来确定h i l l 模型中各向异性参数的标定方法。数值分析表 明，两类模型对单晶体塑性各向异性的描述能力存在着差异。对f c c 单晶材料， h i l l 模型对各向异性的预测能力没有晶体塑性模型细致，晶体塑性模型更能追 踪塑性各向异性的变化，但两种模型对应力应变响应预测的趋势是一致的。2 0 0 4 年，张光、张克实等人【1 2 】研究了不同硬化模型的选取对多晶集合体应力应变响 应的影响，结果显示对多晶集合体应力应变响应的预测有明显影响但对织构演 化的预测影响不大。同年，他们【1 3 】对晶粒数量对多晶集合体初始各向异性的影 响也进行了研究。讨论了三种确定晶体随机取向的方法。计算结果表明：晶粒数 量有限的多晶集合体的应力应变响应仍有一定的各向异性，且随着晶粒数量增 多，多晶集合体的各向异性程度降低；就所包含晶粒数相同的多晶集合体来说， 在确定晶粒随机取向时，选取不同的方法对它的各向异性程度也有一定的影响。 2 0 0 5 年，肖纳敏等人【1 4 】为了定量描述热变形奥氏体在介观尺度上微观变形 的非均匀性，采用了晶体塑性有限元方法( c r y s t a lp l a s t i c i t y f i n i t ee l e m e n t m e t h o d ，c p f e m ) 模拟了c m n 钢热变形奥氏体在不同变形条件下的变形行为， 并将模拟所得到的应力应变曲线与文献测定的应力应变曲线进行了对比，模 拟的曲线和测定的曲线基本一致。同时计算了在晶粒尺度下奥氏体的微观应力 应变分布，结合塑性变形的位错模型，给出了变形储存能的微观分布。通过对 热变形奥氏体的研究发现，即使在外部的均匀变形条件下，无论是在晶粒内部 3 还是晶粒间，材料内部变形都非常不均匀。这种变形不均匀性主要是由晶粒的 初始取向不同、近邻晶粒的取向差、以及变形时滑移系的运动特性不同所引起 的。定量描述了介观尺度上奥氏体变形储存能不均匀分布，为结合介观尺度组 织模拟，实现组织演变的多尺度耦合计算提供了参考。 2 0 0 6 年，郑成武、肖纳敏等人【l5 】采用晶体塑性有限元( c p f e m ) 和元胞自动 机( c a ) 耦合的方法模拟了热变形低碳钢的静态再结晶。c p f e m 的计算结果定 量描述了介观尺度上奥氏体变形储能的不均匀分布，为模拟再结晶的形核和长 大提供了依据，从而在再结晶c a 模型中考虑了不均匀分布对再结晶组织演化 以及动力学的影响。模拟结果显示：变形储能分布不均匀导致再结晶晶粒在不 同位置的形核概率以及形核密度不同，位置分布也不均匀。临界形核储能值较 高时，形核主要集中在原始晶界，在晶内只有少量晶核形成；随着临界形核储 能的降低，晶内形核数量增加，结晶位置分布趋于均匀。增加晶内变形带，使 晶内在结晶形核数量增多，对细化再结晶晶粒有利。 1 2 2 国外晶体塑性塑性理论及其应用的发展 1 ) 晶体塑性有限元模拟制耳形成 制耳往往是由板材中的晶体学织构引起的深冲制品的壁缘呈波浪形现象。 由于制耳的存在，导致材料浪费、产品质量差、成品率低等不良后果。对板材 深冲制耳进行模拟预测和控制具有重大的现实意义。很多研究者运用晶体塑性 有限元方法模拟了制耳及其织构的演化，以达到预测的目的。 b a l a s u b r m a n i a n 和a n a n d 1 6 j 的研究表明晶体塑性理论可以研究并优化深冲 工艺。z h o u 等【l7 】模拟了织构对面心立方材料制耳的影响。i n a l 等人【1 8 】则用一 个六分量屈服函数模拟了制耳的形成过程。 2 ) 晶体塑性有限元模拟组织织构 多晶体金属材料经过塑性变形后导致晶粒取向的分布不均匀，微观上在其 内部形成了明显的织构，宏观上表现出塑性各向异性。晶体塑性有限元法已经 可以成功地模拟预测织构及其组织演化。 s a r m a h 等人【l9 j 提出了一个整合了邻近晶粒相互作用的多晶塑性模型，并模 拟预测了织构的演化。研究结果显示此模型能够提高织构演化的预测能力。 k u m a r 等人【2 0 j 在非欧拉取向空间对晶体织构演化进行了有限元模拟，预测结果 与传统的欧拉角参数相比，非欧拉参数具有很大的优越性。k u m a r 等人【2 l 】运用 新的织构分析工具在r o d r i g u e s 空间对面心立方材料的变形织构进行了计算机 模拟，模拟结果直接解释了理想织构分量分布和取向织构的关系。c h o i 等【2 2 】 对面心立方板料金属深冲过程中理想取向的稳定性和织构演化进行了研究。 d e l a n n a y 等l 2 3 j 定量预测了中等变形冷轧铝板的织构。关于晶体塑性有限元法在 织构模拟方面的文献还有很多，限于篇幅不复赘述。 4 3 ) 模拟在其他方面的应用 晶体塑性有限元同样也广泛应用于模拟材料的性能响应如屈服平面【z 4 。26 。、 微观应力应变分布【2 孔、冷轧条件下形变对再结晶的影响【2 s 、加工件表面粗糙度 研究 2 9 】和成形极限图 3 0 - 3 1 】等方面。 b e a u d o i n 等【2 4 】通过模拟研究了由聚合体加权取向得出的屈服平面。研究结 果表明试验织构与预测织构一致。k a l i d i n d i 和s c h o e n f e l d 2 5 】利用有限元技术提 出了数值均匀过程预测具有不同织构的面心立方多晶体的屈服平面。结果表明， 在确定各向异性应力响应时，屈服平面可能比择优取向更重要。b o u d e a r u 等【j z j 发展了三维的晶体塑性有限元关系并用之分析了有5 0 个晶粒的多晶体的单向 拉伸变形，还同时考察了形变引发的各向异性对成形极限图( f l d s ) 的影响，首 次把c p f e m 直接应用到生产实践上。r a a b e 等【3 3 j 用元胞自动机结合晶体塑性 有限元的方法对冷加工后a l 的初始静态再结晶进行了模拟。 1 3 选题的背景及意义 随着计算机技术和有限元技术的飞速发展，金属薄板冲压成形过程的计算 机模拟技术已经逐渐从试验室走向实际应用。同时随着航天航空及汽车工业的 发展对金属板材性能的要求也越来越高。 晶体塑性理论将塑性变形归结为：在晶体中特定滑移系上的位错运动在该 滑移系中所造成的剪切应变。晶体材料的宏观变形是细观位错滑移运动积累的 结果，位错滑移运动本质上导致了材料塑性变形的各向异性和率相关性，而且 材料微观残余应力的存在导致了材料后继屈服的各向异性。板材呈现各向异性 特性，对于板材的成形性能有很大的影响。例如在筒形件冲压过程中可导致制 耳产生，并且在深冲过程中，在制耳间的波谷或制耳本身上可能产生开裂的危 险，严重地影响板材的深冲成形与使用。而率相关晶体塑性有限元模型可以从 细观和宏观尺度有效地描述以上材料的属性，揭示材料成形过程的细观变形机 制及演化。所以，运用晶体塑性模型对金属板料进行有限元模拟，已成为国内 外有关金属塑性大变形分析及金属板材成形性能研究领域内的前沿课题及发展 方向之一。在新型材料不断出现、成形性能要求愈来愈复杂的背景下，这对于 指导冲压成形工艺、以及优化材料细微观组织等，都有着十分重要的理论和实 际意义。 通过对国内外的研究状况分析可知，晶体塑性模型的应用较多地局限于二 维模拟和拉伸、镦粗等简单成形过程的三维模拟。而基于率相关晶体塑性模型 对复杂成形过程( 材料非线性、几何非线性和边界条件非线性的大塑性变形问 题) 的三维有限元模拟研究则少见报道。本文的选题正是在这样的背景下提出 来的，并得到国家自然科学基金项目“铝合金粘塑性行为的多尺度研究 ( n o 1 0 7 0 2 0 1 8 ) 的资助。 1 4 本文的主要研究内容 本文的主要研究内容如下： 1 基于晶体塑性滑移理论，采用切线系数法，利用f o r t r a n 语言将率相关本构模型嵌 入到大型商用有限元软件a b a q u s 用户子程序u m a t 中，使之能够用于板料的有限元 模拟： 2 对5 a 0 2 铝合金进行了单向拉伸实验，得到了5 a 0 2 铝合金的力学性能数据，并进 行了晶体塑性有限元模拟，然后将模拟结果与实验结果进行比较分析，验证模拟结果 的准确性； 3 对5 a 0 2 铝合金进行了筒形件拉深成形实验，得到了5 a 0 2 铝合金的极限拉深系数 以及典型缺陷一制耳的形成，并用晶体塑性有限元法进行了模拟，成功的模拟出了制 耳现象和极限拉深系数值，然后将模拟结果与实验结果进行比较分析，验证模拟结果 的准确性； 4 模拟分析深冲过程中最大最小厚度的变化； 5 考虑板料与模具之间的摩擦和压边圈的压边力，对深冲过程中制耳轮廓高度和侧 壁厚度变化的影响规律。 1 5 论文的章节安排 本文的章节安排如下，共有五个章节： 第一章，绪论：总结了国内外专家学者在晶体塑性领域研究工作的现状， 分析了他们在研究过程的优缺点，为本文的研究工作提供了思路。并给出了本 文的研究内容、选题背景及意义。 第二章，晶体塑性本构理论与数值算法：全面介绍了晶体塑性本构理论的 基本概念，建立了本文所采用的基于晶体塑性运动学描述的率相关本构关系， 给出了相应的数值算法，并说明了利用商用有限元软件a b a q u s 用户子程序 u m a t 实现本文采用的率相关本构关系的方法。基于有限变形晶体塑性滑移理 论，采用切线系数法，利用f o r t r a n 语言编制了率相关晶体塑性滑移本构模型的 有限元子程序，并将其植入到大型商用有限元软件a b a q u s 的用户子程序 u m a t 中，使之能够用于有限元模拟。 第三章，铝合金塑性变形的实验机理研究：通过对铝合金的拉深实验，得 到了a l 合金的成形极限值( 极限拉深系数) 以及其典型缺陷制耳的形成，分 析摩擦系数和压边力对拉深拉深件制耳轮廓高度和拉深件壁厚的影响，以及测 定了极限拉深系数，为后文采用晶体塑性有限元法模拟奠定了实验基础。 第四章，铝合金塑性变形的数值模拟：本章通过使用自己编写的u m a t 子 程序，采用基于率相关晶体塑性有限元法对5 a 0 2 铝合金板材进行了深冲过程 6 进行了数值模拟，并将模拟结果与实验进行了比较，分析了深冲过程中等效应 力、应变的分布情况以及筒形件整体厚度、侧壁厚度的变化情况，还考虑了摩 擦和压边力对深冲制耳的影响。并迸一步模拟预测了极限拉伸系数值，并与实 验值进行对比，结果证实c p f e m 模拟深冲过程是可行的和准确的。 第五章，全文总结及展望：对全文的研究工作进行了总结，并提出后续研 究的方向。 7 第二章晶体塑性本构理论与数值方法 2 1 引言 本构关系是材料的固有属性，是材料变形过程中必须遵循的客观规律。将 这些规律以数学公式的形式表达出来，并通过这些公式的计算从而预测材料变 形的实际行为，就是材料的本构建模过程。在这中间，对材料变形的客观规律 的正确认识是首要的，同时，将这些正确认识公式化和模型化从而得到可靠的 本构模型则是该过程的关键和核心技术。我们采用的晶体塑性理论是从考虑晶 体材料变形的微观物理机制入手，构建了描述材料宏观变形的细观本构理论框 架。 2 2 晶体塑性理论基础 晶体塑性变形行为自上世纪3 0 年代直至6 0 年代一直是学者们感兴趣的问 题，早期研究的成果主要体现在晶体塑性变形的几何学与晶体学特征上，在这 方面整理出大量的实验规律，为塑性变形的位错机制奠定了基础。正是这些关 于塑性变形和强度理论的研究导致了位错概念的提出和位错理论的发展，位错 理论的成熟，又反过来促进了晶体塑性理论的发展。现代塑性本构关系的描述 较之数十年以前，有了质的提高。 晶体塑性理论的早期研究，应当首推t a y l o r 、o r o w a n 、s c h m i d 和p o l a n y i 的开创性工作。这些先行者已经开始意识到塑性变形和晶体结构的密切关系， 它们的研究工作明确无误地提醒以后的研究者，研究晶体的塑性变形必须从晶 体的微观结构入手，决不能仅停留在宏观层次上。 金属材料在塑性变形过程中往往呈现出各向异性，这主要是由微观晶体的 结构特性所致。与此同时，金属材料的塑性变形常常伴随着应变局部化、塑性 损伤等现象，经历了大应变塑性变形的物体，内部的应力场、应变场也呈现出 明显的非均匀性。宏观唯象模型及其有限元方法，虽然可以迅速、充分地预测 金属塑性变形的宏观力学行为，但却不能很好的预测金属变形过程中微结构的 演化以及由此所引起的力学性能的变化。晶体塑性理论从塑性变形的物理本质 出发，考察与宏观变形相关的微观变形机理，能够更准确地反映材料变形的本 质，因此有必要基于晶体塑性理论对金属材料在塑性变形过程中的各向异性和 非均匀性进行模拟和预测。 2 2 1 晶体塑性变形几何学与运动学 晶体的变形可分解成弹性和塑性两部分，弹性部分由于晶格的畸变( 包括晶 格的刚性转动) 造成，而塑性部分则是由于晶体中的位错在特定的结晶学平面沿 一定的方向滑移造成的。晶格畸变的宏观表现就是弹性变形，可以用弹性力学 8 的方法处理。难点在位错滑移的处理：一方面，位错的分布是离散的，无法用 “连续”介质的概念来描述；另一方面，位错滑移在晶粒内部产生的是间断的 位移，这种不连续的位移也不可能用连续的位移梯度来描述。为了客服这个困 难，力学界提出了一种均匀滑移模型来描述晶体中位错的运动。因为晶粒内部 位错的数量是巨大的，从宏观角度来看滑移在晶粒内部是均匀的，从而就可以 用连续介质的场变量变形梯度张量来描述滑移的宏观效应，这种模型被称为均 匀滑移模型。有了这个模型，就可以用连续介质力学的方法来描述晶体中位错 的运动。 首先考虑单晶体的变形，按照乘法分解，总的变形梯度f 表示为 f = f 8 f pd e t f 8 0( 2 1 ) f 8 = 尺8 u 。( 2 2 ) 式中，f p 是塑性变形梯度，假设塑性变形过程中体积不变，则d e t f p = 1 ，f 8 是 弹性变形梯度，r 8 表示点阵旋转张量，u 。表示晶格畸变的右伸长张量。 韧始构形 象卸裁焉构形 ( 或当前构形 卸载后构形 ( 或等倾中间构形) 图2 - 1 晶体弹塑性变形的乘法分解 如图2 1 所示，当晶体滑移时，晶格矢量并不发生变化，而当晶格畸变时， 晶格矢量将发生伸长和转动。因此先介绍下面两个公式 m 口= f 8 m a ( 2 3 ) 卉口= 刀口 f 8 r 1 ( 2 - 4 ) 式中，当前构形中第口个滑移系滑移方向上的单位向量用m 口来表示，滑移面上 的单位法向量用i t 口表示，晶格畸变后滑移方向上的向量将变为磊口，滑移面上 9 的法向量变为五口。 考虑( 2 1 ) 式，晶体当前的速度梯度张量用上表示，可分解为弹性和塑性两 部分 l = 户e f 。一1 + f e f p f p l f 。一1 ( 2 5 ) 不难证明，得 声f p 。1 = y 多口朋口。刀口 ( 2 - 6 ) z 一 式中，7 口表示第口滑移系引起的切应变率。 把将公式( 2 6 ) 代入到公式( 2 5 ) 得到 l = l 。+ r ( 2 7 ) l e = 户e f “( 2 8 ) 口= d 尸+ w 尸= 羔乡口茄口q 磊口 ( 2 9 ) 其中r 是弹性变形速度梯度张量，p 是塑性变形速度梯度张量。 三可分解为对称部分( d ，变形率张量) 和反对称部分( 形，旋转率张量) 变形率张量d d=d8+dp( 2 10 ) d p = 芝p 口夕口 ( 2 11 ) p 口= 吉( r 矗口。矗口+ n a9r ；l 口) ( 2 - 1 2 ) 旋转率张量形 w = w 。+ 矽p w k 形口尹口 一 玎，口= 言( r 矗口 n + 口一矗口。r h 口) 这些都是晶体变形与运动学的基本公式， 形率之间的联系。 1 0 ( 2 1 3 ) ( 2 1 4 ) ( 2 15 ) 它们建立了滑移剪切率和宏观变 2 2 2 率相关单晶体塑性本构理论 假设晶体的弹性性质不受滑移变形的影响，根据h i l l 和r i c e 3 3 1 的描述， 可采用如下单晶体弹性本构关系 ；。：；一we r + f w 。= e m t ：d 。 ( 2 16 ) 式中；。表示以中间构形为基准状态( 基于w 。) 的k i r c h h o f f 应力张量f 的刚性 导数( j a u m a n n 导数) ，；是k i r c h h o f f 应力张量f 的物质导数，e m t 是晶体弹 性模量的四阶张量( e l a s t i cm o d u l u st e n s o r ) 。 当单晶体进一步发生塑性变形后，则基于形的j a u m a n n 导数f v ( 应力率) 表示为 ；= ；一w f + f w ( 2 1 7 ) 分别由公式( 2 16 ) 和公式( 2 1 7 ) 定义的应力率之差为 ；e f v = 羔口乡口 ( 2 - 1 8 ) 其中 口= w 口f f 形口 联立式( 2 1 6 ) 和( 2 1 8 ) 得到单晶本构方程为 ；：善。一( 善e 一；)f2 f 。一( f 。一r ) = e m t ：d 8 一口夕口 = e m t ：( d d 尸) 一窆g 乡口 将公式( 2 1 1 ) 代入上式，得： ；：e m t ：d 3 一夕尺a 户 f2： 一 夕ky “ 式中的尺口即为 r 。= e m t ：p 。+ b 公式( 2 - 2 0 ) 的本构方程中多口的表达式尚未确定， 表示出来之后，此本构方程才算完整。 下面求乡口。 根据p e i r c e 等人3 4 1 提出的率相关滑移模型： ( 2 1 9 ) ( 2 2 0 ) ( 2 2 1 ) 用应力张量和应变张量把乡口 胪罗口s s n ) ( 2 - 2 2 ) s g nc f 口，= 二。 耋：三吕( 2 - 2 3 ) 式中，s g n ( ) 是符号函数，表示取正负的意思，尹席是切应变率，多。是参考切 应变率，g 口是硬化函数，r 口是滑移系口( 口值为1 - 1 2 ) 上的分切应力，m 是 滑移系应变率敏感指数。 为了计算方便认为硬化函数g 口仅是滑移总变形量，的函数 g 口= g 口( ，) ( 2 - 2 4 ) y227 口i ( 2 - 2 5 ) 当y = o 时，g 口指定为f o ，而在多滑移系中，单晶应变硬化规律为 羔锄l 少 p = l ( 2 2 6 ) 式中，硬化模量矩阵k 为y 的函数，其公式3 5 琊1 为 k = 乃( 7 ) = 办。s e c 乃2 i ；与f ( 不对口求和) ( 2 - 2 7 ) 矗船= q 乃( y ) ( 口)( 2 - 2 8 ) 式中：h o 是初始硬化率，f o 是屈服切应力，f 。是流动应力饱和值。q 是硬化比 率，即潜硬化和自硬化效应的比值，一般取l q 1 ，迭代照常进行，若p n e w d t l ，则必须终止迭代，放 弃当前的增量步，尝试采用一个较小的增量步，新的增量步取 p n e w d t d t i m e 。经过多次试验，在对各状态变量做跟踪观察之后，发现对 时间步长最敏感的变量为滑移应变率，如果滑移应变率的数值能稳定下来，其 它状态变量一般都能稳定下来。根据这个经验，我们决定用滑移应变率来控制 迭代过程。 各晶体的初始取向角通过状态矩阵q 指定，q 实际上就代表各矢量或张量 从参考状态到初始状态的旋转角。参考状态下的滑移方向矢量、滑移面法向矢 量、正交晶体的三个弹性常数都作为不变值由p r o p s 命令输入。滑移系八渭 移系硬度g 口、k i r c h h o f f 应力张量丁和旋转张量f 。在积分点上的值均作为材料 状态参量存储在u m a t 的状态变量中，t = 0 时f 。反映晶格的初始取向，即 f 8 ( f ：o ) = q 。 2 4 2 本构关系的实现 应力率用k i r c h h o f f 应力张量f 的j a u m a n n 导数表示，即 1 6 ；= e m t ：d 羔r 口户 ( 2 4 4 ) 上式中的e m t 是旋转以后的晶体弹性模量张量。 根据上节式( 2 4 3 ) ，利用中心差值，得到滑移应变速率的表达式为 尹口= 口+ f 口：d ( 2 4 5 ) 代入本构关系中得 ；= 脚一r 口p 】：d 一r 口口 ( 2 4 6 ) 上式右端求得的应力率乘以时间步长就得到应力增量，上增量步的应力 旋转到当前积分点的取向中，与应力增量相加，用公式表示就是 f 似= i ，。y r + a t e m t 一r 口严】：d 一缸尺口夕口 ( 2 - 4 7 ) 式中y 即为a b a q u s 提供的数组d r o t 。 在设计子程序时，将这样求得的总滑移量7 、滑移系硬度g 口等作为内部变 量储存在u m a t 中，刚度矩阵d d s d d e 乘以应变增量d