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

《材料成形有限元法基础：材料成形有限元法基础》由会员分享，可在线阅读，更多相关《材料成形有限元法基础：材料成形有限元法基础（323页珍藏版）》请在金锄头文库上搜索。

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

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

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

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

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

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

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

8、插值函数比较容易推导单元平衡方程，特别是易于进行微分和积分运算。随着多项式函数阶次的增加，可以提高有限元法的计算精度。从理论上说，无限提高多项式的阶数，可以求得系统的精确解。单元模型构造方法 n整体坐标系法n局部坐标系法 nLagrange插值方法nHermite插值方法单元模型构造方法n2节点线单元12 oxu1u2x1x2ux1. 假设插值多项式2. 利用节点值求 a0 和 a1 单元模型构造方法3. 代入a0 和 a1，得插值多项式 u(x)4. 按u1 和 u2合并同类项，设 l = x2- x1单元模型构造方法n关键 如何构造插值多项式 u ?二维问题三维问题，如何构造插值多项式?n

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

10、均匀状态的物理意义就是刚体位移，一阶导数的均匀状态对应的是常应变状态。插值多项式收敛性条件n协调单元 满足插值多项式收敛性条件和的单元 n完备单元 满足插值多项式收敛性条件的单元ncr 阶连续性 场函数的第r阶导数是连续的 插值多项式收敛性条件n非协调单元与部分协调单元 对于一般固体力学问题来说，协调性要求单元在变形时，相邻单元之间不应引起开裂、重叠或其它不连续现象。例如，梁、板、壳等单元，在单元边界不但要求位移是连续的，而且其一阶导数也必须是连续的。板、壳单元位移函数沿单元边界的法向导数（转角）的连续性一般比较难实现，因此出现了许多不完全满足协调性要求的“非协调单元”或“部分协调单元”，有时

11、它们的精度也很好。 插值多项式选择条件 n插值多项式应该尽可能满足其收敛性条件（收敛性）n由插值多项式所确定的场函数变化应该与局部坐标系的选择无关（各向同性） n假设的插值多项式系数的数量应该等于单元的节点数（解的唯一性） 选择条件插值多项式选择条件n深入分析由收敛性条件可知，插值多项式中必须含有常数项（刚体位移项），高阶项的次数必须依次增加，不允许有跳跃插值多项式选择条件由选择条件可知，插值多项式函数在所有自由度方向上要满足各向同性性，这样就不会随局部坐标系变化而改变了 n深入分析插值多项式选择条件n深入分析选择条件是为了能由单元节点值唯一确定插值多项式 4节点四边形的插值多项式应该是 插值

12、多项式系数i (i = 0,1,2,3) 也是4个 单元模型构造整体坐标系法n基本思想 针对弹性体有限元网格建立一个统一的坐标系，每个单元的插值多项式都在这个坐标系上建立 单元模型构造整体坐标系法n2节点线单元12 oxu1u2x1x2ux1. 假设插值多项式2. 利用节点值求 a0 和 a1 单元模型构造整体坐标系法3. 代入a0 和 a1，得插值多项式 u(x)4. 按u1 和 u2合并同类项，设 l = x2- x1单元模型构造整体坐标系法N1 和 N2 称为单元的形函数；N 称为单元的形函数矩阵；ue 称为单元节点位移向量。 n2节点线的单元形函数单元模型构造整体坐标系法n二维3节点三

13、角形单元 建立整体坐标系oxy 单元模型构造整体坐标系法1. 假设插值多项式2. 首先，利用节点值求 0 、 1 和 2 n二维3节点三角形单元 单元模型构造整体坐标系法A为单元面积单元模型构造整体坐标系法3. 将 0 、 1 和 2 代入插值多项式，按u1、u2、u3合并同类项单元模型构造整体坐标系法4. 同理可得单元模型构造整体坐标系法5. 单元插值多项式为单元模型构造整体坐标系法6. 单元插值多项式写成矩阵形式（常用）单元模型构造整体坐标系法7. 单元插值多项式的另一种矩阵形式（不常用）单元模型构造整体坐标系法n4节点四面体单元单元模型构造整体坐标系法1. 假设插值多项式2. 插值多项式

14、为单元模型构造整体坐标系法（i=1,2,3,4） 循环轮换脚标1、2、3、4，相应可以得到a2，b2 ， c2 ， d2 、 a3 ， b3 ， c3 ， d3 、 a4 ， b4 ， c4 ， d4 单元模型构造整体坐标系法3. 单元插值多项式写成矩阵形式（常用）单元模型构造整体坐标系法4. 单元插值多项式另一种矩阵形式（不常用）单元模型构造整体坐标系法n从理论上讲，整体坐标系法可以求任意单元的形函数，但计算过程太复杂n只能求一维2节点线单元、二维3节点三角形单元和三维4节点四面体单元3种简单单元的形函数n复杂的或二次以上的单元必须采用局部坐标系法求n位移场 u 是形函数 Ni 的线性组合，

