材料成形有限元法基础：材料成形有限元法基础

1、材料成形有限元法基础教材董湘怀. 材料成形计算机模拟. 北京：机械工业出版社，2001. 讲课：12学时主要参考书王勖成，邵 敏. 有限单元法基本原理与数值方法. 北京：清华大学出版社，1996. R.D.库克著，程耿东等译. 有限元分析的概念和应用科学出版社.考试与成绩课堂作业与期中闭卷考试课堂作业: 20分（4次，每次5分）期中闭卷考试：80分有限元法基础有限元发展过程60年代 70年代 80年代 90年代 2000年以后有限元应用制造业/机械/电子/航空航天有限元发展方向大规模/高精度/高速度/智能化有限元法的基本思想基本思想 1）将连续的求解系统离散为一组由节点相互联在一起的单元组合体

2、 2）在每个单元内假设近似函数来分片表示系统的求解场函数 有限元法的基本思想有限元法的基本思想有限元法的基本思想有限元法的基本思想有限元法的基本思想离散为单元网格的冲压件仍然要保证是一个连续体，单元与单元之间没有裂缝、不能重叠，所有单元通过单元节点相互关联着板料无论产生多大的塑性变形，单元与单元之间依然不会产生裂缝、交叉和重叠，关联单元的节点也不能脱开有限元法的基本思想不合格单元单元裂缝单元重叠有限元法的基本思想变形前后单元之间都是连续的变形前的网格变形后的网格有限元法的基本思想基本思想通过在单元内假设不同的插值函数，建立不同的单元模型，适应各种各样的变形模式和受力模式XFXF有限元法的基本思

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

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

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

6、六面体线性单元20节点六面体二次单元常用单元模型准三维空间单元桁架单元一维2节点线单元+单元局部随体坐标系 为什么要建立单元局部随体坐标系 ？简化分析问题的复杂程度。在局部坐标系中，空间桁架的每根杆每变成了一维2节点线单元常用单元模型准三维空间单元框架单元三维梁单元+一维2节点线单元+单元局部随体坐标系 两端都是刚性联结 可以要承受拉压、弯曲、扭转3种变形模式 框架单元的特点常用单元模型准三维空间单元板单元薄板单元中厚板单元弯曲和横向剪切2种变形模式抵抗板的变形如果板很薄，忽略横向剪切抗力，认为抵抗载荷的主要因素是弯矩常用单元模型准三维空间单元壳单元 抵抗拉压变形的二维单元+板单元+单元局部随

7、体坐标系。适合于薄壳单元和中厚壳单元从几何上分为薄壳单元和中厚壳单元 组合单元常用单元模型准三维空间单元 壳理论单元 由空间壳理论严格构造的壳单元。适合于薄壳单元和中厚壳单元 退化单元 由三维实体单元退化成的壳单元。只适合于中厚壳单元 单元模型构造 有限元法的基本思想 通过单元分片近似，在每个单元内假设近似函数来分片表示系统的场函数 选择近似函数简单、实用的原则在有限元法中，近似函数称为插值函数 单元模型构造插值函数 一般都采用多项式函数，主要原因是: 采用多项式插值函数比较容易推导单元平衡方程，特别是易于进行微分和积分运算。随着多项式函数阶次的增加，可以提高有限元法的计算精度。从理论上说，无

8、限提高多项式的阶数，可以求得系统的精确解。单元模型构造方法 整体坐标系法局部坐标系法 Lagrange插值方法Hermite插值方法单元模型构造方法2节点线单元12 oxu1u2x1x2ux1. 假设插值多项式2. 利用节点值求 a0 和 a1 单元模型构造方法3. 代入a0 和 a1，得插值多项式 u(x)4. 按u1 和 u2合并同类项，设 l = x2- x1单元模型构造方法关键 如何构造插值多项式 u ?二维问题三维问题，如何构造插值多项式?收敛性条件 在单元内，场函数必须是连续的； 当单元无限缩小时，插值多项式应该能够描述场函数及其各阶导数（直至最高阶导数）的均匀状态； 场函数及其偏

