有限元法基础

1、有限元法基础华中科技大学材料成形与模具技术国家重点实验室教材n董湘怀. 材料成形计算机模拟，第2版. 北京：机械工业出版社，2005 n讲课：24学时主要参考书n王勖成. 有限单元法. 北京：清华大学出版社，2003 nT. R. Chandrupatla, A.D. Belegundu著，曾攀译. 工程中的有限元方法. 北京：清华大学出版社，2006n曾攀. 有限元分析及应用. 北京：清华大学出版社，2004主要参考书n谢水生，李雷. 金属塑性成形的有限元模拟技术及应用. 北京：科学出版社，2008n李尚健. 金属塑性成形过程模拟. 北京：机械工业出版社，1999 n钟志华，李光耀. 薄板冲

2、压成型过程的计算机仿真与应用. 北京：北京理工大学出版社，1998材料成形问题的复杂性n非线性 几何非线性 物理非线性 状态非线性n模具和工艺的复杂性 复杂几何形状 多工序n多物理场耦合 变形 热传导n材料组织性能变化 相变、再结晶 织构 损伤数值模拟方法的优越性n经济、快速、优化、并行n 结果详尽 应力、应变、温度 组织性能变化n虚拟、灵活板料成形数值模拟实例板料成形数值模拟实例有限元法概述n有限元法的发展过程 1960年，Clough求解平面弹性问题，第一次提出“有限单元法”的名称n有限元法的工程应用 杆件结构问题 平板弯曲问题 壳体问题 热传导问题 动力学问题 流固耦合问题 接触碰撞问题

3、有限元法概述n有限元分析软件 NASTRAN，SAP， MARC，ANSYS，ABAQUS， JIFEX， LS-DYNA3D n有限元法的发展方向 本构模型，单元模型，数值分析方案，CAD/CAE/CAM集成 有限元法的基本思想n 化整为零n 积零为整有限元法的基本思想n基本思想 1）将连续的求解系统离散为一组由节点相互联在一起的单元组合体 2）在每个单元内假设近似函数来分片表示系统的求解场函数 有限元法的基本思想n有限元法分类 1）位移法：基于最小势能原理或虚功原理 2）力法： 基于最小余能原理 3）杂交法：基于修正余能原理 4）混合法：基于Reissner变分原理 有限元法的基本思想n位

4、移法基本过程 1）离散化过程 3）约束处理过程 2）单元平衡方程组装过程 5）应变、应力回代过程 4）方程组求解过程 离散化过程n最小势能原理 P V A G 弹性体 弹性体的势能peipWW 为弹性体变形后所具有的内能 iW为弹性体所受的外力功 eWdVVT21VAdVdAGuPuTT离散化过程 为弹性体的应变 为弹性体的应力 u为弹性体的可容位移 弹性体处于平衡状态时，其势能应为最小 P0TTTVAVdVdAdVGuPu0离散化过程n单元插值关系 eNuu n单元几何关系 n单元本构关系 Lu DeN为单元形函数矩阵 L为单元几何微分算子 为单元弹性矩阵 eD0TTTvavPdvdadvG

