材料加工过程中的数值模拟.ppt

1、材料加工过程中的数值模拟,有限元法与板料成形模拟,柳玉起,华中科技大学模具技术国家重点实验室,有限元法基础,有限元发展过程 60年代 70年代 80年代 90年代 2000年以后 有限元应用 制造业/机械/电子/航空航天 有限元发展方向 大规模/高精度/高速度/智能化,有限元法的基本思想,基本思想,1）将连续的求解系统离散为一组由节点相互联在一起的单元组合体,2）在每个单元内假设近似函数来分片表示系统的求解场函数,有限元法的基本思想,有限元法的基本思想,有限元法的基本思想,有限元法的基本思想,有限元法的基本思想,基本思想 通过在单元内假设不同的插值函数，建立不同的单元模型，适应各种各样的变形模

2、式和受力模式,有限元在材料成形中的应用,有限元法的基本思想,有限元法分类,1）位移法：基于最小势能原理或虚功原理,2）力法： 基于最小余能原理,3）杂交法：基于修正余能原理,4）混合法：基于Reissner变分原理,有限元法的基本思想,位移法基本过程,1）离散化过程,3）约束处理过程,2）单元平衡方程组装过程,5）应变、应力回代过程,4）方程组求解过程,离散化过程,最小势能原理,弹性体的势能,为弹性体变形后所具有的内能,为弹性体所受的外力功,离散化过程, 为弹性体的应变, 为弹性体的应力,u为弹性体的可容位移,弹性体处于平衡状态时，其势能应为最小,0,离散化过程,单元插值关系,单元几何关系,单

3、元本构关系,N为单元形函数矩阵,L为单元几何微分算子,为单元弹性矩阵,单元节点自由度向量,离散化过程,B 称为应变矩阵,单元平衡方程或单元刚度方程,k 称为单元刚度矩阵,f 称为单元载荷向量,单元刚度矩阵的特性,对称性,奇异性,主元恒正且对角占优,离散化过程,线弹性问题几何方程三维问题,三维问题,线弹性问题几何方程二维问题,二维问题,平面应力和平面应变状态,线弹性问题几何方程二维问题,二维问题,轴对称状态,线弹性问题几何方程一维问题,一维问题,线弹性问题本构方程三维问题,三维问题,E为弹性模量；为泊松比,线弹性问题本构方程平面应力,二维问题,平面应力状态,线弹性问题本构方程平面应力,平面应力状

4、态,线弹性问题本构方程平面应变,二维问题,平面应变状态,线弹性问题本构方程平面应变,平面应变状态,线弹性问题本构方程轴对称,二维问题,轴对称状态,线弹性问题本构方程轴对称,二维问题,轴对称状态,线弹性问题本构方程轴对称,轴对称状态,线弹性问题本构方程一维问题,一维问题,常用单元模型,单元模型,插值关系,一一对应,单元类型,一维单元、二维单元、三维单元,等参单元、超参单元、次参单元,常用单元模型,一维单元,2节点线单元,3节点线单元,梁单元,常用单元模型,二维单元,3节点三角形线性单元,6节点三角形二次单元,常用单元模型,二维单元,10节点三角形三次单元,4节点四边形双线性单元,常用单元模型,二

5、维单元,8节点四边形二次单元,12节点四边形三次单元,常用单元模型,三维单元,4节点四面体线性单元,10节点四面体二次单元,常用单元模型,三维单元,8节点六面体线性单元,20节点六面体二次单元,常用单元模型,准三维空间单元,桁架单元,一维2节点线单元+单元局部随体坐标系,为什么要建立单元局部随体坐标系 ？,简化分析问题的复杂程度。 在局部坐标系中，空间桁架的每根杆每变成了一维2节点线单元,常用单元模型,准三维空间单元,框架单元,三维梁单元+一维2节点线单元+单元局部随体坐标系,两端都是刚性联结,可以要承受拉压、弯曲、扭转3种变形模式,框架单元的特点,常用单元模型,准三维空间单元,板单元,薄板单

6、元,中厚板单元,弯曲和横向剪切2种变形模式抵抗板的变形,如果板很薄，忽略横向剪切抗力，认为抵抗载荷的主要因素是弯矩,常用单元模型,准三维空间单元,壳单元,抵抗拉压变形的二维单元+板单元+单元局部随体坐标系。适合于薄壳单元和中厚壳单元,从几何上分为薄壳单元和中厚壳单元, 组合单元,常用单元模型,准三维空间单元, 壳理论单元,由空间壳理论严格构造的壳单元。适合于薄壳单元和中厚壳单元, 退化单元,由三维实体单元退化成的壳单元。只适合于中厚壳单元,单元模型构造,有限元法的基本思想,通过单元分片近似，在每个单元内假设近似函数来分片表示系统的场函数,选择近似函数,简单、实用的原则 在有限元法中，近似函数称

7、为插值函数,单元模型构造,插值函数,一般都采用多项式函数，主要原因是:,采用多项式插值函数比较容易推导单元平衡方程，特别是易于进行微分和积分运算。,随着多项式函数阶次的增加，可以提高有限元法的计算精度。从理论上说，无限提高多项式的阶数，可以求得系统的精确解。,单元模型构造方法,整体坐标系法 局部坐标系法,Lagrange插值方法 Hermite插值方法,单元模型构造方法,2节点线单元,1,2,u1,u2,x1,x2,u,x,1. 假设插值多项式,2. 利用节点值求 a0 和 a1,单元模型构造方法,3. 代入a0 和 a1，得插值多项式 u(x),4. 按u1 和 u2合并同类项，设 l =

