三角形单元有限元弹性力学

1、结构有限元法结构有限元法第第1章章 三角形常应变单元的有限元法三角形常应变单元的有限元法第第2章章 有限元程序设计与分析软件有限元程序设计与分析软件第第3章章 平面问题高阶单元的有限元法平面问题高阶单元的有限元法第第4章章 空间实体的有限元法空间实体的有限元法第第5章章 杆系结构的有限元法杆系结构的有限元法第第6章章 板壳问题的有限元法板壳问题的有限元法第第7章章 结构动力问题的有限元法？结构动力问题的有限元法？第第8章章 弹塑性问题的有限元法弹塑性问题的有限元法结结 构构 有有 限限 元元 分分 析析 第第1 1章章 三角形单元的有限元法三角形单元的有限元法1.1 有限元法的基本思想有限元法

2、的基本思想 有限元法在有限元法在20世纪世纪50年代起源于飞机结构的矩阵年代起源于飞机结构的矩阵分析，其分析，其基本思想基本思想是用有限个离散单元的集合体代是用有限个离散单元的集合体代替原连续体，采用能量原理研究单元及其离散集合替原连续体，采用能量原理研究单元及其离散集合体的平衡，以计算机为工具进行结构数值分析。它体的平衡，以计算机为工具进行结构数值分析。它避免了经典弹性力学获得连续解的困难（建立和求避免了经典弹性力学获得连续解的困难（建立和求解偏微分方程），使大型、复杂结构的计算容易地解偏微分方程），使大型、复杂结构的计算容易地在计算机上完成，应用十分广泛。在计算机上完成，应用十分广泛。AN

3、SYS, SAP2K 把整体结构离散为有限个单元，研究单元的平把整体结构离散为有限个单元，研究单元的平衡和变形协调；再把这有限个离散单元集合还原成衡和变形协调；再把这有限个离散单元集合还原成结构，研究离散结构的平衡和变形协调。结构，研究离散结构的平衡和变形协调。 划分的单元大小和数目根据计算精度和计算机能划分的单元大小和数目根据计算精度和计算机能力来确定。力来确定。12345678910P5764 56345678 弹性悬臂板弹性悬臂板剖分与集合剖分与集合单元、节点需编号单元、节点需编号有限元法主要优点：有限元法主要优点：（1） 概念浅显，容易掌握。（离散、插值、能概念浅显，容易掌握。（离散、

4、插值、能量原理、数学分析）量原理、数学分析）（2）适用性强，应用范围广，几乎适用于所有）适用性强，应用范围广，几乎适用于所有连续体和场问题的分析。（结构、热、流体、连续体和场问题的分析。（结构、热、流体、电磁场和声学等问题）电磁场和声学等问题）（3）计算规格化（采用矩阵表示），便于计算）计算规格化（采用矩阵表示），便于计算机编程。机编程。1.1.1 有限元法的分析步骤有限元法的分析步骤 （1）结构）结构离散化：离散化：用点、线或面把结构剖分为有用点、线或面把结构剖分为有限个离散单元体，并在单元指定点设置限个离散单元体，并在单元指定点设置节点节点。研究。研究单元的平衡和变形协调，形成单元平衡方程

5、。单元的平衡和变形协调，形成单元平衡方程。l/2l/2P123 1、F1 2、F2 3、F3 4、F4l/212l/223 1、F1 2、F2 3、F3 4、F4单元的单元的节点上节点上有位移有位移 和力和力F（2）单元集合：把所有离散的有限个单元集合起来）单元集合：把所有离散的有限个单元集合起来代替原结构，形成离散结构节点平衡方程。代替原结构，形成离散结构节点平衡方程。（3）由平衡方程求解得节点位移和计算单元应力。）由平衡方程求解得节点位移和计算单元应力。 1、F1 2、F2 3、F3 4、F4l/212l/223 1、F1 2、F2 3、F3 4、F4l/2l/2P1231.1.2 有限元

6、法分析思路流程有限元法分析思路流程解综合方程解综合方程K= P求结构节点位移求结构节点位移计算结构内力和应力计算结构内力和应力系统分析系统分析(把单元刚度矩阵集合成结构刚度矩阵把单元刚度矩阵集合成结构刚度矩阵K形成等价节点荷载形成等价节点荷载P )离散（剖分）结构离散（剖分）结构为若干单元为若干单元单元分析单元分析(建立单元刚度矩阵建立单元刚度矩阵ke形成单元等价节点力形成单元等价节点力)（1-1）Tsysxsysxsqqqqq2、单元内任意点的、单元内任意点的体积力体积力列阵列阵 qV （1-2）TVyVxVyVxVqqqqq1、单元表面或边界上任意点的、单元表面或边界上任意点的表面力表面力

