第一章有限元法基本原理

1、农业机械有限元软件方法农业机械有限元软件方法吉林大学生物与农业工程学院吉林大学生物与农业工程学院韩志武韩志武 第一章第一章 有限元法基本原理有限元法基本原理1.11.1有限元法方法及其历史有限元法方法及其历史1.1 1.1 有限元法简介及其历史有限元法简介及其历史v精确解：少数方程性质简单，形状规则；精确解：少数方程性质简单，形状规则； v复杂问题解：简化假设复杂问题解：简化假设+ +数值解法数值解法 v有限差分法：网格，用差分方程近似微分方程有限差分法：网格，用差分方程近似微分方程流体应用；流体应用； v其他方法：配点法、最小二乘法、其他方法：配点法、最小二乘法、GalerkinGalerk

2、in法、力矩法、里法、力矩法、里兹法；兹法；v共同点：在整个求解域上假设近似函数。共同点：在整个求解域上假设近似函数。 有限元的发展简史有限元的发展简史有限元方法的提出有限元方法的提出现代有限元现代有限元 19431943年，年， Courant Courant应用定义在三角形应用定义在三角形区域上的分片连续区域上的分片连续函数和最小位能原函数和最小位能原理相结合，求解理相结合，求解St.VenantSt.Venant 特殊特殊问题。问题。19561956年，年， Turner,Turner,Clough,Clough,刚架刚架位移法推广应位移法推广应用于弹性力学用于弹性力学问题，分析飞问题，

3、分析飞机结构。机结构。 19601960年，年，CloughClough处理处理平面弹性问题，平面弹性问题，第一次提出第一次提出“有限单元法有限单元法”。19631964,19631964,BesselingBesseling, ,Melosh,JonesMelosh,Jones证明有限元法是证明有限元法是基于变分原理基于变分原理的里兹法，的里兹法，确认了有限元确认了有限元法是处理连续法是处理连续介质问题的介质问题的一种普遍方法。一种普遍方法。有限元法假定的近似函数不是在全求解域而有限元法假定的近似函数不是在全求解域而是在单元上规定的，而且事先不要求满足任是在单元上规定的，而且事先不要求满足任

4、何边界条件，可以处理很复杂的连续介质问何边界条件，可以处理很复杂的连续介质问题。题。四十多年来，弹性力学平面问题扩展到空四十多年来，弹性力学平面问题扩展到空间问题，板壳问题，由静力平衡问题扩展间问题，板壳问题，由静力平衡问题扩展到稳定问题，动力问题和波动问题。分析到稳定问题，动力问题和波动问题。分析对象从弹性材料扩展到塑性，粘弹性，粘对象从弹性材料扩展到塑性，粘弹性，粘塑性和复合材料，从固体力学扩展到流体塑性和复合材料，从固体力学扩展到流体力学，传热学等连续介质力学领域。力学，传热学等连续介质力学领域。变分法建立变分法建立有限元方程有限元方程与经典里兹与经典里兹法的主要法的主要区别区别有限元法

5、有限元法的应用的应用有限元法的基本解题思路：有限元法的基本解题思路： 将连续求解区域离散一组有限个且按一定方将连续求解区域离散一组有限个且按一定方式相互联结在一起的单元的组合体，利用在每一个式相互联结在一起的单元的组合体，利用在每一个单元内假设的近似函数来分片的表示全求解域上待单元内假设的近似函数来分片的表示全求解域上待求的未知场函数。单元内的近似函数由插值函数表求的未知场函数。单元内的近似函数由插值函数表达，未知场函数及其导数在各个结点上的数值就成达，未知场函数及其导数在各个结点上的数值就成为新的未知量，从而使一个连续的无限自由度问题为新的未知量，从而使一个连续的无限自由度问题变成离散的有限

6、自由度问题，插值求整个解。变成离散的有限自由度问题，插值求整个解。 参考书目：参考书目： v有限元法及其在锻压工程中的应用有限元法及其在锻压工程中的应用吕丽萍主编，吕丽萍主编，西北工业大学出版社西北工业大学出版社 v弹性和塑性力学中的有限单元法弹性和塑性力学中的有限单元法丁皓江等主编，丁皓江等主编，机械工业出版社机械工业出版社v有限元分析的基本方法及工程应用有限元分析的基本方法及工程应用周昌玉、贺周昌玉、贺小华小华 编著，化学工业出版社编著，化学工业出版社1.2 1.2 位移函数与形状函数位移函数与形状函数1 1、坐标系、坐标系以杆单元为例：以杆单元为例： X YZi jiu iUiv iVi