8、x2- x1,单元模型构造方法,关键,如何构造插值多项式 u ?,二维问题三维问题，如何构造插值多项式?,收敛性条件 在单元内，场函数必须是连续的； 当单元无限缩小时，插值多项式应该能够描述场函数及其各阶导数（直至最高阶导数）的均匀状态； 场函数及其偏导数（直至比最高阶导数低一阶的各阶导数）在各单元边界必须连续。,插值多项式收敛性条件,收敛：当单元逐渐缩小时，如果插值多项式满足收敛性条件，则数值解将收敛于精确解,插值多项式收敛性条件,多项式形式本身就是连续的，因此是自动满足的。 均匀状态是场函数最基本的变化状态，插值多项式应该能够描述这种状态。类似地，当单元逐渐缩小时，场函数的各阶导数在单元内

9、将趋于常数，因此也要求插值多项式具备这种表达能力。 从固体力学角度来看，插值多项式零阶导数所描述的场函数均匀状态的物理意义就是刚体位移，一阶导数的均匀状态对应的是常应变状态。,插值多项式收敛性条件,协调单元 满足插值多项式收敛性条件和的单元 完备单元 满足插值多项式收敛性条件的单元 cr 阶连续性 场函数的第r阶导数是连续的,插值多项式收敛性条件,非协调单元与部分协调单元 对于一般固体力学问题来说，协调性要求单元在变形时，相邻单元之间不应引起开裂、重叠或其它不连续现象。例如，梁、板、壳等单元，在单元边界不但要求位移是连续的，而且其一阶导数也必须是连续的。板、壳单元位移函数沿单元边界的法向导数（

10、转角）的连续性一般比较难实现，因此出现了许多不完全满足协调性要求的“非协调单元”或“部分协调单元”，有时它们的精度也很好。,插值多项式选择条件,插值多项式应该尽可能满足其收敛性条件（收敛性） 由插值多项式所确定的场函数变化应该与局部坐标系的选择无关（各向同性） 假设的插值多项式系数的数量应该等于单元的节点数（解的唯一性）,选择条件,插值多项式选择条件,深入分析,由收敛性条件可知，插值多项式中必须含有常数项（刚体位移项），高阶项的次数必须依次增加，不允许有跳跃,插值多项式选择条件,由选择条件可知，插值多项式函数在所有自由度方向上要满足各向同性性，这样就不会随局部坐标系变化而改变了,深入分析,插值

11、多项式选择条件,深入分析,选择条件是为了能由单元节点值唯一确定插值多项式,4节点四边形的插值多项式应该是,插值多项式系数i (i = 0,1,2,3) 也是4个,单元模型构造整体坐标系法,基本思想,针对弹性体有限元网格建立一个统一的坐标系，每个单元的插值多项式都在这个坐标系上建立,单元模型构造整体坐标系法,2节点线单元,1,2,u1,u2,x1,x2,u,x,1. 假设插值多项式,2. 利用节点值求 a0 和 a1,单元模型构造整体坐标系法,3. 代入a0 和 a1，得插值多项式 u(x),4. 按u1 和 u2合并同类项，设 l = x2- x1,单元模型构造整体坐标系法,N1 和 N2 称

12、为单元的形函数； N 称为单元的形函数矩阵； ue 称为单元节点位移向量。,2节点线的单元形函数,单元模型构造整体坐标系法,二维3节点三角形单元,建立整体坐标系oxy,单元模型构造整体坐标系法,1. 假设插值多项式,2. 首先，利用节点值求 0 、 1 和 2,二维3节点三角形单元,单元模型构造整体坐标系法,A为单元面积,单元模型构造整体坐标系法,3. 将 0 、 1 和 2 代入插值多项式，按u1、u2、u3合并同类项,单元模型构造整体坐标系法,4. 同理可得,单元模型构造整体坐标系法,5. 单元插值多项式为,单元模型构造整体坐标系法,6. 单元插值多项式写成矩阵形式（常用）,单元模型构造整

13、体坐标系法,7. 单元插值多项式的另一种矩阵形式（不常用）,单元模型构造整体坐标系法,4节点四面体单元,单元模型构造整体坐标系法,1. 假设插值多项式,2. 插值多项式为,单元模型构造整体坐标系法,（i=1,2,3,4）,循环轮换脚标1、2、3、4，相应可以得到a2，b2 ， c2 ， d2 、 a3 ， b3 ， c3 ， d3 、 a4 ， b4 ， c4 ， d4,单元模型构造整体坐标系法,3. 单元插值多项式写成矩阵形式（常用）,单元模型构造整体坐标系法,4. 单元插值多项式另一种矩阵形式（不常用）,单元模型构造整体坐标系法,从理论上讲，整体坐标系法可以求任意单元的形函数，但计算过程太

14、复杂 只能求一维2节点线单元、二维3节点三角形单元和三维4节点四面体单元3种简单单元的形函数 复杂的或二次以上的单元必须采用局部坐标系法求 位移场 u 是形函数 Ni 的线性组合，因此形函数Ni同样具有插值多项式的特性,单元刚度矩阵2节点线单元,一维2节点线单元,单元插值关系,单元几何关系,单元本构关系,N=N1 N2,De=E,单元刚度矩阵2节点线单元,单元刚度矩阵,A为单元截面积；l为单元长度,矩阵B,单元刚度矩阵三角形单元,二维3角形单元,单元插值关系,单元刚度矩阵三角形单元,单元几何关系,单元刚度矩阵三角形单元,单元本构关系,平面应力问题,单元刚度矩阵三角形单元,矩阵B,单元刚度矩阵三

15、角形单元,单元刚度矩阵,h为单元厚度,k为对称的6*6常数矩阵,A为单元面积,单元模型构造整体坐标系法,单元形函数的特性,正规性：单元形函数之和等于1。,正交性：形函数在本节点的值等于1，在其它节点的值等于0。,单元模型构造局部坐标系法,为什么要采用局部坐标系法?,在单元上假设一种局部坐标系，确定局部坐标系的度量，并在单元节点上标出局部坐标值,根据插值多项式选择条件假设形函数多项式,利用单元形函数的特性(正交性)求单元的形函数,一般需要如下几个过程,单元模型构造局部坐标系法,一维单元长度坐标,定义自然坐标L1、L2,L1、L2也叫长度坐标,单元模型构造局部坐标系法,长度坐标L1、L2特性,L1

