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

1、刚塑性有限元法及其在轧制中的应用刚塑性有限元法及其在轧制中的应用轧制技术及连轧自动化国家重点实验室徐建忠学习的主要内容学习的主要内容l刚塑性有限元的基本概念和基本理论；刚塑性有限元的基本概念和基本理论；l刚塑性有限元相关技术问题的处理方法；刚塑性有限元相关技术问题的处理方法；l求解轧制过程的刚塑性有限元程序。1.1.刚塑性有限元法的基本理论刚塑性有限元法的基本理论 1.1 1.1 有限元法的基本概念有限元法的基本概念 有限元法有限元法：把工件划分成有限结点相连接的单元，以结点上的速度(位移)作为未知量，利用最小能原理求解相应的方程组确定此未知量，按结点速度与单元内部应变以及单元内部应力之间的关

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

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

4、； (4) 材料的体积不变； (5) 不考虑加工硬化，忽略变形抗力对变形速度的敏感性。1.3.1.1 1.3.1.1 理想刚塑性材料模型的应力应变关系理想刚塑性材料模型的应力应变关系图1-1 理想刚塑性材料的应力应变关系1.3.1.2 1.3.1.2 理想刚塑性材料模型的特点理想刚塑性材料模型的特点l 只要等效应力达到一恒定数值，材料便发生屈服，而且材料在整个变形过程中屈服应力不再发生变化。l 采用理想刚塑性材料模型进行能量积分时，可以把等效应力做为常数提到积分号之外，从而使积分过程得到简化。1.3.1.3 1.3.1.3 采用该模型进行采用该模型进行FEMFEM求解应该注意的问题求解应该注意

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

6、积不变； (5) 考虑加工硬化和变形抗力对变形速度的敏感性。 1.3.2.1 1.3.2.1 刚塑性硬化材料的应力应变关系刚塑性硬化材料的应力应变关系 图1-2 刚塑性硬化材料的应力应变关系1.3.2.2 1.3.2.2 刚塑性硬化材料的变形抗力刚塑性硬化材料的变形抗力 对于刚塑性硬化材料来说，当材料的化学成份和物理状态一定时，通常把变形抗力表示成变形温度、变形速度和变形程度的函数：( , , ) sf T 0nsa()nmsbc(1-1)(1-2)(1-3)1.3.2.3 1.3.2.3 Mises流动法则流动法则 理想刚塑性材料和刚塑性硬化材料都假设材料是不可压缩的，根据Mises流动法则

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

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

9、ees屈服曲面平面1122 , 23图1-3 刚塑性可压缩材料主应力空间的屈服曲面(2) (2) 可压缩参数对屈服条件的影响可压缩参数对屈服条件的影响 在平面上(1 + 2 + 3 =3m=0 ),椭球体与Mises屈服圆柱相切，即刚塑性可压缩材料与Mises屈服条件一致。在平面以外(m0)，刚塑性可压缩材料比理想刚塑性材料容易进入屈服状态。可压缩参数g值越小，椭球体长轴延长，越接近于Mises屈服条件。当考虑加工硬化时，椭球体的体积将随着加工硬化而膨胀。 (3) (3) 塑性势和变形速度分量塑性势和变形速度分量 假设刚塑性可压缩材料的屈服函数为塑性势 ：1222222221()()()()6

10、()2ijxyyzzxxyyzzxmFg(1-9)()ijijijFdd根据塑性势的定义：()ijijijFd则(1-10)dtdddt 两边同除以时间增量，令(4) (4) 变形速度分量变形速度分量 利用(1-9)式和(1-10)式可直接求出应变速度分量：2(1)92(1)92(1)9222xxmyymzzmxyxyyzyzzxzxgdgdgdddd(1-11)(5) (5) 与与 之间关系之间关系 下面以 为例推导 与 之间关系dd()ijxxFddd1222222221()()()()6()2ijxyyzzxxyzmFg1222222222()()()312()()()6()22223(

11、)33 2 ()2 ()3 xyzxmijxxyyzzxxyyzzxmxyzmxxyzmijijxgFgggFF由于2233(1)392 ()2 ()mmxmijijggFF (1-12)32 ()ijddF dd(5) (5) 与与 之间关系之间关系 (6) (6) 的简化形式的简化形式 根据单位体积的塑性变形功率 ：dxxyyzzxyxyyzyzzxzxdw (1-13)把(2-11)式代入(2-13)式整理可得：212921292129xxmxyxyyymyzyzzzmzxzxgddgddgdd所以 32d (1-14)222222222222222222222222()1()92 2(

12、)3192 2()331 333xyzxyyzzxmxyxxyzxyyzzxmxyzxyyzzxmmxydgdgdgd 2222222222222222222222236()921 333()6()231 2222226()231 ()(3zxyyzzxmmxyzxyzxyyzzxmxyzxyyzzxxyyzzxmxygdgdgd 2222222)()6()22 3yzzxxyyzzxmgd(7) (7) 刚塑性可压缩材料刚塑性可压缩材料 与与 g g 和和 关系关系 22333vxyzxyzmmg dgdg(1-15)当 一定时，g 越小，体积变形速度小；当g = 0时，此时刚塑性可压缩材料