9、导数（直至比最高阶导数低一阶的各阶导数）在各单元边界必须连续。 插值多项式收敛性条件 收敛：当单元逐渐缩小时，如果插值多项式满足收敛性条件，则数值解将收敛于精确解 插值多项式收敛性条件多项式形式本身就是连续的，因此是自动满足的。均匀状态是场函数最基本的变化状态，插值多项式应该能够描述这种状态。类似地，当单元逐渐缩小时，场函数的各阶导数在单元内将趋于常数，因此也要求插值多项式具备这种表达能力。从固体力学角度来看，插值多项式零阶导数所描述的场函数均匀状态的物理意义就是刚体位移，一阶导数的均匀状态对应的是常应变状态。插值多项式收敛性条件协调单元 满足插值多项式收敛性条件和的单元 完备单元 满足插值多

10、项式收敛性条件的单元cr 阶连续性 场函数的第r阶导数是连续的 插值多项式收敛性条件非协调单元与部分协调单元 对于一般固体力学问题来说，协调性要求单元在变形时，相邻单元之间不应引起开裂、重叠或其它不连续现象。例如，梁、板、壳等单元，在单元边界不但要求位移是连续的，而且其一阶导数也必须是连续的。板、壳单元位移函数沿单元边界的法向导数（转角）的连续性一般比较难实现，因此出现了许多不完全满足协调性要求的“非协调单元”或“部分协调单元”，有时它们的精度也很好。 插值多项式选择条件 插值多项式应该尽可能满足其收敛性条件（收敛性）由插值多项式所确定的场函数变化应该与局部坐标系的选择无关（各向同性） 假设的

11、插值多项式系数的数量应该等于单元的节点数（解的唯一性） 选择条件插值多项式选择条件深入分析由收敛性条件可知，插值多项式中必须含有常数项（刚体位移项），高阶项的次数必须依次增加，不允许有跳跃插值多项式选择条件由选择条件可知，插值多项式函数在所有自由度方向上要满足各向同性性，这样就不会随局部坐标系变化而改变了 深入分析插值多项式选择条件深入分析选择条件是为了能由单元节点值唯一确定插值多项式 4节点四边形的插值多项式应该是 插值多项式系数i (i = 0,1,2,3) 也是4个 单元模型构造整体坐标系法基本思想 针对弹性体有限元网格建立一个统一的坐标系，每个单元的插值多项式都在这个坐标系上建立 单元

12、模型构造整体坐标系法2节点线单元12 oxu1u2x1x2ux1. 假设插值多项式2. 利用节点值求 a0 和 a1 单元模型构造整体坐标系法3. 代入a0 和 a1，得插值多项式 u(x)4. 按u1 和 u2合并同类项，设 l = x2- x1单元模型构造整体坐标系法N1 和 N2 称为单元的形函数；N 称为单元的形函数矩阵；ue 称为单元节点位移向量。 2节点线的单元形函数单元模型构造整体坐标系法二维3节点三角形单元 建立整体坐标系oxy 单元模型构造整体坐标系法1. 假设插值多项式2. 首先，利用节点值求 0 、 1 和 2 二维3节点三角形单元 单元模型构造整体坐标系法A为单元面积单

13、元模型构造整体坐标系法3. 将 0 、 1 和 2 代入插值多项式，按u1、u2、u3合并同类项单元模型构造整体坐标系法4. 同理可得单元模型构造整体坐标系法5. 单元插值多项式为单元模型构造整体坐标系法6. 单元插值多项式写成矩阵形式（常用）单元模型构造整体坐标系法7. 单元插值多项式的另一种矩阵形式（不常用）单元模型构造整体坐标系法4节点四面体单元单元模型构造整体坐标系法1. 假设插值多项式2. 插值多项式为单元模型构造整体坐标系法（i=1,2,3,4） 循环轮换脚标1、2、3、4，相应可以得到a2，b2 ， c2 ， d2 、 a3 ， b3 ， c3 ， d3 、 a4 ， b4 ，