16、和L2不是相互独立的,L1和L2与形函数的关系,长度坐标系与整体坐标系之间的变换关系,单元模型构造局部坐标系法,一维二次单元,建立长度坐标系的目的是为了求一维高阶单元的形函数,一维3节点单元,单元模型构造局部坐标系法,假设单元形函数Ni为,（i=1,2,3）,根据Ni的正交性，N1在1节点处的值等于1，在2、3节点处的值等于0,单元模型构造局部坐标系法,得到方程组,得到形函数N1,同理得到N2和N3,单元模型构造局部坐标系法,说明,计算过程表明，采用长度坐标系求2次单元形函数比整体坐标系法要简单得多,由于长度坐标Li本身就含有常数项和一次项， Ni完全满足插值多项式选择条件求要,采用这种方法还

17、可以求更高阶单元的形函数,单元模型构造局部坐标系法,说明,由于长度坐标L1和L2不是相互独立的，形函数多项式的假设就会出现多种形式，只要它们满足插值多项式选择条件求要即可,（i=1,2,3）,（i=1,2,3）,（i=1,2,3）,单元模型构造局部坐标系法,一维三次单元,一维4节点单元,单元模型构造局部坐标系法,形函数形式假设为,根据形函数Ni的正交性，可分别求得单元的形函数,单元模型构造局部坐标系法,二维单元面积坐标,定义面积坐标L1、L2 、L3,单元模型构造局部坐标系法,面积坐标L1、L2 、L3特性,L1 、 L2和L3不是相互独立的,L1 、 L2和L3与形函数的关系,单元模型构造局

18、部坐标系法,面积坐标系与整体坐标系之间的转换关系,单元模型构造局部坐标系法,二维6节点三角形单元,标注节点面积坐标,单元模型构造局部坐标系法,假设单元形函数,（i=1,2,3,4,5,6）,根据Ni的正交性，求得单元形函数分别为,单元模型构造局部坐标系法,二维10节点三角形单元,标注节点面积坐标,单元模型构造局部坐标系法,假设单元形函数,根据Ni的正交性，求得单元形函数分别为,单元模型构造局部坐标系法,三维单元体积坐标,定义体积坐标Li,单元模型构造局部坐标系法,体积坐标L1、L2 、L3 、L4特性,L1 、 L2 、 L3和L4不是相互独立的,L1 、 L2 、 L3和L4与形函数的关系,

19、单元模型构造局部坐标系法,体积坐标系与整体坐标系之间的转换关系,单元模型构造局部坐标系法,三维10节点四面体单元,标注节点面积坐标,单元模型构造局部坐标系法,假设单元形函数,根据Ni的正交性，求得单元形函数分别为,单元模型构造局部坐标系法,小结 长度坐标只适合求一维线单元的形函数 2节点线性单元、3节点二次单元 面积坐标只适合求二维三角形单元的形函数 3节点线性单元、6节点二次单元、10节点二次单元 体积坐标只适合求三维四面体单元的形函数 4节点线性单元、10节点二次单元,单元模型构造局部坐标系法,正规自然坐标,求单元的形函数时，最常用的是采用正规自然坐标,正规自然坐标系是一种正规化的曲线坐标

20、系,单元模型构造局部坐标系法,一维情况,一维2节点线性单元,假设单元形函数,正规自然坐标,单元模型构造局部坐标系法,根据Ni的正交性,求得0和1,单元模型构造局部坐标系法,得形函数N1为,同理得形函数N2为,单元模型构造局部坐标系法,一维3节点二次单元,假设单元形函数,正规自然坐标,单元模型构造局部坐标系法,根据Ni的正交性，求得单元的形函数分别为,采用正规自然坐标方法可求更高阶单元的形函数,单元模型构造局部坐标系法,二维情况,二维4节点四边形单元,正规自然坐标,单元模型构造局部坐标系法,假设单元形函数,根据Ni的正交性，求得单元的形函数分别为,单元模型构造局部坐标系法,二维8节点四边形 二次

21、单元,正规自然坐标,单元模型构造局部坐标系法,假设单元形函数,根据Ni的正交性，求得单元的形函数分别为,单元模型构造局部坐标系法,二维12节点四边形 三次单元,正规自然坐标,单元模型构造局部坐标系法,假设单元形函数,根据Ni的正交性，求得单元的形函数,单元模型构造局部坐标系法,三维情况,三维8节点六面体单元,正规自然坐标,单元模型构造局部坐标系法,假设单元形函数,根据Ni的正交性，求得单元的形函数,单元模型构造局部坐标系法,三维20节点六面体 二次单元,正规自然坐标,单元模型构造局部坐标系法,假设单元形函数,根据Ni的正交性，求得单元的形函数,单元模型等参单元,等参单元,单元内任意一点的位移u

22、与单元节点位移ue之间的关系为,一般单元坐标的插值关系也采用与位移插值关系相同的变换关系即单元内任意一点的坐标x与单元节点坐标xe之间的关系为,单元模型等参单元,等参单元,凡是几何形状和位移场采用同阶同参数插值关系来描述的单元，称为等参单元,前面介绍的所有单元都属于等参单元,在描述单元的几何形状和位移场时，并不一定非采用同阶插值关系,单元模型等参单元,等参单元,3节点三角形等参单元,单元模型等参单元,超参单元,如果几何形状插值函数的阶数高于位移场插值函数的阶数，称为超参单元,次参单元,如果几何形状插值函数的阶数低于位移场插值函数的阶数，称为次参单元,单元刚度矩阵2节点线单元,单元插值关系,单元