7、w iWjujUjv jVjw jW 和和 是是 两点沿方向的位移分量。两点沿方向的位移分量。 和和 是是 两点沿方向的节点力分量。两点沿方向的节点力分量。统一规定：和坐标轴正向一致的为正。统一规定：和坐标轴正向一致的为正。 iuiviwjujvjwijiUiViWjUjVjWijv内位移：杆单元在节点力的作用下所产生的内位移：杆单元在节点力的作用下所产生的位移称为内位移。位移称为内位移。v位移函数：描绘内位移的函数。位移函数：描绘内位移的函数。)(iiUu)(jjUuijlX材料力学：仅受轴向力作用的杆，其中各点的位移是材料力学：仅受轴向力作用的杆，其中各点的位移是沿杆的轴线按线性规律变化，

8、即：沿杆的轴线按线性规律变化，即： 为杆单元位移函数。为杆单元位移函数。 其中：其中： ， 待定常数，由待定常数，由 ，节点的位移确定。，节点的位移确定。用矩阵表示为用矩阵表示为 ： xaaxu211 . 1 1a2aij aQaaxxuf2112 . 1式中：式中：1 1和和 为基底函数，为基底函数， 为基底函数矩阵。为基底函数矩阵。 ， 为为单元的广义位移，单元的广义位移， 为广义位移列阵。为广义位移列阵。由单元的边界条件，确定广义位移：由单元的边界条件，确定广义位移： ， ; , ; , 。 代入代入 式：式： x Q1a2a a0 x iuu0lx ju lu1 . 1laauauji

9、211矩阵形式为：矩阵形式为：式中：式中： 为单元的节点位移为单元的节点位移 。 acaaluujie21101 ejiuu则：则： 为变换矩阵为变换矩阵 。式中：式中： ace3 . 1 c eca14 . 1 111010111clll则：则： 代入代入 式：式：121011iijijuuauuualll5 . 1iua 1juuaij21 . 1设：设： ， 代入上式，得：代入上式，得：矩阵形式：矩阵形式： jiijiulxulxxluuuxu1lxNi1lxNj jjiiuNuNxua6 . 1 jijiuuNNxub6 . 1即：即： 为形状函数矩阵。为形状函数矩阵。由式由式 可知：

10、当可知：当 ， 时，时， ； 当当 ， 时，时， 。 形状函数的力学含义：当单元的一个节点位移为单形状函数的力学含义：当单元的一个节点位移为单位值，其他节点的位移为零时，单元内位移的分布规律。位值，其他节点的位移为零时，单元内位移的分布规律。 eNfc6 . 1 Na6 . 11iu0ju iNxu0iu1ju jNxu数学意义：如果说结构被有限个自然节点离散化为有数学意义：如果说结构被有限个自然节点离散化为有限个单元的集合，实现了结构模型的离散化，那么限个单元的集合，实现了结构模型的离散化，那么形状函数则完成了数学模型。形状函数则完成了数学模型。将式（将式（1.41.4）代入式（）代入式（1

11、.21.2）得：）得：由式（由式（1.61.6） eCQxu1 eNxuf则：则： 位移函数或形状函数的选择是有限元分位移函数或形状函数的选择是有限元分析的关键，位移函数选择的优劣，会直接影析的关键，位移函数选择的优劣，会直接影响到解的收敛性及解的精确度。响到解的收敛性及解的精确度。 1CQN7 . 11.3 1.3 单元应力和应变单元应力和应变位移函数位移函数几何方程（几何方程（应变）应变）物理方程（物理方程（应力）应力）杆单元的几何方程为：杆单元的几何方程为：dxdux eeexllxdxdlxdxdNdxd1111简化为：简化为： 为几何矩阵。为几何矩阵。 杆单元的物理方程为：杆单元的物

12、理方程为：或或 ： 为弹性矩阵为弹性矩阵 。 eB8 . 1 111lB9 . 1 BxxE D D将式（将式（1.81.8）代入上式，得：）代入上式，得： 为应力矩阵为应力矩阵 。对于杆单元：对于杆单元： eeSBD10. 1 BDS 11. 1 S 11 11 1ESDBEll 为弹性矩阵，对于杆单元为弹性矩阵，对于杆单元 ，是，是1X1 1X1 阶矩阵。阶矩阵。 （1.81.8）、（）、（1.101.10）是两个常数公式。）是两个常数公式。 D ED 1.4 1.4 虚功原理虚功原理设有一受外力作用的物体，如下图所示：设有一受外力作用的物体，如下图所示：ABCijiUjUiVjViWjW

13、XYZO 节点外力：节点外力： ， 和和 ； 节点外力：节点外力： ， 和和 。外力用外力用 表示；内力用表示；内力用 表示。表示。iiUiViWjjUjVjW F TjjjiiiWVUWVUF Tzxyzxyzyx设物体在外力设物体在外力 和内力和内力 以及边界固定点以及边界固定点A A、B B、C C处支反力作用下处于平衡状态。处支反力作用下处于平衡状态。假设物体发生了虚位移：假设物体发生了虚位移： ， ， ， ， ， ；由虚位移产生的虚应变：由虚位移产生的虚应变： ， ， ， ， ， 。 F *iu*iv*iw*ju*jv*jw*x*y*z*xy*yz*zx发生虚位移时，外力在虚位移上所

14、做的虚功是：发生虚位移时，外力在虚位移上所做的虚功是：应力在虚应变上的虚功是：应力在虚应变上的虚功是：整个物体的虚应变能是：整个物体的虚应变能是： FwWvVuUwWvVuUTjjjjjjiiiiii* Tzxzxyzyzxyxyzzyyxx* dxdydzT*虚功原理：如果在虚位移发生之前，物体处于虚功原理：如果在虚位移发生之前，物体处于平衡状态，那么在虚位移发生时，外力所做平衡状态，那么在虚位移发生时，外力所做的虚功等于物体的虚应变能。即：的虚功等于物体的虚应变能。即： 称为弹性体的虚功方程。称为弹性体的虚功方程。 VTTdxdydzF*12. 11.5 1.5 单元刚度矩阵单元刚度矩阵假

15、设杆单元假设杆单元L L：其中：其中： , , , , 为参数。为参数。liiUujjUuAEl杆单元应力杆单元应力- -应变关系为：应变关系为： 则：则： 由力的平衡条件：由力的平衡条件：则：则： luuEEAUijjijjuulEAU13. 10jiUUijiuulEAU14. 1（1.131.13）和（）和（1.141.14）用矩阵表示：）用矩阵表示：即：即： 为单元为单元 的节点力向量的节点力向量 ； 为单元为单元 的节点位移向量。的节点位移向量。 jijjjiijiijijiuuKKKKuulEAUU111115. 1 eeeKF16. 1e jieuue 为单元为单元 的刚度矩阵。

16、的刚度矩阵。刚度系数：单位节点位移分量所引起的节点力分量。刚度系数：单位节点位移分量所引起的节点力分量。单元：节点力为单元：节点力为 ，内力为，内力为 ，节点虚位移，节点虚位移 ，单元虚位移单元虚位移 ，单元虚应变，单元虚应变 。 1111lEAKKKKKjjjiijiie17. 1e eF e* ef* e*单元的外力虚功：单元的外力虚功： 单元的虚应变能：单元的虚应变能： 虚功原理虚功原理 ： eTeFV* VTedVU*UV将将 代入上式，得：代入上式，得： eeBDB eTTeeTedVBDBF* 0*eTeTedVBDBF由虚位移的任意性：由虚位移的任意性：与（与（1.161.16）

17、式对比，得：）式对比，得：称为单元的刚度矩阵，简称单刚。称为单元的刚度矩阵，简称单刚。 0eTedVBDBF dVBDBKTe1.6 整体刚度矩阵的集成整体刚度矩阵的集成整体结构整体结构离散离散有限个单元有限个单元单元分析单元分析单元特性集合单元特性集合台阶轴：台阶轴： 台阶轴台阶轴 离散化模型离散化模型 单元模型单元模型1122233 1P3P a b c11uP11uP33uP33uP2u22uP22uP a b c则单元（则单元（1 1）的刚度方程，考虑（）的刚度方程，考虑（1.151.15）式：）式：单元（单元（2 2）的刚度方程：）的刚度方程： 1211122121112111121

18、1uuKKKKPP 121221112112121121111111uKuKPuKuKP 23222332322232222322uuKKKKPP台阶轴在外力台阶轴在外力 和和 的作用下发生变形时，在节点处的作用下发生变形时，在节点处的变形必须是连续的，即：的变形必须是连续的，即： 。三个节点的位移：三个节点的位移： ， ， 232332223223232232222222uKuKPuKuKP1P3P 2212uu 111uu 22122uuu 233uu 台阶轴处于平衡状态，各个节点的力是平衡的。台阶轴处于平衡状态，各个节点的力是平衡的。v节点节点1 1：外力：外力 ，节点力，节点力 ；v节

19、点节点2 2：外力为：外力为0 0，节点力，节点力 ；v节点节点3 3：外力：外力 ，节点力，节点力 。1P3P 02212 PP 11P 23P列出各节点的力的平衡方程式：列出各节点的力的平衡方程式：在节点在节点1 1处：处： 在节点在节点2 2处：处： 在节点在节点3 3处：处： 211211111211211111111uKuKuKuKPP 3223222212211212322322222122111212322322222121221112122120uKuKKuKuKuKKuKuKuKuKuKPP 32332232323322232233uKuKuKuKPP用矩阵表示：用矩阵表示：

20、简写为：简写为： 其中：其中： 为外载荷列阵；为外载荷列阵； 为节点力列阵。为节点力列阵。 3212332322232221221211121112322121131000uuuKKKKKKKKPPPPPP KFP19. 1 TPPP310 TPPPPF23221211 为整体刚度矩阵。为整体刚度矩阵。 为整体位移列阵。为整体位移列阵。 整体刚度方程是结构的力的平衡方程式，其中每整体刚度方程是结构的力的平衡方程式，其中每一行都表示了在一个节点上的力的平衡方程式。一行都表示了在一个节点上的力的平衡方程式。 工程应用中，采用叠加原理集成整体刚度矩阵。工程应用中，采用叠加原理集成整体刚度矩阵。 23

21、323222322212212111211100KKKKKKKKK20. 1 Tuuu321以台阶轴为例讲述叠加原理：以台阶轴为例讲述叠加原理： 首先列出单元的刚度方程：首先列出单元的刚度方程： e1:e1: e2:e2: 12111221211121111211uuKKKKPP 23222332322232222322uuKKKKPP将将2X22X2阶的单元刚度矩阵扩充为阶的单元刚度矩阵扩充为3X33X3阶的贡献矩阵：阶的贡献矩阵： 3211221121121111211000000uuuKKKKPP 3212332322232222322000000uuuKKKKPP将贡献矩阵叠加：将贡献

22、矩阵叠加： 在实际工程计算中，单元数目往往有上百个，无需将在实际工程计算中，单元数目往往有上百个，无需将每个单元刚度矩阵都扩充为贡献矩阵，而是按刚度系数的下每个单元刚度矩阵都扩充为贡献矩阵，而是按刚度系数的下标直接加到整体刚度矩阵中去。标直接加到整体刚度矩阵中去。 3212332322232221221211121112322121100uuuKKKKKKKKPPPP KF 23323222322212212111211100KKKKKKKKK整体刚度矩阵的性质：对称性： 。在计算中，只存贮矩阵的上三角部分或下三角部分。稀疏性：绝大多数元素都是零，这是因为和某一个节点相关的节点数一般不会超过9

23、个。整体网格分的越细，则 的稀疏性越突出，利用这个特点可设法只存贮 中的非零元素，从而节省存贮容量。 带形分布规律：非零元素分别在以主对角线为中心的带形区域内。 在包括对角元素在内的半个带形区域中，每行具有的元素数叫半带宽（d）。 d=(相邻节点号的最大差值1)*2 jiijKK K K 若结构有若结构有N N个节点，则迭加后的整体刚度矩阵为个节点，则迭加后的整体刚度矩阵为2N2N* *2N2N阶。阶。 例如，某结构节点数例如，某结构节点数N=55N=55，节点最大号码差，节点最大号码差D=6D=6，则半带宽则半带宽: d=(6+1): d=(6+1)* *2=142=14。 总自由度数总自由

24、度数: 2N=2: 2N=2* *55=11055=110。 半带宽存贮量：半带宽存贮量：2N2N* *d=1540d=1540 整体刚度矩阵上三角部分存贮量：整体刚度矩阵上三角部分存贮量： 1/2 1/2* *2N2N* *(2N+1)=6105(2N+1)=6105 奇异性：整体结构在无约束的条件下作刚体运动，奇异性：整体结构在无约束的条件下作刚体运动，必须处理边界条件。必须处理边界条件。 111123212321111221000000000000aaaaaaaaaaaaaaaaaaaaaaaannnnnnnnnndN2N21.7 处理边界条件处理边界条件 节点节点3 3固定，则固定，则 。x12311uP2u3u03u 32123323222322212212111211132100uuuKKKKKKKKPPP引入边界条件：引入边界条件：节点节点1 1处有外载荷处有外载荷 ，节点，节点2 2处外载荷为零，节点处外载荷为零，节点3 3处处位移为零，代入刚度方程：位