13、的屈服条件变为Mises屈服条件，应力应变关系变为Levy-Mises流动法则，刚塑性可压缩材料模型变为不可压缩的理想刚塑性或刚塑性硬化材料模型。从这个意义来讲，刚塑性不可压缩材料可以看作刚塑性可压缩材料在g=0时的特例，因此更具有普遍性。 0v 根据(1-11)和(1-14)式可得体积变形速率与可压缩参数和静水压力关系 ：mV1222222221()()()6()2xyyzzxxyyzzxmgm(8) (8) 应力应变速度关系应力应变速度关系 把(1-14)式代入变形速度分量(1-11)式21921392219222xxmyymvzzmvmmxyxyyzyzzxzxgdgddgdggddd

14、(8) (8) 应力应变速度关系应力应变速度关系 2213922 13922 139212 39xxmvxvxxvgggggg0.5/2 3xyxyxyd如果采用张量形式，并引入克罗内克尔符号，可得刚塑性可压缩材料的应力应变速度关系：212()39ijijijvg(1-16)(9) (9) 等效应变速度等效应变速度 把应变速度关系(1-16)式 代入单位体积的塑性变形功率 (1-13)式：212()39ijijijvg212 =39212 39212 3922 33xxyyzzxyxyyzyzzxzxxxvyyvzzvxyxyyzyxggg 23zxzx 12222222221()()()6(

15、)9xyyzzxxyyzzxvg(1-17) 这样，利用(1-16)式就可以从变形速度分量中直接求出应力场。因此，当从运动许可速度场出发以节点速度为未知量求解塑性变形过程时，能够简便地求出变形和力能参数。2222222222222222222124()()()()3932124 ()()3932 3()6(9xyzvxyzvxyzxyyzzxxyzVVxyyzzxxyzVxygg222222222221)21 ()()()6()9yzzxVxyyzzxxyyzzxVgg(10) (10) 刚塑性可压缩材料的流动法则刚塑性可压缩材料的流动法则 根据体积变形速率与可压缩参数和静水压力关系 (1-1

16、5)式可得： 将上式代入变形速度分量(2-11)式，可得：29mmgd(1-18)219221993 2xmmxmxmxmxmmxmgdgdgdd(10) (10) 刚塑性可压缩材料的流动法则刚塑性可压缩材料的流动法则 上式表明，刚塑性可压缩材料的流动法则是偏差应变速度分量与偏差应力分量成正比。与Levy-Mises流动法则一致，只不过是体积变形速率不等于零。32ymxyyzxmzmzxxmymzmxyyzzxd(1-19)(11) (11) 刚塑性可压缩材料的特点刚塑性可压缩材料的特点 刚塑性可压缩材料与刚塑性硬化材料的主要差别是，放松了体积不变条件的约束，即假设屈服与静水压力有关。体积变化

17、率取决于静水压力，当求出材料的屈服应力、等效应变速率和给定材料的体积可压缩参数后，可以直接从速度场求得应力场，所得结果与体积不可压缩条件的解十分接近，而且计算过程得到简化。(12) (12) 可压缩参数可压缩参数g g对体积变形的影响对体积变形的影响 将式(1-15)两端乘以时间增量dt，可得：vmg(1-20)以单向压缩应力状态为例，其应力分量及等效应力为：0 0 0 0 03xyxyyzzxzm (1-21)9/39zmgg所以体积变形为：(1-23)(1-22)9vgg(12)(12) 可压缩参数可压缩参数g对体积变形程度的影响对体积变形程度的影响 利用上式计算不同等效应变下，可压缩参数

18、g对体积变形程度的影响 表1-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.010.0001。1.4 1.4 刚塑性材料的基本方程刚塑性材料的基本方程 刚塑性材料发生塑性变形时，

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

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