15、因此形函数Ni同样具有插值多项式的特性单元刚度矩阵2节点线单元n一维2节点线单元n单元插值关系 n单元几何关系 n单元本构关系 N=N1 N2 De=E单元刚度矩阵2节点线单元n单元刚度矩阵A为单元截面积；l为单元长度n矩阵B单元刚度矩阵三角形单元n二维3角形单元n单元插值关系 单元刚度矩阵三角形单元n单元几何关系 单元刚度矩阵三角形单元n单元本构关系 平面应力问题单元刚度矩阵三角形单元n矩阵B单元刚度矩阵三角形单元n单元刚度矩阵h为单元厚度k为对称的6*6常数矩阵A为单元面积单元模型构造整体坐标系法n单元形函数的特性正规性：单元形函数之和等于1。 正交性：形函数在本节点的值等于1，在其它节点

16、的值等于0。 单元模型构造局部坐标系法n为什么要采用局部坐标系法?在单元上假设一种局部坐标系，确定局部坐标系的度量，并在单元节点上标出局部坐标值根据插值多项式选择条件假设形函数多项式利用单元形函数的特性(正交性)求单元的形函数n一般需要如下几个过程 单元模型构造局部坐标系法n一维单元长度坐标 定义自然坐标L1、L2L1、L2也叫长度坐标 单元模型构造局部坐标系法n长度坐标L1、L2特性L1和L2不是相互独立的 L1和L2与形函数的关系 长度坐标系与整体坐标系之间的变换关系 单元模型构造局部坐标系法n一维二次单元 建立长度坐标系的目的是为了求一维高阶单元的形函数 一维3节点单元 单元模型构造局部

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

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

19、元模型构造局部坐标系法n二维10节点三角形单元 标注节点面积坐标单元模型构造局部坐标系法假设单元形函数根据Ni的正交性，求得单元形函数分别为 单元模型构造局部坐标系法n三维单元体积坐标 定义体积坐标Li单元模型构造局部坐标系法n体积坐标L1、L2 、L3 、L4特性L1 、 L2 、 L3和L4不是相互独立的 L1 、 L2 、 L3和L4与形函数的关系 单元模型构造局部坐标系法n体积坐标系与整体坐标系之间的转换关系 单元模型构造局部坐标系法n三维10节点四面体单元 标注节点面积坐标单元模型构造局部坐标系法假设单元形函数根据Ni的正交性，求得单元形函数分别为 单元模型构造局部坐标系法n小结n长

20、度坐标只适合求一维线单元的形函数 2节点线性单元、3节点二次单元n面积坐标只适合求二维三角形单元的形函数 3节点线性单元、6节点二次单元、10节点二次单元n体积坐标只适合求三维四面体单元的形函数 4节点线性单元、10节点二次单元单元模型构造局部坐标系法n正规自然坐标 求单元的形函数时，最常用的是采用正规自然坐标正规自然坐标系是一种正规化的曲线坐标系 单元模型构造局部坐标系法n一维情况 v一维2节点线性单元 假设单元形函数正规自然坐标单元模型构造局部坐标系法根据Ni的正交性求得0和1 单元模型构造局部坐标系法得形函数N1为 同理得形函数N2为 单元模型构造局部坐标系法v一维3节点二次单元 假设单

21、元形函数正规自然坐标单元模型构造局部坐标系法根据Ni的正交性，求得单元的形函数分别为 采用正规自然坐标方法可求更高阶单元的形函数 单元模型构造局部坐标系法n二维情况 v二维4节点四边形单元正规自然坐标单元模型构造局部坐标系法假设单元形函数根据Ni的正交性，求得单元的形函数分别为 单元模型构造局部坐标系法v二维8节点四边形 二次单元正规自然坐标单元模型构造局部坐标系法假设单元形函数根据Ni的正交性，求得单元的形函数分别为 单元模型构造局部坐标系法v二维12节点四边形 三次单元正规自然坐标单元模型构造局部坐标系法假设单元形函数根据Ni的正交性，求得单元的形函数 单元模型构造局部坐标系法n三维情况

22、v三维8节点六面体单元正规自然坐标单元模型构造局部坐标系法假设单元形函数根据Ni的正交性，求得单元的形函数 单元模型构造局部坐标系法v三维20节点六面体 二次单元正规自然坐标单元模型构造局部坐标系法假设单元形函数根据Ni的正交性，求得单元的形函数 单元模型等参单元n等参单元 单元内任意一点的位移u与单元节点位移ue之间的关系为 一般单元坐标的插值关系也采用与位移插值关系相同的变换关系即单元内任意一点的坐标x与单元节点坐标xe之间的关系为 单元模型等参单元n等参单元凡是几何形状和位移场采用同阶同参数插值关系来描述的单元，称为等参单元 前面介绍的所有单元都属于等参单元 在描述单元的几何形状和位移场