5、uPu0)()()(TTTTTTvevaeeeedvdadvGNuPNuuBDBu0TTTvvaeedvdadvGNPNuBDBeu为单元节点自由度向量离散化过程0)()()(TTTTTTvevaeeeedvdadvGNuPNuuBDBu0TTTvvaeedvdadvGNPNuBDB0TTTvavPdvdadvGuPuB 称为应变矩阵 LNB fku e单元平衡方程或单元刚度方程 k 称为单元刚度矩阵 vedvBDBkT f 称为单元载荷向量 vadvdaGNPNfTT单元刚度矩阵的特性 n对称性 n奇异性 n主元恒正且对角占优 离散化过程线弹性问题几何方程三维问题 Luwvuxzyzxyzy

6、xzuxwywzvxvyuzwyvxuzxyzxyzzyyxx000000000三维问题Lu Luwvuxzyzxyzyxzuxwywzvxvyuzwyvxuzxyzxyzzyyxx000000000线弹性问题几何方程二维问题 Luvuxyyxyvxuyvxuxyyyxx00二维问题平面应力和平面应变状态 线弹性问题几何方程二维问题 二维问题轴对称状态 Luwurzzrrrwzuzwrururzzzrr0010Luwvuxzyzxyzyxzuxwywzvxvyuzwyvxuzxyzxyzzyyxx000000000线弹性问题几何方程一维问题 一维问题 LuuxxuxxLuvuxyyxyvxuy

7、vxuxyyyxx00线弹性问题本构方程三维问题 三维问题De221000000221000000221000000100010001)21)(1 (EeDE为弹性模量；为泊松比 线弹性问题本构方程平面应力 二维问题平面应力状态 00yzxz000yzxzzzzxyzxyzzyyxx000 xyyyxxxyyyxxzxyzxyzzyyxxxyzzyyxx线弹性问题本构方程平面应力 zxyzxyzzyyxxezxyzxyzzyyxxD2100010112EeDxyyyxxxyyyxxE2100010112平面应力状态 00000 xyzzyyxxexyyyxxD线弹性问题本构方程平面应变 二维问

8、题平面应变状态 000yzxzzz00yzxzzxyzxyzzyyxx00 xyzzyyxxxyzzyyxxzxyzxyzzyyxxxyyyxx00000 xyyyxxexyzzyyxxD线弹性问题本构方程平面应变 zxyzxyzzyyxxezxyzxyzzyyxxDxyyyxxxyyyxxE221000101)21)(1 (平面应变状态 221000101)21)(1 (EeD线弹性问题本构方程轴对称 二维问题轴对称状态 00yzxz00yzxzzxyzxyzzyyxx00 xyzzyyxxxyzzyyxxzxyzxyzzyyxxxyzzyyxx线弹性问题本构方程轴对称 二维问题轴对称状态

9、00zr00zrzrzrzzrrzrzzrr00zrzzrrzrzrzzrrzrzzrr线弹性问题本构方程轴对称 轴对称状态 zxyzxyzzyyxxezxyzxyzzyyxxDzrzrzzrrezrzrzzrrDzrzzrrezrzzrr0000D221000010101)21)(1 (EeDzrzzrrzrzzrrE221000010101)21)(1 (线弹性问题本构方程一维问题 一维问题xxxxEEeD常用单元模型 n单元模型插值关系一一对应n单元类型一维单元、二维单元、三维单元等参单元、超参单元、次参单元常用单元模型n一维单元 2节点线单元 1 2 1 2 3 1 2 3节点线单元梁

10、单元常用单元模型n二维单元3节点三角形线性单元 1 2 3 6节点三角形二次单元 1 2 3 5 6 4 常用单元模型n二维单元10节点三角形三次单元 10 1 2 3 6 9 4 5 7 8 4节点四边形双线性单元 1 2 3 4 常用单元模型n二维单元8节点四边形二次单元12节点四边形三次单元 1 2 3 4 8 7 6 5 8 7 1 2 3 4 11 9 5 6 10 12 常用单元模型n三维单元4节点四面体线性单元10节点四面体二次单元 1 2 3 4 1 10 9 8 4 7 2 3 6 5 常用单元模型n三维单元8节点六面体线性单元20节点六面体二次单元 8 6 1 2 3 4

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

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

13、别是易于进行微分和积分运算。随着多项式函数阶次的增加，可以提高有限元法的计算精度。从理论上说，无限提高多项式的阶数，可以求得系统的精确解。单元模型构造方法 n整体坐标系法n局部坐标系法 nLagrange插值方法nHermite插值方法n1）Lagrange插值方法 单元模型构造基本思想与方法 n对于有n个节点的单元，如果它的节点参数中只含有场函数的节点值ui，i=1,2,n ，则单元内某一确定自由度方向的场函数u可插值为 1niiiuN un式中，Ni是插值函数，它有以下性质 单元模型构造基本思想与方法 1(,)1ijjjijniiN xyzNn式中，ij是kronecker符号，xj、yj

14、、zj 是单元节点坐标值。 n上式也反映了插值函数的一个共同特性：插值函数Ni在i节点处的值等于1，在其他节点处的值于0，也称为插值函数的正交性。 单元模型构造基本思想与方法 n插值函数的和等于1，也称为插值函数的正规性。 n2）Hermite插值方法 单元模型构造基本思想与方法 n如果希望在单元间的公共节点上还保持场函数导数的连续性，则节点参数中还应包含场函数导数的节点值。这时可以采用Hermite多项式作为单元的插值函数。例如，对于只有两个节点的一维单元，形函数Ni采用Hermite多项式插值为 单元模型构造基本思想与方法 221(0)11d( )( )( )diiiiiiuuHuH或者4

15、1( )( )iiiuNQ式中(0)(0)1122(1)(1)3142( ),( )( ),( )NHNHNHNH11223412,dd,ddQu QuuuQQnHermite多项式具有以下性质 单元模型构造基本思想与方法 (0)(0)(1)(1)d( )(), 0dd( )()0, djjiijijiijijHHHHn上式在两个节点处最高保持场函数的一阶导数连续性，这种多项式称为一阶Hermite多项式。零阶Hermite多项式就退化为Lagrange多项式。进一步推广，在节点处保持至场函数的n阶导数连续性，就称为n阶Hermite多项式。 单元模型构造基本思想与方法 221(0)11d(

16、)( )( )diiiiiiuuHuH单元模型构造方法n2节点线单元12 oxu1u2x1x2ux1. 假设插值多项式xaaxu10)(2. 利用节点值求 a0 和 a1 21021101xaauxaau12121xxuua1212210 xxxuxua单元模型构造方法3. 代入a0 和 a1，得插值多项式 u(x)xxxuuxxxuxuxu1212121221)(4. 按u1 和 u2合并同类项，设 l = x2- x1212121122112)(uuNNuulxxlxxulxxulxxxu单元模型构造方法n关键 如何构造插值多项式 u ?二维问题三维问题，如何构造插值多项式?n收敛性条件

17、在单元内，场函数必须是连续的； 完备性：插值多项式的阶次必须由低到高依次增加，不能出现跳跃现象； 协调性：各单元边界必须连续，单元边界不能出现开裂现象。 插值多项式收敛性条件 n收敛：当单元逐渐缩小时，如果插值多项式满足收敛性条件，则数值解将收敛于精确解 插值多项式收敛性条件n协调单元 满足插值多项式收敛性条件和的单元 n完备单元 满足插值多项式收敛性条件的单元ncr 阶连续性 插值多项式的第r阶导数是连续的 插值多项式收敛性条件n非协调单元与部分协调单元 对于一般固体力学问题来说，协调性要求单元在变形时，相邻单元之间不应引起开裂、重叠或其它不连续现象。例如，梁、板、壳等单元，在单元边界不但要

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

19、项式选择条件由选择条件可知，插值多项式函数在所有自由度方向上要满足各向同性性，这样就不会随局部坐标系变化而改变了 n深入分析xyyxy, xu3210)(23210)(xyxy, xu插值多项式选择条件n深入分析选择条件是为了能由单元节点值唯一确定插值多项式 4节点四边形的插值多项式应该是 xyyxy, xu3210)(插值多项式系数i (i = 0,1,2,3) 也是4个 单元模型构造整体坐标系法n基本思想 针对弹性体有限元网格建立一个统一的坐标系，每个单元的插值多项式都在这个坐标系上建立 y x o 弹性体 单元模型构造整体坐标系法n2节点线单元12 oxu1u2x1x2ux1. 假设插值

20、多项式xaaxu10)(2. 利用节点值求 a0 和 a1 21021101xaauxaau12121xxuua1212210 xxxuxua单元模型构造整体坐标系法3. 代入a0 和 a1，得插值多项式 u(x)xxxuuxxxuxuxu1212121221)(4. 按u1 和 u2合并同类项，设 l = x2- x1212121122112)(uuNNuulxxlxxulxxulxxxu单元模型构造整体坐标系法lxxN21Nu21212211)(uuNNuNuNxulxxN12N1 和 N2 称为单元的形函数；N 称为单元的形函数矩阵；u 称为单元节点位移向量。 n2节点线单元的形函数单元

21、模型构造整体坐标系法n二维3节点三角形单元 (x1 , y1) (u1 , v1) 1 y x o v 3 2 (x2 , y2) (x3 , y3) u (u3 , v3) (u2 , v2) (u , v) (x , y) 建立整体坐标系oxy 单元模型构造整体坐标系法1. 假设插值多项式yxyxu210),(2. 首先，利用节点值求 0 、 1 和 2 n二维3节点三角形单元 323103222102121101yxuyxuyxuyxyxv210),(单元模型构造整体坐标系法)(21)(21)(21332211233221113322110ucucucAubububAuauauaA332

22、21111121yxyxyxA A为单元面积单元模型构造整体坐标系法12332123132ax yx ybyycxx23113231213ax yx ybyycxx31221312321ax yx ybyycxx3. 将 0 、 1 和 2 代入插值多项式，按u1、u2、u3合并同类项332211),(uNuNuNyxu单元模型构造整体坐标系法)(21)(21)(21333322221111ycxbaANycxbaANycxbaAN4. 同理可得332211),(vNvNvNyxv单元模型构造整体坐标系法5. 单元插值多项式为332211332211),(),(vNvNvNyxvuNuNuNy

23、xu)(21)(21)(21333322221111ycxbaANycxbaANycxbaAN单元模型构造整体坐标系法332211321321000000vuvuvuNNNNNNvu6. 单元插值多项式写成矩阵形式（常用）单元模型构造整体坐标系法321321321321000000vvvuuuNNNNNNvu7. 单元插值多项式的另一种矩阵形式（不常用）单元模型构造整体坐标系法n4节点四面体单元 2 1 (u1 , v1 , w1) x y z 4 3 (x1 , y1 , z1) (x3 , y3 , z3) (x4 , y4 , z4) (u3 , v3 , w3) (u2 , v2 ,

24、w2) (x2 , y2 , z2) (u4 , v4 , w4) (x , y , z) (u , v , w) 单元模型构造整体坐标系法zyxzyxwzyxzyxvzyxzyxu421042104210),(),(),(1. 假设插值多项式443322114433221144332211),(),(),(wNwNwNwNzyxwvNvNvNvNzyxvuNuNuNuNzyxu2. 插值多项式为单元模型构造整体坐标系法)(61zdycxbaVNiiiii（i=1,2,3,4） 4443332221zyxzyxzyxa 4433221111zyzyzyb4433221111zxzxzxc111

25、4433221yxyxyxd 444333222111111161zyxzyxzyxzyxV 循环轮换脚标1、2、3、4，相应可以得到a2，b2 ， c2 ， d2 、 a3 ， b3 ， c3 ， d3 、 a4 ， b4 ， c4 ， d4 单元模型构造整体坐标系法444333222111432143214321000000000000000000000000wvuwvuwvuwvuNNNNNNNNNNNNwvu3. 单元插值多项式写成矩阵形式（常用）单元模型构造整体坐标系法4. 单元插值多项式另一种矩阵形式（不常用）432143214321432143214321000000000000

26、000000000000wwwwvvvvuuuuNNNNNNNNNNNNwvu单元模型构造整体坐标系法n从理论上讲，整体坐标系法可以求任意单元的形函数，但计算过程太复杂n只能求一维2节点线单元、二维3节点三角形单元和三维4节点四面体单元3种简单单元的形函数n复杂的或二次以上的单元必须采用局部坐标系法求n位移场 u 是形函数 Ni 的线性组合，因此形函数Ni同样具有插值多项式的特性单元模型构造整体坐标系法n单元形函数的特性正规性：单元形函数之和等于1。 正交性：形函数在本节点的值等于1，在其它节点的值等于0。 单元刚度矩阵2节点线单元n一维2节点线单元n单元插值关系 eNuu n单元几何关系 n

27、单元本构关系 Lu DeN=N1 N2 De=E21uueuxLlxxNlxxN/ )(/ )(1221单元刚度矩阵2节点线单元n单元刚度矩阵vedvBDBkTvedvBDBkT x yzedxdydzBDBTxedxABDBTA为单元截面积；l为单元长度n矩阵BLNB LNB /21xNxN/1/1llBDBeAlT1111lAE/1/1/1/1llEllAl单元刚度矩阵三角形单元n二维3角形单元n单元插值关系 eNuu T332211vuvuvueu321321000000NNNNNNN)(21)(21)(21333322221111ycxbaANycxbaANycxbaAN单元刚度矩阵三

28、角形单元n单元几何关系 Lu Luvuxyyxxyyyxx00单元刚度矩阵三角形单元n单元本构关系 DexyyyxxxyyyxxE2100010112平面应力问题单元刚度矩阵三角形单元n矩阵BLNB LNB xNyNxNyNxNyNyNyNyNxNxNxN33221132132100000033221132132100000021bcbcbccccbbbA单元刚度矩阵三角形单元n单元刚度矩阵vedvBDBkTvedvBDBkT x yzedxdydzBDBTx yedxdyhBDBTh为单元厚度BDBehATk为对称的6*6常数矩阵A为单元面积单元模型等参单元n等参单元 单元内任意一点的位移u

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

30、的阶数，称为超参单元 n次参单元如果几何形状插值函数的阶数低于位移场插值函数的阶数，称为次参单元 单元平衡方程组装过程 n为什么要组装 ? 2 F 1 3 消除内力n组装的原则是什么 ? 单元自由度与结构自由度对应单元平衡方程组装过程 2 F 1 3 U3U4U2U1U5U6654321UUUUUUU结构自由度向量U单元平衡方程组装过程 3 U611 U2U12u1u2U5u3u443214321UUUUuuuueu65214321UUUUuuuueu2 1 U3U4U2U11u1u2u3u42单元平衡方程组装过程1413121143211441431421411341331321311241

31、2312212111411311211143214321ffffuuuukkkkkkkkkkkkkkkk2 1 U3U4U2U11u1u2u3u42单元平衡方程组装过程 3 U611 U2U12u1u2U5u3u424232221432124424324224123423323223122422322222121421321221165216521ffffuuuukkkkkkkkkkkkkkkk单元平衡方程组装过程n单元矩阵扩阶11111112131411112122232411113132333411114142434412345610020030040050000006000000kkkk

32、kkkkkkkkkkkk2 1 U3U4U2U11u1u2u3u42单元平衡方程组装过程 3 U611 U2U12u1u2U5u3u4222211121314222221222324222231323334222241424344 1 2 3 4 5 610020030000004000000500600kkkkkkkkkkkkkkkkn单元矩阵扩阶单元平衡方程组装过程n组装总刚12121122111112121314131412121122212122222324232411113132333411114142434422223132333422224142434412 3456123004

33、00500600kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk单元平衡方程组装过程n总体平衡方程1212112211111121213141314121211222212122222324232411113313233341111441424344222231323334522224142434461 2 345612300400500600UkkkkkkkkUkkkkkkkkUkkkkUkkkkkkkkUkkkkU1211122213142324ffffffff单元平衡方程组装过程00000000000000000000006543216543211413121165432

34、1144143142141134133132131124123122121114113112111ffffUUUUUUkkkkkkkkkkkkkkkk组装单元单元平衡方程组装过程242314132212211165432124424324224123423323223114414314214113413313213122422312412322212222112121421311411321211221111100000000654321654321ffffffffUUUUUUkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk再组装单元FKU 总体刚度方程 K 称为总体刚度矩阵

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

36、VVVVUUUU121110约束处理过程约束处理方法n位移约束处理方法 赋0赋1法 乘大数法 约束处理过程赋0赋1法n位移约束处理例子FKU 654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKCkUU34关联约束方程 约束处理过程赋0赋1法n赋0赋1法6543216543216665646362615655545352514645444342413635343332312625242

37、32221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK646545444343242141654321666564636261565554535251464544434241363534333231262524232221161514131211000000CKFCKFCKFCKFCKFCKFUUUUUUKKkKKKKKKkKKKKKKkKKKKKKkKKKKKKkKKKKKKkKKKK约束处理过程赋0赋1法n有6个方程，5个未知数，如果约束方程可以消除有限元平衡方程组的奇异性，则取任意5个方程联立求解，都会得到方程组

38、的唯一一组解。 n系数矩阵由原来的对称的变成了非对称的，这对于大规模有限元方程组求解是十分不利的，采用相同的求解方法，在求解时间和矩阵存贮容量方面都增加了一倍。 约束处理过程赋0赋1法n如何解决 ?可以发现, 约束处理过程实际上是在系数矩阵上做了一次列初等变换。 为了保持平衡方程组系数矩阵的对称性，对变换后的系数矩阵再做一次相同的行初等变换。 约束处理过程赋0赋1法6465454443432421416543216665646362615655545352514636453544433433423241312625242322211615141312110)(000010000)(00CKFC

39、KFCKFkCKFCKFCKFUUUUUUKKkKKKKKKkKKKKkKKkKKkKKkkKKkKKkKKKKkKKKKKKkKKKK具体做法：第4行乘以系数k加到第3行，并去掉第4行。为了保持系数矩阵的阶数，将第4行的所有元素赋0，在其对角线元素位置赋1。即所谓赋0赋1法。 646545444444343242141654321666564636261565554535251464544434241463645354443343342324131262524232221161514131211)(0000)(00CKFCKFCKFCKFkCKFCKFCKFUUUUUUKKkKKKKKKkK

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

41、5655545352514636453524443343342324131262524232221161514131211)(0000100)(00CKFCKFCkCCKFkCKFCKFCKFUUUUUUKKkKKKKKKkKKKKkkKKkKKkkkKKkkKKkKKkKKKKkKKKKKKkKKKK约束处理过程赋0赋1法n基本原理 利用初等变换对求解方程组进行相同的行列变换，既保证方程组解不会改变，又可以保持方程组系数矩阵的对称性。 在进行初等变换时，只要保证对方程组系数矩阵做相同的行列变换，就可以保持方程组系数矩阵的对称性。 约束处理过程赋0赋1法n固定和强制位移约束条件处理 通过关联位

42、移约束方法简化。如果k = 0，则退化成强制位移边界条件即 U4=C 如果k =C= 0，则退化成固定位移边界条件即 U4=0 6465454443432421416543216665646362615655545352514636453524443343342324131262524232221161514131211)(0000100)(00CKFCKFCkCCKFkCKFCKFCKFUUUUUUKKkKKKKKKkKKKKkkKKkKKkkkKKkkKKkKKkKKKKkKKKKKKkKKKK约束处理过程赋0赋1法n强制位移约束条件处理 U4=C646545343242141654321

43、6665636261565553525136353332312625232221161513121100001000000CKFCKFCCKFCKFCKFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKK约束处理过程赋0赋1法n固定位移约束条件处理 U4=06543216665636261565553525136353332312625232221161513121100001000000UUUUUUKKKKKKKKKKKKKKKKKKKKKKKKK653210FFFFF646545343242141CKFCKFCCKFCKFCKF约束处理过程赋0赋1法n固定位移和强制位移约束处理

44、后的系数矩阵是相同的，只是简单地将方程组系数矩阵中要约束自由度的行列分别赋0，对角线元素赋1，这也是赋0赋1法的由来。约束处理过程赋0赋1法n在方程组载荷右端项的处理方法上两者是不同，处理固定位移边界条件时，只要将对应自由度的载荷赋0即可。处理强制位移边界条件时，要在方程组系数矩阵未赋0赋1法前，先将对应自由度的列乘以系数C减到载荷右端项，再将对应自由度的载荷位置赋C。约束处理过程乘大数法n乘大数法n基本原理 利用矩阵的初等变换不改变方程组解的思想。 约束处理过程乘大数法65432165432166656463626156555453525146454443424136353433323126

45、2524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKCkUU34关联约束方程 约束处理过程乘大数法CUkU43关联约束方程 ACAUAkU43A是一个大数，是系数矩阵中对角线元素K44的1010倍量级以上 为什么要乘以大数A ？放大位移约束方程的优势约束处理过程乘大数法关联约束方程 ACAUAkU43将关联约束加到有限元平衡方程对应自由度行，第3行或第4行，这里取第4行6543216543216665646362615655545352514645444342413635343332312625242322

46、21161514131211FFACFFFFUUUUUUKKKKKKKKKKKKKKAKAkKKKKKKKKKKKKKKKKKKKKK约束处理过程乘大数法n如果关联约束方程可以消除有限元平衡方程的奇异性，约束后的方程组就存在唯一的一组解。n约束后的方程组的系数矩阵是非对称的。n利用初等变换方法将系数矩阵变换成对称的 约束处理过程乘大数法关联约束方程 ACAUAkU43再乘以系数-k AkCAkUUAk432加到约束后方程组的第3行 约束处理过程乘大数法654321654321666564636261565554535251464544434241363534333231262524232221

47、161514131211FFACFFFFUUUUUUKKKKKKKKKKKKKKAKAkKKKKKKKKKKKKKKKKKKKKK6543216543216665646362615655545352514645444342413635342333231262524232221161514131211FFACFAkCFFFUUUUUUKKKKKKKKKKKKKKAKAkKKKKKAkKAkKKKKKKKKKKKKKKK约束处理过程乘大数法n强制位移边界条件 )0(4kCU约束后的方程组简化为 6543216543216665646362615655545352514645444342413635

48、34333231262524232221161514131211FFACFFFFUUUUUUKKKKKKKKKKKKKKAKKKKKKKKKKKKKKKKKKKKKK6543216543216665646362615655545352514645444342413635342333231262524232221161514131211FFACFAkCFFFUUUUUUKKKKKKKKKKKKKKAKAkKKKKKAkKAkKKKKKKKKKKKKKKK约束处理过程乘大数法n固定位移边界条件 k = 0，C = 0 约束后的方程组简化为 65432165432166656463626156555

49、4535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKAKKKKKKKKKKKKKKKKKKKKKK654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFACFFFFUUUUUUKKKKKKKKKKKKKKAKKKKKKKKKKKKKKKKKKKKKK约束处理过程乘大数法n固定位移和强制位移边界条件的乘大数约束处理相对比较简单，而且它们的系数矩阵约束后是相同的，只是简单地

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

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

52、4434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFfFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK6544321654321666564636261565554535251464544434241363534333231262524232221161

53、514131211FFkUFFFFUUUUUUKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK654321654321666564636261565554535251464544434241363534333231262524232221161514131211FFFFFFUUUUUUKKKKKKKKKKKKKKkKKKKKKKKKKKKKKKKKKKKKK方程组求解过程特点n方程组求解是有限元计算过程中很重要的一部分，在有限元法的发展过程中，有限元方程的求解效率一直是其应用的最大瓶颈之一 n有限元方程组的特点： 有限元方程组的系数矩阵具有对称、稀疏、带状分布以及正定、

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

55、元最大节点号之差+1)*节点自由度数 n带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相差很大 方程组求解过程带宽n带宽是网格节点标注方法直接决定的，不同标注方法带宽可能相差很大 1 2 3 4 5 6 9 10 7 8 11 12 14 13 16 15 20 19 18 17 21 22 23 24 28 26 27 25 1 2 3 4 5 6 9 10 7 8 11 12 14 13 16 15 20 19 18 17 21 22 23 24 28 26 27 25 1 2 3 4 5 6 9 10 7 8 11 12 14 13 16 15 20 19 18 17 21 22

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

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

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

59、单元应变和应力回代求解 通过求解有限元平衡方程得到有限元节点位移后，就可以进行系统的刚度校核。如果所分析问题要进行强度校核，就要回代求解单元的应变和应力。n由插值关系和几何关系可得单元应变，再通过本构关系得到单元应力eeBuLNuDe板料成形有限元法分类n弹塑性有限元法n刚塑性有限元法弹塑性刚塑性板料成形有限元法分类n静力隐式有限元法n静力显式有限元法n动力显式有限元法n逆算法(一步成形)板料成形有限元法单元模型n单元模型n 薄膜单元n 薄壳单元n 中厚壳单元n 等效弯曲单元板料成形有限元法单元平衡方程n单元平衡方程n单元插值关系(局部坐标系) eNuu n单元几何关系(局部坐标系) n单元本

60、构关系(局部坐标系) Lu De0TTTvavPdvdadvGuPun最小势能原理(局部坐标系) n坐标变换关系 eeTUu 板料成形有限元法单元平衡方程0TTTvavPdvdadvGuPu0)()()(TTTTTTvevaeeeedvdadvGNuPNuuBDBu代入局部坐标系的单元3个关系式, 得代入坐标转换关系式, 得0)()()(TTTTTTTTTvevaeeeedvdadvGNTUPNTUTUBDBTU0TTTTTTvvaeedvdadvGNTPNTTUBDBT0)()()(TTTTTTvvaeedvdadvGNTPNTTUBDBTeeefTTUkTTTeeeFUK板料成形有限元法单