21、的基本原理l 刚塑性材料的变分原理是刚塑性有限元法的理论基础，变分原理通过能量积分把偏微分方程组的求解问题变成了求泛函极值问题，从而为各种实际问题的求解提供了一种新方法。l 材料模型不同，变分原理的形式也不相同。根据附加条件的情况，变分原理可分为一般变分原理、不完全广义变分原理和完全广义变分原理。1.5.1 Mapkov原理原理 (第一变分原理第一变分原理) 设刚塑性变形体的体积为V，表面积为S0。在Sp上给定表面力 ，在Sv上给定速度 vi，在满足变形的几何条件、速度边界条件和体积不变条件的一切运动许可速度场中，真实的速度场使泛函 取极小值。ip23PsijijiiVSdVp v ds (1

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

23、束条件的称为不完全广义变分原理，全部约束条件同时引入的称为完全广义变分原理。 (1) 不完全广义变分原理不完全广义变分原理 设定初始速度场时，几何条件和速度边界条件容易满足，而体积不可压缩条件则不易满足，所以把体积不变条件通过拉格朗日乘子引入能量泛函中，得到新能量泛函： 在满足几何方程和速度边界条件的一切速度场 中，真实的速度场使上述能量泛函取驻值。 ivpiiijijVSVdVp v dsdV (1-33)(2) 完全广义变分原理完全广义变分原理 在一切位移速度、应变速度和应力函数中，使能量泛函 取驻值的vi、 和 必为刚塑性材料的正确解。ijij, ()1 ()2pviiijijiiiVS

24、VSijiji jj iVdVp v dsdVvv dsvvdV(1-34)第一变分原理与广义变分原理的差别第一变分原理与广义变分原理的差别 第一变分原理要求能量泛函取最小值，而广义变分原理仅要求能量泛函取驻值； 第一变分原理要求速度场满足运动许可条件，静力许可条件是通过变分过程近似满足的，而广义变分原理的速度场不受任何约束，所有方程均由变分过程近似满足； 第一变分原理比广义变分原理所得到结果精度高，但前者初始速度场的设定比后者难。 1.6 速率敏感材料的总体能量泛函速率敏感材料的总体能量泛函 求解轧制问题时，总体能量泛函的表达式为： 式中： 塑性变形功率； 接触表面摩擦功率； 外张力功率，前

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

26、fggvssvmddsvkk ds (1-39)1.6.3 外张力功率外张力功率 式中： 作用有外张力的表面； 张应力，轧件受拉时取正值，受推时取负值； 相应表面处的位移速度。isv ds(1-40)siv1.6.4 速度不连续面上的剪切功率速度不连续面上的剪切功率 式中： 速度不连续面； 屈服剪应力； 速度不连续面上的速度不连续量。(1-41)ksssvkssssv ds1.6.5 速率敏感材料轧制过程总体能量泛函速率敏感材料轧制过程总体能量泛函 在轧制过程中，边界上的外力为接触表面上的摩擦应力、变形区端部的张应力和速度不连续面上的剪切应力，相应的总体能量泛函为： 对于简单轧制过程来说，总体

27、能量泛函可简化为： 221 1 ()13sgiksVssssmdVvkk dsvdsv dsm(1-42)221 1 ()13sgVsmdVvkk dsm(1-43)1.7 刚塑性有限元的求解途径刚塑性有限元的求解途径l设定运动许可速度场。l建立总体能量泛函，并把泛函表示成速度的函数。l利用数学上的极值理论求解泛函的驻值或最小值。l利用几何方程和本构方程确定应变速度场和应力场。l通过接触表面应力积分求出总轧制力和平均单位压力，利用速度场求出轧件宽展、前滑和轧件侧面形状参数。1.7.1 轧制变形区的有限元离散化轧制变形区的有限元离散化 (1) (1) 研究目标的选择研究目标的选择zox根据变形特

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

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

30、中：nm单元总数； nl 接触表面单元数目。ee221 1111 ()13enmnmnlsgSLeeemdSvkk dLm(1-46)(3) 接触表面的速度边界条件接触表面的速度边界条件 轧制过程接触表面上的节点必须满足速度边界条件： tgyxvv (1-47)图1-5 轧制特征角(4) 总体能量泛函总体能量泛函 由于接触表面节点的两个速度分量只有一个是未知量 ，同时对于 x 轴的节点垂直方向的速度已知(为零)，除去节点已知速度分量之后，把节点未知速度分量统一用 表示，总体能量泛函表示成未知速度分量的函数：式中：nx未知量总数。xviv12( ,)inxf v vvv(1-48)1.7.3 总

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

32、eTeeekkkkvvvvv211111 eenmnmnmeeeTeeeekkkvvvv (1-51)(1-52)(2) (2) 阻尼因子阻尼因子 采用(1-52)式求解时，先给定一个初始速度场 ，然后用迭代法按上式求出 ，直到相邻两次迭代速度近似解的偏差充分小时为止。迭代时，为了防止发散一般采用公式： vv1 kkvvv(1-53)01式中： 阻尼因子， 。2. 刚塑性有限元的基本公式刚塑性有限元的基本公式 从第一章可以看出，刚塑性有限元求解过程的关键内容是计算单元的能量泛函计算单元的能量泛函。为此，需要解决以下几个方面问题：l设定初始速度场初始速度场；l设定单元的速度分布函数速度分布函数插

33、值函数；l建立单元变形速度单元变形速度与节点速度节点速度的关系；l计算单元能量泛函单元能量泛函；l求单元能量泛函的一阶偏导数即梯度梯度；l求单元能量泛函的二阶偏导数即Hessian矩阵矩阵。 2. 1 单元坐标及速度插值、形状函数单元坐标及速度插值、形状函数 单元内任意一点的速度及坐标可以用单元节点的速度及坐标进行线性线性插值计算。 对速度及坐标使用相同的线性插值函数相同的线性插值函数，该函数只与单元的形状和节点的配置有关，因此，称之为形状函数或形函数形函数，与此相应的单元称为线性等参单元线性等参单元。(1)(1) 四边形四边形4节点线形等参单元节点线形等参单元 为了便于插值插值与积分积分，单

34、元内部采用局部坐标系 ，整体坐标系的任意四边形单元映射到局部坐标系后就变为正方形单元，单元节点编号及坐标系的映射变换如图2-1所示。 42y13xO3124O( , ) 图2-1 单元节点编号及坐标系的映射变换(2)(2) 单元局部坐标与整体坐标的变换单元局部坐标与整体坐标的变换 ( , ) , 14kkxyk (, )kkint(1)/ 2)12341231-1-1 11( 1)-1 1-11( 1)kkkk (2-1)(2-2)41 122334414112233441( , )( , )kkkkkkxN xN xN xN xNxyN yN yN yN yNy 整体坐标系中4个节点坐标节点