23、时，并不一定非采用同阶插值关系 单元模型等参单元n等参单元3节点三角形等参单元 单元模型等参单元n超参单元如果几何形状插值函数的阶数高于位移场插值函数的阶数，称为超参单元 n次参单元如果几何形状插值函数的阶数低于位移场插值函数的阶数，称为次参单元 单元刚度矩阵2节点线单元n单元插值关系 n单元几何关系 n单元本构关系 N=N1 N2 De=En一维2节点线单元单元刚度矩阵2节点线单元n矩阵B单元刚度矩阵2节点线单元l为单元长度单元刚度矩阵2节点线单元n单元刚度矩阵A为单元截面积；l为单元长度单元刚度矩阵四边形单元n二维4边形单元n单元插值关系 单元刚度矩阵四边形单元n单元几何关系 单元刚度矩阵

24、四边形单元n单元本构关系 平面应力问题单元刚度矩阵四边形单元n矩阵B由于Ni是关于自然坐标r和s的函数，因此不能直接偏导单元刚度矩阵四边形单元将Ni(i=1,2,3,4)，看成是关于整体坐标x和y的隐函数，利用链锁规则求偏导单元刚度矩阵四边形单元n雅可比矩阵 J利用等参单元关系 单元刚度矩阵四边形单元雅可比矩阵J中的4个元素分别表示为 单元刚度矩阵四边形单元B进一步表示为 单元刚度矩阵由于bi和ci是关于自然坐标r、s的函数，所以矩阵B和k也是它们的函数，一般要采用高斯积分方法计算单刚k 单元平衡方程组装过程 n为什么要组装 ? 消除内力n组装的原则是什么 ? 单元自由度与结构自由度对应单元平

25、衡方程组装过程 2 F 1 3 U3U4U2U1U5U6结构自由度向量U单元平衡方程组装过程 2 1 U3U4U2U11u1u2u3u4 3 U611 U2U12u1u2U5u3u42单元平衡方程组装过程2 1 U3U4U2U11u1u2u3u42单元平衡方程组装过程组装单元单元平衡方程组装过程 3 U611 U2U12u1u2U5u3u4单元平衡方程组装过程再组装单元总体刚度方程 K 称为总体刚度矩阵 U 称为位移向量 F 称为载荷向量 总体刚度矩阵K的特性 n对称性 n奇异性 n稀疏性 n非零元素带状分布 约束处理过程 n为什么要约束处理 ?总体平衡方程组是奇异的消除无限制的刚体运动 使总

26、体平衡方程组存在唯一一组解约束处理过程边界条件n 边界条件分类 力(载荷)边界条件位移边界条件 集中载荷力 表面分布力 自重力热交换引起的温度载荷 固定位移约束 强制位移约束 关联位移约束 约束处理过程模型简化xy约束处理过程模型简化yx xy约束处理过程约束方程123456789101112yx约束处理过程约束处理方法n位移约束处理方法 赋0赋1法 乘大数法 约束处理过程赋0赋1法n位移约束处理例子关联约束方程 约束处理过程赋0赋1法n赋0赋1法约束处理过程赋0赋1法n有6个方程，5个未知数，如果约束方程可以消除有限元平衡方程组的奇异性，则取任意5个方程联立求解，都会得到方程组的唯一一组解。

27、 n系数矩阵由原来的对称的变成了非对称的，这对于大规模有限元方程组求解是十分不利的，采用相同的求解方法，在求解时间和矩阵存贮容量方面都增加了一倍。 约束处理过程赋0赋1法n如何解决 ?可以发现, 约束处理过程实际上是在系数矩阵上做了一次列初等变换。 为了保持平衡方程组系数矩阵的对称性，对变换后的系数矩阵再做一次相同的行初等变换。 约束处理过程赋0赋1法具体做法：第4行乘以系数k加到第3行，并去掉第4行。为了保持系数矩阵的阶数，将第4行的所有元素赋0，在其对角线元素位置赋1。即所谓赋0赋1法。 约束处理过程赋0赋1法n经过初等变换，方程组的系数矩阵仍然保持对称性 n初等变换不会改变方程组的解 n

28、约束后的方程组可以求得5个未知数 n通过关联约束方程回代求解U4约束处理过程赋0赋1法n进一步改进方法 将关联约束方程直接引进约束后的方程组）中，使U4也直接参与方程组求解。 用关联约束取代约束后方程组的赋0赋1行（第4行），然后再做对称化处理。即将取代后的一行方程（第4行）乘以k加到第3行，则系数矩阵仍然保持对称性。约束处理过程赋0赋1法约束处理过程赋0赋1法n基本原理 利用初等变换对求解方程组进行相同的行列变换，既保证方程组解不会改变，又可以保持方程组系数矩阵的对称性。 在进行初等变换时，只要保证对方程组系数矩阵做相同的行列变换，就可以保持方程组系数矩阵的对称性。 约束处理过程赋0赋1法n

