刚塑性有限元法及其在轧制中的应用

1、刚塑性有限元法及其在轧制中的应用刚塑性有限元法及其在轧制中的应用轧制技术及连轧自动化国家重点实验室1.1.学习目的和要求学习目的和要求l了解现代轧钢生产和轧制技术的发展概况；l了解现代轧制理论研究的基本任务；l掌握刚塑性有限元的基本概念；l掌握刚塑性有限元的基本理论；l掌握刚塑性有限元的基本方法；2.2.学习的主要内容学习的主要内容l刚塑性有限元的基本概念和基本理论；l刚塑性有限元相关技术问题的处理方法；l求解轧制过程的刚塑性有限元程序。3.3.本课程的基础和相关知识本课程的基础和相关知识l 现代塑性加工力学现代塑性加工力学 基本方程、变分原理、有限元基础知识；l 工程数学工程数学 矩阵分析、

2、优化方法、数值分析；l计算机基础知识计算机基础知识 操作系统、FORTRAN语言和FORTRAN4.0编程软件。4.4.讲课和学习方法讲课和学习方法l 课堂讲授课堂讲授 基本概念、基本理论、基本方法 程序剖析；l 课外自学课外自学 消化理解、阅读程序；l上机实践上机实践 调试程序1.1.绪绪 论论 1.1 1.1 现代轧制理论研究的发展概况现代轧制理论研究的发展概况1.1.1 1.1.1 现代轧钢生产的发展现代轧钢生产的发展1.1.2 1.1.2 轧制技术的发展轧制技术的发展1.1.3 1.1.3 现代轧制理论研究的基本任务现代轧制理论研究的基本任务1.2 1.2 轧制理论数值方法轧制理论数值

3、方法1.2.1 1.2.1 初等理论中的数值方法初等理论中的数值方法1.2.2 1.2.2 滑移线理论及其数值解法滑移线理论及其数值解法1.2.3 1.2.3 能量法及其数值解法能量法及其数值解法1.2.4 1.2.4 弹塑性有限元法弹塑性有限元法1.1 1.1 现代轧制理论研究的发展概况现代轧制理论研究的发展概况 20世纪世纪60年代前年代前，轧钢生产过程手工手工操作和使用单体单体设备。 轧制理论主要解决问题轧制力、力矩、功率、宽展和前滑等参数的近似近似计算。 主要进展提出卡尔曼和奥罗万方程，采用一些假设条件推导出轧制力和宽展等公式，逐步形成了以工程法为核心的传统轧制理论体系。 20世纪世纪

4、60年代以后年代以后，随着轧钢生产轧钢生产和轧制技术轧制技术的飞跃发展和用用户对产品质量户对产品质量要求的日益提高，以计算机计算机为工具，以现代数值现代数值分析方法分析方法的为特征的现代轧制理论现代轧制理论得到了迅速发展。1.1.1 1.1.1 现代轧钢生产的发展现代轧钢生产的发展l 20世纪世纪5070年代年代发展趋势是大型化、高速化和连续化发展趋势是大型化、高速化和连续化 1960年前建立的热带钢轧机，辊身范围11202490mm，年生产能力100200万吨，带钢卷重614吨，最大精轧速度为1012m/s，技术进步是将AGC应用于精轧机； 20世纪6070年代，轧机向现代化技术方面发展，同

5、时连铸技术发展成熟。大型连铸坯、步进式加热炉、大型化的粗轧机、7机架精轧机组、AGC、升速轧制、层流冷却技术以及轧制过程计算机控制的全面应用。l 20世纪世纪80年代以后年代以后轧钢生产主要向提高产品质量、降低消轧钢生产主要向提高产品质量、降低消耗、优化轧制过程、开发新钢材和新品种方向发展。耗、优化轧制过程、开发新钢材和新品种方向发展。 板形、厚度及超级钢我国轧钢生产的发展我国轧钢生产的发展 1957年鞍钢第一套2800/1700mm半连续式板带钢轧机 到目前为止，辊身长度在1422mm以上的热轧宽带钢轧机8套、薄板坯连铸连轧带钢轧机10余套。 武钢、本钢1700mm3/4连续式热带钢轧机各一

6、套 宝钢2050mm3/4连续式热带钢轧机 攀钢1450mm半连续式热带钢轧机 太钢1549 mm半连续式热带钢轧机 梅钢1422mm全连续式热带钢轧机 宝钢1580、鞍钢1780mm半连续式热带钢轧机各一套 珠钢1500、邯钢1900和包钢1750薄板坯连铸连轧机各一套。1.1.2 1.1.2 轧制技术的发展轧制技术的发展轧钢生产的发展促进了轧制技术的进步l 连铸技术l 连铸直接轧制技术(CC-DR)l 连铸热装直接轧制技术(CC-HCR)l AGC、AFC、ATCl SFR及无头轧制技术l ISP及CSP薄板坯连铸连轧技术 1.1.3 1.1.3 现代轧制理论研究的基本任务现代轧制理论研究

7、的基本任务求解轧制变形区各种分布量，如应力场、应变场、速度场和温度场等，为板形板厚控制和型钢孔型设计提供理论基础。对轧制过程中工具及工件的温度温度与变形变形进行综合研究，为钢的高精度轧制及轧机的高精度控制服务。对轧件不均匀变形及轧件头尾不稳定变形过程的理论研究，为提高产品质量和成材率、进一步优化轧制规程服务。提高轧制过程参数的理论解析精度，建立和完善建立和完善控制轧制过程的数学模型。开展轧制过程热力学及冶金学参数的综合研究，对轧制过程的变形温度、变形程度、金属的微观组织及产品的最终性能进行综合模拟，实现实现根据产品使用进行钢材成份及轧制过程的预设计。1.2.1 1.2.1 初等理论中的数值方法

