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

刚塑性有限元法及其在轧制中的应用,轧制技术及连轧自动化国家重点实验室,1.学习目的和要求,了解现代轧钢生产和轧制技术的发展概况；了解现代轧制理论研究的基本任务；掌握刚塑性有限元的基本概念；掌握刚塑性有限元的基本理论；掌握刚塑性有限元的基本方法；,2.学习的主要内容,刚塑性有限元的基本概念和基本理论；刚塑性有限元相关技术问题的处理方法；求解轧制过程的刚塑性有限元程序。,3.本课程的基础和相关知识,现代塑性加工力学基本方程、变分原理、有限元基础知识；工程数学矩阵分析、优化方法、数值分析；计算机基础知识操作系统、FORTRAN语言和FORTRAN4.0编程软件。,4.讲课和学习方法,课堂讲授基本概念、基本理论、基本方法程序剖析；课外自学消化理解、阅读程序；上机实践调试程序,1.绪论,1.1现代轧制理论研究的发展概况1.1.1现代轧钢生产的发展1.1.2轧制技术的发展1.1.3现代轧制理论研究的基本任务1.2轧制理论数值方法1.2.1初等理论中的数值方法1.2.2滑移线理论及其数值解法1.2.3能量法及其数值解法1.2.4弹塑性有限元法,1.1现代轧制理论研究的发展概况,20世纪60年代前，轧钢生产过程手工操作和使用单体设备。轧制理论主要解决问题轧制力、力矩、功率、宽展和前滑等参数的近似计算。主要进展提出卡尔曼和奥罗万方程，采用一些假设条件推导出轧制力和宽展等公式，逐步形成了以工程法为核心的传统轧制理论体系。20世纪60年代以后，随着轧钢生产和轧制技术的飞跃发展和用户对产品质量要求的日益提高，以计算机为工具，以现代数值分析方法的为特征的现代轧制理论得到了迅速发展。,1.1.1现代轧钢生产的发展,20世纪5070年代发展趋势是大型化、高速化和连续化1960年前建立的热带钢轧机，辊身范围11202490mm，年生产能力100200万吨，带钢卷重614吨，最大精轧速度为1012m/s，技术进步是将AGC应用于精轧机；20世纪6070年代，轧机向现代化技术方面发展，同时连铸技术发展成熟。大型连铸坯、步进式加热炉、大型化的粗轧机、7机架精轧机组、AGC、升速轧制、层流冷却技术以及轧制过程计算机控制的全面应用。20世纪80年代以后轧钢生产主要向提高产品质量、降低消耗、优化轧制过程、开发新钢材和新品种方向发展。板形、厚度及超级钢,我国轧钢生产的发展,1957年鞍钢第一套2800/1700mm半连续式板带钢轧机到目前为止，辊身长度在1422mm以上的热轧宽带钢轧机8套、薄板坯连铸连轧带钢轧机10余套。武钢、本钢1700mm3/4连续式热带钢轧机各一套宝钢2050mm3/4连续式热带钢轧机攀钢1450mm半连续式热带钢轧机太钢1549mm半连续式热带钢轧机梅钢1422mm全连续式热带钢轧机宝钢1580、鞍钢1780mm半连续式热带钢轧机各一套珠钢1500、邯钢1900和包钢1750薄板坯连铸连轧机各一套。,1.1.2轧制技术的发展,轧钢生产的发展促进了轧制技术的进步连铸技术连铸直接轧制技术(CC-DR)连铸热装直接轧制技术(CC-HCR)AGC、AFC、ATCSFR及无头轧制技术ISP及CSP薄板坯连铸连轧技术,1.1.3现代轧制理论研究的基本任务,求解轧制变形区各种分布量，如应力场、应变场、速度场和温度场等，为板形板厚控制和型钢孔型设计提供理论基础。对轧制过程中工具及工件的温度与变形进行综合研究，为钢的高精度轧制及轧机的高精度控制服务。对轧件不均匀变形及轧件头尾不稳定变形过程的理论研究，为提高产品质量和成材率、进一步优化轧制规程服务。提高轧制过程参数的理论解析精度，建立和完善控制轧制过程的数学模型。开展轧制过程热力学及冶金学参数的综合研究，对轧制过程的变形温度、变形程度、金属的微观组织及产品的最终性能进行综合模拟，实现根据产品使用进行钢材成份及轧制过程的预设计。,1.2.1初等理论中的数值方法,采用有限差分方法求解卡尔曼或奥罗万方程。基本思想：在变形区内取微元体，建立力平衡微分方程，然后在变形区内进行差分网格划分，在已知边界条件下，采用差分方法求解微元体上的力平衡方程。特点：能够定性地得出变形区中的轧制力和金属流动规律，但计算精度有待于进一步提高。,1.2.2滑移线理论及其数值解法,滑移线法：把轧制过程变形区划分为一系列由滑移线族组成的滑移线网络，每条滑移线均为达到屈服切应力k，根据Henky应力方程可以确定变形区的应力场。近年来，利用计算机可以形成金属成型变形区的滑移线网络，并计算相应的滑移线场。特点：滑移线法只能处理理想刚塑性体平面变形或轴对称变形问题，对三维变形问题、温度和材料性质参数分布不均问题是无能为力的。,1.2.3能量法及其数值解法,能量法的基础是刚塑性材料的变分原理。基本思想：给定边界条件设定含有待定参数的运动许可速度场或静力许可应力场建立相应能量泛函使其最小化确定待定参数得到真实的速度场由塑性力学基本关系求出变形及力能参数得到变形区内的应变场。优点：能量法可以求解三维变形问题，直接得出变形功率、转矩和由速度场决定的宽展、前滑;缺点：由于不能直接得出静水压力，所以不能直接得出应力分布。此外，能量法也难以处理温度、变形抗力等不均匀分布的问题。,1.2.4弹塑性有限元法,弹塑性有限元法分析金属成型时采用弹塑性材料本构关系，考虑变形的历史相关性，在求解时需要采用增量加载，在每一个加载步中，只能有少数单元从弹性状态进入塑性状态，以便减小计算误差，因此，所需计算机的容量较大、计算时间长。优点：不仅可以求解塑性区的扩展、应力、应变分布，而且可以有效地处理卸载问题，计算残余应力、残余应变分布。缺点：存在积累误差，计算机容量较大，计算时间长。,2.刚塑性有限元法的基本理论,2.1有限元法的基本概念有限元法：把工件划分成有限结点相连接的单元，以结点上的速度(位移)作为未知量，利用最小能原理求解相应的方程组确定此未知量，按结点速度与单元内部应变以及单元内部应力之间的关系确定各单元的应力、应变分布。,2.2刚塑性有限元法及基本思想,用有限元方法分析金属塑性成型过程时，采用刚塑性材料本构模型进行求解就是刚塑性有限元法。基本思想：从刚塑性材料的变分原理出发，按有限元模式把能耗率泛函表示为节点速度的非线性函数，利用数学上的最优化理论得出满足极值条件的最优解，即使总能耗率取最小值的运动许可速度场，根据塑性力学的基本关系和本构方程得出应变速度场、应力场以及变形和力能参数。,2.3刚塑性材料模型,金属成型过程中，材料变形的物理过程非常复杂，为了便于数学处理，必须做出一些假设，把变形中的某些过程理想化。用刚塑性有限元法分析材料变形问题时，材料满足下列基本假设：(1)材料均质各向同性；(2)忽略材料的弹性变形，不计体积力与惯性力；(3)材料的变形流动服从Levy-Mises流动理论；(4)材料的体积不变或微可压缩。,2.3.1理想刚塑性材料模型,理想刚塑性材料模型的基本假设如下：(1)材料均质各向同性；(2)忽略材料的弹性变形，不计体积力与惯性力；(3)材料的变形流动服从Levy-Mises流动理论；(4)材料的体积不变；(5)不考虑加工硬化，忽略变形抗力对变形速度的敏感性。,2.3.1.1理想刚塑性材料模型的应力应变关系,图2-1理想刚塑性材料的应力应变关系,2.3.1.2理想刚塑性材料模型的特点,只要等效应力达到一恒定数值，材料便发生屈服，而且材料在整个变形过程中屈服应力不再发生变化。采用理想刚塑性材料模型进行能量积分时，可以把等效应力做为常数提到积分号之外，从而使积分过程得到简化。,2.3.1.3采用该模型进行FEM求解应该注意的问题,在轧制变形区中，由于轧件各点的温度、变形速度和变形程度的不同，屈服应力相差很大，在整个变形区内采用理想刚塑性材料模型必然会给计算结果带来误差。因此，用有限元法求解时，把变形区划分成足够多的单元，这样可以认为每个单元内的温度、变形速度和变形程度相同，在每个单元内采用理想刚塑性材料模型，不同单元采用不同的屈服应力，这样处理才能得到比较接近实际的结果。,2.3.2刚塑性硬化材料模型,刚塑性硬化材料基本假设如下：(1)材料均质各向同性；(2)忽略材料的弹性变形，不计体积力与惯性力；(3)材料的变形流动服从Levy-Mises流动理论；(4)材料的体积不变；(5)考虑加工硬化和变形抗力对变形速度的敏感性。,2.3.2.1刚塑性硬化材料的应力应变关系,图2-2刚塑性硬化材料的应力应变关系,2.3.2.2刚塑性硬化材料的变形抗力,对于刚塑性硬化材料来说，当材料的化学成份和物理状态一定时，通常把变形抗力表示成变形温度、变形速度和变形程度的函数：,2.3.2.3Mises流动法则,理想刚塑性材料和刚塑性硬化材料都假设材料是不可压缩的，根据Mises流动法则，变形速度分量与偏差应力分量成正比，即这种材料的变形速度场与偏差应力场一一对应，但由于静水压力是不确定的，所以当以速度为未知量进行求解时，不能直接求得应力场。而体积可压缩材料模型可巧妙地解决这一问题。,(2-4),2.3.3刚塑性可压缩材料模型,刚塑性可压缩材料的基本假设如下：(1)材料均质各向同性；(2)忽略材料的弹性变形，不计体积力与惯性力；(3)材料的变形流动服从Levy-Mises流动理论；(4)材料的体积微可压缩；(5)考虑加工硬化和变形抗力对变形速度的敏感性。,(1)刚塑性可压缩材料的屈服条件,刚塑性可压缩材料与刚塑性硬化材料的主要差别是放松了体积不变条件的约束，即假设屈服与静水压力有关，屈服条件不仅取决于偏差应力的二次不变量，也取决于应力的一次不变量。,(2-5),(2-6),刚塑性可压缩材料的屈服条件,在主轴条件下，屈服条件：,(2-7),(2-8),Misees屈服曲面,平面,图2-3刚塑性可压缩材料主应力空间的屈服曲面,可压缩参数对屈服条件的影响,在平面上(1+2+3=3m=0),椭球体与Mises屈服圆柱相切，即刚塑性可压缩材料与Mises屈服条件一致。在平面以外(m0)，刚塑性可压缩材料比理想刚塑性材料容易进入屈服状态。可压缩参数g值越小，椭球体长轴延长，越接近于Mises屈服条件。当考虑加工硬化时，椭球体的体积将随着加工硬化而膨胀。,数量场和矢量场,数量场：对于空间或部分空间的任意一点M，都有一个确定的数值f(M)与之对应，则称在这个空间或部分空间上确定了一个数量场。该数量场可用数值函数f(M)来确定。例如温度场就是一个数量场。矢量场：对于空间或部分空间的任意一点M，都有一个确定的矢量f(M)与之对应，则称在这个空间或部分空间上确定了一个矢量场。该矢量场可用矢量函数f(M)来确定。例如速度场、应力场均是矢量场。,势和势场,梯度：设函数u=f(x,y,z)在空间或部分空间具有一阶连续偏导数，则对空间或部分空间的每一个点P(x,y,z)都可以确定一个矢量：这个矢量就称作函数u=f(x,y,z)在点P(x,y,z)的梯度。利用场的概念，该矢量函数在空间或部分空间确定了一个矢量场，即梯度场。它是由数值函数产生的，称数值函数为这个矢量场的势，矢量函数则称为有势场或势场。,(2)塑性势和变形速度分量,假设刚塑性可压缩材料的屈服函数为塑性势：,(2-9),根据塑性势的定义：,变形速度分量,利用(2-9)式和(2-10)式可直接求出应变速度分量：,(2-11),与之间关系,下面以为例推导与之间关系,与之间关系,可得：,(2-12),的简化形式,根据单位体积的塑性变形功率：,(2-13),把(2-11)式代入(2-13)式整理可得：,所以,(2-14),(3)应力应变速度关系,把(2-14)式代入(2-11)式,刚塑性可压缩材料与g和关系,(2-15),当一定时，g越小，体积变形速度小；当g=0时，此时刚塑性可压缩材料的屈服条件变为Mises屈服条件，应力应变关系变为Levy-Mises流动法则，刚塑性可压缩材料模型变为不可压缩的理想刚塑性或刚塑性硬化材料模型。从这个意义来讲，刚塑性不可压缩材料可以看作刚塑性可压缩材料在g=0时的特例，因此更具有普遍性。,根据(2-11)和(2-14)式可得体积变形速率与可压缩参数和静水压力关系：,(3)应力应变速度关系,如果采用张量形式，并引入克罗内克尔符号，可得刚塑性可压缩材料的应力应变速度关系：,(2-16),(4)等效应变速度,把(2-16)代入(2-13)式：,(2-17),这样，利用(2-16)式就可以从变形速度分量中直接求出应力场。因此，当从运动许可速度场出发以节点速度为未知量求解塑性变形过程时，能够简便地求出变形和力能参数。,(5)刚塑性可压缩材料的流动法则,根据(2-15)式可得：将上式代入(2-11)式，可得：,(2-18),(5)刚塑性可压缩材料的流动法则,上式表明，刚塑性可压缩材料的流动法则是偏差应变速度分量与偏差应力分量成正比。与Levy-Mises流动法则一致，只不过是体积变形速率不等于零。,(2-19),(6)刚塑性可压缩材料的特点,刚塑性可压缩材料与刚塑性硬化材料的主要差别是，放松了体积不变条件的约束，即假设屈服与静水压力有关。体积变化率取决于静水压力，当求出材料的屈服应力、等效应变速率和给定材料的体积可压缩参数后，可以直接从速度场求得应力场，所得结果与体积不可压缩条件的解十分接近，而且计算过程得到简化。,(7)可压缩参数g对体积变形的影响,将式(2-15)两端乘以时间增量dt，可得：,(2-20),以单向压缩应力状态为例，其应力分量及等效应力为：,(2-21),所以体积变形为：,(2-23),(2-22),(7)可压缩参数g对体积变形程度的影响,利用上式计算不同等效应变下，可压缩参数g对体积变形程度的影响表2-1可压缩参数g对体积变形程度的影响,从上表可见，在1030%的等效应变条件下，可压缩参数g取0.010.0001之间时，材料的体积变化不超过0.1%，可以近似满足体积不可压缩条件。因此，求解金属材料的轧制过程，可压缩参数g可取0.010.0001。,2.4刚塑性材料的基本方程,刚塑性材料发生塑性变形时，由表面积S所围成的体积V中，应力、速度和应变速率应满足下列基本方程:1.力平衡微分方程,(2-24),2.几何方程,(2-25),2.4刚塑性材料的基本方程,2.4刚塑性材料的基本方程,3.本构关系,(2-26),2.4刚塑性材料的基本方程,4.屈服准则5.体积不可压缩条件6.边界条件,(2-27),(2-28),(2-29),(2-30),2.4刚塑性材料的基本方程,式中：、等效应变速率和等效应力；偏差应力；屈服剪切应力；外力表面单位法线矢量的方向余弦。对于刚塑性材料发生塑性变形时，需要对上述基本方程进行联立求解来确定变形区的速度、应变速率和应力。由于直接求解这些偏微分方程组非常困难，因此，人们寻求其它求解途径，即利用变分原理求解。,2.5刚塑性有限元的基本原理,刚塑性材料的变分原理是刚塑性有限元法的理论基础，变分原理通过能量积分把偏微分方程组的求解问题变成了求泛函极值问题，从而为各种实际问题的求解提供了一种新方法。材料模型不同，变分原理的形式也不相同。根据附加条件的情况，变分原理可分为一般变分原理、不完全广义变分原理和完全广义变分原理。,2.5.1Mapkov原理(第一变分原理),设刚塑性变形体的体积为V，表面积为S0。在Sp上给定表面力，在Sv上给定速度vi，在满足变形的几何条件、速度边界条件和体积不变条件的一切许可速度场中，真实的速度场使泛函取极小值。,(2-31),2.5.2刚塑性可压缩材料的变分原理,在满足变形的几何条件、速度边界条件和体积微可压缩条件的一切许可速度场中，真实的速度场使泛函取极小值。,(2-32),2.5.3刚塑性材料的广义变分原理,在第一变分原理中，运动许可速度场必须满足几何条件、速度边界条件和体积不可压缩条件。在处理实际问题时，有些条件比较容易满足，而有的条件不易满足，因此引用拉格朗日乘子把这些约束条件全部或部分引入总体能量泛函中，通过泛函变分使这些约束条件得到满足。引入约束条件后，变分原理的表达式要发生变化，统称为广义变分原理。引入部分约束条件的称为不完全广义变分原理，全部约束条件同时引入的称为完全广义变分原理。,(1)不完全广义变分原理,设定初始速度场时，几何条件和速度边界条件容易满足，而体积不可压缩条件则不易满足，所以把体积不变条件通过拉格朗日乘子引入能量泛函中，得到新能量泛函：在满足几何方程和速度边界条件的一切速度场中，真实的速度场使上述能量泛函取驻值。,(2-33),(2)完全广义变分原理,在一切位移速度、应变速度和应力函数中，使能量泛函取驻值的vi、和必为刚塑性材料的正确解。,(2-34),第一变分原理与广义变分原理的差别,第一变分原理要求能量泛函取最小值，而广义变分原理仅要求能量泛函取驻值；第一变分原理要求速度场满足运动许可条件，静力许可条件是通过变分过程近似满足的，而广义变分原理的速度场不受任何约束，所有方程均由变分过程近似满足；第一变分原理比广义变分原理所得到结果精度高，但前者初始速度场的设定比后者难。,2.6速率敏感材料的总体能量泛函,求解轧制问题时，总体能量泛函的表达式为：式中：塑性变形功率；接触表面摩擦功率；外张力功率，前张力取负号，后张力取正号；速度不连续面上的剪切功率。,(2-35),2.6.1内部塑性变形功率,对于速率敏感性材料或在高温下成型的金属,单位体积的塑性变形功率：速率敏感材料的内部塑性变形功率为：,(2-37),(2-38),2.6.2接触表面上的摩擦功率,轧制过程中，轧辊与轧件接触表面上存在中性面，为了避免因摩擦力变向出现第二类奇异点而导致能量泛函不收敛的问题，可采用摩擦应力模型：摩擦应力是相对滑动速度vg的函数，因此摩擦功率为：,(2-38),(2-39),2.6.3外张力功率,式中：作用有外张力的表面；张应力，轧件受拉时取正值，受推时取负值；相应表面处的位移速度。,(2-40),2.6.4速度不连续面上的剪切功率,式中：速度不连续面；屈服剪应力；速度不连续面上的速度不连续量。,(2-41),速率敏感材料轧制过程总体能量泛函,在轧制过程中，边界上的外力为接触表面上的摩擦应力、变形区端部的张应力和速度不连续面上的剪切应力，相应的总体能量泛函为：对于简单轧制过程来说，总体能量泛函可简化为：,(2-42),(2-43),2.7刚塑性有限元的求解途径,设定运动许可速度场。建立总体能量泛函，并把泛函表示成速度的函数。利用数学上的极值理论求解泛函的驻值或最小值。利用几何方程和本构方程确定应变速度场和应力场。通过接触表面应力积分求出总轧制力和平均单位压力，利用速度场求出轧件宽展、前滑和轧件侧面形状参数。,2.7.1轧制变形区的有限元离散化,(1)研究目标的选择,根据变形特征确定研究目标。利用对称面上一个方向速度为零的特点，可使问题简化。如果变形过程的几何条件、物理条件在变形区内部及边界上对称，那么其真实解也必然对称。,图2-2板带轧制过程的对称面,(2)单元类型的选取,研究目标确定之后，便可着手在所选区域内设置节点、划分单元。除了速度已知边界之外，每个节点都有一组速度未知量，对于2维问题节点速度为,3维问题为，所以设置的节点越多，求解时未知数越多。单元类型的选择，可根据所处理问题的特点来确定，求解平面变形问题，常用四边形单元，求解3维变形问题常用立方体单元。,2.7.2总能量泛函的离散化,把轧制变形区划分为有限个单元之后，总能量泛函便可由单元能量泛函迭加求和来得到。下面以简单轧制过程为例进行说明。,(1)单元能量泛函,对于平面变形条件下的简单轧制过程，设单元的面积为Se，外力已知边界为Le，则单元能量泛函为：当单元的速度插值函数设定之后，单元内的便可确定。在给定材料屈服应力和摩擦边界条件下，单元能量泛函实际上是单元节点速度的函数：,(2-44),(2-45),(2)总能量泛函,根据(2-44)式对所研究区域全部单元的能量泛函求和，便可得到总能量泛函：式中：nm单元总数；nl接触表面单元数目。,(2-46),(3)接触表面的速度边界条件,轧制过程接触表面上的节点必须满足速度边界条件：,(2-47),图2-3轧制特征角,(4)总体能量泛函,由于接触表面节点的两个速度分量只有一个是未知量，同时对于x轴的节点垂直方向的速度已知(为零)，除去节点已知速度分量之后，把节点未知速度分量统一用表示，总体能量泛函表示成未知速度分量的函数：式中：nx未知量总数。,(2-48),2.7.3总能量泛函的最小化,用刚塑性有限元法求解时，根据变分原理应求总体能量泛函的极小值。从式(2-48)可知，总体能量泛函是节点速度矢量的多元函数。根据多元函数求极值的条件，则有：由于变分是任意的，必须有下列方程组成立：,(2-49),(2-50),(1)线性化处理,采用Newton-Raphson方法线性化求解时，假定泛函连续并存在各阶导数，则借助于泰勒级数在任一点展开。在第次迭代时，泛函可在第次迭代求得的值上展开，略去二阶以上的高阶微量，从而得到以为未知量的线性方程组：,(2-51),(2-52),(2)阻尼因子,采用(2-52)式求解时，先给定一个初始速度场，然后用迭代法按上式求出，直到相邻两次迭代速度近似解的偏差充分小时为止。迭代时，为了防止发散一般采用公式：,(2-53),式中：阻尼因子，。,2.8刚塑性有限元的基本公式,从上节可以看出，刚塑性有限元求解过程的一个基本步骤是计算单元的能量泛函。为此，需要解决以下几个方面问题：设定初始速度场；设定单元的速度分布函数插值函数；建立单元变形速度与节点速度的关系；计算单元能量泛函；求单元能量泛函的一阶偏导数即梯度；求单元能量泛函的二阶偏导数即Hessian矩阵。,2.8.1单元坐标及速度插值、形状函数,单元内任意一点的速度及坐标可以用单元节点的速度及坐标进行插值计算。对速度及坐标的插值使用相同的线性插值函数，该函数只与单元的形状和节点的配置有关，因此，称之为形状函数或形函数，与此相应的单元称为线性等参单元。,(1)四边形线形等参单元,为了便于插值与积分，单元内部采用局部坐标系，整体坐标系的任意四边形单元映射到局部坐标系后就变为正方形单元，单元节点编号及坐标系的映射变换如图2-4所示。,图2-4单元节点编号及坐标系的映射变换,(2)单元局部坐标与整体坐标的变换,整体坐标系中4个节点坐标为，局部节点坐标如下：单元局部坐标与整体坐标的变换关系：,(2-54),(2-56),(3)四边形等参单元的形函数及其性质,式(2-56)中称为单元的形函数，具有以下性质：形函数是与坐标插值函数相同形式的双线性插值函数；形函数在节点上，；形函数在节点上，。,(4)形函数表达式的确定,以为例说明其确定过程：,上图局部坐标系中，单元24、34两条边方程可表示为：,(2-57),(4)形函数表达式的确定,根据形函数的性质(3)，在节点2、3、4处的值为零，而直线24、34经过这些节点，所以可取：根据形函数的性质(2)，在节点1处其值为1，把1点的局部坐标值代入上式得：,(2-58),(2-59),(2-60),(5)四边形等参单元的形函数,(2-61),(2-62),(6)单元内任一点的速度与节点速度的关系,(2-63),(2-64),2.8.2单元变形速度与节点速度关系、B矩阵,(1)平面变形的几何方程,(2-68),(2)以节点速度表示的几何方程,(2-69),(3)几何方程的矩阵形式,(2-70),式中：,(3)几何方程的矩阵形式,令,(2-74),(2-71),(2-73),(2-75),式中：单元应变矩阵或B矩阵；典型子矩阵。,2.8.3Jacobi矩阵及其逆矩阵和行列式,单元应变B矩阵中的元素bi及ci是形函数对x及y的偏导数，为了便于积分计算，形函数对x及y的偏导数应该用对及的偏导数来表示。根据复合函数的求导法则：,(2-76),(1)2维问题的Jacobi矩阵,将上式写成矩阵形式：,(2-77),(2-78),(1)2维问题的Jacobi矩阵,(2-79),(2-80),(2)Jacobi矩阵的逆矩阵及行列式,(2-81),(2-82),(2-83),2.8.4高斯积分,对于2维变形问题，建立单元内任意一点变形速度与节点速度关系之后，在给定初始速度场条件下，即可对单元能量泛函进行积分计算。在单元局部坐标系中，单元能量泛函式可表示成：刚塑性有限元中的数值积分常用高斯积分，即按照数学上的规则在单元内选取若干个积分点，用积分点处的函数值与求积系数之积的累加结果近似代替原积分。,(2-84),(1)单元塑性变形功率高斯积分表示,式中：n单元积分点个数，Hi、Hj求积系数。,(2-85),对于2维问题采用线性单元时，n=2,Hi=Hj=1，此时单元内有4个高斯积分点，积分点坐标为0.57735027，如图2-5所示。,(2)单元摩擦功率的积分,单元摩擦功率的积分，可采用1维高斯积分近似计算：,(2-86),对于1维问题采用线性单元时，n=2,Hi=1，此时单元内有2个高斯积分点，积分点坐标为0.57735027。,(3)1维问题的Jacobi矩阵行列式,式中：x1、x2为摩擦表面单元两个节点的整体坐标值；N1、N2为1维单元的两个形函数。,(2-87),(2-88),(4)单元的能量泛函,采用高斯积分后，2维变形简单轧制过程的单元能量泛函为：式中：分别为单元高斯积分点处的等效应力、等效应变及Jacobi矩阵行列式的值。,(2-90),2.8.5矩阵分析中的公式,(1)数值函数对向量变量的偏导数,设是给定的向量，是向量变量，且，则有,(2-a),(1)数值函数对向量变量的偏导数,设是给定的矩阵，是向量变量，且，则有,(2-b),(2-c),(2)向量函数对向量变量的偏导数,设是给定的矩阵，是向量变量，且，则有,(2-d),2.8.6能量泛函的一阶偏导数(梯度),(1)平面变形问题单元的等效应变速度列阵平面变形问题单元的变形速度列阵,(2-91),(2-92),式中的常数矩阵和向量,(2-94),(2-93),(1)2维平面变形条件下的等效应变,(2-95),(1)2维平面变形条件下的等效应变,上式可以进一步写成更为简洁的形式：,(2-96),(2-97),等效应变不用(2-96)式简洁形式？,因为采用体积可压缩方法进行求解时，如果等效应变采用(2-96)式简洁形式，则要求单元内每一个高斯积分点处的体积变形速度都很小，这样的约束条件过于苛刻，经常出现迭代计算发散现象，导致有限元数值求解过程无法进行。,(1)2维平面变形条件下的等效应变,为了排除每一个高斯积分点处体积变形速度都很小这一严格约束，单元的体积变形速度取各高斯积分点处体积变形速度的平均值，即：,(2-98),(2-100),(2-99),(2)等效应变速度对速度向量的一阶偏导数,等效应变速度(2-100)式对速度向量的一阶偏导数：,(2-101),(2-102),(3)单元梯度,把变形抗力模型代入单元塑性变形功率表达式，并对单元节点速度向量求一阶偏导数：,(2-103),(4)单元摩擦功率对速度向量的一阶偏导数,接触表面任意一点的速度边界条件：,接触表面任意一点的相对滑动速度：,(2-104),单元内任意一点的速度与节点速度的关系,对于2维轧制过程的接触表面单元而言，单元节点上只有一个未知量，因此，单元内任意一点的速度与节点速度的关系如下：,(2-105),(2-106),对节点速度向量的一阶偏导数,把(2-105)式代入(2-104)式可得相对滑动速度：,(2-107),(2-108),(2-109),单元能量泛函的一阶偏导数即梯度,接触表面单元摩擦功率的一阶偏导数：单元能量泛函的一阶偏导数即梯度：,(2-110),(2-111),2.8.7能量泛函的二阶偏导数(Hessian矩阵),2.8.7能量泛函的二阶偏导数(Hessian矩阵),(2-112),2.8.8总体能量泛函能量泛函的梯度矩阵和Hessian矩阵,根据迭加原理对所有单元的梯度矩阵和Hessian矩阵求和，可得到总体能量泛函的梯度和Hessian矩阵：总体能量泛函的梯度和Hessian矩阵中的元素是按节点未知速度分量的顺序进行排列的，因此，上述求和需要按整体节点编号进行。,(2-113),2.93维问题的刚塑性有限元基本公式,严格来讲，各种轧制过程均是3维变形过程。虽然传统的轧制理论认为平辊轧制可近似看成平面变形，这在计算机轧制力、力矩和前滑等积分量或均值量时尚可，但要从理论上解决近年来板材生产中提出的板形控制、横向厚差控制以及组织性能控制等问题，不可避免地要涉及到板材轧制过程中的横向金属流动、辊形曲线和板凸度问题以及变形区内轧件的温度和速度分布问题。因此，研究3维有限元法具有更大的实际意义。,2.9.1空间八节点六面体等参单元,对于三维轧制问题，一般采用空间八节点六面体等参单元。这种单元在整体坐标系中是任意六面体，映射到局部坐标系之后变为边长为2的立方体单元。单元局部节点坐标及编号如图2-6所示。,单元局部节点坐标,表2-2单元局部节点坐标,(1)单元的形函数,求解3维轧制问题时坐标系选取原则：整体坐标系中的x和局部坐标系中表示轧制方向，y和表示轧件宽度方向，z和表示轧件的厚度方向。单元的形函数如下：,(2-114),式中：单元内任意一点的局部座标；单元节点的局部座标。,(2)整体座标与局部座标的变换关系,(2-115),式中：已知节点的整体座标。,(3)单元速度与节点速度的关系,(2-116),式中：已知节点的速度分量。,(2-117),(2-118),(2-119),(2-120),式中：单元形函数矩阵；单元内任意一点的速度列阵；单元节点速度列阵。,(4)单元应变速度矩阵,(4)单元应变速度矩阵,(4)单元应变速度矩阵,(2-121),(2-122),(2-123),(2-124),(5)形函数对整体坐标的偏导数,(6)Jacobi矩阵,(2-125),(7)形函数对整体坐标的偏导数,(2-126),(8)形函数对局部坐标的偏导数,(2-127),2.9.2单元塑性变形功率的一阶偏导数及二阶偏导数,采用刚塑性体积可压缩材料模型求解时，考虑材料的变形速率敏感性，将单元的塑性变形功率用下式表示：把变形抗力模型代入上式，求一阶偏导：,(2-128),(2-129),(1)等效应变速度,采用体积可压缩法求解时，为了排除单元内每一高斯点的体积变形速度都很小的严格约束，应取单元体积变形速度的平均值，即：,(2-130),(2-131),(2-133),(2-132),(2)等效应变速度的一阶偏导数,(2-134),(2-135),(3)单元塑性变形功率的一阶偏导数,把(2-134)、(2-135)式代入(2-129)式，可得单元塑性变形功率的一阶偏导数：,(2-136),(4)单元塑性变形功率的二阶偏导数,(4)单元塑性变形功率的二阶偏导数,(2-137),(2-138),2.9.3单元摩擦功率的一阶偏导数及二阶偏导数,单元的摩擦功率,对上式求一阶偏导数:,(2-139),(2-140),(1)三维变形接触摩擦表面的相对滑动速度,接触摩擦表面上的速度边界条件:,(2-141),接触摩擦表面的相对滑动速度：,(2-143),(2-142),(2)接触表面单元速度与节点速度的关系,(2-148),(2-147),(2-146),(2-145),(2-144),(3)接触摩擦表面的相对滑动速度,(2-151),把(2-144)、(2-149)和(2-150)式代入(2-143)式可得：,(2-149),(2-150),对节点速度列阵求偏导,(2-154),(2-152),(2-153),(4)摩擦功率的一阶偏导数,把(2-152)式代入(2-140)式可得单元摩擦功率的一阶偏导数：,(2-155),(5)单元摩擦功率的二阶偏导数,(5)单元摩擦功率的二阶偏导数,(2-156),(6)单元能量泛函的梯度和Hessian矩阵,把单元的塑性变形功率及摩擦变形功率的一阶偏导数(2-136)、(2-155)及二阶偏导数(2-137)、(2-156)求和，可得单元能量泛函的梯度及Hessian矩阵：,(2-157),3.刚塑性有限元基本问题及其处理方法,单元节点编号初始速度场的设定奇异点的处理刚塑性有限元分析过程的收敛准则刚性区与塑性区分界的处理轧件侧面形状和宽展的确定接触表面速度边界条件的处理,3.1单元编号和节点总体编号,刚塑性有限元在求解过程中，总体刚度矩阵即Hessian矩阵是采用一维压缩存贮方式进行存贮的，该矩阵一般称之为优化矩阵。在有限元网格划分之后，单元编号和节点总体编号虽然不能改变总体刚度矩阵的规模、节点未知速度的总数，但是却能改变总体刚度矩阵的带宽而影响优化矩阵的规模，从而给计算机存贮容量和计算时间带来很大的影响。,(1)单元总体节点编号的原则,同一单元最大节点编号与单元最小节点编号的差值最小，使总体刚度矩阵的带宽最小。较好节点编号方法是从坐标原点开始，首先沿单元数目最少的方向增加单元和节点的编号，然后沿单元次少的方向进行编号，最后沿单元数目最多的方向进行编号。,(3-1),(2)节点编号方式对总刚带宽的影响,图(a)同一单元节点编号的极差为：35+31+2-1=19；相应总体刚度矩阵的半带宽为：3(19+1)=60。图(b)同一单元节点编号的极差为：105+101+2-1=61；相应总体刚度矩阵的半带宽为：3(61+1)=186。,3.2初始速度场的设定,设定轧制过程初始速度场的方法有以下几种：1.初等方法设定的均匀初始速度场；2.用逐步细化网格法设定初始速度场；3.函数法设定近似的初始速度场。,3.2.1用初等方法设定速度场,对于平面变形轧制过程，初始速度场设定采用假设：平断面假设，即同一垂直横断面各点水平速度vx相同；厚度方向vy呈线性分布；通过任意垂直横断面的金属秒流量相等。,(1)中性面处轧件的半高度,秒流量与中性面的位置有关，中性角：,(3-2),式中：咬入角；f摩擦系数。中性面处轧件的半高度：,(3-3),(2)金属秒流量,中性面处轧件的水平方向运动速度：,(3-4),金属秒流量为：,(3-5),(3)各节点水平速度的设定,根据秒流量相等条件，接触表面节点i处的水平速度：,根据平断面假设，可得：,(3-6),(3-7),(3-8),(4)各节点垂直速度的设定,根据速度边界条件，可得接触表面节点i处的垂直速度：根据厚度方向速度呈线性分布及对称面上节点垂直速度的特点，可得内部节点垂直速度：,(3-9),(3-10),显然，按照(3-7)至(3-10)式设定的初始速度场满足速度边界条件和体积不可压缩条件。,3.2.2用逐步细化网格法设定初始速度场,基本思想：先用少数几个单元，以均匀速度场作为初始速度场进行计算，可以很快得到收敛计算结果。单元细化，可把原来一个单元划分成2、4、或8个单元，新增节点可利用原来单元相邻两个节点速度插值求得。将此速度场作为细化后网格的初始速度场进行计算。该细化过程可重复进行，直到得到令人满意的网格密度为止。,关键技术,图3-3有限元网格细化示意图,利用细化网格的方法设定初始速度场时，在预处理程序中要有相应的子程序来处理新老两种网格中节点编号的对应关系，在此基础上正确处理新产生节点的坐标和速度插值关系。,3.2.3G函数法设定初始速度场,采用G函数法设定初始速度场是森谦一朗提出的。基本思想是寻找一个与总体能量泛函相近的泛函G，泛函G应便于求导，并能在G=0时得到一个线性方程组。可见，G函数的选择十分关键，函数可选择如下形式：,(3-11),G函数的偏导数,(3