29、固定和强制位移约束条件处理 通过关联位移约束方法简化。如果k = 0，则退化成强制位移边界条件即 U4=C 如果k =C= 0，则退化成固定位移边界条件即 U4=0 约束处理过程赋0赋1法n强制位移约束条件处理 U4=C约束处理过程赋0赋1法n固定位移约束条件处理 U4=0约束处理过程赋0赋1法n固定位移和强制位移约束处理后的系数矩阵是相同的，只是简单地将方程组系数矩阵中要约束自由度的行列分别赋0，对角线元素赋1，这也是赋0赋1法的由来。约束处理过程赋0赋1法n在方程组载荷右端项的处理方法上两者是不同，处理固定位移边界条件时，只要将对应自由度的载荷赋0即可。处理强制位移边界条件时，要在方程组系

30、数矩阵未赋0赋1法前，先将对应自由度的列乘以系数C减到载荷右端项，再将对应自由度的载荷位置赋C。约束处理过程乘大数法n乘大数法n基本原理 利用矩阵的初等变换不改变方程组解的思想。 约束处理过程乘大数法关联约束方程 约束处理过程乘大数法关联约束方程 A是一个大数，是系数矩阵中对角线元素K44的1010倍量级以上 为什么要乘以大数A ？放大位移约束方程的优势约束处理过程乘大数法关联约束方程 将关联约束加到有限元平衡方程对应自由度行，第3行或第4行，这里取第4行约束处理过程乘大数法n如果关联约束方程可以消除有限元平衡方程的奇异性，约束后的方程组就存在唯一的一组解。n约束后的方程组的系数矩阵是非对称的

31、。n利用初等变换方法将系数矩阵变换成对称的 约束处理过程乘大数法关联约束方程 再乘以系数-k 加到约束后方程组的第3行 约束处理过程乘大数法约束处理过程乘大数法n强制位移边界条件 约束后的方程组简化为 约束处理过程乘大数法n固定位移边界条件 k = 0，C = 0 约束后的方程组简化为 约束处理过程乘大数法n固定位移和强制位移边界条件的乘大数约束处理相对比较简单，而且它们的系数矩阵约束后是相同的，只是简单地将方程组系数矩阵中要约束自由度的对角线元素加上一个相对大数A即可 n乘大数法的叫法并不十分准确，应该叫加大数法更贴切 n乘大数和加大数的效果是一样的约束处理过程两种方法比较n赋0赋1法在约束

32、处理过程中是严格精确的，而乘大数法是一种近似约束处理方法，它的精度取决于所乘大数A值 两种方法都可以消除有限元平衡方程的奇异性，得到符合实际边界条件的唯一一组解。但两种方法还是有很大的区别 约束处理过程两种方法比较n采用乘大数法约束处理后的有限元平衡方程在求解时可能造成解的失真，大数A值越大可能解的偏差会越大，而赋0赋1法就不会出现类似的问题，它在约束过程和求解过程都是精确的n乘大数法相对于赋0赋1法在约束处理过程上简单一些 约束处理过程两种方法比较n赋0赋1法实际上是将关联位移约束方程代入到有限元平衡方程中的，是代入法。而乘大数是将占绝对优势的关联位移约束方程合并到有限元平衡方程中的，是罚方

33、法，计算误差来自于合并过程，计算精度取决于关联位移约束方程的优势大小n商业软件中，位移边界条件的约束处理都采用赋0赋1法，乘大数很少被采用主要原因是它是一种近似方法，而且大数的大小也不好确定，有时还会造成求解失败 约束处理过程弹簧单元假设柔性弹簧 kOXYU4f f = kU4k约束处理过程弹簧单元弹簧约束方程 f = kU4方程组求解过程特点n方程组求解是有限元计算过程中很重要的一部分，在有限元法的发展过程中，有限元方程的求解效率一直是其应用的最大瓶颈之一 n有限元方程组的特点： 有限元方程组的系数矩阵具有对称、稀疏、带状分布以及正定、主元占优。有效地利用这些特点，以减少系数矩阵的存贮量，提

34、高方程组求解效率 方程组求解过程分类比较n线性方程组的解法主要分两大类: 直接解法：以高斯消去法基础，以等带宽或变带宽方式存贮系数矩阵内元素，对于求解规模比较大的问题，要存贮的元素非常巨大。 迭代解法：只需要存贮系数矩阵中非零元素，存贮量很小，一般是变带宽存贮量的20%或更少，有些算法的求解效率也非常高，适合求解大规模线性方程组。但是这种解法对接近病态的方程组很难保证收敛性。 方程组求解过程带宽定义n有限元方程组系数矩阵是稀疏的、非零元素呈带状分布，带宽就是它的宽度，带宽的大小是由系统有限元网格的节点号排序决定的，具体求法是 带宽=(单元最大节点号之差+1)*节点自由度数 n带宽是网格节点标注