7、列阵列阵 qs ijmxyijmxyqVqs1.2 基本力学量矩阵表示基本力学量矩阵表示图图1-1ijmxyuv3、单元内任意点的位移列阵、单元内任意点的位移列阵 f Tuf（1-3） 4、单元内任意点的应变列阵、单元内任意点的应变列阵 Txyyx（1-4）ijmxy5、单元内任意点的应力列阵、单元内任意点的应力列阵 Txyyx（1-5）6、几何方程几何方程Txyuyxu（1-6）xvyuyvxuxyyx,将上式代入式将上式代入式(1-4），），ijmxyTxyyx（1-4）7、物理方程矩阵式、物理方程矩阵式xyyxxyyxE21001112称对（1-7）式中式中 E、 弹性模量、泊松比。弹性

8、模量、泊松比。上式可简写为上式可简写为D（1-8）其中其中 对于弹性力学的平面应力问题，对于弹性力学的平面应力问题，物理方程物理方程的矩阵形的矩阵形式可表示为：式可表示为：)(12yxxE21001112称对ED(1-9）矩阵矩阵D称为弹性矩阵。称为弹性矩阵。对于平面应变问题，将式对于平面应变问题，将式(1-9）中的中的E换为换为 ， 换为换为 。21E1D(1-8） 各种类型结构的弹性物理方程都可用式各种类型结构的弹性物理方程都可用式(1-8）描）描述。但结构类型不同，力学性态述。但结构类型不同，力学性态 (应力分量、应变分应力分量、应变分量量)有区别，有区别， 弹性矩阵弹性矩阵D的体积和元

9、素是不同的。的体积和元素是不同的。1.3 位移函数和形函数位移函数和形函数 1、位移函数概念、位移函数概念 由于有限元法采用能量原理进行单元分析，因而由于有限元法采用能量原理进行单元分析，因而必须事先设定位移函数。必须事先设定位移函数。 “位移函数位移函数”也称也称 “位移位移模式模式”，是，是单元内部位移变化的数学表达式，设为坐单元内部位移变化的数学表达式，设为坐标的函数标的函数。一般而论，位移函数选取会影响甚至严重影响一般而论，位移函数选取会影响甚至严重影响计算结果的精度。在弹性力学中，恰当选取位移函数计算结果的精度。在弹性力学中，恰当选取位移函数不是一件容易的事情；但不是一件容易的事情；

10、但在有限元中，当单元划分得在有限元中，当单元划分得足够小时，把位移函数设定为简单的多项式就可以获足够小时，把位移函数设定为简单的多项式就可以获得相当好的精确度。得相当好的精确度。这正是有限单元法具有的重要优这正是有限单元法具有的重要优势之一。势之一。 不同类型结构会有不同的位移函数。这里，仍不同类型结构会有不同的位移函数。这里，仍以平面问题三角形单元（图以平面问题三角形单元（图1-2）为例，说明设定位）为例，说明设定位移函数的有关问题。移函数的有关问题。 图图1-2是一个三节点三角形是一个三节点三角形单元，其节点单元，其节点i、j、m按按逆时针逆时针方向排列。每个节点位移在单方向排列。每个节点

11、位移在单元平面内有两个分量：元平面内有两个分量：),(mjiuTiii（1-10） 一个三角形单元有一个三角形单元有3个节点（以个节点（以 i、j、m为为 序），序），共有共有6个节点位移分量。其个节点位移分量。其单元位移或单元节点位移单元位移或单元节点位移列阵列阵为：为：图图1-2ijmuiujumvivjvmxy2、位移函数设定、位移函数设定 本问题选位移函数（单元中任意一点的位移与节点本问题选位移函数（单元中任意一点的位移与节点位移的关系）为简单多项式：位移的关系）为简单多项式：yaxaayaxaau654321(1-12）式中：式中：a1、a2、a6待定常数，由单元位移的待定常数，由单

12、元位移的6个分量确定。个分量确定。a1、a4代表刚体位移，代表刚体位移，a2、 a3 、 a5 、 a6 代表单元中的常应变，而且，位移函数是代表单元中的常应变，而且，位移函数是连续函数。连续函数。 Tmmjjiimjiuuu(1-11）ijmuiujumvivjvmxyuv625352,aaaaayvaxuxyyxu选取位移函数应考虑的问题选取位移函数应考虑的问题 （1）位移函数的个数位移函数的个数等于等于单元中任意一点的位移分量个数。本单元中单元中任意一点的位移分量个数。本单元中有有u和和v，与此相应，有，与此相应，有2个位移函数；个位移函数； （3）位移函数中待定常数个数位移函数中待定常

13、数个数 待定常数个数应等于待定常数个数应等于单元节点自由度总数单元节点自由度总数，以，以便用单元节点位移确定位移函数中的待定常数。本便用单元节点位移确定位移函数中的待定常数。本单元有单元有6个节点自由度，两个位移函数中共包含个节点自由度，两个位移函数中共包含6个个待定常数。待定常数。（2）位移函数是坐标的函数位移函数是坐标的函数 本单元的坐标系为：本单元的坐标系为：x、y； （4）位移函数中必须包含单元的刚体位移。位移函数中必须包含单元的刚体位移。 （5）位移函数中必须包含单元的常应变。位移函数中必须包含单元的常应变。 （6）位移函数在单元内要连续。相邻单元间要尽位移函数在单元内要连续。相邻单

14、元间要尽 量协调。量协调。 条件（条件（4）、（）、（5）构成单元的）构成单元的完备性完备性准则。准则。 条件（条件（6）是单元的位移）是单元的位移协调性协调性条件。条件。 理论和实践都已证明，完备性准则是有限元解收理论和实践都已证明，完备性准则是有限元解收敛于真实解的必要条件。单元的位移协调条件构成有敛于真实解的必要条件。单元的位移协调条件构成有限元解收敛于真实解的充分条件。限元解收敛于真实解的充分条件。 容易证明，三角形三节点常应变单元满足以上必容易证明，三角形三节点常应变单元满足以上必要与充分条件。要与充分条件。 （7）位移函数的形式位移函数的形式 一般选为完全多项式。为实现（一般选为完

15、全多项式。为实现（4）（6）的要）的要求，根据求，根据Pascal三角形由低阶到高阶按顺序、对称地三角形由低阶到高阶按顺序、对称地选取；多项式的项数等于（或稍大于）单元节点自由选取；多项式的项数等于（或稍大于）单元节点自由度数。度数。4322343223221yxyyxyxxyxyyxxyxyxyxyaxaayaxaau654321例：平面应力矩形板被划分为若干三角形单元。例：平面应力矩形板被划分为若干三角形单元。 位移函数中包含了单元的常应变。位移函数中包含了单元的常应变。 xvyuyvxuxyyx,(a2, a6, a3+a5 ) 位移函数中包含了单元的刚体位移。位移函数中包含了单元的刚体

16、位移。 (a1, a4 )254136对任一单元，如对任一单元，如单元，取位移函数：单元，取位移函数：、单元的位移函数都是单元的位移函数都是yaxaayaxaau654321可以看出：可以看出： 位移函数在单元内是连续的；位移函数在单元内是连续的； 以以、的边界的边界26为例为例2562635623xyuu6u2uu6u2两条直线上有两个点重合，此两条直线必全重合。两条直线上有两个点重合，此两条直线必全重合。位移函数在单元之间的边界上也连续吗？是。位移函数在单元之间的边界上也连续吗？是。3、形函数、形函数 形函数是用单元节点位移分量来描述位移函数的形函数是用单元节点位移分量来描述位移函数的插值

17、函数插值函数。iiiiiiyaxaayaxaau654321jjjjjjyaxaayaxaau654321mmmmmmyaxaayaxaau654321(1-13） （1）形函数确定）形函数确定 现在，通过单元节点位移确定位移函数中的待定现在，通过单元节点位移确定位移函数中的待定常数常数a1、a2、a6 。设节点。设节点i、j、m的坐标分别为的坐标分别为（xi、yi）、（）、（ xj、yj ）、（）、（ xm、ym ），节点位移分别），节点位移分别为（为（ui、vi）、）、 （uj、vj） 、 （um、vm）。将它们）。将它们代入式（代入式（1-12），有），有)121 (654321yaxa

18、ayaxaau从式从式(1-13）左边）左边3个方程中解出待定系数个方程中解出待定系数a1、a2、a3为为 mmmjjjiiiyxuyxuyxuAa211mmjjiiyuyuyuAa111212mmjjiiuxuxuxAa111213(1-14） 式中，式中， A为三角形单元的面积，有为三角形单元的面积，有 mmjjiiyxyxyxA11121(1-15） 特别指出：特别指出：为使求得面积的值为正值，本单元节点号为使求得面积的值为正值，本单元节点号的次序必须是的次序必须是逆时针逆时针转向，如图所示。至于将哪个节转向，如图所示。至于将哪个节点作为起始节点点作为起始节点i，则没有关系。，则没有关系

19、。 将式将式(1-14）代入式）代入式(1-12）的第一式，整理后得）的第一式，整理后得)()()(21mmmmjjjjiiiiuycxbauycxbauycxbaAu同理同理)()()(21mmmmjjjjiiiiycxbaycxbaycxbaAijmxy（2）（1）（7）)()()(21mmmmjjjjiiiiuycxbauycxbauycxbaAu)()()(21mmmmjjjjiiiiycxbaycxbaycxbaA(1-16）式中式中 ),(mjijmmjiyxyxamjiyybmjixxc(1-17） ijm式式(1-17）中（）中（i、j、m）意指：按）意指：按i、j、m依次轮换

20、下依次轮换下标，可得到标，可得到aj、bj、cjam、bm、cm。后面出现类似情。后面出现类似情况时，照此推理。式况时，照此推理。式(1-17）表明：）表明： aj、bj、cjam、bm、cm是单元三个节点坐标的函数。是单元三个节点坐标的函数。)()()(21mmmmjjjjiiiiuycxbauycxbauycxbaAu)()()(21mmmmjjjjiiiiycxbaycxbaycxbaA(1-16）令令 )(21ycxbaANiiii),(mji(1-18） 位移模式位移模式(1-16）可以简写为）可以简写为(1-19） mmjjiimmjjiiNNNuNuNuNu 式式(1-19）中的

21、）中的Ni、Nj、Nm是坐标的函数，反应了是坐标的函数，反应了单元的位移形态，称为单元位移函数的形函数。数学单元的位移形态，称为单元位移函数的形函数。数学上它反应了节点位移对单元内任一点位移的插值，又上它反应了节点位移对单元内任一点位移的插值，又称插值函数。称插值函数。 )()()(21mmmmjjjjiiiiuycxbauycxbauycxbaAu)()()(21mmmmjjjjiiiiycxbaycxbaycxbaA(1-16）用形函数把式用形函数把式(1-16）写成矩阵，有）写成矩阵，有mmjjiimjimjivuvuvuNNNNNNvu000000缩写为缩写为Nf(1-20） 形函数是

22、有限单元法中的一个重要函数，它具形函数是有限单元法中的一个重要函数，它具有以下性质：有以下性质：N为形函数矩阵，写成分块形式：为形函数矩阵，写成分块形式：mjiNNNN (1-21）其中子矩阵其中子矩阵),(00mjiINNNNiiii(1-22）I是是22的单位矩阵。的单位矩阵。 （2）形函数性质）形函数性质性质性质1 形函数形函数Ni在节点在节点i上的值等于上的值等于1，在其它节点，在其它节点 上的值等于上的值等于0。对于本单元，有。对于本单元，有 0),(0),(1),(mmijjiiiiyxNyxNyxN（i、j、m）性质性质2 在单元中任一点，所有形函数之和等于在单元中任一点，所有形

23、函数之和等于1。对。对 于本单元，有于本单元，有1),(),(),(yxNyxNyxNmjimmjjiimmjjiiNNNuNuNuNuxyN（i，j，m）Ni =1ijm图图1-3？公式证明、和利用iiiiiiicbaycxbaAN)(21xyN（I，j，m）Ni =1ijmNj =1ijmNm =1ijmNi =1ijmNj =1Nm =1图图1-4也可利用行列式代数余子式与某行或列元素也可利用行列式代数余子式与某行或列元素乘积的性质（等于行列式值或乘积的性质（等于行列式值或0）证明。）证明。性质性质3 在三角形单元的边界在三角形单元的边界ij上任一点（上任一点（x，y），有），有 0),

24、(),(1),(yxNxxxxyxNxxxxyxNmijijijiixxixjxyNi(xi，yi)j (xj，yj)m (xm，ym)Ni(x、y)1jijixxxxyxN1),(ijijijiijijixxxxxxxxxxxxxxyxN1),(证证图图1-50)(21)(21),()() 3(imimmmmmmiimmycxbaAycxbaAyxNyxxcbyij方程：1),(),(),()2(yxNyxNyxNmji（1）性质性质4 形函数在单元上的面积分和在边界上的线积分形函数在单元上的面积分和在边界上的线积分公式为公式为 ijdlNAdxdyNijiAi213(1-23）式中式中 为

25、为 边的长度。边的长度。 ijij1.4 单元应变和应力单元应变和应力 根据几何方程根据几何方程(1-6）和位移函数）和位移函数(1-16）可以求）可以求得单元应变。得单元应变。 1、单元应变、单元应变xvyuyvxuxyyx(1-6）对位移函数（式对位移函数（式(1-16）mmjjiimmmmjjjjiiiixyyxuuubccbbccbbccbA00000021(1-24）)()()(21mmmmjjjjiiiiuycxbauycxbauycxbaAu)()()(21mmmmjjjjiiiiycxbaycxbaycxbaA(1-16）求导后代入式求导后代入式(1-6），得到应变和节点位移的

26、关系式。），得到应变和节点位移的关系式。上式简写一般式：上式简写一般式： B(1-25）式中，式中， B单元应变矩阵。单元应变矩阵。 对本问题，维数为对本问题，维数为36。它的分块形式为：。它的分块形式为：mjiBBBB 子矩阵子矩阵 ),(0021mjibccbABiiiii(1-26） 由于由于 与与x、y无关，都是常量，因此无关，都是常量，因此B矩阵也是常量。单元中任一点的应变分量是矩阵也是常量。单元中任一点的应变分量是B矩阵矩阵与单元位移的乘积，因而也都是常量。因此，这种单元与单元位移的乘积，因而也都是常量。因此，这种单元被称为常应变单元。被称为常应变单元。mmjjiicbcbcbA,

27、 2、单元应力、单元应力 将式将式(1-25）代入物理方程式）代入物理方程式(1-8），得），得 单元应力单元应力DBD(1-27）也可写为也可写为 S(1-28）其中：其中：S称为称为单元应力矩阵单元应力矩阵，并有，并有(1-29）633363BDS 这里，这里，D是是33矩阵，矩阵，B是是36矩阵，因此矩阵，因此S也是也是36矩阵。它可写为分块形式矩阵。它可写为分块形式 BmjiSSSS (1-30）将弹性矩阵（式将弹性矩阵（式(1-9） 和应变矩阵（式和应变矩阵（式(1-26）代）代入，得子矩阵入，得子矩阵Si由式由式(1-29）iiBDS ),(2121)1 (22mjibcccbbA

28、ESiiiiiii(1-31）式式(1-31）是平面应力的结果。对于平面应变问题，）是平面应力的结果。对于平面应变问题，只要将上式中的只要将上式中的E换成换成 ， 换成换成 即得。即得。21E1),()1 (221)1 (22111)21)(1 (2)1 (mjibccbcbAESiiiiiii(1-32） 由于同一单元中的由于同一单元中的D、B矩阵都是常数矩阵，矩阵都是常数矩阵，所以所以S矩阵也是常数矩阵。也就是说，三角形三节矩阵也是常数矩阵。也就是说，三角形三节点单元内的应力分量也是常量。点单元内的应力分量也是常量。 当然，相邻单元的当然，相邻单元的bi、ci(i，j，m)一般不完全相同，

29、一般不完全相同，因而具有不同的应力，这就造成在相邻单元的公共因而具有不同的应力，这就造成在相邻单元的公共边上存在着边上存在着应力突变现象应力突变现象。但是随着网格的细分，。但是随着网格的细分，这种突变将会迅速减小，收敛于平衡被满足。这种突变将会迅速减小，收敛于平衡被满足。1.5 单元平衡方程单元平衡方程 1、 单元应变能单元应变能 对于平面应力问题中的三角形单元，设单元厚度对于平面应力问题中的三角形单元，设单元厚度为为h 。ATAxyxyyyxxhdxdyhdxdyU21)(21将式将式(1-25）和）和(1-8）代入上式进行矩阵运算，并注）代入上式进行矩阵运算，并注意到弹性矩阵意到弹性矩阵D

30、的对称性，有的对称性，有AThdxdyDU21 ATThdxdyBDB21应变能应变能 U为为ijmxyh TTTTDD)(B(1-25)D(1-8）由于由于和和T是常量，提到积分号外，上式可写成是常量，提到积分号外，上式可写成 21ATThdxdyBDBU引入矩阵符号引入矩阵符号k，且有，且有AThdxdyBDBk(1-33a）式式(1-33a）是针对平面问题三角形单元推出的。注意）是针对平面问题三角形单元推出的。注意到其中到其中hdxdy的实质是任意的微体积的实质是任意的微体积dv，于是得计算，于是得计算k的的一般式一般式。 dvBDBkvT(1-33） 式式(1-33）不仅适合于平面问题

31、三角形单元，）不仅适合于平面问题三角形单元，也是计算各种类型单元也是计算各种类型单元k的一般式。的一般式。 ATThdxdyBDBU21dv 1.6节中将明确节中将明确k的力学意义是单元刚度矩阵。的力学意义是单元刚度矩阵。式式(1-33）便是计算单元刚度矩阵的基本矩阵式。它）便是计算单元刚度矩阵的基本矩阵式。它适合于各种类型的单元。适合于各种类型的单元。单元应变能写成单元应变能写成 21kUT(1-34） 2、 单元外力势能单元外力势能 单元受到的外力一般包括单元受到的外力一般包括体积力体积力、表面力表面力和和集中力集中力。自重属于体积力范畴。表面力指作用在单元表面的分自重属于体积力范畴。表面

32、力指作用在单元表面的分布载荷，如风力、压力，以及相邻单元互相作用的内布载荷，如风力、压力，以及相邻单元互相作用的内力等。力等。 dvBDBkvT(1-33） （1） 体积力势能体积力势能 单位体积中的体积力单位体积中的体积力如式如式(1-35）所示。）所示。单元上体积力具有的势能单元上体积力具有的势能Vv为为AVTvhdxdyqfV(1-35）TVyVxVqqqijmxyqVxqVyijmxyuv注意到式注意到式(1-20）AVTTAVTvhdxdyqNhdxdyqNV)(有有Nf(1-20） （2） 表面力势能表面力势能 面积力虽然包括单元之间公共边上互相作用的分布面积力虽然包括单元之间公共

33、边上互相作用的分布力，但它们属于结构内力，成对出现，集合时互相抵力，但它们属于结构内力，成对出现，集合时互相抵消，在结构整体分析时可以不加考虑，因此单元分析消，在结构整体分析时可以不加考虑，因此单元分析时也就不予考虑。时也就不予考虑。 现在，只考虑弹性体边界上的表面力，它只在部现在，只考虑弹性体边界上的表面力，它只在部分单元上形成表面力（右下图）。设分单元上形成表面力（右下图）。设边界面上单位面边界面上单位面积受到的表面力积受到的表面力如下式：如下式：dAqNhdlqfVlSTTlSTSl单元边界长度单元边界长度h单元厚度单元厚度A表面力作用面积表面力作用面积Tsysxsysxsqqqqq q

34、s qs 沿厚度均匀分布，沿厚度均匀分布，则单元表面力的势能则单元表面力的势能Vs为为 （3） 集中力势能集中力势能 当结构受到集中力时，通常在划分单元网格时就当结构受到集中力时，通常在划分单元网格时就把集中力的作用点设置为节点。于是单元集中力把集中力的作用点设置为节点。于是单元集中力 Pc 的势能的势能Vc为为CTCPVp p/2C （4）总势能总势能CSVVVVV把把(1-35）式中原括符内的部分用列阵）式中原括符内的部分用列阵 Fd 代替，代替， 综合以上诸式，单元外力的总势能综合以上诸式，单元外力的总势能V为为 lCSTAVTTphdlqNhdxdyqN(1-35） Fd 具有和具有和

35、相同的行、列数。则相同的行、列数。则dTFV(1-36） 由单元的应变能由单元的应变能U(1-34）和外力势能）和外力势能V(1-36），），可得单元的总势能可得单元的总势能 21dTTFkVU(1-37）将式将式(1-37）代入，）代入，0根据弹性力学最小势能原理：结构处于根据弹性力学最小势能原理：结构处于稳定平衡稳定平衡的的必要和充分条件是必要和充分条件是总势能有极小值总势能有极小值。3、单元平衡方程、单元平衡方程于是有，于是有，21kUT(1-34）dTFV(1-36）式式(1-38）是从能量原理导出的单元平衡方程。这个方）是从能量原理导出的单元平衡方程。这个方程表达了单元力与单元位移之

36、间的关系。其中，程表达了单元力与单元位移之间的关系。其中， Fd 和和单元节点力单元节点力 F 具有相同的意义。具有相同的意义。 dFk(1-38）即得单元平衡方程即得单元平衡方程 1.6 单元刚度矩阵单元刚度矩阵 平衡方程平衡方程(1-38）中的矩阵）中的矩阵k是单元力和单元是单元力和单元位移关系间的系数矩阵，代表了单元的刚度特性，位移关系间的系数矩阵，代表了单元的刚度特性，称为单元刚度矩阵。单元刚度矩阵的体积为称为单元刚度矩阵。单元刚度矩阵的体积为nj nj， nj 是单元位移总数。其一般计算公式为：是单元位移总数。其一般计算公式为： 1、一般计算公式、一般计算公式它与单元应变矩阵它与单元

37、应变矩阵B和弹性矩阵和弹性矩阵D有关。有关。 dvBDBkvT(1-33） 对于平面应力三角形单元，应变矩阵对于平面应力三角形单元，应变矩阵B是常数是常数矩阵，同时弹性矩阵矩阵，同时弹性矩阵D也是常数矩阵，于是式也是常数矩阵，于是式(1-33）可以化简为可以化简为 式中式中A表示三角形单元的面积。表示三角形单元的面积。h是单元厚度。是单元厚度。 2、平面问题三角形单元刚度矩阵、平面问题三角形单元刚度矩阵（1）平面应力三角形单元）平面应力三角形单元 hASBhABDBkTT(1-39） 将式将式(1-9）和）和(1-26）代入上式，）代入上式，即得平面应力三角形单元刚度矩阵。写成分块形即得平面应

38、力三角形单元刚度矩阵。写成分块形式，有式，有mmmjmijmjjjiimijiikkkkkkkkkk(1-40）21001112称对ED(1-9）),(0021mjibccbABiiiii(1-26）式式(1-40)中子矩阵中子矩阵krs为为22矩阵，有矩阵，有 (1-41）（2）平面应变三角形单元）平面应变三角形单元对于平面应变问题，须将上式中的对于平面应变问题，须将上式中的E换为换为 ， 换为换为 ，于是有，于是有21E1srsrsrsrsrsrsrsrsTrrsbbcccbbcbccbccbbAEhhABDBk21212121)1 ( 42),(mjisr，组合见式，组合见式(1-40）

39、其中，其中，bi(j,m)、ci(j,m)是形函数式是形函数式(1-16）中的系数）中的系数(式式2-17)。(1-42） 平面问题的单元刚度矩阵平面问题的单元刚度矩阵k不随单元（或坐标轴）不随单元（或坐标轴）的平行移动或作的平行移动或作n 角度（角度（n为整数）的转动而改变。为整数）的转动而改变。 由公式由公式(1-41）、）、(1-42）知，）知，krs矩阵和其中的矩阵和其中的br、cr 、 bs、cs （r、s=i、j、m）有关。）有关。 单元平移时，单元平移时， bi、ci不变。不变。srsrsrsrsrsrsrsrrsbbcccbbcbccbccbbAhEk)1 ( 221)1 (

40、2211)1 ( 2211)1 ( 221)21)(1 ( 4)1 (),(mjisr，组合见式，组合见式(1-40）(3) 三角形单元刚度矩阵与坐标系无关三角形单元刚度矩阵与坐标系无关ijmxyo 单元转动时，单元转动时， bi、ci不变。不变。 当单元旋转时，各节点的编号保持不变。如图当单元旋转时，各节点的编号保持不变。如图1-7所示，图所示，图a所示的单元旋转所示的单元旋转 时，到达图时，到达图b所示位置。所示位置。mjiyybmjixxc(1-17） ijmyjymijm图图1-7xyo(b)xyo(a)jim可以证明，这两种情形的可以证明，这两种情形的k是相同的。是相同的。 其实，推

41、演公式其实，推演公式(1-40）、）、(1-41）、）、(1-42）时并）时并没有规定坐标系的方位，当坐标系旋转任意角度时，没有规定坐标系的方位，当坐标系旋转任意角度时，也不影响刚度矩阵的结果。因此，也不影响刚度矩阵的结果。因此，平面问题的单元刚平面问题的单元刚度矩阵可以认为是结构坐标系中的单元刚度矩阵，没度矩阵可以认为是结构坐标系中的单元刚度矩阵，没有坐标变换问题。有坐标变换问题。jjjjjjjjjjjjnjinjinnjninnnjnjjjijjinijiiiinjinjiFFFFFkkkkkkkkkkkkkkkkkkkkkkkkk212121212122222211111211(1-38

42、） （1）单元刚度矩阵中每个元素有明确的物理意义）单元刚度矩阵中每个元素有明确的物理意义例如，例如，kij表示单元第表示单元第j个自由度产生单位位移个自由度产生单位位移( j=1)，其他自由度固定（其他自由度固定（=0）时，在第）时，在第i个自由度产生的节个自由度产生的节点力点力Fi。主对角线上元素主对角线上元素kii(i=1,nj)恒为正值。恒为正值。3、单元刚度矩阵性质、单元刚度矩阵性质（2）k的每一行或每一列元素之和为零的每一行或每一列元素之和为零F1 =0F2 =0F3=0Fi=0Fj =0Fnj =0rst11injinjjijiiiiiFkkkkk2211以上式中第以上式中第i行为

43、例，行为例，当所有节点沿当所有节点沿x向或向或y向向都产生单位位移时，都产生单位位移时，单元作平动运动，无应变，单元作平动运动，无应变，也无应力。也无应力。则有：则有：121nj021injijiiiikkkkk即：即：k的每一行元素之和为零。根据对称性，每一列元的每一行元素之和为零。根据对称性，每一列元素之和也为零。素之和也为零。rstxy图图1-6iinjijiiFkkkk21（3）k是对称矩阵是对称矩阵 由由k各元素各元素的表达式，可知的表达式，可知k具有对称性。具有对称性。 jjjjjjjjjjnnjninnnjnjjjijjinijiiiinjinjikkkkkkkkkkkkkkkk

44、kkkkkkkkkk21212122222211111211njnj对于主对角线元素对称。对称表达式：对于主对角线元素对称。对称表达式：kij = kji证明证明 kij表示当单元位移中第表示当单元位移中第j个元素为个元素为1( j=1)其余元其余元素为零时，引起的单元力中的第素为零时，引起的单元力中的第i个节点力个节点力Fi kji表示当单元位移中第表示当单元位移中第i个元素为个元素为1( i=1)其余元其余元素为零时，引起的单元力中的第素为零时，引起的单元力中的第j个节点力个节点力Fj第第 i自由度自由度 第第 j自由度自由度位移位移 i=1 j=1力力Fi = kijFj = kji虚功

45、虚功Fi i = kijFj j = kji由虚功原理，得由虚功原理，得 kij = kji（4）单元刚度矩阵是奇异矩阵）单元刚度矩阵是奇异矩阵 即即k的行列式为零（由行列式性质）的行列式为零（由行列式性质） 。 单元刚度矩阵是在单元处于平衡状态的前提下得单元刚度矩阵是在单元处于平衡状态的前提下得出的。单元作为分离体看待，作用在它上面的外力出的。单元作为分离体看待，作用在它上面的外力（单元力）必定是平衡力系。然而，研究单元平衡时（单元力）必定是平衡力系。然而，研究单元平衡时没有引入约束。没有引入约束。承受平衡力系作用的无约束单元，其承受平衡力系作用的无约束单元，其变形是确定的，但位移不是确定的

46、。变形是确定的，但位移不是确定的。所以出现性质所以出现性质（3）中的）中的“平动问题平动问题”，即单元可以发生任意的刚，即单元可以发生任意的刚体运动。从数学上讲，方程体运动。从数学上讲，方程(1-28）的解不是唯一的或）的解不是唯一的或不能确定的。由此，单元刚度矩阵一定是奇异的。不能确定的。由此，单元刚度矩阵一定是奇异的。（5）单元刚度矩阵是常量矩阵）单元刚度矩阵是常量矩阵单元力和单元位移成线性关系是基于弹性理论的结果。单元力和单元位移成线性关系是基于弹性理论的结果。4、例：平面应力直角三角形单元刚度矩阵、例：平面应力直角三角形单元刚度矩阵 图图1-8示出一平面应力直角三角形单元，直角边示出一

47、平面应力直角三角形单元，直角边长分别为长分别为a、b，厚度为，厚度为h，弹性模量为，弹性模量为E，泊松比为，泊松比为 ，计算单元刚度矩阵。计算单元刚度矩阵。图图1-8ijmabxy 第一步：计算第一步：计算bi、ci和单元和单元 面积面积A。 图图1-8mjiyybmjixxc（1-17） ijmabxyXi(j,m)Yi(j,m)bi(j,m)ci(j,m)ia0b0j0b0am00-b-a表表2-1 单元节点坐标和单元节点坐标和bi、ci值（值（i、j、m）参数参数节点节点单元面积单元面积: A=ab/2 计算步骤计算步骤 第二步：求子矩阵第二步：求子矩阵 由式由式(1-41），算得），算

48、得 2222100)1 (2bbabEhkii0210)1 (22abababEhkij其他从略。其他从略。 第三步：形成第三步：形成k将将kii等按式等按式(1-40）组集成）组集成k 。(1-43a） 222222222222222222121212121212121002121021210212102121000)1 ( 2baababaabbabababababaabbaabaababaaabbababbabbabbabEhk 2i-1 2i 2j-1 2j 2m-1 2m2i-12i2j-12j2m-12m红色号码红色号码是单元位移（是单元位移（ 1、 2、）在结构中对应的）在结构中

49、对应的节点位移的序号。节点位移的序号。ijmijmi、j、m表示单元中表示单元中3个节点在结构系统中的编号。个节点在结构系统中的编号。 当当a=b时，即等腰直角三角形单元，有时，即等腰直角三角形单元，有 (1-43b）2321121212123212111100212102121021210212101001)1 (22Ehk i j mijm1.7 等价节点力等价节点力 从前面单元分析可以看出：单元平衡所用到的从前面单元分析可以看出：单元平衡所用到的的量均要属于节点的量，如单元位移、单元力。载的量均要属于节点的量，如单元位移、单元力。载荷亦应如此，必须将体积力、表面力转化到节点上荷亦应如此，

50、必须将体积力、表面力转化到节点上去，成为等价节点力（载荷）。在第去，成为等价节点力（载荷）。在第2.5节中已经得节中已经得到了公式到了公式(1-35）和）和(1-36） 。这里，这里， Fd 就是体积力、表面力和集中力之和的总等就是体积力、表面力和集中力之和的总等价节点力。价节点力。CSVVVVVlCSTAVTTphdlqNhdxdyqN(1-35）dTFV(1-36）cSVdpFFF(1-44）把总等价节点力把总等价节点力 Fd 分解成体积力、表面力和集中分解成体积力、表面力和集中力的等价节点力之和，有力的等价节点力之和，有 FV 单元上体积力的等价节点力单元上体积力的等价节点力 FS 单元

51、上表面力的等价节点力单元上表面力的等价节点力 pC 单元上节点上的集中力单元上节点上的集中力注意到式注意到式(1-35），得体积力等价节点力计算公式：），得体积力等价节点力计算公式：表面力的等价节点力计算公式：表面力的等价节点力计算公式：AVTVhdxdyqNF(1-45）hdlqNFSTlS(1-46） 1、体积力的等价节点力、体积力的等价节点力 2、表面力的等价节点力、表面力的等价节点力 3、等价节点力计算举例、等价节点力计算举例（1）单元自重）单元自重 图图1-9所示平面应力三角形单元，单元厚度为所示平面应力三角形单元，单元厚度为h。单元单位体积自重为单元单位体积自重为 ，自重指向，自重

52、指向y轴的负方向。轴的负方向。 PvixPviyPvjxPvjyPvmxPvmyAVTVhdxdyqNF(1-45）0Vq 计算式计算式mjiNNNN (1-21）图图1-9xyijm-0000000dxdyNNNNNNhFAmmjjiiV注意到形函数的性质注意到形函数的性质4：3AdxdyNAi(1-23）得自重荷载的等价节点力得自重荷载的等价节点力 00INNNNiiii(1-22）（i,j,m）根据体积力和式根据体积力和式(1-45）、）、(1-21）、）、(1-22），得），得TTVhAAAAAAAhF101010310300330033003(1-47）上式表明：自重载荷的等价节点力

53、为单元重量的上式表明：自重载荷的等价节点力为单元重量的1/3。 （2）均布面力）均布面力ijm图图1-10 xyqs单元边界上作用了均匀的分布力，单元边界上作用了均匀的分布力，如图如图1-10所示，其集度为所示，其集度为 qs 。 SySxSqqq hdlqNFSTlS(1-46）mjiNNNN (1-21）根据式根据式(1-46）、）、(1-21）和）和(1-22） 计算式计算式SySxlTmmjjiiSqqdlNNNNNNhF000000注意到形函数性质注意到形函数性质4 ：(1-23）ijij idlN21得得TSySxSySxSqqqqhlF002(1-48）00INNNNiiii(1

54、-22）均匀分布力的等价节点力为均匀分布力的等价节点力为 式式(1-48）表明：在）表明：在ij边上受均布面力的平面问题边上受均布面力的平面问题三角形单元，其等价节点力等于将均布面力合力之半三角形单元，其等价节点力等于将均布面力合力之半简单地简化到简单地简化到i、j节点上，方向与分布力方向相同。节点上，方向与分布力方向相同。m节点上为零。节点上为零。TSySxSySxSqqqqhlF002(1-48）ijmxyqsxFs1Fs3ijmxyqsyFs2Fs4（3）线性分布面力）线性分布面力ijm图图1-11xys 表面力集度在表面力集度在i点为点为qsx qsyT，而在而在j点为点为0。设坐标轴

55、。设坐标轴s的原点取在的原点取在j点，沿点，沿ji为正向，为正向， 。 lssij , 0ij边上任一点的面力集度边上任一点的面力集度 qs lsqlsqqsysxssqsiqsijm图图1-12xysl在在ij边上有：边上有：lsNilsNj10mN将将 qs 和上式代入式和上式代入式(1-46），有），有由形函数的性质由形函数的性质3：0),(),(1),(yxNxxxxyxNxxxxyxNmijijijiiSySxTqqllllh000060063003TSySxSySxqqqqhl00313132322(1-49） dsqlsqlslslslslshFsysxTlS0000100100

56、式式(1-49）表明：）表明：ij边受线性分布面力：边受线性分布面力： i点为点为qsx, qsyT，j点为点为0时，其等价节点力可将总载荷的时，其等价节点力可将总载荷的2/3分配给分配给i点，点，1/3分分配给配给j点，点，m点为零得出。点为零得出。 xyijmqsiqs体积力和表面力向节点的移置体积力和表面力向节点的移置符合静力等效原理的前提条件符合静力等效原理的前提条件是：线性位移模式。是：线性位移模式。1.7 1.7 系统分析系统分析1.7.1 坐标系坐标系研究各离散单元集合成整体结构，集合整体结构的研究各离散单元集合成整体结构，集合整体结构的平衡和变形协调，建立整体结构平衡方程。平衡

57、和变形协调，建立整体结构平衡方程。 单元分析时采用的坐标系成为局部坐标或单元坐单元分析时采用的坐标系成为局部坐标或单元坐标（单元刚度矩阵的通用性）。而结构系统分析时，标（单元刚度矩阵的通用性）。而结构系统分析时，必须在统一的坐标系内进行（各力学量才能叠加），必须在统一的坐标系内进行（各力学量才能叠加），称为称为“结构坐标结构坐标”或或“整体坐标整体坐标”，如图，如图1-13所示。所示。单元坐标系下，单元位移、单元力、单元刚度矩阵单元坐标系下，单元位移、单元力、单元刚度矩阵表示为：表示为： eeekF 、整体坐标系下，单元位移、单元力、单元刚度矩阵整体坐标系下，单元位移、单元力、单元刚度矩阵表示

58、为：表示为： kF 、XYXYPP图（图（1-13）(a) 平面桁架平面桁架(杆件单元杆件单元) 悬臂深梁悬臂深梁 (平面三角形单平面三角形单元元) xyxyxyxy如何从单元坐标转化为结构坐标将在第如何从单元坐标转化为结构坐标将在第4章中讨论。章中讨论。1.7.2 整体刚度矩阵 假设整体结构被划分为假设整体结构被划分为ne个单元和个单元和n个节点，在个节点，在整体坐标系下，对于每个单元均有：整体坐标系下，对于每个单元均有：eeneXnnennnFk11111111:对于整体结构，有Fk 将上述这些方程集合起来（整体坐标下叠加），将上述这些方程集合起来（整体坐标下叠加），便可得到整个结构的平衡

59、方程。为此，需要将便可得到整个结构的平衡方程。为此，需要将k、F体积膨胀，分别扩大为体积膨胀，分别扩大为n1n1、n11和和n11的矩阵才能相加。膨胀后，原有节点号对应位的矩阵才能相加。膨胀后，原有节点号对应位置的元素不变，而其它元素均为零。置的元素不变，而其它元素均为零。PK 组装方法：组装方法：建立一个体积为建立一个体积为n1n1的方阵，按的方阵，按单元序号依次把结构坐标单元刚度矩阵的元素放入单元序号依次把结构坐标单元刚度矩阵的元素放入该方阵中。该方阵中。 放入方法：放入方法：（1）按单元节点编码）按单元节点编码对号入座对号入座； （2）同位置元素累加同位置元素累加。eenennnennF

60、PkK11111111, , 式中：式中：K为整体刚度矩阵，为整体刚度矩阵，为整体节点位移为整体节点位移列阵；列阵；P为整体等价节点荷载列阵。如下：为整体等价节点荷载列阵。如下：（1-50）000000000000000022mmmjmijmjjjiimijiinnkkkkkkkkkkijmijm例：平面三角单元例：平面三角单元双行双列双行双列1.7.3 结构刚度矩阵特性结构刚度矩阵特性1、结构刚度矩阵元素的力学意义、结构刚度矩阵元素的力学意义 把方程把方程(1-50）写开，）写开，=1=0=0=0=0=011111111111113211321321321321333333231222232