8、初等理论中的数值方法 采用有限差分方法求解卡尔曼或奥罗万方程。 基本思想：在变形区内取微元体，建立力平衡微分方程，然后在变形区内进行差分网格划分，在已知边界条件下，采用差分方法求解微元体上的力平衡方程。 特点：能够定性地得出变形区中的轧制力和金属流动规律，但计算精度有待于进一步提高。1.2.2 1.2.2 滑移线理论及其数值解法滑移线理论及其数值解法 滑移线法：把轧制过程变形区划分为一系列由滑移线族组成的滑移线网络，每条滑移线均为达到屈服切应力k，根据Henky应力方程可以确定变形区的应力场。 近年来，利用计算机可以形成金属成型变形区的滑移线网络，并计算相应的滑移线场。 特点：滑移线法只能处理

9、理想刚塑性体平面变形或轴对称变形问题，对三维变形问题、温度和材料性质参数分布不均问题是无能为力的。 1.2.3 1.2.3 能量法及其数值解法能量法及其数值解法 能量法的基础是刚塑性材料的变分原理。 基本思想：给定边界条件设定含有待定参数的运动许可速度场或静力许可应力场建立相应能量泛函使其最小化确定待定参数得到真实的速度场由塑性力学基本关系求出变形及力能参数得到变形区内的应变场。 优点：能量法可以求解三维变形问题，直接得出变形功率、转矩和由速度场决定的宽展、前滑; 缺点：由于不能直接得出静水压力，所以不能直接得出应力分布。此外，能量法也难以处理温度、变形抗力等不均匀分布的问题。 1.2.4 1

10、.2.4 弹塑性有限元法弹塑性有限元法 弹塑性有限元法分析金属成型时采用弹塑性材料本构关系，考虑变形的历史相关性，在求解时需要采用增量加载，在每一个加载步中，只能有少数单元从弹性状态进入塑性状态，以便减小计算误差，因此，所需计算机的容量较大、计算时间长。 优点：不仅可以求解塑性区的扩展、应力、应变分布，而且可以有效地处理卸载问题，计算残余应力、残余应变分布。 缺点：存在积累误差，计算机容量较大，计算时间长。2.2.刚塑性有限元法的基本理论刚塑性有限元法的基本理论 2.1 2.1 有限元法的基本概念有限元法的基本概念 有限元法有限元法：把工件划分成有限结点相连接的单元，以结点上的速度(位移)作为

11、未知量，利用最小能原理求解相应的方程组确定此未知量，按结点速度与单元内部应变以及单元内部应力之间的关系确定各单元的应力、应变分布。2.2 2.2 刚塑性有限元法及基本思想刚塑性有限元法及基本思想 用有限元方法分析金属塑性成型过程时，采用刚塑性用有限元方法分析金属塑性成型过程时，采用刚塑性材料本构模型进行求解材料本构模型进行求解就是刚塑性有限元法。 基本思想基本思想：从刚塑性材料的变分原理变分原理出发，按有限元模式把能耗率泛函能耗率泛函表示为节点速度的非线性函数，利用数学上的最优化理论最优化理论得出满足极值条件的最优解，即使总能耗率取最小值的运动许可速度场运动许可速度场，根据塑性力学的基本关系和

12、本构方程得出应变速度场、应力场以及变形和力能参数。 2.3 2.3 刚塑性材料模型刚塑性材料模型 金属成型过程中，材料变形的物理过程非常复杂，为了便于数学处理，必须做出一些假设，把变形中的某些过程理想化。用刚塑性有限元法分析材料变形问题时，材料满足下列基本假设： (1) 材料均质各向同性； (2) 忽略材料的弹性变形，不计体积力与惯性力； (3) 材料的变形流动服从Levy-Mises流动理论； (4) 材料的体积不变或微可压缩。 2.3.1 2.3.1 理想刚塑性材料模型理想刚塑性材料模型 理想刚塑性材料模型的基本假设如下： (1) 材料均质各向同性； (2) 忽略材料的弹性变形，不计体积力

13、与惯性力； (3) 材料的变形流动服从Levy-Mises流动理论； (4) 材料的体积不变； (5) 不考虑加工硬化，忽略变形抗力对变形速度的敏感性。2.3.1.1 2.3.1.1 理想刚塑性材料模型的应力应变关系理想刚塑性材料模型的应力应变关系图2-1 理想刚塑性材料的应力应变关系2.3.1.2 2.3.1.2 理想刚塑性材料模型的特点理想刚塑性材料模型的特点l 只要等效应力达到一恒定数值，材料便发生屈服，而且材料在整个变形过程中屈服应力不再发生变化。l 采用理想刚塑性材料模型进行能量积分时，可以把等效应力做为常数提到积分号之外，从而使积分过程得到简化。2.3.1.3 2.3.1.3 采用

14、该模型进行采用该模型进行FEMFEM求解应该注意的问题求解应该注意的问题 在轧制变形区中，由于轧件各点的温度、变形速度和变形程度的不同，屈服应力相差很大，在整个变形区内采用理想刚塑性材料模型必然会给计算结果带来误差。因此，用有限元法求解时，把变形区划分成足够多的单元，这样可以认为每个单元内的温度、变形速度和变形程度相同，在每个单元内采用理想刚塑性材料模型，不同单元采用不同的屈服应力，这样处理才能得到比较接近实际的结果。2.3.2 2.3.2 刚塑性硬化材料模型刚塑性硬化材料模型 刚塑性硬化材料基本假设如下： (1) 材料均质各向同性； (2) 忽略材料的弹性变形，不计体积力与惯性力； (3)