35、方法直接决定的，不同标注方法带宽可能相差很大 方程组求解过程带宽n带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相差很大 方程组求解过程带宽n所示四边形网格的三种节点号标注方法，每个节点是2个自由度n结构的带宽分别是12，18，56，相差很大，其中12和56之间相差近5倍，这就意味着系数矩阵的存贮量也是相差5倍，因此，对于大规模复杂系统的节点号优化是十分必要的 方程组求解过程系数矩阵存贮 n系数矩阵存贮 如果节点号排序优化的比较好，系数矩阵的存贮量就会减少很多。根据系数矩阵的对称性，一般都是按半带宽存贮。n系数矩阵存贮的方法 二维等带宽存贮 一维变带宽存贮 方程组求解过程二维等带宽存贮

36、 n二维等带宽存贮 方程组求解过程二维等带宽存贮n二维等带宽存贮消除了最大带宽以外的全部零元素，节省了系数矩阵元素的存贮量。但是由于取最大带宽为存贮范围，因此不能排除在带宽内的大量零元素。当系数矩阵的各行带宽变化不大时，适合采用二维等带宽存贮，方程组求解过程中系数矩阵元素的寻址也比较方便，求解效率较高。n当出现局部带宽特别大的情况时，采用二维等带宽存贮时，将由于局部带宽过大而使整体系数矩阵的存贮大大增加。 方程组求解过程一维变带宽存贮 n一维变带宽存贮 一维变带宽存贮方法就是把变化的带宽内的元素按一定的顺序存贮在一个一维数组中。由于它不按最大带宽存贮，因此比二维等带宽存贮更节省内存。按照解法可

37、分为按行一维变带宽存贮和按列一维变带宽存贮。 n按行一维变带宽存贮 方程组求解过程一维变带宽存贮 辅助的寻址数组M n一维变带宽存贮是最节省内存的一种方法，但是由于要借助于寻址数组寻找系数矩阵元素的位置，相对二维等带宽存贮方法来说要复杂一些，而且在程序实现时也要复杂得多，方程组求解过程中也要消耗一些数组寻址时间。因此，在选用存贮方法时要权衡二者的利弊，统盘考虑。一般当带宽变化不大，计算机内存允许时，采用二维等带宽存贮方法是比较合适的。 方程组求解过程一维变带宽存贮 方程组求解过程求解方法n方程组求解方法 高斯消去法 三角分解法 雅可比（Jacobi）迭代法 高斯-赛德尔（Gauss-Seide

38、l）迭代法 共轭梯度迭代法 预优共轭梯度迭代法应变、应力回代过程 n单元应变和应力回代求解 通过求解有限元平衡方程得到有限元节点位移后，就可以进行系统的刚度校核。如果所分析问题要进行强度校核，就要回代求解单元的应变和应力。n由插值关系和几何关系可得单元应变，再通过本构关系得到单元应力数值积分 n为什么要进行数值积分?2节点线单元、3节点三角形单元和4节点四面体单元3种单元的单元刚度矩阵是常数矩阵，不需要再进行数值积分运算。 除了这3种单元外，一般其它单元的刚度矩阵都是积分变量的函数，要采用数值积分方法进行计算。 数值积分主要方法n常用的单元面内数值积分方法主要有： Hammer积分 Gauss

39、积分 数值积分Hammer积分 nHammer积分 在三角形单元和四面体单元中，自然坐标是面积坐标和体积坐标。采用这些坐标建立的单元形函数，其单元刚度矩阵的一般形式为 二维 数值积分Hammer积分三维 一维单元在实际有限元应用中，一般都采用正规自然坐标系法建立单元的形函数。它采用Gauss积分方案 数值积分Hammer积分nHammer积分运算 三角形单元的Hammer积分表示为 数值积分Hammer积分四面体单元的Hammer积分表示为 数值积分Hammer积分数值积分Hammer积分数值积分Gauss 积分 nGauss 积分 采用正规自然坐标确定形函数的单元，其单元刚度矩阵的一般形式为

40、 一维 二维 三维 数值积分Gauss 积分n这些积分形式在积分限上与高斯积分完全一致，因此高斯积分方案被广泛应用于此类积分形式中。它们的具体数值形式为 一维 二维 数值积分Gauss 积分三维 数值积分Gauss 积分数值积分积分阶次选择 n数值积分的阶次选择 求解单元平衡方程时，绝大多数情况要采用数值积分方法，如何选择数值积分的阶次将直接影响计算精度和计算量。如果积分阶次选择不当，有时甚至会导致计算失败 数值积分积分阶次选择n选择积分阶次的原则主要依据以下两点： 积分精度 积分阶次n与被积分多项式的阶次m有直接关系。一般来说，有限元应用的经验公式 积分项有两个应变矩阵B相乘，因此m一定是偶