23、几何关系,单元本构关系,N=N1 N2,De=E,一维2节点线单元,单元刚度矩阵2节点线单元,矩阵B,单元刚度矩阵2节点线单元,l为单元长度,单元刚度矩阵2节点线单元,单元刚度矩阵,A为单元截面积；l为单元长度,单元刚度矩阵四边形单元,二维4边形单元,单元插值关系,单元刚度矩阵四边形单元,单元几何关系,单元刚度矩阵四边形单元,单元本构关系,平面应力问题,单元刚度矩阵四边形单元,矩阵B,由于Ni是关于自然坐标r和s的函数，因此不能直接偏导,单元刚度矩阵四边形单元,将Ni(i=1,2,3,4)，看成是关于整体坐标x和y的隐函数，利用链锁规则求偏导,单元刚度矩阵四边形单元,雅可比矩阵 J,利用等参单

24、元关系,单元刚度矩阵四边形单元,雅可比矩阵J中的4个元素分别表示为,单元刚度矩阵四边形单元,B进一步表示为,单元刚度矩阵,由于bi和ci是关于自然坐标r、s的函数，所以矩阵B和k也是它们的函数，一般要采用高斯积分方法计算单刚k,单元平衡方程组装过程,为什么要组装 ?,消除内力,组装的原则是什么 ?,单元自由度与结构自由度对应,单元平衡方程组装过程,2,F,1,3,U3,U4,U2,U1,U5,U6,结构自由度向量U,单元平衡方程组装过程,2,单元平衡方程组装过程,2,单元平衡方程组装过程,单元平衡方程组装过程,单元平衡方程组装过程,总体刚度方程,K 称为总体刚度矩阵,U 称为位移向量,F 称为

25、载荷向量,总体刚度矩阵K的特性,对称性,奇异性,稀疏性,非零元素带状分布,约束处理过程,为什么要约束处理 ?,总体平衡方程组是奇异的,消除无限制的刚体运动,使总体平衡方程组存在唯一一组解,约束处理过程边界条件,边界条件分类,力(载荷)边界条件,位移边界条件,集中载荷力,表面分布力,自重力,热交换引起的温度载荷,固定位移约束,强制位移约束,关联位移约束,约束处理过程模型简化,约束处理过程模型简化,约束处理过程约束方程,1,2,3,4,5,6,7,8,9,10,11,12,约束处理过程约束处理方法,位移约束处理方法,赋0赋1法,乘大数法,约束处理过程赋0赋1法,位移约束处理例子,关联约束方程,约束

26、处理过程赋0赋1法,赋0赋1法,约束处理过程赋0赋1法,有6个方程，5个未知数，如果约束方程可以消除有限元平衡方程组的奇异性，则取任意5个方程联立求解，都会得到方程组的唯一一组解。,系数矩阵由原来的对称的变成了非对称的，这对于大规模有限元方程组求解是十分不利的，采用相同的求解方法，在求解时间和矩阵存贮容量方面都增加了一倍。,约束处理过程赋0赋1法,如何解决 ?,可以发现, 约束处理过程实际上是在系数矩阵上做了一次列初等变换。,为了保持平衡方程组系数矩阵的对称性，对变换后的系数矩阵再做一次相同的行初等变换。,约束处理过程赋0赋1法,具体做法： 第4行乘以系数k加到第3行，并去掉第4行。 为了保持

27、系数矩阵的阶数，将第4行的所有元素赋0，在其对角线元素位置赋1。 即所谓赋0赋1法。,约束处理过程赋0赋1法,经过初等变换，方程组的系数矩阵仍然保持对称性 初等变换不会改变方程组的解 约束后的方程组可以求得5个未知数 通过关联约束方程回代求解U4,约束处理过程赋0赋1法,进一步改进方法 将关联约束方程直接引进约束后的方程组）中，使U4也直接参与方程组求解。,用关联约束取代约束后方程组的赋0赋1行（第4行），然后再做对称化处理。即将取代后的一行方程（第4行）乘以k加到第3行，则系数矩阵仍然保持对称性。,约束处理过程赋0赋1法,约束处理过程赋0赋1法,基本原理 利用初等变换对求解方程组进行相同的行

28、列变换，既保证方程组解不会改变，又可以保持方程组系数矩阵的对称性。 在进行初等变换时，只要保证对方程组系数矩阵做相同的行列变换，就可以保持方程组系数矩阵的对称性。,约束处理过程赋0赋1法,固定和强制位移约束条件处理 通过关联位移约束方法简化。,如果k = 0，则退化成强制位移边界条件 即 U4=C,如果k =C= 0，则退化成固定位移边界条件 即 U4=0,约束处理过程赋0赋1法,强制位移约束条件处理 U4=C,约束处理过程赋0赋1法,固定位移约束条件处理 U4=0,约束处理过程赋0赋1法,固定位移和强制位移约束处理后的系数矩阵是相同的，只是简单地将方程组系数矩阵中要约束自由度的行列分别赋0，

29、对角线元素赋1，这也是赋0赋1法的由来。,约束处理过程赋0赋1法,在方程组载荷右端项的处理方法上两者是不同，处理固定位移边界条件时，只要将对应自由度的载荷赋0即可。处理强制位移边界条件时，要在方程组系数矩阵未赋0赋1法前，先将对应自由度的列乘以系数C减到载荷右端项，再将对应自由度的载荷位置赋C。,约束处理过程乘大数法,乘大数法,基本原理 利用矩阵的初等变换不改变方程组解的思想。,约束处理过程乘大数法,关联约束方程,约束处理过程乘大数法,关联约束方程,A是一个大数，是系数矩阵中对角线元素K44的1010倍量级以上,为什么要乘以大数A ？,放大位移约束方程的优势,约束处理过程乘大数法,关联约束方程

