数值模拟：第五讲平面问题(二)——三角形单元分析

1、5.1 离散化要解决的问题离散化要解决的问题如何进行离散结构的求解？如何进行离散结构的求解？如何得到小单元的刚度特性？如何得到小单元的刚度特性？为什么离散化结构的解能作为原连续问题的近似解？为什么离散化结构的解能作为原连续问题的近似解？ 研究上述问题就涉及到有限元法的基本原理和基本理论。研究上述问题就涉及到有限元法的基本原理和基本理论。5 5 三节点三角形单元解平面问题三节点三角形单元解平面问题 目 标: 掌握平面问题 简单三角形单元位移模式。5.2 5.2 三节点三角形单元的特性分析三节点三角形单元的特性分析5.2.1 单元作为分析对象单元作为分析对象按前面结构矩阵位移法分析思想，要求解按前

2、面结构矩阵位移法分析思想，要求解平面问题的有限元离散结构，需要知道单平面问题的有限元离散结构，需要知道单元（三角形薄片）在节点自由度上受力时元（三角形薄片）在节点自由度上受力时的弹性特性或刚度特性。这是一个新问题，的弹性特性或刚度特性。这是一个新问题，一个特殊的弹性力学问题。一个特殊的弹性力学问题。 三角形顶点设为节点，其三角形顶点设为节点，其局部编号局部编号为为l，m，n（逆时针）。逆时针）。每节点有总体坐标每节点有总体坐标x，y方向两个待求位移分量：方向两个待求位移分量：u，v。单元共单元共有有6个位移分量个位移分量6个自由度。个自由度。下面研究有限元法中特有的求解该特下面研究有限元法中特

3、有的求解该特殊弹性力学问题的方法。殊弹性力学问题的方法。 有限元离散结构受力平衡后，取出一个典型三节点三角形单元有限元离散结构受力平衡后，取出一个典型三节点三角形单元e。 单元平衡时要在节点处受到节点力（节点对单元的作用单元平衡时要在节点处受到节点力（节点对单元的作用力），每节点有力），每节点有2个节点力分量，单元有个节点力分量，单元有6个节点力分量。个节点力分量。 （2） 单元节点力列阵单元节点力列阵 Tnnmmllnmlevuvuvu Tynxnymxmylxlnmlepppppppppp下面要研究的问题是下面要研究的问题是该三角形薄片弹性体在保持平衡时所受节该三角形薄片弹性体在保持平衡时