41、数，则积分阶数n等于0.5、1.5、2.5、 数值积分积分阶次选择n常用单元的积分阶次选择 一维单元 一般都采用正规自然坐标系法得到的形函数 在单元平衡方程中雅可比矩阵中虽然也含有自然坐标，但是它只是单刚的一个系数，只对单刚中的每个元素的大小有相同的影响，不会改变单刚的特性 数值积分积分阶次选择2节点线单元只能取高斯积分点n=13节点线单元可以取n=1或n=24节点线单元可以取n=2或n=32节点线单元：m=0，n=0.53节点线单元：m=2，n=1.54节点线单元：m=4，n=2.5按经验公式计算：实际应用计算：数值积分积分阶次选择在有限元法中，把3节点单元取n=1以及4节点单元取n=2的积

42、分方案称为减缩积分，而3节点单元取n=2以及4节点单元取n=3的积分方案称为完全积分 实际数值结果表明，有时减缩积分方案会带来很大的计算误差，产生零能模式(沙漏模式) 正常积分方案的有时计算结果也会偏小，产生(剪切)闭锁现象 数值积分积分阶次选择造成这些现象的原因有很多，例如，单元形状、单元相对大小、单元受力状况、分析问题的类型等等。为了避免零能模式和闭锁现象的发生，一般采用减缩积分加阻尼矩阵方法。采用减缩积分方案时，对每个节点施加一个柔性弹簧，通过弹簧的阻尼增加刚度矩阵的稳定性，阻止零能模式的发生。但是弹簧的刚性系数越大，计算误差就越大，因此弹簧系数的选择也有一定的困难 数值积分积分阶次选择

43、三角形单元 3节点三角形单元按经验公式计算的积分阶次n=0.5，实际计算时只能取n=1 结果造成计算结果偏硬，有时会产生闭锁现象，实际有限元计算时也证明了这一点 数值积分积分阶次选择6节点三角形单元积分阶次是比较精确的 按经验公式计算的积分阶次应该取n=1.5，在单元面内应该是3个积分点 但是，这并不意味着单元的精度就比较高，因为单元的精度是由插值多项式本身决定的 数值积分积分阶次选择n四边形单元 4节点四边形单元：m=2，n=1.58节点四边形单元：m=4，n=2.512节点四边形单元：m=6，n=3.5按经验公式计算：4节点四边形单元可以取 n=1 或 n=28节点四边形单元可以取 n=2

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

45、分方案 n=220节点六面体单元取减缩积分方案 n=2或完全积分方案 n=3有限元基本方法总结n位移法基本过程 1）离散化过程 3）约束处理过程 2）单元平衡方程组装过程 5）应变、应力回代过程 4）方程组求解过程 有限元基本方法总结 为弹性体的应变 为弹性体的应力 u为弹性体的可容位移 弹性体处于平衡状态时，其势能应为最小 P 为表面分布力 G 为弹性体的体力 A 为弹性体的表面积 V 为弹性体的体积 有限元基本方法总结n单元插值关系 n单元几何关系 n单元本构关系 N为单元形函数矩阵 L为单元几何微分算子 为单元弹性矩阵 单元节点自由度向量有限元基本方法总结B 称为应变矩阵 单元平衡方程或

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

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

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

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

50、总结n对称性 n奇异性 n稀疏性 n非零元素带状分布 有限元基本方法总结n约束处理 总体平衡方程组是奇异的消除无限制的刚体运动 使总体平衡方程组存在唯一一组解有限元基本方法总结n位移约束处理方法 赋0赋1法 乘大数法 有限元基本方法总结n赋0赋1法的基本原理 利用初等变换对求解方程组进行相同的行列变换，既保证方程组解不会改变，又可以保持方程组系数矩阵的对称性。 在进行初等变换时，只要保证对方程组系数矩阵做相同的行列变换，就可以保持方程组系数矩阵的对称性。 有限元基本方法总结n乘大数法的基本原理 利用矩阵的初等变换不改变方程组解的思想 放大位移约束方程的优势 大数A是系数矩阵中对角线元素的101

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

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

53、项有两个应变矩阵B相乘，因此m一定是偶数，则积分阶数n等于0.5、1.5、2.5、 有限元基本方法总结减缩积分：低于经验公式的阶次正常积分：高于经验公式的阶次实际数值结果表明，有时减缩积分方案会带来很大的计算误差，产生零能模式 正常积分方案的有时计算结果也会偏小，产生闭锁现象 有限元基本方法总结造成这些现象的原因有很多，例如，单元形状、单元相对大小、单元受力状况、分析问题的类型等等。为了避免零能模式和闭锁现象的发生，一般采用减缩积分加阻尼矩阵方法。采用减缩积分方案时，对每个节点施加一个柔性弹簧，通过弹簧的阻尼增加刚度矩阵的稳定性，阻止零能模式的发生。但是弹簧的刚性系数越大，计算误差就越大，因此