30、,将关联约束加到有限元平衡方程对应自由度行，第3行或第4行，这里取第4行,约束处理过程乘大数法,如果关联约束方程可以消除有限元平衡方程的奇异性，约束后的方程组就存在唯一的一组解。 约束后的方程组的系数矩阵是非对称的。 利用初等变换方法将系数矩阵变换成对称的,约束处理过程乘大数法,关联约束方程,再乘以系数-k,加到约束后方程组的第3行,约束处理过程乘大数法,约束处理过程乘大数法,强制位移边界条件,约束后的方程组简化为,约束处理过程乘大数法,固定位移边界条件,k = 0，C = 0,约束后的方程组简化为,约束处理过程乘大数法,固定位移和强制位移边界条件的乘大数约束处理相对比较简单，而且它们的系数矩

31、阵约束后是相同的，只是简单地将方程组系数矩阵中要约束自由度的对角线元素加上一个相对大数A即可 乘大数法的叫法并不十分准确，应该叫加大数法更贴切 乘大数和加大数的效果是一样的,约束处理过程两种方法比较,赋0赋1法在约束处理过程中是严格精确的，而乘大数法是一种近似约束处理方法，它的精度取决于所乘大数A值,两种方法都可以消除有限元平衡方程的奇异性，得到符合实际边界条件的唯一一组解。但两种方法还是有很大的区别,约束处理过程两种方法比较,采用乘大数法约束处理后的有限元平衡方程在求解时可能造成解的失真，大数A值越大可能解的偏差会越大，而赋0赋1法就不会出现类似的问题，它在约束过程和求解过程都是精确的 乘大

32、数法相对于赋0赋1法在约束处理过程上简单一些,约束处理过程两种方法比较,赋0赋1法实际上是将关联位移约束方程代入到有限元平衡方程中的，是代入法。而乘大数是将占绝对优势的关联位移约束方程合并到有限元平衡方程中的，是罚方法，计算误差来自于合并过程，计算精度取决于关联位移约束方程的优势大小 商业软件中，位移边界条件的约束处理都采用赋0赋1法，乘大数很少被采用主要原因是它是一种近似方法，而且大数的大小也不好确定，有时还会造成求解失败,约束处理过程弹簧单元,假设柔性弹簧 k,f = kU4,约束处理过程弹簧单元,弹簧约束方程 f = kU4,方程组求解过程特点,方程组求解有限元计算过程中很重要的一部分，

33、在有限元法的发展过程中，有限元方程的求解效率一直是其应用的最大瓶颈之一 有限元方程组的特点： 有限元方程组的系数矩阵具有对称、稀疏、带状分布以及正定、主元占优。有效地利用这些特点，以减少系数矩阵的存贮量，提高方程组求解效率,方程组求解过程分类比较,线性方程组的解法主要分两大类: 直接解法：以高斯消去法基础，以等带宽或变带宽方式存贮系数矩阵内元素，对于求解规模比较大的问题，要存贮的元素非常巨大。 迭代解法：只需要存贮系数矩阵中非零元素，存贮量很小，一般是变带宽存贮量的20%或更少，有些算法的求解效率也非常高，适合求解大规模线性方程组。但是这种解法对接近病态的方程组很难保证收敛性。,方程组求解过程

34、带宽定义,有限元方程组系数矩阵是稀疏的、非零元素呈带状分布，带宽就是它的宽度，带宽的大小是由系统有限元网格的节点号排序决定的，具体求法是,带宽=(单元最大节点号之差+1)*节点自由度数,带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相关很大,方程组求解过程带宽,带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相关很大,方程组求解过程带宽,所示四边形网格的三种节点号标注方法，每个节点是2个自由度 结构的带宽分别是12，18，56，相差很大，其中12和56之间相差近5倍，这就意味着系数矩阵的存贮量也是相差5倍，因此，对于大规模复杂系统的节点号优化是十分必要的,方程组求解过程系数矩阵存

35、贮,系数矩阵存贮 如果节点号排序优化的比较好，系数矩阵的存贮量就会减少很多。根据系数矩阵的对称性，一般都是按半带宽存贮。 系数矩阵存贮的方法 二维等带宽存贮 一维变带宽存贮,方程组求解过程二维等带宽存贮,二维等带宽存贮,方程组求解过程二维等带宽存贮,二维等带宽存贮消除了最大带宽以外的全部零元素，节省了系数矩阵元素的存贮量。但是由于取最大带宽为存贮范围，因此不能排除在带宽内的大量零元素。当系数矩阵的各行带宽变化不大时，适合采用二维等带宽存贮，方程组求解过程中系数矩阵元素的寻址也比较方便，求解效率较高。 当出现局部带宽特别大的情况时，采用二维等带宽存贮时，将由于局部带宽过大而使整体系数矩阵的存贮大

36、大增加。,方程组求解过程一维变带宽存贮,一维变带宽存贮 一维变带宽存贮方法就是把变化的带宽内的元素按一定的顺序存贮在一个一维数组中。由于它不按最大带宽存贮，因此比二维等带宽存贮更节省内存。按照解法可分为按行一维变带宽存贮和按列一维变带宽存贮。,按行一维变带宽存贮,方程组求解过程一维变带宽存贮,辅助的寻址数组M,一维变带宽存贮是最节省内存的一种方法，但是由于要借助于寻址数组寻找系数矩阵元素的位置，相对二维等带宽存贮方法来说要复杂一些，而且在程序实现时也要复杂得多，方程组求解过程中也要消耗一些数组寻址时间。因此，在选用存贮方法时要权衡二者的利弊，统盘考虑。一般当带宽变化不大，计算机内存允许时，采用

37、二维等带宽存贮方法是比较合适的。,方程组求解过程一维变带宽存贮,方程组求解过程求解方法,方程组求解方法 高斯消去法 三角分解法 雅可比（Jacobi）迭代法 高斯-赛德尔（Gauss-Seidel）迭代法 共轭梯度迭代法 预优共轭梯度迭代法,应变、应力回代过程,单元应变和应力回代求解 通过求解有限元平衡方程得到有限元节点位移后，就可以进行系统的刚度校核。如果所分析问题要进行强度校核，就要回代求解单元的应变和应力。 由插值关系和几何关系可得单元应变，再通过本构关系得到单元应力,数值积分,为什么要进行数值积分?,2节点线单元、3节点三角形单元和4节点四面体单元3种单元的单元刚度矩阵是常数矩阵，不需