4、所受节点力和节点位移的关系点力和节点位移的关系。（1）单元节点位移列阵）单元节点位移列阵5.2.2 单元位移模式单元位移模式按弹性力学位移法求近似解的思路，位移作为基本未知量时，需按弹性力学位移法求近似解的思路，位移作为基本未知量时，需要要对单元上位移的分布作出假设对单元上位移的分布作出假设，即构造含待定参量的简单位移即构造含待定参量的简单位移函数函数位移模式位移模式。yaxaayxvyaxaayxu654321),(),(为待定系数，称为广义坐标为待定系数，称为广义坐标。61 aa通常用多项式函数作位移模式，通常用多项式函数作位移模式，对三节点三角形单元对三节点三角形单元，有，有6个待定个待

5、定节点位移分量，所以单元上的位移函数只能是节点位移分量，所以单元上的位移函数只能是含含6个待定系数的完个待定系数的完全一次多项式全一次多项式：3211aaayxu321111aaayxyxyxuuunnmmllnml654111aaayxyxyxvvvnnmmllnml位移多项式写成矩阵形式位移多项式写成矩阵形式：代入各节点取值条件后代入各节点取值条件后：坐标取节点值坐标取节点值6541aaayxv坐标取节点值坐标取节点值 由于有限元法中未知量是节点位移，所以上面单由于有限元法中未知量是节点位移，所以上面单元位移模式需要转换为以节点位移分量为待定参量的元位移模式需要转换为以节点位移分量为待定参

6、量的形式。过程如下形式。过程如下：321111aaayxyxyxuuunnmmllnmlnmlnmlnmlnmlnmlnnmmlluuucccbbbaaauuuyxyxyxaaa211111321其中：其中：nnmmllyxyxyx1112为三角形面积由第一组方程求解由第一组方程求解:31aa节点坐标行列式节点坐标行列式分别解出分别解出6个待定系数个待定系数：654111aaayxyxyxvvvnnmmllnmlnmlnmlnmlnmlnmlnnmmllvvvcccbbbaaavvvyxyxyxaaa211111654为常数。元素的代数余子式，均个，行第，的第分别是节点坐标行列式321)(,n

7、mlkkcbakkk由第二组方程求解由第二组方程求解:64aa上面求出的待定系数上面求出的待定系数 代回位移多项式，得到代回位移多项式，得到：61 aanmlnmlnmlnmluuucccbbbaaayxaaayxu2111321nmliiinnmmllvNvNvNvNv,nmlnnnmmmllluuuycxbaycxbaycxba21nmlnmluuuNNNnmliiinnmmlluNuNuNuN,至此，单元位移模式已转换为节点位移的插值形式。至此，单元位移模式已转换为节点位移的插值形式。上式中上式中：)(21ycxbaNiiii），（nmli 位移插值基函数，称为位移插值基函数，称为形状函

8、数（形函数）形状函数（形函数） enmlnmlNNNNNNNvuvu000000合并后用矩阵表示为：和 N 称为形函数矩阵，是对单元节点位移进行插值得到单元位移称为形函数矩阵，是对单元节点位移进行插值得到单元位移分布函数的转换矩阵分布函数的转换矩阵。 用节点位移插值表示单元位移模式是有限元法中除了离用节点位移插值表示单元位移模式是有限元法中除了离散化之外最具代表性，最重要的步骤散化之外最具代表性，最重要的步骤!5.2.3 形函数及其性质形函数及其性质对于单元位移模式对于单元位移模式：nnmmllnnmmllvNvNvNvuNuNuNu)(21ycxbaNiiii），（nmli 假设假设：01n

9、mluuu得到得到：lNyxu),(形函数。对应一个形函数。单元每个节点节点的形状函数，简称因此称为的位移分布，节点固定不动时，单元节点发生单位位移，是lnmlNl,显然，显然，形函数决定了单元上位移分布的形态形函数决定了单元上位移分布的形态。事实上，单元位。事实上，单元位移模式就是所有形函数的线性组合。移模式就是所有形函数的线性组合。一个单元的位移模式决定了该单元描述局部位移场的能力，决定一个单元的位移模式决定了该单元描述局部位移场的能力，决定求解的精度、收敛性等，而形函数是最重要的因素。求解的精度、收敛性等，而形函数是最重要的因素。 针对三节点三角形单元，可以导出单元形函数针对三节点三角形

10、单元，可以导出单元形函数的下列性质。的下列性质。性质性质1：单元上某节点的形函数在该节点的值为单元上某节点的形函数在该节点的值为1，在其它节点，在其它节点 的值为零的值为零。性质性质2：单元上所有形函数之和等于单元上所有形函数之和等于1。0),(0),(1),(nnlmmllllyxNyxNyxN)(nml，1nmlNNN（简单三角形单元的形函数只有（简单三角形单元的形函数只有2个独立）个独立）性质性质3（推论）（推论）：简单三角形单元的形函数在边界上的性质。简单三角形单元的形函数在边界上的性质。 某节点的形函数在该点邻边上呈线性分布，取值在某节点的形函数在该点邻边上呈线性分布，取值在01之间

11、，之间，在该点对边上值为零在该点对边上值为零。简单三角形单元形函数的几何意义简单三角形单元形函数的几何意义v 由形函数表达式和性质由形函数表达式和性质1可画出下列形函数几何图形。可画出下列形函数几何图形。v根据位移模式表达式及其形函数的性质，可以推断出两个相邻根据位移模式表达式及其形函数的性质，可以推断出两个相邻三角形单元上位移分布形状和公共边界上位移的情况：三角形单元上位移分布形状和公共边界上位移的情况：nnmmllnnmmllvNvNvNvuNuNuNu两个单元上位移线性连续分布，各单两个单元上位移线性连续分布，各单元在公共边界上位移线性分布，数值元在公共边界上位移线性分布，数值相同相同边

12、界位移协调！边界位移协调！v 由图形几何性质可以推断简单三角形单元形函数的下列结论：由图形几何性质可以推断简单三角形单元形函数的下列结论：1）三角形形心上三角形形心上：31nmlNNNlmllmlldsNdsNlmlm212）形函数在边界上的积分形函数在边界上的积分：3）形函数在单元上的积分形函数在单元上的积分：3),(AdxdyyxNi),(nmli 5.2.4 单元应变和应力单元应变和应力已知节点位移插值形式的单元位移模式已知节点位移插值形式的单元位移模式： eNvu代入平面问题几何方程（应变代入平面问题几何方程（应变位移关系）得到单元应变：位移关系）得到单元应变： eenmlennmml

13、lnmlnmlBBBBxNyNxNyNxNyNyNyNyNxNxNxN000000 enmlnmlxyyxNNNNNNxyyxvuxyyx0000000000上式简写为上式简写为： eB nnmmllnmlnmlbcbcbccccbbbB00000021常应变三角形单元。常数，故该单元又称为应变在单元上是是常数矩阵，该单元的所以该单元的应变矩阵关的常数，均为与单元节点坐标有、由于上式中),(nmlicbii 为应变矩为应变矩阵，其一个阵，其一个子块的计算子块的计算式为式为：B),(00nmlixNyNyNxNBiiiii该式建立了用单元节点位移表该式建立了用单元节点位移表达单元上应变分布的关系

14、。达单元上应变分布的关系。对简单三角形单元，应变矩阵为对简单三角形单元，应变矩阵为：把单元应变代入平面问题物理方程，即得到单元应力把单元应变代入平面问题物理方程，即得到单元应力： Dxyyx eS 平面问题弹性矩阵单元应力矩阵DBDS eBv应力矩阵的分块形式应力矩阵的分块形式： nmlnmlSSSBBBDS iiBDS ),(nmli eeSBDv对于平面应力问题，应力矩阵的子块为： iiiiiiibccbcbES2121)1 (22),(nmli v由于均质材料的弹性系数为常数，应力矩阵也是常数矩阵，即该单元上应力分布也是常数。对平面应变问题的应力矩阵，只要把上式中弹性系数作相应变换：21

15、EEu1u 平面问题简单三角形单元应力的讨论平面问题简单三角形单元应力的讨论单元上应力和应变均为常数单元上应力和应变均为常数，其大小与单元几何和节点位移，其大小与单元几何和节点位移有关，一般各单元应力值不相同，应力在单元之间不连续，有关，一般各单元应力值不相同，应力在单元之间不连续，这是这是有限元解的近似性的反映有限元解的近似性的反映。一般把单元上这个常应力值。一般把单元上这个常应力值作为单元中心的应力值较合理，有较高的精度。作为单元中心的应力值较合理，有较高的精度。理论和数值试验均可证明，用该单元求解问题时，误差随单元理论和数值试验均可证明，用该单元求解问题时，误差随单元尺寸的减小和单元数目

16、的增加而减小，这就是尺寸的减小和单元数目的增加而减小，这就是有限元解的收敛有限元解的收敛性性。事实上，当单元变得越来越小时，结构在一个小单元区域。事实上，当单元变得越来越小时，结构在一个小单元区域上的应力趋于均匀，而大量小单元拼接在一起就可能很精确地上的应力趋于均匀，而大量小单元拼接在一起就可能很精确地模拟结构中变化的应力。模拟结构中变化的应力。5.2.5 单元刚度方程与刚度矩阵单元刚度方程与刚度矩阵前面已经得到了用单元节点前面已经得到了用单元节点位移表示的单元内部位移、位移表示的单元内部位移、应变、应力的分布：应变、应力的分布： Tnnmmllevuvuvu Tynxnymxmylxlepp

17、ppppp eNvu eBD eB结构受力平衡时，每个单元在节点力结构受力平衡时，每个单元在节点力 作作用下也保持平衡，节点产生相应位移用下也保持平衡，节点产生相应位移 。 e ep下面下面用弹性体虚位移原理（虚功原理）建立单元的平衡关系用弹性体虚位移原理（虚功原理）建立单元的平衡关系。v设单元节点虚位移：设单元节点虚位移： Tnnmmllevuvuvu*则相应有单元内部虚位移则相应有单元内部虚位移： eNvu*相应单元虚应变为：相应单元虚应变为： eB*v外力虚功计算外力虚功计算： eTeynnxnnymmxmmyllxllppvpupvpupvpuW)(*对一个单元而言，外力即单元节点力对

18、一个单元而言，外力即单元节点力。v按虚功原理，虚位移下单元的外力虚功等于虚应变能按虚功原理，虚位移下单元的外力虚功等于虚应变能。v虚应变能计算虚应变能计算单元虚应变能即单元内应力在虚应变上所作虚功。单元虚应变能即单元内应力在虚应变上所作虚功。 dVdVUTVxyxyyyxVxee*)( eB* dVBDBdVUeTVeVTee)(* eBD eVeTTedVBDB)(* eVTTeedVBDB)(*v由虚功原理得到由虚功原理得到：根据上面推导，有：根据上面推导，有： eTepW)(* eVTTeedVBDBU)(*外力虚功外力虚功：虚应变能虚应变能：UW eVTTeeTeedVBDBp)()(

19、*考虑到虚位移考虑到虚位移 的任意性，由上式立即得到下列方程：的任意性，由上式立即得到下列方程：e* eVTeedVBDBp上式简写成： eeekp其中单元刚度矩阵： eVTedVBDBk该方程建立了单元节点力和节点位移之间的关系，为单元刚度方程。至此，我们已经得到了关于弹性实体小单元的力学特性方程。由于弹性矩阵 是对称的，由上式可以看出，单元刚度矩阵 是对称方阵。对3节点三角形单元，单元刚度矩阵是66 方阵。 D ek上面单元刚度方程和单元刚度矩阵的推导虽然在平面问题的简单三角形单元下得到，但推导的过程和公式的形式具有一般性，推导原理和结果适用一般的连续体力学问题的所有单元。提提 示：示：单

20、元刚度方程和单元刚度矩阵的建立是单元分析的核心内容。一般情况下，单元应变矩阵是坐标的函数矩阵，所以单元刚度矩阵的计算需要进行积分运算。所建立的单元刚度矩阵反映了一般弹性体小单元近似的弹性性质，是单元特性的核心。 eeekp单元刚度矩阵的计算单元刚度矩阵的计算弹性力学平面问题的单元刚度矩阵通式弹性力学平面问题的单元刚度矩阵通式： eVTedVBDBk（单元刚度矩阵通式） ehdxdyBDBkTe（平面问题单刚通式）单元厚度（常数）h由于应变矩阵是常数矩阵，所以由平面问题刚度矩阵通式得到：平面问题简单三角形单元刚度矩阵计算 nmlTnTmTlTeBBBDBBBhhBDBk nnnmnlmnmmml

21、lnlmllekkkkkkkkkk33分块形式对平面应力问题，上述刚度矩阵的一个22子矩阵为： srsrsrsrsrsrsrsrsTrersbbcccbbcbccbccbbEhhBDBk21212121)1 (42）（nmlsr, 单元刚度方程的分块形式 eeekp nnnmnlmnmmmllnlmllekkkkkkkkkk nmle nmlepppp简单三角形单元刚度方程按节点分块形式：nmlnnnmnlmnmmmllnlmllnmlkkkkkkkkkppp展开的单元刚度方程： ),(nmlikkkpninmimlili每行各节点刚度子块反映了节点位移对节点力的贡献。5.2.6 简单三角形单

22、元分析计算实例简单三角形单元分析计算实例(1)(1) 形函数和形函数矩阵前面针对一维杆单元和平面问题简单三角形单元引入的单元形函数，应变、应力矩阵，单元刚度矩阵等是结构有限元分析中最基本、最重要的概念。dd厚度h下面通过分析计算一个具体的单元加以说明。按公式计算各形函数按公式计算各形函数)(21ycxbaNiiii），（nmli 200101011112dddyxyxyxnnmmlldddddcccbbbaaanmlnmlnml00002代数余子式矩阵：按形函数公式，得到三个形函数：计算节点坐标行列式：dd厚度h容易验证其满足形函数性质dydxNdyNdxNnml1按形函数性质和几何意义，直接

23、写出形函数表达式按形函数性质和几何意义，直接写出形函数表达式dd厚度h形函数 的图形是一个平面，其方程是：lNdxNl上述两个函数定义在三角形区域的部分就是形函数。同理，直接写出形函数 图形的方程为：mNdyNm由形函数的性质，得到第三个形函数：dydxNNNmln11显然，与公式计算得到的结果相同！3）形函数矩阵形函数矩阵根据计算所得单元形函数，写出单元的形函数矩阵： dydxdydxdydxdydxN10000100(2) (2) 应变矩阵应变矩阵由公式计算单元应变矩阵： 1101101010000100011dB),(00nmlixNyNyNxNBiiiii nnmmllnmlnmlbcbcbccccbbbB00000021dddddcccbbbaaanmlnmlnml00002(3) (3) 应力矩阵应力矩阵 212102121011001001)1 (2dESdddddcccbbbaaanmlnmlnml00002 iiiiiiibccbcbES2121)1 (22),(nmli 对平面应力问题：(4) (4) 单元刚度矩阵计算单元刚度矩阵计算dddddcccbbbaaanmlnmlnml00002 srsrsrsrsrsrsrsrsTrersbbcccbbcbccbccbbEhhBDBk21212121)1 ( 42）（nmlsr, 平