54、弹簧系数的选择也有一定的困难 有限元基本方法总结n采用乘大数法约束处理后的有限元平衡方程在求解时可能造成解的失真，大数A值越大可能解的偏差会越大，而赋0赋1法就不会出现类似的问题，它在约束过程和求解过程都是精确的n乘大数法相对于赋0赋1法在约束处理过程上简单一些 板料成形有限元法分类n静力隐式有限元法n静力显式有限元法n动力显式有限元法n逆算法(一步成形)板料成形有限元法单元模型n单元模型n 薄膜单元n 薄壳单元n 中厚壳单元n 等效弯曲单元板料成形有限元法薄膜单元n薄膜单元是由二维三角形单元或四边形单元构造的空间板壳单元n薄膜单元适用于板料在变形过程中主要以拉伸和压缩变形为主，局部弯曲变形对

55、整个成形问题不产生大的影响 n液压胀形、半球冲头胀形等一类问题 板料成形有限元法薄膜单元n三角形薄膜单元 板料成形有限元法薄膜单元n局部坐标系以1节点为坐标原点，x轴与单元1、2边重合并指向2节点，z轴与单元法向量nz平行。单位法向量nz为 板料成形有限元法薄膜单元n随体局部坐标系oxyz与空间整体坐标系OXYZ之间坐标转换矩阵 为 板料成形有限元法薄膜单元n单元局部坐标自由度向量ue与整体坐标自由度向量Ue的变换关系为 板料成形有限元法薄膜单元n四边形薄膜单元 板料成形有限元法薄膜单元n局部坐标系以1节点为坐标原点，x轴与单元1、2边重合并指向2节点，z轴与单元法向量nz平行。单位法向量nz

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

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

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

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

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

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

62、力收敛准则为 为第n迭代步的失衡力 失衡力收敛准则当系统处于失稳状态时，失衡力的微小变化将引起位移增量的很大偏差，此时不能采用失衡力收敛准则 板料成形有限元法收敛准则能量收敛准则 常用能量收敛准则为 位移收敛准则和失衡力收敛准则都有一定的缺陷，相对而言能量收敛准则是比较好的，因为它同时控制位移增量和失衡力。能量收敛准则是把每次迭代后的内能增量与初始内能增量相比较，即 板料成形有限元法主要方法n板料成形有限元方法主要有：n基于Green应变和第二类Kirchhoff应力能量共轭的虚功原理，以初始时刻为参考构形的全量拉格朗日（TL）有限元方法n基于Green应变和第二类Kirchhoff应力能量共

63、轭的虚功原理，以当前时刻为参考构形的修正拉格朗日（UL）有限元方法板料成形有限元法主要方法n基于Lagrange（第一类Kirchhoff）应力与速度对物质坐标偏导数能量共轭的虚功原理，以当前时刻为参考构形的虚功率增量型有限元方法n基于动力学原理和中心差分方法建立的动力显式有限元方法n基于全量应变空间，全量形变理论和虚功原理的有限元逆算法板料成形有限元法TL法nTL法 TL法是取初始时刻构形作为参考构形，在所有的时间步长内的计算都参照时刻t0 = 0构形来定义。则t+t时刻的虚功方程为 板料成形有限元法TL法n单元插值关系 n单元本构关系 为时间增量步t内第二类Kirchhoff应力增量 为时

64、间增量步t内Green应变增量 为切线本构矩阵 板料成形有限元法TL法n单元几何关系t 时刻 t+t 时刻 板料成形有限元法TL法n线性化的TL法有限元平衡方程 KT为t时刻的单元切线刚度矩阵 Ks为初应力单元刚度矩阵 板料成形有限元法UL法nUL法 UL法是取当前时刻构形作为参考构形，在所有的时间步长内的计算都参照当前时刻t构形来定义。则t+t时刻的虚功方程可以表示为 板料成形有限元法UL法n单元插值关系 n单元本构关系 为时间增量步t内第二类Kirchhoff应力增量 为时间增量步t内Green应变增量 为切线本构矩阵 板料成形有限元法UL法n单元几何关系t 时刻 t+t 时刻 板料成形有

65、限元法UL法n线性化的UL法有限元平衡方程 Ks为初应力单元刚度矩阵 Rs为初应力节点力向量 板料成形有限元法虚功率增率法n虚功原理 由于Lagrange（第一类Kirchhoff）应力与速度对物质坐标偏导数是能量共轭的，可得时刻现实构形，以初始时刻为参考构形的弹塑性大变形虚功方程为 板料成形有限元法虚功率增率法n单元插值关系 n单元本构关系 单元内任意一点的速度v与单元节点速度向量ve的关系为 板料成形有限元法虚功率增率法n单元几何关系变形率张量dij为 为速度梯度张量 以当前时刻为参考构形时，变形率张量dij等于应变率张量 n线性化的虚功率增率法有限元平衡方程 板料成形有限元法虚功率增率法