38、要再进行数值积分运算。,除了这3种单元外，一般其它单元的刚度矩阵都是积分变量的函数，要采用数值积分方法进行计算。,数值积分主要方法,常用的单元面内数值积分方法主要有： Hammer积分 Gauss积分,数值积分Hammer积分,Hammer积分,在三角形单元和四面体单元中，自然坐标是面积坐标和体积坐标。采用这些坐标建立的单元形函数，其单元刚度矩阵的一般形式为,二维,数值积分Hammer积分,三维,一维单元在实际有限元应用中，一般都采用正规自然坐标系法建立单元的形函数。它采用Gauss积分方案,数值积分Hammer积分,Hammer积分运算,三角形单元的Hammer积分表示为,数值积分Hamme

39、r积分,四面体单元的Hammer积分表示为,数值积分Hammer积分,数值积分Hammer积分,数值积分Gauss 积分,Gauss 积分 采用正规自然坐标确定形函数的单元，其单元刚度矩阵的一般形式为,一维,二维,三维,数值积分Gauss 积分,这些积分形式在积分限上与高斯积分完全一致，因此高斯积分方案被广泛应用于此类积分形式中。它们的具体数值形式为,一维,二维,数值积分Gauss 积分,三维,数值积分Gauss 积分,数值积分积分阶次选择,数值积分的阶次选择,求解单元平衡方程时，绝大多数情况要采用数值积分方法，如何选择数值积分的阶次将直接影响计算精度和计算量。如果积分阶次选择不当，有时甚至会

40、导致计算失败,数值积分积分阶次选择,选择积分阶次的原则主要依据以下两点：,积分精度,积分阶次n与被积分多项式的阶次m有直接关系。一般来说，有限元应用的经验公式,积分项有两个应变矩阵B相乘，因此m一定是偶数，则积分阶数n等于0.5、1.5、2.5、,数值积分积分阶次选择,常用单元的积分阶次选择,一维单元,一般都采用正规自然坐标系法得到的形函数,在单元平衡方程中雅可比矩阵中虽然也含有自然坐标，但是它只是单刚的一个系数，只对单刚中的每个元素的大小有相同的影响，不会改变单刚的特性,数值积分积分阶次选择,2节点线单元只能取高斯积分点n=1 3节点线单元可以取n=1或n=2 4节点线单元可以取n=2或n=

41、3,2节点线单元：m=0，n=0.5 3节点线单元：m=2，n=1.5 4节点线单元：m=4，n=2.5,按经验公式计算：,实际应用计算：,数值积分积分阶次选择,在有限元法中，把3节点单元取n=1以及4节点单元取n=2的积分方案称为减缩积分，而3节点单元取n=2以及4节点单元取n=3的积分方案称为正常积分,实际数值结果表明，有时减缩积分方案会带来很大的计算误差，产生零能模式,正常积分方案的有时计算结果也会偏小，产生闭锁现象,数值积分积分阶次选择,造成这些现象的原因有很多，例如，单元形状、单元相对大小、单元受力状况、分析问题的类型等等。为了避免零能模式和闭锁现象的发生，一般采用减缩积分加阻尼矩阵

42、方法。采用减缩积分方案时，对每个节点施加一个柔性弹簧，通过弹簧的阻尼增加刚度矩阵的稳定性，阻止零能模式的发生。但是弹簧的刚性系数越大，计算误差就越大，因此弹簧系数的选择也有一定的困难,数值积分积分阶次选择,三角形单元,3节点三角形单元按经验公式计算的积分阶次n=0.5，实际计算时只能取n=1,结果造成计算结果偏硬，有时会产生闭锁现象，实际有限元计算时也证明了这一点,数值积分积分阶次选择,6节点三角形单元积分阶次是比较精确的,按经验公式计算的积分阶次应该取n=1.5，在单元面内应该是3个积分点,但是，这并不意味着单元的精度就比较高，因为单元的精度是由插值多项式本身决定的,数值积分积分阶次选择,四

43、边形单元,4节点四边形单元：m=2，n=1.5 8节点四边形单元：m=4，n=2.5 12节点四边形单元：m=6，n=3.5,按经验公式计算：,4节点四边形单元可以取 n=1 或 n=2 8节点四边形单元可以取 n=2 或 n=3 12节点四边形单元可以取 n=3 或 n=4,实际应用计算：,数值积分积分阶次选择,这些单元在数值积分时，同样会象一维单元一样，出现零能模式或闭锁现象,正常积分方案,4节点取 n=1 8节点取 n=2 12节点取 n=3,减缩积分方案,4节点取 n=2 8节点取 n=3 12节点取 n=4,数值积分积分阶次选择,四面体单元,四面体单元与三角形单元相似，它们都没有减缩

44、积分方案。4节点单元取 n=1，但计算结果偏硬，有时会产生闭锁现象。其它高阶单元积分阶次是比较精确的，但是单元的计算精度都比较低,数值积分积分阶次选择,六面体单元,六面体单元与四边形单元相似,8节点六面体单元取减缩积分方案 n=2 或正常积分方案 n=3,20节点六面体单元取减缩积分方案 n=3或正常积分方案 n=4,有限元基本方法总结,位移法基本过程,1）离散化过程,3）约束处理过程,2）单元平衡方程组装过程,5）应变、应力回代过程,4）方程组求解过程,有限元基本方法总结, 为弹性体的应变, 为弹性体的应力,u为弹性体的可容位移,弹性体处于平衡状态时，其势能应为最小,P 为表面分布力,G 为