15、材料的变形流动服从Levy-Mises流动理论； (4) 材料的体积不变； (5) 考虑加工硬化和变形抗力对变形速度的敏感性。 2.3.2.1 2.3.2.1 刚塑性硬化材料的应力应变关系刚塑性硬化材料的应力应变关系 图2-2 刚塑性硬化材料的应力应变关系2.3.2.2 2.3.2.2 刚塑性硬化材料的变形抗力刚塑性硬化材料的变形抗力 对于刚塑性硬化材料来说，当材料的化学成份和物理状态一定时，通常把变形抗力表示成变形温度、变形速度和变形程度的函数：( , , ) sf T 0nsa()nmsbc(2-1)(2-2)(2-3)2.3.2.3 2.3.2.3 Mises流动法则流动法则 理想刚塑性

16、材料和刚塑性硬化材料都假设材料是不可压缩的，根据Mises流动法则，变形速度分量与偏差应力分量成正比，即 这种材料的变形速度场与偏差应力场一一对应，但由于静水压力是不确定的，所以当以速度为未知量进行求解时，不能直接求得应力场。而体积可压缩材料模型可巧妙地解决这一问题。 yxyyzxzxzxyzxyyzzxd(2-4)2.3.3 2.3.3 刚塑性可压缩材料模型刚塑性可压缩材料模型 刚塑性可压缩材料的基本假设如下： (1) 材料均质各向同性； (2) 忽略材料的弹性变形，不计体积力与惯性力； (3) 材料的变形流动服从Levy-Mises流动理论； (4) 材料的体积微可压缩 ； (5) 考虑加

17、工硬化和变形抗力对变形速度的敏感性。(1) (1) 刚塑性可压缩材料的屈服条件刚塑性可压缩材料的屈服条件 刚塑性可压缩材料与刚塑性硬化材料的主要差别是放松了体积不变条件的约束，即假设屈服与静水压力有关，屈服条件不仅取决于偏差应力的二次不变量，也取决于应力的一次不变量。122222222223()1 =()()()6()6xyzxyzmxyyzzxxyyzzxxyyzzxxyyzzxJJ (2-5)(2-6)刚塑性可压缩材料的屈服条件刚塑性可压缩材料的屈服条件 在主轴条件下，屈服条件：1222221223311()()() 2mg(2-7)122222222122211()()()6()2 39

18、xyyzzxxyyzzxmggJJ(2-8)Misees屈服曲面平面1122 , 23图2-3 刚塑性可压缩材料主应力空间的屈服曲面可压缩参数可压缩参数对对屈服条件的影响屈服条件的影响 在平面上(1 + 2 + 3 =3m=0 ),椭球体与Mises屈服圆柱相切，即刚塑性可压缩材料与Mises屈服条件一致。在平面以外(m0)，刚塑性可压缩材料比理想刚塑性材料容易进入屈服状态。可压缩参数g值越小，椭球体长轴延长，越接近于Mises屈服条件。当考虑加工硬化时，椭球体的体积将随着加工硬化而膨胀。 数量场和矢量场数量场和矢量场l 数量场：对于空间或部分空间的任意一点M，都有一个确定的数值f(M)与之对

19、应，则称在这个空间或部分空间上确定了一个数量场。该数量场可用数值函数f(M)来确定。例如温度场就是一个数量场。l 矢量场：对于空间或部分空间的任意一点M，都有一个确定的矢量f(M)与之对应，则称在这个空间或部分空间上确定了一个矢量场。该矢量场可用矢量函数f(M)来确定。例如速度场、应力场均是矢量场。 势和势场势和势场l 梯度：设函数u = f ( x, y, z )在空间或部分空间具有一阶连续偏导数，则对空间或部分空间的每一个点P ( x, y, z )都可以确定一个矢量： 这个矢量就称作函数u = f ( x, y, z )在点P ( x, y, z )的梯度。l 利用场的概念，该矢量函数在