14、c4 ， d4 单元模型构造整体坐标系法3. 单元插值多项式写成矩阵形式（常用）单元模型构造整体坐标系法4. 单元插值多项式另一种矩阵形式（不常用）单元模型构造整体坐标系法从理论上讲，整体坐标系法可以求任意单元的形函数，但计算过程太复杂只能求一维2节点线单元、二维3节点三角形单元和三维4节点四面体单元3种简单单元的形函数复杂的或二次以上的单元必须采用局部坐标系法求位移场 u 是形函数 Ni 的线性组合，因此形函数Ni同样具有插值多项式的特性单元刚度矩阵2节点线单元一维2节点线单元单元插值关系 单元几何关系 单元本构关系 N=N1 N2 De=E单元刚度矩阵2节点线单元单元刚度矩阵A为单元截面积

15、；l为单元长度矩阵B单元刚度矩阵三角形单元二维3角形单元单元插值关系 单元刚度矩阵三角形单元单元几何关系 单元刚度矩阵三角形单元单元本构关系 平面应力问题单元刚度矩阵三角形单元矩阵B单元刚度矩阵三角形单元单元刚度矩阵h为单元厚度k为对称的6*6常数矩阵A为单元面积单元模型构造整体坐标系法单元形函数的特性正规性：单元形函数之和等于1。 正交性：形函数在本节点的值等于1，在其它节点的值等于0。 单元模型构造局部坐标系法为什么要采用局部坐标系法?在单元上假设一种局部坐标系，确定局部坐标系的度量，并在单元节点上标出局部坐标值根据插值多项式选择条件假设形函数多项式利用单元形函数的特性(正交性)求单元的形

16、函数一般需要如下几个过程 单元模型构造局部坐标系法一维单元长度坐标 定义自然坐标L1、L2L1、L2也叫长度坐标 单元模型构造局部坐标系法长度坐标L1、L2特性L1和L2不是相互独立的 L1和L2与形函数的关系 长度坐标系与整体坐标系之间的变换关系 单元模型构造局部坐标系法一维二次单元 建立长度坐标系的目的是为了求一维高阶单元的形函数 一维3节点单元 单元模型构造局部坐标系法假设单元形函数Ni为 （i=1,2,3） 根据Ni的正交性，N1在1节点处的值等于1，在2、3节点处的值等于0 单元模型构造局部坐标系法得到方程组 得到形函数N1 同理得到N2和N3 单元模型构造局部坐标系法说明计算过程表

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

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

19、L3 、L4特性L1 、 L2 、 L3和L4不是相互独立的 L1 、 L2 、 L3和L4与形函数的关系 单元模型构造局部坐标系法体积坐标系与整体坐标系之间的转换关系 单元模型构造局部坐标系法三维10节点四面体单元 标注节点面积坐标单元模型构造局部坐标系法假设单元形函数根据Ni的正交性，求得单元形函数分别为 单元模型构造局部坐标系法小结长度坐标只适合求一维线单元的形函数 2节点线性单元、3节点二次单元面积坐标只适合求二维三角形单元的形函数 3节点线性单元、6节点二次单元、10节点二次单元体积坐标只适合求三维四面体单元的形函数 4节点线性单元、10节点二次单元单元模型构造局部坐标系法正规自然坐

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

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

22、一点的位移u与单元节点位移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利用等参单元关系 单元刚度矩阵四边形单元雅可比矩阵J中的

24、4个元素分别表示为 单元刚度矩阵四边形单元B进一步表示为 单元刚度矩阵由于bi和ci是关于自然坐标r、s的函数，所以矩阵B和k也是它们的函数，一般要采用高斯积分方法计算单刚k 单元平衡方程组装过程 为什么要组装 ? 消除内力组装的原则是什么 ? 单元自由度与结构自由度对应单元平衡方程组装过程 2 F 1 3 U3U4U2U1U5U6结构自由度向量U单元平衡方程组装过程 2 1 U3U4U2U11u1u2u3u4 3 U611 U2U12u1u2U5u3u42单元平衡方程组装过程2 1 U3U4U2U11u1u2u3u42单元平衡方程组装过程组装单元单元平衡方程组装过程 3 U611 U2U12