45、弹性体的体力,A 为弹性体的表面积,V 为弹性体的体积,有限元基本方法总结,单元插值关系,单元几何关系,单元本构关系,N为单元形函数矩阵,L为单元几何微分算子,为单元弹性矩阵,单元节点自由度向量,有限元基本方法总结,B 称为应变矩阵,单元平衡方程或单元刚度方程,k 称为单元刚度矩阵,f 称为单元载荷向量,单元刚度矩阵的特性,对称性,奇异性,主元恒正且对角占优,有限元基本方法总结,单元模型构造,整体坐标系法 局部坐标系法 采用多项式函数 易于进行微分和积分运算 随着多项式函数阶次的增加，可以提高有限元法的计算精度,有限元基本方法总结,收敛性条件 在单元内，场函数必须是连续的； 当单元无限缩小时，

46、插值多项式应该能够描述场函数及其各阶导数（直至最高阶导数）的均匀状态； 场函数及其偏导数（直至比最高阶导数低一阶的各阶导数）在各单元边界必须连续。,有限元基本方法总结,均匀状态是场函数最基本的变化状态，当单元逐渐缩小时，场函数的各阶导数在单元内将趋于常数，因此也要求插值多项式具备这种表达能力。 从固体力学角度来看，插值多项式零阶导数所描述的场函数均匀状态的物理意义就是刚体位移，一阶导数的均匀状态对应的是常应变状态。,有限元基本方法总结,协调单元 满足插值多项式收敛性条件和的单元 完备单元 满足插值多项式收敛性条件的单元 cr 阶连续性 场函数的第r阶导数是连续的,有限元基本方法总结,插值多项式

47、应该尽可能满足其收敛性条件（收敛性） 由插值多项式所确定的场函数变化应该与局部坐标系的选择无关（各向同性） 假设的插值多项式系数的数量应该等于单元的节点数（解的唯一性）,选择条件,有限元基本方法总结,由收敛性条件可知，插值多项式中必须含有常数项（刚体位移项），高阶项的次数必须依次增加，不允许有跳跃 由选择条件可知，插值多项式函数在所有自由度方向上要满足各向同性性，这样就不会随局部坐标系变化而改变了 选择条件是为了能由单元节点值唯一确定插值多项式,有限元基本方法总结,从理论上讲，整体坐标系法可以求任意单元的形函数，但计算过程太复杂 只能求一维2节点线单元、二维3节点三角形单元和三维4节点四面体单

48、元3种线性单元的形函数 复杂的或二次以上的单元必须采用局部坐标系法求 位移场 u 是形函数 Ni 的线性组合，因此形函数Ni同样具有插值多项式的特性,有限元基本方法总结,单元形函数的特性,正规性：单元形函数之和等于1。,正交性：形函数在本节点的值等于1，在其它节点的值等于0。,有限元基本方法总结,局部坐标系法 长度坐标 一维线单元 面积坐标 二维三角形单元 体积坐标 三维四面体单元 正规的自然坐标 一维线单元、二维四边形单元、 三维六面体单元,有限元基本方法总结,等参单元,单元内任意一点的位移u与单元节点位移ue之间的关系为,一般单元坐标的插值关系也采用与位移插值关系相同的变换关系即单元内任意

49、一点的坐标x与单元节点坐标xe之间的关系为,有限元基本方法总结,等参单元,凡是几何形状和位移场采用同阶同参数插值关系来描述的单元，称为等参单元,前面介绍的所有单元都属于等参单元,在描述单元的几何形状和位移场时，并不一定非采用同阶插值关系,有限元基本方法总结,对称性,奇异性,稀疏性,非零元素带状分布,有限元基本方法总结,约束处理,总体平衡方程组是奇异的,消除无限制的刚体运动,使总体平衡方程组存在唯一一组解,有限元基本方法总结,位移约束处理方法,赋0赋1法,乘大数法,有限元基本方法总结,赋0赋1法的基本原理 利用初等变换对求解方程组进行相同的行列变换，既保证方程组解不会改变，又可以保持方程组系数矩

50、阵的对称性。 在进行初等变换时，只要保证对方程组系数矩阵做相同的行列变换，就可以保持方程组系数矩阵的对称性。,有限元基本方法总结,乘大数法的基本原理 利用矩阵的初等变换不改变方程组解的思想 放大位移约束方程的优势 大数A是系数矩阵中对角线元素的1010倍量级以上,有限元基本方法总结,赋0赋1法在约束处理过程中是严格精确的，而乘大数法是一种近似约束处理方法，它的精度取决于所乘大数A值,两种方法都可以消除有限元平衡方程的奇异性，得到符合实际边界条件的唯一一组解。但两种方法还是有很大的区别,有限元基本方法总结,方程组求解有限元计算过程中很重要的一部分，在有限元法的发展过程中，有限元方程的求解效率一直

51、是其应用的最大瓶颈之一 有限元方程组的特点： 有限元方程组的系数矩阵具有对称、稀疏、带状分布以及正定、主元占优。有效地利用这些特点，以减少系数矩阵的存贮量，提高方程组求解效率,有限元基本方法总结,有限元方程组系数矩阵是稀疏的、非零元素呈带状分布，带宽就是它的宽度，带宽的大小是由系统有限元网格的节点号排序决定的，具体求法是,带宽=(单元最大节点号之差+1)*节点自由度数,带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相关很大,有限元基本方法总结,系数矩阵存贮 如果节点号排序优化的比较好，系数矩阵的存贮量就会减少很多。根据系数矩阵的对称性，一般都是按半带宽存贮。 系数矩阵存贮的方法 二维等

52、带宽存贮 一维变带宽存贮,有限元基本方法总结,常用的单元面内数值积分方法主要有： Hammer积分 Gauss积分,有限元基本方法总结,选择积分阶次的原则：,积分精度,积分阶次n与被积分多项式的阶次m有直接关系。一般来说，有限元应用的经验公式,积分项有两个应变矩阵B相乘，因此m一定是偶数，则积分阶数n等于0.5、1.5、2.5、,有限元基本方法总结,减缩积分：低于经验公式的阶次 正常积分：高于经验公式的阶次,实际数值结果表明，有时减缩积分方案会带来很大的计算误差，产生零能模式,正常积分方案的有时计算结果也会偏小，产生闭锁现象,有限元基本方法总结,造成这些现象的原因有很多，例如，单元形状、单元相