20、空间或部分空间确定了一个矢量场，即梯度场。它是由数值函数产生的，称数值函数为这个矢量场的势，矢量函数则称为有势场或势场。 kzfjyfixfzyxf),( grad(2) (2) 塑性势和变形速度分量塑性势和变形速度分量 假设刚塑性可压缩材料的屈服函数为塑性势 ：1222222221()()()()6()2ijxyyzzxxyyzzxmFg(2-9)()ijijijFdd根据塑性势的定义：()ijijijFd则(2-10)dtdddt 两边同除以时间增量，令变形速度分量变形速度分量 利用(2-9)式和(2-10)式可直接求出应变速度分量：2(1)92(1)92(1)9222xxmyymzzmx

21、yxyyzyzzxzxgdgdgdddd(2-11) 与与 之间关系之间关系 下面以 为例推导 与 之间关系dd()ijxxFddd1222222221()()()()6()2ijxyyzzxxyzmFg1222222222()()()312()()()6()22223()33 2 ()2 ()3 xyzxmijxxyyzzxxyyzzxmxyzmxxyzmijijxgFgggFF由于2233(1)392 ()2 ()mmxmijijggFF 与与 之间关系之间关系 可得：(2-12)32 ()ijddF dd的简化形式的简化形式 根据单位体积的塑性变形功率 ：dxxyyzzxyxyyzyzz

22、xzxdw (2-13)把(2-11)式代入(2-13)式整理可得：212921292129xxmxyxyyymyzyzzzmzxzxgddgddgdd所以 32d (2-14)222222222222222222222222()1()92 2()3192 2()331 333xyzxyyzzxmxyxxyzxyyzzxmxyzxyyzzxmmxydgdgdgd 2222222222222222222222236()921 333()6()231 2222226()231 ()(3zxyyzzxmmxyzxyzxyyzzxmxyzxyyzzxxyyzzxmxygdgdgd 2222222)()

23、6()22 3yzzxxyyzzxmgd(3) (3) 应力应变速度关系应力应变速度关系 把(2-14)式代入(2-11)式21921392219222xxmyymvzzmvmmxyxyyzyzzxzxgdgddgdggddd 刚塑性可压缩材料刚塑性可压缩材料 与与 g 和和 关系关系 22333vxyzxyzmmg dgdg(2-15)当 一定时，g 越小，体积变形速度小；当g = 0时，此时刚塑性可压缩材料的屈服条件变为Mises屈服条件，应力应变关系变为Levy-Mises流动法则，刚塑性可压缩材料模型变为不可压缩的理想刚塑性或刚塑性硬化材料模型。从这个意义来讲，刚塑性不可压缩材料可以看

24、作刚塑性可压缩材料在g=0时的特例，因此更具有普遍性。 根据(2-11)和(2-14)式可得体积变形速率与可压缩参数和静水压力关系 ：mV1222222221()()()6()2xyyzzxxyyzzxmgm(3) (3) 应力应变速度关系应力应变速度关系 2213922 13922 139212 39xxmvxvxxvgggggg0.5/2 3xyxyxyd如果采用张量形式，并引入克罗内克尔符号，可得刚塑性可压缩材料的应力应变速度关系：212()39ijijijvg(2-16)(4) (4) 等效应变速度等效应变速度 把(2-16) 代入(2-13)式：212()39ijijijvg212

25、=39212 39212 3922 33xxyyzzxyxyyzyzzxzxxxvyyvzzvxyxyyzyxggg 23zxzx 12222222221()()()6()9xyyzzxxyyzzxvg(2-17) 这样，利用(2-16)式就可以从变形速度分量中直接求出应力场。因此，当从运动许可速度场出发以节点速度为未知量求解塑性变形过程时，能够简便地求出变形和力能参数。2222222222222222222124()()()()3932124 ()()3932 3()6(9xyzvxyzvxyzxyyzzxxyzVVxyyzzxxyzVxygg 222222222221)21 ()()()6

26、()9yzzxVxyyzzxxyyzzxVgg(5) (5) 刚塑性可压缩材料的流动法则刚塑性可压缩材料的流动法则 根据(2-15)式可得： 将上式代入(2-11)式，可得：v239mmgd(2-18)219221993 2xmmxmxmxmxmmxmgdgdgdd(5) (5) 刚塑性可压缩材料的流动法则刚塑性可压缩材料的流动法则 上式表明，刚塑性可压缩材料的流动法则是偏差应变速度分量与偏差应力分量成正比。与Levy-Mises流动法则一致，只不过是体积变形速率不等于零。32ymxyyzxmzmzxxmymzmxyyzzxd(2-19)(6) (6) 刚塑性可压缩材料的特点刚塑性可压缩材料的

27、特点 刚塑性可压缩材料与刚塑性硬化材料的主要差别是，放松了体积不变条件的约束，即假设屈服与静水压力有关。体积变化率取决于静水压力，当求出材料的屈服应力、等效应变速率和给定材料的体积可压缩参数后，可以直接从速度场求得应力场，所得结果与体积不可压缩条件的解十分接近，而且计算过程得到简化。(7) (7) 可压缩参数可压缩参数g g对体积变形的影响对体积变形的影响 将式(2-15)两端乘以时间增量dt，可得：vmg(2-20)以单向压缩应力状态为例，其应力分量及等效应力为：0 0 0 0 03xyxyyzzxzm (2-21)9/39zmgg所以体积变形为：(2-23)(2-22)9vgg(7)(7)

28、 可压缩参数可压缩参数g对体积变形程度的影响对体积变形程度的影响 利用上式计算不同等效应变下，可压缩参数g对体积变形程度的影响 表2-1 可压缩参数g对体积变形程度的影响 g10.050.010.0010.00010.13.160.1660.03330.003330.000330.26.320.3320.06660.006670.000670.39.490.4990.09990.010000.00100 从上表可见，在1030%的等效应变条件下，可压缩参数g取0.010.0001之间时，材料的体积变化不超过0.1%，可以近似满足体积不可压缩条件。因此，求解金属材料的轧制过程，可压缩参数g可取0

29、.010.0001。2.4 2.4 刚塑性材料的基本方程刚塑性材料的基本方程 刚塑性材料发生塑性变形时，由表面积 S 所围成的体积 V 中，应力 、速度 和应变速率 应满足下列基本方程: 1.力平衡微分方程(2-24)ijivij,0000 xyxxzyxyyzij jzyzxzxyzxyzxyz2.几何方程(2-25),1211()2212yxxxxyyyziji jj iyyzxzzzzxvvvxyxvvvvvxzyvvvxxz2.4 2.4 刚塑性材料的基本方程刚塑性材料的基本方程2.4 2.4 刚塑性材料的基本方程刚塑性材料的基本方程3.本构关系(2-26)33223332223322

30、xxmxyxyijijyymyzyzzzmzxzxs2.4 2.4 刚塑性材料的基本方程刚塑性材料的基本方程 4.屈服准则 5.体积不可压缩条件 6.边界条件上在给定速度表面viiSvv ijjipnpS在给定外力表面上(2-27)(2-28)(2-29)(2-30)223ijijs skk0vijijxyz 2.4 2.4 刚塑性材料的基本方程刚塑性材料的基本方程式中： 、 等效应变速率和等效应力； 偏差应力； 屈服剪切应力； 外力表面单位法线矢量的方向余弦。 对于刚塑性材料发生塑性变形时，需要对上述基本方程进行联立求解来确定变形区的速度、应变速率和应力。由于直接求解这些偏微分方程组非常困难

31、，因此，人们寻求其它求解途径，即利用变分原理求解。 ijskjn2.5 刚塑性有限元的基本原理刚塑性有限元的基本原理l 刚塑性材料的变分原理是刚塑性有限元法的理论基础，变分原理通过能量积分把偏微分方程组的求解问题变成了求泛函极值问题，从而为各种实际问题的求解提供了一种新方法。l 材料模型不同，变分原理的形式也不相同。根据附加条件的情况，变分原理可分为一般变分原理、不完全广义变分原理和完全广义变分原理。2.5.1 Mapkov原理原理 (第一变分原理第一变分原理) 设刚塑性变形体的体积为V，表面积为S0。在Sp上给定表面力 ，在Sv上给定速度 vi，在满足变形的几何条件、速度边界条件和体积不变条

32、件的一切许可速度场中，真实的速度场使泛函 取极小值。ip23PsijijiiVSdVp v ds (2-31)2.5.2 刚塑性可压缩材料的变分原理刚塑性可压缩材料的变分原理 在满足变形的几何条件、速度边界条件和体积微可压缩条件的一切许可速度场 中，真实的速度场使泛函 取极小值。 ivpiiVSdVp v ds(2-32)2.5.3 刚塑性材料的广义变分原理刚塑性材料的广义变分原理 在第一变分原理中，运动许可速度场必须满足几何条件、速度边界条件和体积不可压缩条件。在处理实际问题时，有些条件比较容易满足，而有的条件不易满足，因此引用拉格朗日乘子把这些约束条件全部或部分引入总体能量泛函中，通过泛函

33、变分使这些约束条件得到满足。引入约束条件后，变分原理的表达式要发生变化，统称为广义变分原理。引入部分约束条件的称为不完全广义变分原理，全部约束条件同时引入的称为完全广义变分原理。 (1) 不完全广义变分原理不完全广义变分原理 设定初始速度场时，几何条件和速度边界条件容易满足，而体积不可压缩条件则不易满足，所以把体积不变条件通过拉格朗日乘子引入能量泛函中，得到新能量泛函： 在满足几何方程和速度边界条件的一切速度场 中，真实的速度场使上述能量泛函取驻值。 ivpiiijijVSVdVp v dsdV (2-33)(2) 完全广义变分原理完全广义变分原理 在一切位移速度、应变速度和应力函数中，使能量

34、泛函 取驻值的vi、 和 必为刚塑性材料的正确解。ijij, ()1 ()2pviiijijiiiVSVSijiji jj iVdVp v dsdVvv dsvvdV(2-34)第一变分原理与广义变分原理的差别第一变分原理与广义变分原理的差别 第一变分原理要求能量泛函取最小值，而广义变分原理仅要求能量泛函取驻值； 第一变分原理要求速度场满足运动许可条件，静力许可条件是通过变分过程近似满足的，而广义变分原理的速度场不受任何约束，所有方程均由变分过程近似满足； 第一变分原理比广义变分原理所得到结果精度高，但前者初始速度场的设定比后者难。 2.6 速率敏感材料的总体能量泛函速率敏感材料的总体能量泛函

35、 求解轧制问题时，总体能量泛函的表达式为： 式中： 塑性变形功率； 接触表面摩擦功率； 外张力功率，前张力取负号，后张力取正号； 速度不连续面上的剪切功率。pfs(2-35)pfS2.6.1 内部塑性变形功率内部塑性变形功率 对于速率敏感性材料或在高温下成型的金属 ,单位体积的塑性变形功率： 速率敏感材料的内部塑性变形功率为： 001() 1nmijEdadm (2-37)1 1pVdVm (2-38)2.6.2 接触表面上的摩擦功率接触表面上的摩擦功率 轧制过程中，轧辊与轧件接触表面上存在中性面，为了避免因摩擦力变向出现第二类奇异点而导致能量泛函不收敛的问题，可采用摩擦应力模型： 摩擦应力是

36、相对滑动速度vg的函数，因此摩擦功率为： 1223gsfgvmvk(2-38)2210()()3gsffggvssvmddsvkk ds (2-39)2.6.3 外张力功率外张力功率 式中： 作用有外张力的表面； 张应力，轧件受拉时取正值，受推时取负值； 相应表面处的位移速度。isv ds(2-40)siv2.6.4 速度不连续面上的剪切功率速度不连续面上的剪切功率 式中： 速度不连续面； 屈服剪应力； 速度不连续面上的速度不连续量。(2-41)ksssvkssssv ds速率敏感材料轧制过程总体能量泛函速率敏感材料轧制过程总体能量泛函 在轧制过程中，边界上的外力为接触表面上的摩擦应力、变形区

37、端部的张应力和速度不连续面上的剪切应力，相应的总体能量泛函为： 对于简单轧制过程来说，总体能量泛函可简化为： 221 1 ()13sgiksVssssmdVvkk dsvdsv dsm(2-42)221 1 ()13sgVsmdVvkk dsm(2-43)2.7 刚塑性有限元的求解途径刚塑性有限元的求解途径l设定运动许可速度场。l建立总体能量泛函，并把泛函表示成速度的函数。l利用数学上的极值理论求解泛函的驻值或最小值。l利用几何方程和本构方程确定应变速度场和应力场。l通过接触表面应力积分求出总轧制力和平均单位压力，利用速度场求出轧件宽展、前滑和轧件侧面形状参数。2.7.1 轧制变形区的有限元离

38、散化轧制变形区的有限元离散化 (1) (1) 研究目标的选择研究目标的选择zox根据变形特征确定研究目标。利用对称面上一个方向速度为零的特点，可使问题简化。如果变形过程的几何条件、物理条件在变形区内部及边界上对称，那么其真实解也必然对称。图2-2 板带轧制过程的对称面 (2) 单元类型的选取单元类型的选取 研究目标确定之后，便可着手在所选区域内设置节点、划分单元。除了速度已知边界之外，每个节点都有一组速度未知量，对于2维问题节点速度为 ,3维问题为 ，所以设置的节点越多，求解时未知数越多。单元类型的选择，可根据所处理问题的特点来确定，求解平面变形问题，常用四边形单元，求解3维变形问题常用立方体

39、单元。 yxvv 、zyxvvv、2.7.2 总能量泛函的离散化总能量泛函的离散化 把轧制变形区划分为有限个单元之后，总能量泛函便可由单元能量泛函迭加求和来得到。下面以简单轧制过程为例进行说明。(1) 单元能量泛函单元能量泛函 对于平面变形条件下的简单轧制过程，设单元的面积为Se，外力已知边界为Le，则单元能量泛函为： 当单元的速度插值函数设定之后，单元内的 便可确定。在给定材料屈服应力和摩擦边界条件下，单元能量泛函实际上是单元节点速度的函数： ee221 1 ()13esgSLmdSvkk dLm(2-44)ij 、e112244(,)xyxyxyf vvvvvv(2-45)(2) 总能量泛

40、函总能量泛函 根据(2-44)式对所研究区域全部单元的能量泛函求和，便可得到总能量泛函：式中：nm单元总数； nl 接触表面单元数目。ee221 1111 ()13enmnmnlsgSLeeemdSvkk dLm(2-46)(3) 接触表面的速度边界条件接触表面的速度边界条件 轧制过程接触表面上的节点必须满足速度边界条件： tgyxvv (2-47)图2-3 轧制特征角(4) 总体能量泛函总体能量泛函 由于接触表面节点的两个速度分量只有一个是未知量 ，同时对于 x 轴的节点垂直方向的速度已知(为零)，除去节点已知速度分量之后，把节点未知速度分量统一用 表示，总体能量泛函表示成未知速度分量的函数

41、：式中：nx未知量总数。xviv12( ,)inxf v vvv(2-48)2.7.3 总能量泛函的最小化总能量泛函的最小化 用刚塑性有限元法求解时，根据变分原理应求总体能量泛函的极小值。从式(2-48)可知，总体能量泛函是节点速度矢量的多元函数。根据多元函数求极值的条件，则有： 由于变分是任意的，必须有下列方程组成立：1 0 enxeeevvvv(2-49)10 enxeevv(2-50)(1) (1) 线性化处理线性化处理 采用Newton-Raphson方法线性化求解时，假定泛函连续并存在各阶导数，则借助于泰勒级数在任一点展开。在第次迭代时，泛函可在第次迭代求得的值上展开，略去二阶以上的

42、高阶微量，从而得到以为未知量的线性方程组：211111 0 eenmnmnmeeeeTeeekkkkvvvvv211111 eenmnmnmeeeTeeeekkkvvvv (2-51)(2-52)(2) (2) 阻尼因子阻尼因子 采用(2-52)式求解时，先给定一个初始速度场 ，然后用迭代法按上式求出 ，直到相邻两次迭代速度近似解的偏差充分小时为止。迭代时，为了防止发散一般采用公式： vv1 kkvvv(2-53)01式中： 阻尼因子， 。2.8 刚塑性有限元的基本公式刚塑性有限元的基本公式 从上节可以看出，刚塑性有限元求解过程的一个基本步骤是计算单元的能量泛函计算单元的能量泛函。为此，需要解

43、决以下几个方面问题：l设定初始速度场；l设定单元的速度分布函数插值函数；l建立单元变形速度与节点速度的关系；l计算单元能量泛函；l求单元能量泛函的一阶偏导数即梯度；l求单元能量泛函的二阶偏导数即Hessian矩阵。 2.8.1 单元坐标及速度插值、形状函数单元坐标及速度插值、形状函数 单元内任意一点的速度及坐标可以用单元节点的速度及坐标进行插值计算。对速度及坐标的插值使用相同的线性插值函数，该函数只与单元的形状和节点的配置有关，因此，称之为形状函数或形函数，与此相应的单元称为线性等参单元。(1)(1) 四边形线形等参单元四边形线形等参单元 为了便于插值与积分，单元内部采用局部坐标系 ，整体坐标

44、系的任意四边形单元映射到局部坐标系后就变为正方形单元，单元节点编号及坐标系的映射变换如图2-4所示。 42y13xO3124O( , ) 图2-4 单元节点编号及坐标系的映射变换(2)(2) 单元局部坐标与整体坐标的变换单元局部坐标与整体坐标的变换 整体坐标系中4个节点坐标为 ，局部节点坐标 如下： 单元局部坐标 与整体坐标 的变换关系：( , ) , 14kkxyk (, )kkint(1)/ 2)12341231-1-1 11( 1)-1 1-11( 1)kkkk (2-54),(yx( , ) (2-56)41 122334414112233441( , )( , )kkkkkkxN x

45、N xN xN xNxyN yN yN yN yNy (3)(3) 四边形等参单元的形函数及其性质四边形等参单元的形函数及其性质 式(2-56)中 称为单元的形函数，具有以下性质：l形函数 是与坐标插值函数相同形式的双线性插值函数；l形函数 在节点上， ；l形函数 在节点上， 。 ( , )kN ( , )kN ( , )kN ( , )kN (,)1kkkN (,)0 , kllNlk (4)(4)形形函数表达式的确定函数表达式的确定 以 为例说明其确定过程： 1( , )N 42y13xO3124O上图局部坐标系中，单元24、34两条边方程可表示为： 1-01-0(2-57)(4)(4)形

46、形函数表达式的确定函数表达式的确定 根据形函数的性质(3)， 在节点2、3、4处的值为零，而直线24、34经过这些节点，所以可取： 根据形函数的性质(2)， 在节点1处其值为1，把1点的局部坐标值代入上式得：1( , )N 1( , )(1)(1)NC (2-58)1( , )N 11(1 1)(1 1)4CC(2-59)11111( , )(1)(1)(1)(1)44N (2-60)(5)(5)四边形等参单元的形函数四边形等参单元的形函数 22233344411( , )(1)(1)(1)(1)4411( , )(1)(1)(1)(1)4411( , )(1)(1)(1)(1)44NNN 1

47、( , )(1)(1) (1 4)4kkkNk (2-61)(2-62)(6)(6) 单元内任一点的速度与节点速度的关系单元内任一点的速度与节点速度的关系 4141( , )( , )xkxkkykykkvNvvNv (2-63) evNv(2-64) xyvvv 4321432100000000NNNNNNNNN Tyxyxyxyxevvvvvvvvv443322112.8.2 单元变形速度与节点速度关系、单元变形速度与节点速度关系、B矩阵矩阵 (1) 平面变形的几何方程平面变形的几何方程xxyyyxxyvxvyvvyx(2-68)(2) 以节点速度表示的几何方程 4414141411( ,

48、 )( , )()ixxxixixkxkykiyyyiiykykkyiixxyyixixyiNvvxxvNvvNvyyvNvvNNvvvxyyx (2-69)(3) 几何方程的矩阵形式 11212342123431122334434400000000 xyxxyyxxyyxyvvvbbbbvccccvcbcbcbcbvvv (2-70)4) , 3 , 2 , 1( , , iyNcxNbiiii式中：(3) 几何方程的矩阵形式 Txyxy4321BBBBB iiiiibccbB00 eBv令 (2-74)(2-71)(2-73)(2-75)BiB式中： 单元应变矩阵或B矩阵； 典型子矩阵。2

49、.8.3 Jacobi矩阵及其逆矩阵和行列式矩阵及其逆矩阵和行列式 单元应变B矩阵中的元素bi及ci是形函数 对x及y的偏导数，为了便于积分计算，形函数 对x及y的偏导数应该用 对 及 的偏导数来表示。根据复合函数的求导法则： ( , )iN ( , )iN ( , )iN TiiiiiTiiiiiNNNNNxyxyxyxyNNNNNxyxyxyxy(2-76)(1) 2维问题的维问题的Jacobi矩阵矩阵 将上式写成矩阵形式：iiiiNxyNxNNxyy(2-77) xyJxy(2-78)(1) 2维问题的维问题的Jacobi矩阵矩阵 1( , )(1)(1)4( , )1(1) 4 (1

50、4)( , )1(1)4kkkiiiiiiNNiN (2-79)444111444111( , ) ( , )iiiiiiiiiiiiiiiiiiNNxyxyxNxJxyNNyNyxy (2-80)(2) Jacobi矩阵的逆矩阵及行列式矩阵的逆矩阵及行列式 1iiiiiiiiNNxyNNxxJNNNxyNyy(2-81)1 1 yxJyxJ(2-82)xyxyJ(2-83)2.8.4 高斯积分高斯积分 对于2维变形问题，建立单元内任意一点变形速度与节点速度关系之后，在给定初始速度场条件下，即可对单元能量泛函进行积分计算。在单元局部坐标系中，单元能量泛函式可表示成： 刚塑性有限元中的数值积分常

51、用高斯积高斯积分，即按照数学上的规则在单元内选取若干个积分点，用积分点处的函在单元内选取若干个积分点，用积分点处的函数值与求积系数之积的累加结果近似代替原积分数值与求积系数之积的累加结果近似代替原积分。 1112211111 ()13esgfmJ d dvkk Jdm (2-84)(1)单元塑性变形功率单元塑性变形功率高斯积分表示高斯积分表示 式中：n单元积分点个数，Hi、Hj求积系数。 41111111111 111nnepijijiJ d dH HJJmmm (2-85)G13124OG3 G2G4 对于2维问题采用线性单元时，n=2, Hi=Hj =1，此时单元内有4个高斯积分点，积分点

52、坐标为0.57735027，如图2-5所示。(2)单元摩擦功率的积分单元摩擦功率的积分 单元摩擦功率的积分，可采用1维高斯积分近似计算：12211221122211()3 ()3 =()3esfgfnsigfisgfimvkk JdmHvkk Jmvkk J(2-86) 对于1维问题采用线性单元时，n=2, Hi=1，此时单元内有2个高斯积分点，积分点坐标为0.57735027。12OG1G2(3) 1维问题的维问题的Jacobi矩阵行列式矩阵行列式 式中：x1、x2为摩擦表面单元两个节点的整体坐标值； N1、N2为1维单元的两个形函数。212121ifiiNNNxJxxx(2-87)121(

53、1)21(1)2NN(2-88)(4) 单元的能量泛函单元的能量泛函 采用高斯积分后，2维变形简单轧制过程的单元能量泛函为： 式中： 分别为单元高斯积分点处的等效应力、等效应变及Jacobi 矩阵行列式的值。42221111()13esgfiimJvkk JmfJJ、 、 和(2-90)2.8.5 矩阵分析中的公式矩阵分析中的公式(1) 数值函数对向量变量的偏导数数值函数对向量变量的偏导数TTTTA XAXX AAX 12TnAaaa12TxXxxx()TTf XA XX A设 是给定的向量， 是向量变量，且 ，则有 (2-a)(1) 数值函数对向量变量的偏导数数值函数对向量变量的偏导数 ()

54、Tf XA XAXXijn nAaT12xXxxx()Tf XX AX设 是给定的矩阵， 是向量变量，且 ，则有() 2Tf XAAAXX如果 ，则 (2-b)(2-c)(2) 向量函数对向量变量的偏导数向量函数对向量变量的偏导数 ijn nAaT12nXxxx()f XAX设 是给定的矩阵， 是向量变量，且 ，则有()()TTTf XAXf XAX(2-d)2.8.6 能量泛函的一阶偏导数能量泛函的一阶偏导数(梯度梯度)(1) 平面变形问题单元的等效应变速度列阵平面变形问题单元的等效应变速度列阵 平面变形问题单元的变形速度列阵 TexyxyB v222411931 xxyyxyvTTTgZC

55、Cg (2-91)(2-92)式中的常数矩阵和向量式中的常数矩阵和向量 4209924 0991003Z 1 1 0C (2-94)(2-93) TZ 222224209924 0991003422419999342241419999393xTxyxyyxyxxyxyxyyxyxxyxyyxyxxyyZ 2xy 1 TTCCg 21 1111 10011TTxxyxyyxyxyxyVCCgggg (1) 2维平面变形条件下的等效应变维平面变形条件下的等效应变 222 4111() 93TexyxyTTxxyyxyvB vZgg 1 TeTTeeTTTeCCvBZ B vvBCC B vg(2-

56、95)(1) 2维平面变形条件下的等效应变维平面变形条件下的等效应变 上式可以进一步写成更为简洁的形式： TeTTeZvBZB v412109912141 0991003TggZZCCggg (2-96)(2-97)等效应变不用等效应变不用(2-96)式简洁形式？式简洁形式？ 因为采用体积可压缩方法进行求解时，如果等效应变采用(2-96)式简洁形式，则要求单元内每一个高斯积分点处的体积变形速度都很小，这样的约束条件过于苛刻，经常出现迭代计算发散现象，导致有限元数值求解过程无法进行。(1) 2维平面变形条件下的等效应变维平面变形条件下的等效应变 为了排除每一个高斯积分点处体积变形速度都很小这一严

57、格约束，单元的体积变形速度取各高斯积分点处体积变形速度的平均值，即：444111111 444eeeVkkkkkkCC BvCBvC B v1 eTTeeTTTevBZ B vvBCC B vg(2-98)(2-100)411 4kkBB(2-99)(2) 等效应变速度对速度向量的一阶偏导数等效应变速度对速度向量的一阶偏导数 等效应变速度(2-100)式对速度向量的一阶偏导数：1 TeTTeePeBZBvBCCBvKvgv1 TTTPKBZBBCC Bg(2-101)(2-102) PK式中：8 8的矩阵。(3) 单元梯度单元梯度 把变形抗力模型代入单元塑性变形功率表达式，并对单元节点速度向量

58、求一阶偏导数：44111144114411111111 nmenmPiiemnnmPeeeiiePeiiaJaJmmaJaJmvvvJJKvv(2-103)(4) 单元摩擦功率对速度向量的一阶偏导数单元摩擦功率对速度向量的一阶偏导数 接触表面任意一点的速度边界条件： vxvytgyxvv 接触表面任意一点的相对滑动速度：secgRxvvv(2-104)单元内任意一点的速度与节点速度的关系单元内任意一点的速度与节点速度的关系 对于2维轧制过程的接触表面单元而言，单元节点上只有一个未知量，因此，单元内任意一点的速度与节点速度的关系如下： efexffvvNv120.5(1- ) 0.5(1) fT

59、efxxNvvv(2-105)(2-106) 对节点速度向量的一阶偏导数对节点速度向量的一阶偏导数gv 把(2-105)式代入(2-104)式可得相对滑动速度： secegRffvvNv(2-107)sec gTfefvNv 2222()22 secgggeeffggTfefvkkvvvvvkvNv (2-108)(2-109)单元能量泛函的一阶偏导数即梯度单元能量泛函的一阶偏导数即梯度 接触表面单元摩擦功率的一阶偏导数： 单元能量泛函的一阶偏导数即梯度：22221111() sec 33egfTssfffeeiiffvkkmmJJNvv 42111 sec3eefePefeeTsPffiiv

60、vmJKvJN(2-110)(2-111)2.8.7 能量泛函的二阶偏导数能量泛函的二阶偏导数(Hessian矩阵矩阵) 1223 (1) (1) (1) nmeTeTnmeTePeTeTePavvmavKvmvvmKv 2.8.7 能量泛函的二阶偏导数能量泛函的二阶偏导数(Hessian矩阵矩阵)412421324221 (1) (1) eePPeieePPPeeTeTiTePeTTeeePPPPeeTiJKvvJKJKvvvvmKvvKvKvJKmvv 由于 (2-112)2.8.8 总体能量泛函能量泛函的梯度矩阵和总体能量泛函能量泛函的梯度矩阵和Hessian矩阵矩阵 根据迭加原理对所有