25、u1u2U5u3u4单元平衡方程组装过程再组装单元总体刚度方程 K 称为总体刚度矩阵 U 称为位移向量 F 称为载荷向量 总体刚度矩阵K的特性 对称性 奇异性 稀疏性 非零元素带状分布 约束处理过程 为什么要约束处理 ?总体平衡方程组是奇异的消除无限制的刚体运动 使总体平衡方程组存在唯一一组解约束处理过程边界条件 边界条件分类 力(载荷)边界条件位移边界条件 集中载荷力 表面分布力 自重力热交换引起的温度载荷 固定位移约束 强制位移约束 关联位移约束 约束处理过程模型简化xy约束处理过程模型简化yxxy约束处理过程约束方程123456789101112yx约束处理过程约束处理方法位移约束处理方

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

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

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

29、矩阵中要约束自由度的行列分别赋0，对角线元素赋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法，乘大数很少被采用主要原因是它是一种近似方法，而且大数的大小也不好确定，有时还会造成求解失败 约束处理过程弹簧单元假设柔性弹簧 kOXYU4f f = kU4k约束处理过程弹簧单元弹簧约束方程 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种单元的单元刚度矩阵是常数矩阵，不需要再进行数值积分运算。 除了这3种单

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

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

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

41、54节点线单元：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个积分点 但是，这并不意味着单元的精度就比较高，因为单元的精度是由插值多项式本身决定的 数值积分积分阶次选择四边形单元 4节点四边形单元：m=2，n=1.58节点四边形单元：m

43、=4，n=2.512节点四边形单元：m=6，n=3.5按经验公式计算：4节点四边形单元可以取 n=1 或 n=28节点四边形单元可以取 n=2 或 n=312节点四边形单元可以取 n=3 或 n=4实际应用计算：数值积分积分阶次选择这些单元在数值积分时，同样会象一维单元一样，出现零能模式或闭锁现象 完全积分方案 4节点取 n=18节点取 n=212节点取 n=3减缩积分方案 4节点取 n=28节点取 n=312节点取 n=4数值积分积分阶次选择四面体单元四面体单元与三角形单元相似，它们都没有减缩积分方案。4节点单元取 n=1，但计算结果偏硬，有时会产生闭锁现象。其它高阶单元积分阶次是比较精确的

44、，但是单元的计算精度都比较低 数值积分积分阶次选择六面体单元 六面体单元与四边形单元相似 8节点六面体单元取减缩积分方案 n=1 或完全积分方案 n=220节点六面体单元取减缩积分方案 n=2或完全积分方案 n=3有限元基本方法总结位移法基本过程 1）离散化过程 3）约束处理过程 2）单元平衡方程组装过程 5）应变、应力回代过程 4）方程组求解过程 有限元基本方法总结 为弹性体的应变 为弹性体的应力 u为弹性体的可容位移 弹性体处于平衡状态时，其势能应为最小 P 为表面分布力 G 为弹性体的体力 A 为弹性体的表面积 V 为弹性体的体积 有限元基本方法总结单元插值关系 单元几何关系 单元本构关

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

46、单元边界必须连续。 有限元基本方法总结均匀状态是场函数最基本的变化状态，当单元逐渐缩小时，场函数的各阶导数在单元内将趋于常数，因此也要求插值多项式具备这种表达能力。从固体力学角度来看，插值多项式零阶导数所描述的场函数均匀状态的物理意义就是刚体位移，一阶导数的均匀状态对应的是常应变状态。有限元基本方法总结协调单元 满足插值多项式收敛性条件和的单元 完备单元 满足插值多项式收敛性条件的单元cr 阶连续性 场函数的第r阶导数是连续的 有限元基本方法总结插值多项式应该尽可能满足其收敛性条件（收敛性）由插值多项式所确定的场函数变化应该与局部坐标系的选择无关（各向同性） 假设的插值多项式系数的数量应该等于

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