35、坐标为 ，局部坐标系中的节点坐标 如下：),(yx( , ) 单元内部任意一点的局部坐标 与整体坐标 的变换关系：(3)(3) 四边形等参单元的形函数及其性质四边形等参单元的形函数及其性质 式(2-2)中 单元的形函数，具有以下性质：l形函数 是与坐标插值函数相同形式的双线性插值函数；l形函数 在节点k上， ；l形函数 在节点l上， 。 ( , )kN ( , )kN ( , )kN ( , )kN (,)1kkkN (,)0 , kllNlk (4)(4)形形函数表达式的确定函数表达式的确定 以 为例说明其确定过程： 1( , )N 42y13xO3124O上图局部坐标系中，单元24、34两

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

37、)(1)4411( , )(1)(1)(1)(1)44NNN 1( , )(1)(1) (1 4)4kkkNk (2-7)(2-8)(6)(6) 单元内任一点的速度与节点速度的关系单元内任一点的速度与节点速度的关系 4141( , )( , )xkxkkykykkvNvvNv (2-9) evNv(2-10) xyvvv 4321432100000000NNNNNNNNN Tyxyxyxyxevvvvvvvvv443322112.2 单元变形速度与节点速度关系、单元变形速度与节点速度关系、B矩阵矩阵 (1) 平面变形的几何方程平面变形的几何方程xxyyyxxyvxvyvvyx(2-11)(2)