66、板料成形有限元法逆算法n虚功方程 在最终状态上建立虚功方程经过单元离散化板料成形有限元法逆算法n单元插值关系 n单元本构关系 Hencky全量形变理论可得本构方程 板料成形有限元法逆算法n单元几何关系是B-1张量的特征值,是第一主应变方向板料成形有限元法逆算法n逆算法有限元平衡方程 板料成形有限元法动力显式算法n板料的运动学微分方程 n系统的虚功方程 板料成形有限元法动力显式算法n单元插值关系 n单元几何关系 板料成形有限元法动力显式算法n线性化有限元方程 M是质量矩阵 C是阻尼矩阵 P是节点外力向量 F是节点内力向量 通常，动力显式积分算法采用集中质量矩阵，既M是一对角矩阵，并取 板料成形有

67、限元法动力显式算法则联立方程组变成（节点数节点*自由度数）个相互独立的方程 n显式时间积分的中心差分算法 板料成形有限元法动力显式算法将节点速度和加速度用差分格式写成 n显式时间积分的中心差分算法 板料成形有限元法动力显式算法得到节点位移和速度的显式计算格式，前提条件是已知前两步的位移和速度 n临界时间步长的确定 动力显式算法几个关键问题由于中心差分算法是条件稳定的，为了保证系统计算的稳定性，对时间增量步长的大小必须加以限制 稳定性条件通常由系统的最高频率 决定 是最高模态中的临界阻尼 满足稳定性条件的时间增量步长可以由膨胀波沿网格中任意单元的最小穿越时间近似得到 n临界时间步长的确定 c 为

68、膨胀波在材料中的传播速度 为第n状态单元e的名义长度 动力显式算法几个关键问题n单元名义长度的确定 对于高阶单元来说，临界时间步长远比低阶单元的临界时间步长来的小。这一事实使得高阶单元对于显式积分算法相当不合适 动力显式算法几个关键问题n模具网格与板料网格剖分 动力显式算法几个关键问题板料网格 模具网格动力显式算法几个关键问题n压机成形速度与数值模拟成形速度 实际压机的成形速度约为0.1m/s0.5m/s 数值模拟时放大压机成形速度，一般假设为4.0m/s8.0m/s放大压机成形速度的目的是为了提高数值模拟的计算速度n动力效应动力显式算法几个关键问题由于放大了模拟成形过程的模具速度，会造成数值

69、模拟过程中的动力效应这种动力效应在复杂压料面压边过程中体现的比较明显减少动力效应的最有效办法是降低模拟过程中模具的速度动力显式算法几个关键问题n动力效应初始平板料 模具装配动力显式算法几个关键问题n动力效应3m/s 10m/s不同压边圈成形速度的动力效应n网格自适应加密技术动力显式算法几个关键问题为了减少计算网格数量，提高冲压成形模拟速度，根据单元和单元间的变形情况自动加密网格密度理论上加密级别是没有限制的，一般实际计算过程中常采用2-6级加密自适应虽然可以大幅度提高求解效率，但也会降低计算精度n网格自适应加密技术动力显式算法几个关键问题 动力显式算法几个关键问题 n网格自适应加密与减密技术虽

70、然上面给出的稳定性准则严格来说是对线性系统而言的，但对于非线性问题也给出了有用的稳定性估计。对线性问题时间步长的80%90%缩小，对于大多数非线性问题保持其系统稳定性是足够的。然而，十分重要的是：在整个计算过程中，要不断的检查能量的平衡问题。任何总能量的增加或损失（5%以上）都将导致失稳。Belytschko指出：常增量时间步不能保持解的稳定性，即使系统的最高频率 不断减小 板料成形有限元法动力显式算法板料成形有限元法算法比较n商业软件n动力显式算法 LS-DYNA3D、DYNAFORM、PAM-STAMPn静力隐式算法 AUTOFORM n静力显式算法 KMAS、ITAES n逆算法 FASTAMP、AUTOFORM、FASTFORM 3D 板料成形有限元法算法比较n算法的优缺点： 速度、求解稳定性 精度、模拟真实性n用户理论基础和工艺知识要求 有限元基本概念、工艺基本知识n软件的易用性 操作设计、界面设计 数模转换、曲面识别、网格剖分 板料成形有限元法发展趋势n高精度n大规模n面向模具设计n与CAD软件无缝集成 (不是现在所说的CAD/CAE/CAM一体化技术)板料成形有限元法发展方向n板料成形数值模拟的发展方向n面向数字化的虚拟制造技术n面向模具设计人员的工艺分析技术n面向专家系统的智能化的优化设计技术