48、法总结单元形函数的特性正规性：单元形函数之和等于1。 正交性：形函数在本节点的值等于1，在其它节点的值等于0。 有限元基本方法总结局部坐标系法 长度坐标 一维线单元 面积坐标 二维三角形单元 体积坐标 三维四面体单元 正规的自然坐标 一维线单元、二维四边形单元、 三维六面体单元有限元基本方法总结等参单元 单元内任意一点的位移u与单元节点位移ue之间的关系为 一般单元坐标的插值关系也采用与位移插值关系相同的变换关系即单元内任意一点的坐标x与单元节点坐标xe之间的关系为 有限元基本方法总结等参单元凡是几何形状和位移场采用同阶同参数插值关系来描述的单元，称为等参单元 前面介绍的所有单元都属于等参单元

49、 在描述单元的几何形状和位移场时，并不一定非采用同阶插值关系 有限元基本方法总结对称性 奇异性 稀疏性 非零元素带状分布 有限元基本方法总结约束处理 总体平衡方程组是奇异的消除无限制的刚体运动 使总体平衡方程组存在唯一一组解有限元基本方法总结位移约束处理方法 赋0赋1法 乘大数法 有限元基本方法总结赋0赋1法的基本原理 利用初等变换对求解方程组进行相同的行列变换，既保证方程组解不会改变，又可以保持方程组系数矩阵的对称性。 在进行初等变换时，只要保证对方程组系数矩阵做相同的行列变换，就可以保持方程组系数矩阵的对称性。 有限元基本方法总结乘大数法的基本原理 利用矩阵的初等变换不改变方程组解的思想

50、放大位移约束方程的优势 大数A是系数矩阵中对角线元素的1010倍量级以上有限元基本方法总结赋0赋1法在约束处理过程中是严格精确的，而乘大数法是一种近似约束处理方法，它的精度取决于所乘大数A值 两种方法都可以消除有限元平衡方程的奇异性，得到符合实际边界条件的唯一一组解。但两种方法还是有很大的区别 有限元基本方法总结方程组求解有限元计算过程中很重要的一部分，在有限元法的发展过程中，有限元方程的求解效率一直是其应用的最大瓶颈之一 有限元方程组的特点： 有限元方程组的系数矩阵具有对称、稀疏、带状分布以及正定、主元占优。有效地利用这些特点，以减少系数矩阵的存贮量，提高方程组求解效率 有限元基本方法总结有

51、限元方程组系数矩阵是稀疏的、非零元素呈带状分布，带宽就是它的宽度，带宽的大小是由系统有限元网格的节点号排序决定的，具体求法是 带宽=(单元最大节点号之差+1)*节点自由度数 带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相关很大 有限元基本方法总结系数矩阵存贮 如果节点号排序优化的比较好，系数矩阵的存贮量就会减少很多。根据系数矩阵的对称性，一般都是按半带宽存贮。系数矩阵存贮的方法 二维等带宽存贮 一维变带宽存贮 有限元基本方法总结常用的单元面内数值积分方法主要有： Hammer积分 Gauss积分 有限元基本方法总结选择积分阶次的原则： 积分精度 积分阶次n与被积分多项式的阶次m有直接

52、关系。一般来说，有限元应用的经验公式 积分项有两个应变矩阵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之间坐标转换矩阵 为 板料成形有限元法薄膜单元单元局部坐标自由度向量ue与整体坐标自由度向量Ue的变换关系为 板料成形有限元法薄膜单元四边形薄膜单元 板料成形有限元法薄膜单元局部坐标系以1节点为坐标原点，x轴与单元1、2边重合并指向2节点，z轴与单元法向量nz平行。单位法向量n

55、z为 板料成形有限元法薄膜单元随体局部坐标系oxyz与空间整体坐标系OXYZ之间坐标转换矩阵 为 板料成形有限元法薄膜单元单元局部坐标自由度向量ue与整体坐标自由度向量Ue的变换关系为 板料成形有限元法单元平衡方程单元平衡方程单元插值关系(局部坐标系) 单元几何关系(局部坐标系) 单元本构关系(局部坐标系) 最小势能原理(局部坐标系) 坐标变换关系 板料成形有限元法单元平衡方程代入局部坐标系的单元3个关系式, 得代入坐标转换关系式, 得板料成形有限元法单元平衡方程ke 为局部坐标系的单元刚度矩阵 f e 为局部坐标系的单元载荷向量 Ke 为整体坐标系的单元刚度矩阵 F e 为整体坐标系的单元载