38、 以节点速度表示的几何方程 4414141411( , )( , )()ixxxixixkxkykiyyyiiykykkyiixxyyixixyiNvvxxvNvvNvyyvNvvNNvvvxyyx (2-12)(3) 几何方程的矩阵形式 11212342123431122334434400000000 xyxxyyxxyyxyvvvbbbbvccccvcbcbcbcbvvv (2-13)4) , 3 , 2 , 1( , , iyNcxNbiiii式中：(3) 几何方程的矩阵形式 Txyxy4321BBBBB iiiiibccbB00 eBv令 (2-16)(2-14)(2-15)(2-17

39、)BiB式中： 单元应变矩阵或B矩阵； 典型子矩阵。2.3 Jacobi矩阵及其逆矩阵和行列式矩阵及其逆矩阵和行列式 单元应变B矩阵中的元素bi及ci是形函数 对x及y的偏导数，为了便于积分计算，形函数 对x及y的偏导数应该用 对 及 的偏导数来表示。根据复合函数的求导法则： ( , )iN ( , )iN ( , )iN TiiiiiTiiiiiNNNNNxyxyxyxyNNNNNxyxyxyxy (2-18)(1) 2维问题的维问题的Jacobi矩阵矩阵 将上式写成矩阵形式：iiiiNxyNxNNxyy(2-19) xyJxy(2-20)(1) 2维问题的维问题的Jacobi矩阵矩阵 1(

40、 , )(1)(1)4( , )1(1) 4 (1 4)( , )1(1)4kkkiiiiiiNNiN (2-21)444111444111( , ) ( , )iiiiiiiiiiiiiiiiiiNNxyxyxNxJxyNNyNyxy (2-22)(2) Jacobi矩阵的逆矩阵及行列式矩阵的逆矩阵及行列式 1iiiiiiiiNNxyNNxxJNNNxyNyy(2-23)1 1 yxJyxJ(2-24)xyxyJ(2-25)2.4 高斯积分高斯积分 对于2维变形问题，建立单元内任意一点变形速度与节点速度关系之后，在给定初始速度场条件下，即可对单元能量泛函进行积分计算。在局部坐标系中，单元能量

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

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

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

44、 矩阵分析中的公式矩阵分析中的公式 12341234123412344538()4538()4538()4538TTTTTTTTTTTAXxxxxf XA XX Axxxxf XA XffffAXXxxxxf XX AffffAxxxxXX(1) 数值函数对向量变量的偏导数数值函数对向量变量的偏导数 ()Tf XA XAXXijn nAaT12xXxxx()Tf XX AX设 是给定的矩阵， 是向量变量，且 ，则有() 2Tf XAAAXX如果 ，则 (2-b)(2-c)(1) 数值函数对向量变量的偏导数数值函数对向量变量的偏导数 1234123411234212343123441234112

45、345383645 135427634325637()34568543 (432)(563TTTAXxxxxxxxxxxxxxxf XX AXxxxxxxxxxxxxxxxxxx 42123431234412347) (3456)(8543)()TTxxxxxxxxxxxxf XffffA XAXXxxxx(2) 向量函数对向量变量的偏导数向量函数对向量变量的偏导数 ijn nAaT12nXxxx()f XAX设 是给定的矩阵， 是向量变量，且 ，则有()()TTTf XAXf XAX(2-d)(2) 向量函数对向量变量的偏导数向量函数对向量变量的偏导数 1234123412341234123

46、4123445383645 1354276345383645()354276345383645()13542763TTAXxxxxxxxxxxxxf XAXxxxxxxxxf XffffxxxxX1234()TTTAf XffffAXxxxx2.6 能量泛函的一阶偏导数能量泛函的一阶偏导数(梯度梯度)(1) 平面变形问题单元的等效应变速度列阵平面变形问题单元的等效应变速度列阵 平面变形问题单元的变形速度列阵 TexyxyB v222411931 xxyyxyvTTTgZCCg (2-32)(2-33)式中的常数矩阵和向量式中的常数矩阵和向量 4209924 0991003Z 1 1 0C (2

47、-35)(2-34) TZ 222224209924 0991003422419999342241419999393xTxyxyyxyxxyxyxyyxyxxyxyyxyxxyyZ 2xy 1 TTCCg 21 1111 10011TTxxyxyyxyxyxyVCCgggg (1) 2维平面变形条件下的等效应变维平面变形条件下的等效应变 222 4111() 93TexyxyTTxxyyxyvB vZgg 1 TeTTeeTTTeCCvBZ B vvBCC B vg(2-36)(1) 2维平面变形条件下的等效应变维平面变形条件下的等效应变 上式可以进一步写成更为简洁的形式： TeTTeZvBZ

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

49、CBvC Bv1 eTTeeTTTevBZ B vvBCC B vg(2-39)(2-41)411 4kkBB(2-40)(2) 等效应变速度对速度向量的一阶偏导数等效应变速度对速度向量的一阶偏导数 等效应变速度(2-41)式对速度向量的一阶偏导数：1 TeTTeePeBZBvBCCBvKvgv1 TTTPKBZBBCC Bg(2-42)(2-43) PK式中：88的矩阵。(3) 单元梯度单元梯度 把变形抗力模型代入单元塑性变形功率表达式，并对单元节点速度向量求一阶偏导数：44111144114411111111 nmenmPiiemnnmPeeeiiePeiiaJaJmmaJaJmvvvJJ

50、Kvv(2-44)(4) 单元摩擦功率对速度向量的一阶偏导数单元摩擦功率对速度向量的一阶偏导数 接触表面任意一点的速度边界条件： tgyxvv 接触表面任意一点的相对滑动速度：secgRxvvv(2-45)vxvy单元内任意一点的速度与节点速度的关系单元内任意一点的速度与节点速度的关系 对于2维轧制过程的接触表面单元而言，单元节点上只有一个未知量，因此，单元内任意一点的速度与节点速度的关系如下： efexffvvNv120.5(1- ) 0.5(1) fTefxxNvvv(2-46)(2-47) 对节点速度向量的一阶偏导数对节点速度向量的一阶偏导数gv 把(2-46)式代入(2-45)式可得相

51、对滑动速度： secegRffvvNv(2-48)sec gTfefvNv 2222()22 secgggeeffggTfefvkkvvvvvkvNv (2-49)(2-50)单元能量泛函的一阶偏导数即梯度单元能量泛函的一阶偏导数即梯度 接触表面单元摩擦功率的一阶偏导数： 单元能量泛函的一阶偏导数即梯度：22221111() sec 33egfTssfffeeiiffvkkmmJJNvv 42111 sec3eefePefeeTsPffiivvmJKvJN(2-51)(2-52)2.7 能量泛函的二阶偏导数能量泛函的二阶偏导数(Hessian矩阵矩阵) 1223 (1) (1) (1) nme

52、TeTnmeTePeTeTePavvmavKvmvvmKv 2.7 能量泛函的二阶偏导数能量泛函的二阶偏导数(Hessian矩阵矩阵)412421324221 (1) (1) eePPeieePPPeeTeTiTePeTTeeePPPPeeTiJKvvJKJKvvvvmKvvKvKvJKmvv 由于 (2-53)2.8 总体能量泛函能量泛函的梯度矩阵和总体能量泛函能量泛函的梯度矩阵和Hessian矩阵矩阵 根据迭加原理对所有单元的梯度矩阵和Hessian矩阵求和，可得到总体能量泛函的梯度和Hessian矩阵： 总体能量泛函的梯度和Hessian矩阵中的元素是按节点未知速度分量的顺序进行排列的，

53、因此，上述求和需要按整体节点编号进行。1221nmellnmell(2-54)3. 刚塑性有限元基本问题及其处理方法刚塑性有限元基本问题及其处理方法 l单元节点编号单元节点编号l初始速度场的设定初始速度场的设定l奇异点的处理奇异点的处理l刚塑性有限元分析过程的收敛准则l刚性区与塑性区分界的处理l轧件侧面形状和宽展的确定l接触表面速度边界条件的处理3.1 单元编号和节点总体编号单元编号和节点总体编号 刚塑性有限元在求解过程中，总体刚度矩阵即Hessian矩阵矩阵是采用一维压缩存贮一维压缩存贮方式进行存贮的，该矩阵一般称之为优化矩阵优化矩阵。在有限元网格划分之后，单元编号和节点总体编号虽然不能改变

54、总体刚度矩阵的规模、节点未知速度的总数，但是却能改变总体刚度矩阵的带宽而影响优化矩带宽而影响优化矩阵的规模阵的规模，从而给计算机存贮容量和计算时间带来很大的影响。单元总体节点编号的原则单元总体节点编号的原则 同一单元最大节点编号与单元最小节点编号的差值最小，使总体刚度矩阵的带宽最小。 较好节点编号方法是从坐标原点开始，首先沿单元数目最少的方向增加单元和节点的编号，然后沿单元次少的方向进行编号，最后沿单元数目最多的方向进行编号。minminmax NN(3-1)xiy21i+1i+22inp- i +1npxyi212ii+2i+1nm3.2 初始速度场的设定初始速度场的设定 设定轧制过程初始速度场的方法有以下几种：设定轧制过程初始速度场的方法有以下几种： 1.初等方法设定的均匀初始速度场；初等方法设定的均匀初始速度场； 2.用逐步