53、对大小、单元受力状况、分析问题的类型等等。为了避免零能模式和闭锁现象的发生，一般采用减缩积分加阻尼矩阵方法。采用减缩积分方案时，对每个节点施加一个柔性弹簧，通过弹簧的阻尼增加刚度矩阵的稳定性，阻止零能模式的发生。但是弹簧的刚性系数越大，计算误差就越大，因此弹簧系数的选择也有一定的困难,有限元基本方法总结,采用乘大数法约束处理后的有限元平衡方程在求解时可能造成解的失真，大数A值越大可能解的偏差会越大，而赋0赋1法就不会出现类似的问题，它在约束过程和求解过程都是精确的 乘大数法相对于赋0赋1法在约束处理过程上简单一些,板料成形有限元法分类,静力隐式有限元法 静力显式有限元法 动力显式有限元法 逆算

54、法(一步成形),板料成形有限元法单元模型,单元模型 薄膜单元 薄壳单元 中厚壳单元 等效弯曲单元,板料成形有限元法薄膜单元,薄膜单元是由二维三角形单元或四边形单元构造的空间板壳单元 薄膜单元适用于板料在变形过程中主要以拉伸和压缩变形为主，局部弯曲变形对整个成形问题不产生大的影响 液压胀形、半球冲头胀形等一类问题,板料成形有限元法薄膜单元,三角形薄膜单元,板料成形有限元法薄膜单元,局部坐标系以1节点为坐标原点，x轴与单元1、2边重合并指向2节点，z轴与单元法向量nz平行。单位法向量nz为,板料成形有限元法薄膜单元,随体局部坐标系oxyz与空间整体坐标系OXYZ之间坐标转换矩阵 为,板料成形有限元

55、法薄膜单元,单元局部坐标自由度向量ue与整体坐标自由度向量Ue的变换关系为,板料成形有限元法薄膜单元,四边形薄膜单元,板料成形有限元法薄膜单元,局部坐标系以1节点为坐标原点，x轴与单元1、2边重合并指向2节点，z轴与单元法向量nz平行。单位法向量nz为,板料成形有限元法薄膜单元,随体局部坐标系oxyz与空间整体坐标系OXYZ之间坐标转换矩阵 为,板料成形有限元法薄膜单元,单元局部坐标自由度向量ue与整体坐标自由度向量Ue的变换关系为,板料成形有限元法单元平衡方程,单元平衡方程,单元插值关系(局部坐标系),单元几何关系(局部坐标系),单元本构关系(局部坐标系),最小势能原理(局部坐标系),坐标变

56、换关系,板料成形有限元法单元平衡方程,代入局部坐标系的单元3个关系式, 得,代入坐标转换关系式, 得,板料成形有限元法单元平衡方程,ke 为局部坐标系的单元刚度矩阵,f e 为局部坐标系的单元载荷向量,Ke 为整体坐标系的单元刚度矩阵,F e 为整体坐标系的单元载荷向量,板料成形有限元法本构方程,在金属塑性大变形有限元分析时经常采用流动理论本构方程，其它本构方程很少采用 基于形变理论或非经典的角点理论本构方程虽然可以比较准确板料失稳后的局部化变形过程，但是板料成形属于强约束过程，对角点本构方程不敏感，而且板料成形也并不十分关心板料失稳后的局部化变形过程 目前板料成形有限元模拟的精度也不高，也没

57、有必要十分强调本构方程的影响,板料成形有限元法本构方程,考虑具有光滑屈服面屈服函数的弹塑性体，假设温度对变形速度的影响很小，可以忽略不计。这样全应变率可以分解为弹性应变率和塑性应变率之和，即,板料成形有限元法本构方程,采用Hooke定律，弹性应变率,板料成形有限元法本构方程,nij是屈服面（应力空间f = 0的曲面）的单位法向量, = 0，应力处于弹性状态, 应力点位于屈服面以内 = 1，应力处于塑性加载状态,应力点位于屈服面上,h 表示当前状态的加工硬化率,板料成形有限元法本构方程,式中,应变率,板料成形有限元法本构方程,如果材料不可压缩,材料本构方程为,板料成形有限元法J2流动理论,J2流

58、动理论,Mises屈服函数为,板料成形有限元法J2流动理论,材料的J2流动本构关系,式中,h 可以由单向拉伸实验确定,板料成形有限元法J2流动理论,E 为弹性模量,Et 为单向拉伸真应力对数应变曲线的切线模量,K 为材料强化系数; n 为材料强化指数,板料成形有限元法屈服函数,Hill正交各向异性函数,一般的，若把各向异性主轴作为随体坐标系的x,y,z轴，则Hill屈服函数可以表示成,F、G、H、L、M、N为各向异性参数，由实验确定,板料成形有限元法屈服函数,在平面应力状态下 Hill各向异性函数,板料成形有限元法迭代解法,非线性方程组迭代解法 板材冲压成形数值模拟是一个强非线性问题，涉及到几何、材料和边界三重非线性。如果采用隐式有限元法就要求解非线性有限元方程组。非线性方程组一般是采用线性化方法，通过一系列线性解逼近非线性解。但是这种方法是有局限性的，而且有时解的漂移误差很大。因此，一般都采用迭代法求解非线性有限元方程组，即Newton-Raphson法,板料成形有限元法迭代解法,Newton-Raphson法,对具有一阶导数连续的函数F(u)在Un处作一阶Taylor展开，并用un表Un，则它在Un处的线性近似公式为,板料成形有限元法迭代解法,非线性方程组在Un附近的近似方程组是一个线性方程组，即,假设,因此 Newton-Raphson法的迭代方程为,板料成形有限元法迭代