56、荷向量 板料成形有限元法本构方程 在金属塑性大变形有限元分析时经常采用流动理论本构方程，其它本构方程很少采用基于形变理论或非经典的角点理论本构方程虽然可以比较准确板料失稳后的局部化变形过程，但是板料成形属于强约束过程，对角点本构方程不敏感，而且板料成形也并不十分关心板料失稳后的局部化变形过程 目前板料成形有限元模拟的精度也不高，也没有必要十分强调本构方程的影响 板料成形有限元法本构方程考虑具有光滑屈服面屈服函数的弹塑性体，假设温度对变形速度的影响很小，可以忽略不计。这样全应变率可以分解为弹性应变率和塑性应变率之和，即 板料成形有限元法本构方程采用Hooke定律，弹性应变率 Piola Kirc

57、hhoff应力的Jaumann速率 用第二表示为 塑性应变率用流动法则和屈服函数 f 表示为板料成形有限元法本构方程nij是屈服面（应力空间f = 0的曲面）的单位法向量 = 0，应力处于弹性状态, 应力点位于屈服面以内 = 1，应力处于塑性加载状态,应力点位于屈服面上 h 表示当前状态的加工硬化率 板料成形有限元法本构方程式中应变率板料成形有限元法本构方程由 与Cauchy应力的Jaumann速率的关系 如果材料不可压缩 材料本构方程为 板料成形有限元法J2流动理论 J2流动理论 Mises屈服函数为 表示应力的偏量板料成形有限元法J2流动理论材料的J2流动本构关系 式中 h 可以由单向拉伸

58、实验确定 板料成形有限元法J2流动理论E 为弹性模量Et 为单向拉伸真应力对数应变曲线的切线模量 K 为材料强化系数; n 为材料强化指数板料成形有限元法屈服函数Hill正交各向异性函数 一般的，若把各向异性主轴作为随体坐标系的x,y,z轴，则Hill屈服函数可以表示成 F、G、H、L、M、N为各向异性参数，由实验确定 板料成形有限元法屈服函数在平面应力状态下 Hill各向异性函数 板料成形有限元法迭代解法 非线性方程组迭代解法 板材冲压成形数值模拟是一个强非线性问题，涉及到几何、材料和边界三重非线性。如果采用隐式有限元法就要求解非线性有限元方程组。非线性方程组一般是采用线性化方法，通过一系列

59、线性解逼近非线性解。但是这种方法是有局限性的，而且有时解的漂移误差很大。因此，一般都采用迭代法求解非线性有限元方程组，即Newton-Raphson法 板料成形有限元法迭代解法Newton-Raphson法 对具有一阶导数连续的函数F(u)在Un处作一阶Taylor展开，并用un表Un，则它在Un处的线性近似公式为 板料成形有限元法迭代解法非线性方程组在Un附近的近似方程组是一个线性方程组，即假设因此 Newton-Raphson法的迭代方程为 板料成形有限元法迭代解法修正的Newton-Raphson法 在Newton-Raphson法中，刚度矩阵是与un有关的，在每个迭代步都要重新计算一次

60、。为了减少计算量，用某一不变的刚度矩阵代替，得到修正的Newton-Raphson法迭代方程 板料成形有限元法迭代解法拟Newton-Raphson法 Newton-Raphson法和修正的Newton-Raphson法都是用切线刚度矩阵进行迭代平衡的。实际计算表明，这种直接迭代法不仅计算量大，而且常常不收敛。因此又提出一种拟Newton-Raphson法，用割线刚度矩阵进行迭代平衡。 板料成形有限元法收敛准则 迭代收敛准则 位移收敛准则 常用位移收敛准则为 当系统含有刚体位移时， 会比较大，此时不适合采用位移收敛准则 板料成形有限元法收敛准则常用失衡力收敛准则为 为第n迭代步的失衡力 失衡力