物理方程矩阵式课件

有限元思路框图有限元思路框图解综合方程解综合方程K=P求结构节点位移求结构节点位移计算结构内力和应力计算结构内力和应力系统分析系统分析(把单元刚度矩阵集合成结构刚度矩阵把单元刚度矩阵集合成结构刚度矩阵K形成等价节点荷载形成等价节点荷载P)离散（剖分）结构离散（剖分）结构为若干单元为若干单元单元分析单元分析(建立单元刚度矩阵建立单元刚度矩阵ke形成单元等价节点力形成单元等价节点力)有限元思路框图解综合方程K=P系统分析离散（1（1）剖分结构时应对单元、节点分别用连续正整数编号。剖分结构时应对单元、节点分别用连续正整数编号。123456789zxyuxuyuz（1）剖分结构时应对单元、节点分别用连续正整数编号。12342（2）从结构中取出单元，进行从结构中取出单元，进行单元分析单元分析52623杆件单元杆件单元板单元板单元（2）从结构中取出单元，进行单元分析52623杆件单元板3第二章第二章 单元分析单元分析 平面问题常应变单元平面问题常应变单元 在用矩阵描述单元各种力学量时，不同性质单元的在用矩阵描述单元各种力学量时，不同性质单元的同一力学量可采用相同的矩阵符号，不同的仅仅是矩阵同一力学量可采用相同的矩阵符号，不同的仅仅是矩阵体积和矩阵元素。体积和矩阵元素。本章主要讲单元分析的一般理论、方法。但为了便本章主要讲单元分析的一般理论、方法。但为了便于理解，以于理解，以平面问题常应变三角形单元平面问题常应变三角形单元为对象进行说明、为对象进行说明、演引。必须指出：尽管说明、演引中具有明显的针对性演引。必须指出：尽管说明、演引中具有明显的针对性（平面问题三角形单元），但原理、方法和主要矩阵公（平面问题三角形单元），但原理、方法和主要矩阵公式都具有普遍性。式都具有普遍性。第二章 单元分析 在用矩阵描述单元各种力学量4单单 元元 分分 析析 的的 内内 容容结点位移(1)单元内部各点位移单元应变单元应力(2)(3)结点力(4)位移协调模式几何方程物理方程平衡方程边界条件单元分析单元刚度矩阵单元刚度矩阵单 元 分 析 的 内 容结点位移(1)单元内部各点位移单元5（2-1）2、单元内任意点的、单元内任意点的体积力体积力列阵列阵 qV（2-2）1、单元表面或边界上任意点的、单元表面或边界上任意点的表面力表面力列阵列阵 qs ijmxyijmxyqVqs2.1 基本力学量矩阵基本力学量矩阵（2-1）2、单元内任意点的体积力列阵qV（2-2）1、6图图2-1ijmxyuv3、单元内任意点的位移列阵、单元内任意点的位移列阵 f（2-3）4、单元内任意点的应变列阵、单元内任意点的应变列阵 （2-4）ijmxy图2-1ijmxyuv3、单元内任意点的位移列阵f（275、单元内任意点的应力列阵、单元内任意点的应力列阵 （2-5）6、几何方程列阵、几何方程列阵（2-6）将上式代入式（将上式代入式（2-4）ijmxy（2-4）5、单元内任意点的应力列阵（2-5）6、几何方程列阵（87、物理方程矩阵式、物理方程矩阵式（2-7）式中式中 E、弹性模量、泊松比。弹性模量、泊松比。上式可简写为上式可简写为（2-8）对对于于弹弹性性力力学学的的平平面面应应力力问问题题，物物理理方方程程的的矩矩阵阵形形式式可可表示为：表示为：7、物理方程矩阵式（2-7）式中 E、弹性模量、泊9（2-9）矩矩阵阵D称称为为弹弹性性矩矩阵阵。式式（2-9）给给出出的的弹弹性性矩矩阵阵D的的矩矩阵元素是按照平面应力问题的物理方程得出的；阵元素是按照平面应力问题的物理方程得出的；对对于于平平面面应应变变问问题题，需需将将式式（2-9）中中的的 E 换换为为 ，换为换为 。（2-8）各种类型结构的弹性物理方程都可用式（各种类型结构的弹性物理方程都可用式（2-8）描述。）描述。但结构类型不同，力学性态但结构类型不同，力学性态(应力分量、应变分量应力分量、应变分量)有区有区别，别，弹性矩阵弹性矩阵D的体积和元素是不同的。的体积和元素是不同的。其中其中:（2-9）矩阵D称为弹性矩阵。式（2-9）给出的弹性矩阵10单单 元元 分分 析析 的的 内内 容容结点位移(1)单元内部各点位移单元应变单元应力(2)(3)结点力(4)位移协调模式几何方程物理方程平衡方程边界条件单元分析单元刚度矩阵单元刚度矩阵？单 元 分 析 的 内 容结点位移(1)单元内部各点位移单元112.2 2.2 位移函数和形函数位移函数和形函数1、位移函数概念、位移函数概念 “位移函数位移函数”也称也称“位移模式位移模式”，是单元内部位移变，是单元内部位移变化的数学表达式，是坐标的函数化的数学表达式，是坐标的函数。有限元法采用能量原理。有限元法采用能量原理进行单元分析，因而必须事先给出（设定）位移函数。一进行单元分析，因而必须事先给出（设定）位移函数。一般而论，般而论，位移函数选取会影响甚至严重影响计算结果的精位移函数选取会影响甚至严重影响计算结果的精度度。弹性力学中，恰当选取位移函数不是一件容易的事情。弹性力学中，恰当选取位移函数不是一件容易的事情。有限单元法中当单元划分得足够小时，把位移函数设定为有限单元法中当单元划分得足够小时，把位移函数设定为简单的多项式也可得到相当精确的结果简单的多项式也可得到相当精确的结果。这正是有限单元。这正是有限单元法具有的重要优势之一。法具有的重要优势之一。2.2 位移函数和形函数1、位移函数概念12 不同类型结构会有不同的位移函数。这里，仍以平面不同类型结构会有不同的位移函数。这里，仍以平面问题三角形单元（图问题三角形单元（图2-2）为例，说明设定位移函数的有）为例，说明设定位移函数的有关问题。关问题。图图2-2是是一一个个三三节节点点三三角角形形单单元元，其其节节点点i、j、m按按逆逆时时针针方方向向排排列列。每每个个节节点点位位移移在在单单元元平平面内有两个分量：面内有两个分量：（2-10）图图2-2ijmuiujumvivjvmxy2、位移函数设定举例、位移函数设定举例 不同类型结构会有不同的位移函数。这里，仍以平13 一个三角形单元有一个三角形单元有3个节点（以个节点（以 i、j、m为为 序），共有序），共有6个节点位移分量。其个节点位移分量。其单元位移或单元节点位移列阵单元位移或单元节点位移列阵为：为：（2-11）一个三角形单元有3个节点（以 i、j、m为 序）14本问题选位移函数为：本问题选位移函数为：（2-12）式中式中:a1、a2、a6待定常数，由单元位移的待定常数，由单元位移的6个分量确定。个分量确定。式（式（2-12）位移函数中，）位移函数中，a1、a4代表刚体位移，代表刚体位移，a2、a3、a5、a6 代表单元中有常应变，且位移函数是连续函数。代表单元中有常应变，且位移函数是连续函数。ijmuiujumvivjvmxyuv本问题选位移函数为：（2-12）式中:a1、a2、a6153、选取位移函数应考虑的问题、选取位移函数应考虑的问题 （1）单元）单元有几个位移函数有几个位移函数 单元中任意一点有几个位移分量就有几个位移函数。单元中任意一点有几个位移分量就有几个位移函数。本单元中有本单元中有u和和v，与此相应，有，与此相应，有2个位移函数；个位移函数；（3）位移函数中待定常数个数位移函数中待定常数个数 待定常数个数应等于单元位移列阵中的位移分量数。以便用单待定常数个数应等于单元位移列阵中的位移分量数。以便用单元位移确定位移函数中的待定常数。本单元位移列阵中有元位移确定位移函数中的待定常数。本单元位移列阵中有6个分量，个分量，为了能把为了能把2个位移函数（个位移函数（u、v）和单元位移的）和单元位移的6个分量联系起来，两个分量联系起来，两个位移函数中包含的待定常数一共应有个位移函数中包含的待定常数一共应有6个。个。（2）位移函数是坐标的函数位移函数是坐标的函数 本单元的坐标系为：本单元的坐标系为：X、Y；3、选取位移函数应考虑的问题 （1）单元有几个位移函数（16（4）位移函数中必须包含单元的刚体位移。位移函数中必须包含单元的刚体位移。（5）位移函数中必须包含单元的常应变。位移函数中必须包含单元的常应变。（6）位移函数在单元内要连续位移函数在单元内要连续;相邻单元间要尽量协调。相邻单元间要尽量协调。条件（条件（4）、（）、（5）构成单元的）构成单元的完备性准则完备性准则，条件，条件（6）是单元的）是单元的协调性条件协调性条件。理论和实践都已证明。理论和实践都已证明:完备完备性准则是有限元解收敛于真实解的必要条件，再加上位移性准则是有限元解收敛于真实解的必要条件，再加上位移协调条件协调条件(充分条件充分条件)才构成有限元解的充要条件。容易证才构成有限元解的充要条件。容易证明，三角形三节点常应变单元满足以上必要与充分条件。明，三角形三节点常应变单元满足以上必要与充分条件。（4）位移函数中必须包含单元的刚体位移。17例：平面应力矩形板被划分为若干三角形单元。例：平面应力矩形板被划分为若干三角形单元。位移函数中包含了单元的常应变。位移函数中包含了单元的常应变。(a2,a6,a3+a5)位移函数中包含了单元的刚体位移。位移函数中包含了单元的刚体位移。(a1,a4)254136对任一单元，如对任一单元，如单元，取位移函数：单元，取位移函数：例：平面应力矩形板被划分为若干三角形单元。位移函数中包含了单18、单元的位移函数都是单元的位移函数都是可以看出可以看出：位移函数在单元内是连续的；位移函数在单元内是连续的；位移函数在单元之间的边界上也连续吗？是。位移函数在单元之间的边界上也连续吗？是。以以、的边界的边界26为例为例:2562635623xyuu6u2uu6u2两条直线上有两个点重合，此两条直线必全重合。两条直线上有两个点重合，此两条直线必全重合。、单元的位移函数都是可以看出：以、的边界26194、形函数、形函数形函数是用单元位移分量来描述位移函数的插值函数形函数是用单元位移分量来描述位移函数的插值函数。（1）形函数定义）形函数定义 现在，通过单元位移确定位移函数中的待定常数现在，通过单元位移确定位移函数中的待定常数a1、a2、a6。设节点。设节点i、j、m的坐标分别为（的坐标分别为（xi、yi）、（）、（xj、yj）、（）、（xm、ym），节点位移分别为（），节点位移分别为（ui、vi）、）、（uj、vj）、（um、vm）。）。将它们代入式（将它们代入式（2-12），有），有:4、形函数形函数是用单元位移分量来描述位移函数的插值函数。（20（2-13）从式（从式（2-13）左边）左边3个方程中解出待定系数个方程中解出待定系数a1、a2、a3为为:（2-13）从式（2-13）左边3个方程中解出待定系数a121（2-14）（2-14）22式中，式中，A为三角形单元的面积，有为三角形单元的面积，有:（2-15）特别指出：特别指出：为使求得面积的值为正值，本单元节点号的次序必须为使求得面积的值为正值，本单元节点号的次序必须是逆时针转向，如图所示。至于将哪个结点作为起始结点是逆时针转向，如图所示。至于将哪个结点作为起始结点i，则没有关系。，则没有关系。ijmxy（2）（1）（3）式中，A为三角形单元的面积，有:（2-15）特别指出：i23将式（将式（2-14）代入式（）代入式（2-12）的第一式，整理后得）的第一式，整理后得同理同理:（2-16）式中式中:（2-17）ijm将式（2-14）代入式（2-12）的第一式，整理后得同理:（24 式（式（2-17）中（）中（i、j、m）意指：按）意指：按i、j、m依次轮换依次轮换下标，可得到下标，可得到aj、bj、cjam、bm、cm。后面出现类似情况。后面出现类似情况时，照此推理。时，照此推理。式（式（2-17）表明：）表明：aj、bj、cjam、bm、cm是单元三是单元三个节点坐标的函数。个节点坐标的函数。式（2-17）中（i、j、m）意指：按i、j25（2-16）令令（2-18）位移模式（位移模式（2-16）可以简写为）可以简写为（2-19）式（式（2-19）中的）中的Ni、Nj、Nm是坐标的函数，反应了单元的是坐标的函数，反应了单元的位移形态，称为单元位移函数的形函数。数学上它反应了位移形态，称为单元位移函数的形函数。数学上它反应了节点位移对单元内任一点位移的插值，又称插值函数。节点位移对单元内任一点位移的插值，又称插值函数。（2-16）令（2-18）位移模式（2-16）可以简写为26形形 函函 数数i,j,m在单元任一点上三个形函数之和等于在单元任一点上三个形函数之和等于1 1第一列与它的代数余子式之和第一列与第二列的代数余子式之和第一列与第三列的代数余子式之和1.三个形函数只有两个是独立的2 当三角形单元的三个结点的位移相等当三角形单元的三个结点的位移相等形 函 数i,j,m在单元任一点上三个形函数之和等于1第27（2-16）用形函数把式（用形函数把式（2-16）写成矩阵，有）写成矩阵，有缩写为缩写为（2-20）形函数是用单元位移分量来描述位移函数的插值函数形函数是用单元位移分量来描述位移函数的插值函数。（2-16）用形函数把式（2-16）写成矩阵，有缩写为（2-28N为形函数矩阵，进一步写成分块形式：为形函数矩阵，进一步写成分块形式：（2-21）其中子矩阵其中子矩阵（2-22）I是是22的单位矩阵。的单位矩阵。下面将会看到，形函数是有限单元法中的一个重要函下面将会看到，形函数是有限单元法中的一个重要函数。了解它的一些基本性质是有益的。数。了解它的一些基本性质是有益的。（2）形函数性质）形函数性质N为形函数矩阵，进一步写成分块形式：（2-21）其中子矩29性质性质1 形函数形函数Ni在节点在节点i上的值等于上的值等于1，在其它节点上的值，在其它节点上的值等于等于0。对于本单元，有。对于本单元，有:（i、j、m）性质1 形函数Ni在节点i上的值等于1，在其它节点上的30 xyN（i，j，m）Ni=1ijm图图2-3xyNNi=1ijm图2-331性质性质2 在单元中任一点，所有形函数之和等于在单元中任一点，所有形函数之和等于1。对于本。对于本单元，有单元，有为什么？为什么？xyN（I，j，m）Ni=1ijmNj=1ijmNm=1ijmNi=1ijmNj=1Nm=1图图2-4性质2 在单元中任一点，所有形函数之和等于1。对于本单32xyN（I，j，m）Ni=1ijmNj=1ijmNm=1ijmNi=1ijmNj=1Nm=1图图2-4上图说明：形函数是线性的（1）单元位移场是线性的；（2）单元位移场与结点位移是协调的；（3）结点位移将影响位移场的数值大小xyNNi=1ijmNj=1ijmNm=1ijmNi 33性质性质3 在三角形单元的边界在三角形单元的边界ij上任一点（上任一点（x，y），有），有:xxixjxyNi(xi，yi)j(xj，yj)m(xm，ym)Ni(x、y)1证图图2-5性质3 在三角形单元的边界ij上任一点（x，y），有:34性质性质4 形函数在单元上的面积分和边界上的线积分公式为形函数在单元上的面积分和边界上的线积分公式为（2-23）式中式中 为为 边的长度。边的长度。相邻单元的位移在公共边上是连续的相邻单元的位移在公共边上是连续的ijpmxxixjxyNi(xi，yi)j(xj，yj)m(xm，ym)Ni(x、y)1Ni=1ijm性质4 形函数在单元上的面积分和边界上的线积分公式为 352.3 单元应变矩阵和应力矩阵单元应变矩阵和应力矩阵 式（2-6）给出了单元内任一点的应变和位移之间关系。（2-6）1、单元应变矩阵、单元应变矩阵2.3 单元应变矩阵和应力矩阵 式（2-6）给出了单元内36对位移函数（式（对位移函数（式（2-16）（2-24）（2-16）求导后代入式（求导后代入式（2-6），得到应变和节点位移的关系式。），得到应变和节点位移的关系式。对位移函数（式（2-16）（2-24）（2-16）求导后代37（2-25）式中式中,B单元应变矩阵单元应变矩阵。对本问题，维数为对本问题，维数为36。它的分块形式为：。它的分块形式为：子矩阵子矩阵:（2-26）由由于于 与与x、y无无关关，都都是是常常量量，因因此此B矩矩阵阵也也是是常常量量。单单元元中中任任一一点点的的应应变变分分量量是是B矩矩阵阵与与单单元元节节点点位位移移的的乘乘积积，因因而而也也都都是是常常量量。因因此此，这这种种单单元元被称为常应变单元被称为常应变单元。（2-25）式中,B单元应变矩阵。对382、单元应力矩阵、单元应力矩阵将式（将式（2-25）代入物理方程式（）代入物理方程式（2-8），得），得（2-8）（2-27）上式也可写为：上式也可写为：（2-28）这这是是单单元元内内任任一一点点应应力力与与单单元元位位移移的的关关系系式式。其其中中S称称为单元应力矩阵为单元应力矩阵，并有：，并有：（2-29）D是是33 弹性矩阵弹性矩阵，B是是36应变矩阵应变矩阵，因此，因此S也是也是36 矩阵。它可写为分块形式矩阵。它可写为分块形式 2、单元应力矩阵将式（2-25）代入物理方程式（2-8），得39（2-30）将将弹弹性性矩矩阵阵（式式（2-9）和和应应变变矩矩阵阵（式式（2-26）代入，得子矩阵代入，得子矩阵Si由式（由式（2-29）得：）得：（2-31）式式（2-31）是是平平面面应应力力问问题题的的结结果果。对对于于平平面面应应变变问问题题，只要将上式中的只要将上式中的E换成换成 ，换成换成 即得。即得。（2-30）将弹性矩阵（式（2-9）和应变矩40（2-32）由由于于同同一一单单元元中中的的D、B矩矩阵阵都都是是常常数数矩矩阵阵，所所以以S矩矩阵阵也也是是常常数数矩矩阵阵。也也就就是是说说，三三角角形形三三节节点点单单元元内内的应力分量也是常量的应力分量也是常量。当当然然，相相邻邻单单元元的的E，A和和bi、ci(i，j，m)一一般般不不完完全全相相同同，因因而而具具有有不不同同的的应应力力，这这就就造造成成在在相相邻邻单单元元的的公公共共边边上上存存在在着着应应力力突突变变现现象象。但但是是随随着着网网格格的的细细分分，这这种种突变将会迅速减小，平衡被满足。突变将会迅速减小，平衡被满足。（2-32）由于同一单元中的D、B矩阵都是41几何关系位移函数几何关系?平衡关系单元刚度矩阵单元刚度矩阵几何关系位移函数几何关系?平衡关系单元刚度矩阵422.4 2.4 单元应变能和外力势能的矩阵表达单元应变能和外力势能的矩阵表达1、单元应变能单元应变能 仍以平面应力问题中的三角形单元说明，设单元厚度为仍以平面应力问题中的三角形单元说明，设单元厚度为h 将将式式（2-25）和和（2-27）代代入入上上式式进进行行矩矩阵阵运运算算，并并注注意意到弹性矩阵到弹性矩阵D的对称性，有的对称性，有应变能应变能 U为：为：ijmxyh2.4 单元应变能和外力势能的矩阵表达1、单元应变能 仍43由于由于和和T是常量，提到积分号外，上式可写成是常量，提到积分号外，上式可写成 引入矩阵符号引入矩阵符号kk，且有：，且有：（2-33a）式（式（2-33a）是针对平面问题三角形单元推出的。注意到其）是针对平面问题三角形单元推出的。注意到其中中hdxdy的实质是任意的微体积的实质是任意的微体积dv，于是得，于是得 k的一般式。的一般式。（2-33）式（式（2-332-33）不仅适合于平面问题三角形单元，也是计）不仅适合于平面问题三角形单元，也是计算各种类型单元算各种类型单元kk的一般式的一般式。由于和T是常量，提到积分号外，上式可写成 引入矩442.6节中将明确节中将明确k的力学意义是单元刚度矩阵。式（的力学意义是单元刚度矩阵。式（2-33）便是计算单元刚度矩阵的基本矩阵式。它适合于各种类型便是计算单元刚度矩阵的基本矩阵式。它适合于各种类型的单元。的单元。单元应变能单元应变能写成写成（2-34）2、单元外力势能单元外力势能 单单元元受受到到的的外外力力一一般般包包括括体体积积力力、表表面面力力和和集集中中力力。自自重重属属于于体体积积力力范范畴畴。表表面面力力指指作作用用在在单单元元表表面面的的分分布布载载荷，如风力、压力，以及相邻单元互相作用的内力等。荷，如风力、压力，以及相邻单元互相作用的内力等。（2-33）2.6节中将明确k的力学意义是单元刚度矩阵。式（2-3345（1）体积力势能体积力势能单位体积中的体积力单位体积中的体积力如式（如式（2-2）所示。）所示。单元上体积力具有的势能单元上体积力具有的势能Vv为为（2-2）ijmxyqVxqVyijmxyuv（1）体积力势能单位体积中的体积力如式（2-2）所示。单46注意到式（注意到式（2-20）有有（2-20）（2）表面力势能表面力势能 面面积积力力虽虽然然包包括括单单元元之之间间公公共共边边上上互互相相作作用用的的分分布布力力，但但它它们们属属于于结结构构内内力力，成成对对出出现现，集集合合时时互互相相抵抵消消，在在结结构构整整体体分分析析时时可可以以不不加加考考虑虑，因因此此单单元元分分析析时时也也就就不不予考虑。予考虑。注意到式（2-20）有（2-20）（2）表面力势能 47 现在，只考虑弹性体边界上的表面力，它只在部分单现在，只考虑弹性体边界上的表面力，它只在部分单元上形成表面力（右下图）。设元上形成表面力（右下图）。设边界单位长度上受到的表边界单位长度上受到的表面力面力如式（如式（2-1）。）。l单元边界长度单元边界长度h单元厚度单元厚度A表面力作用面积表面力作用面积（2-1）qs 则单元表面力的势能则单元表面力的势能Vs为为 现在，只考虑弹性体边界上的表面力，它只在部分48（3）集中力势能集中力势能 当结构受到集中力时，通常在划分单元网格时就把集中当结构受到集中力时，通常在划分单元网格时就把集中力的作用点设置为节点。于是单元集中力力的作用点设置为节点。于是单元集中力 Pc 的势能的势能Vc为为（4）总势能总势能 如果把（如果把（2-35）式中原括号内的部分用列阵）式中原括号内的部分用列阵 Fd 代替，代替，综合以上诸式，单元外力的总势能综合以上诸式，单元外力的总势能V为为（2-35）Fd 具有和具有和相同的行、列数。则：相同的行、列数。则：（3）集中力势能 当结构受到集中力时，通常在划49 由由单单元元的的应应变变能能U（2-34）和和外外力力势势能能V（2-36）,可可得得单元的总势能单元的总势能（2-37）以节点位移为未知量，对总势能取极值问题变成了一以节点位移为未知量，对总势能取极值问题变成了一个多元函数的极值问题。有极值条件个多元函数的极值问题。有极值条件 2.5 能量原理和单元平衡方程能量原理和单元平衡方程（2-36）由单元的应变能U（2-34）和外力势能V（2-36）50式（式（2-38）是从能量原理导出的单元平衡方程。这）是从能量原理导出的单元平衡方程。这个方程表达了单元力与单元位移之间的关系。其中，个方程表达了单元力与单元位移之间的关系。其中，Fd 和单元节点力和单元节点力 F 具有相同的意义。具有相同的意义。（2-38）于是，将式（于是，将式（2-37）代入，即得单元平衡方程）代入，即得单元平衡方程:根据弹性力学能量原理：根据弹性力学能量原理：结构处于稳定平衡的必要和充分结构处于稳定平衡的必要和充分条件是总势能有极小值条件是总势能有极小值。式（2-38）是从能量原理导出的单元平衡方程。这个方程表达了512.6 单元刚度矩阵单元刚度矩阵 平衡方程（平衡方程（2-38）中的矩阵）中的矩阵k是单元力和单元位移关是单元力和单元位移关系间的系数矩阵，代表了单元的刚度特性，称为单元刚度系间的系数矩阵，代表了单元的刚度特性，称为单元刚度矩阵。单元刚度矩阵的体积为矩阵。单元刚度矩阵的体积为nj nj，nj 是单元位移总数。是单元位移总数。1、计算单元刚度矩阵的一般公式、计算单元刚度矩阵的一般公式 计算各类单元的单元刚度矩阵可用式（计算各类单元的单元刚度矩阵可用式（2-33）执行。）执行。它与单元应变矩阵它与单元应变矩阵B和弹性矩阵和弹性矩阵D有关。有关。（2-33）2.6 单元刚度矩阵 平衡方程（2-38）中52 对于平面应力三角形单元，应变矩阵对于平面应力三角形单元，应变矩阵B是常数矩阵，是常数矩阵，同时弹性矩阵同时弹性矩阵D也是常数矩阵，于是式（也是常数矩阵，于是式（2-33）可以化）可以化简为简为 式中式中A表示三角形单元的面积。表示三角形单元的面积。2 2、平面问题三角形单元刚度矩阵、平面问题三角形单元刚度矩阵（1）平面应力三角形单元）平面应力三角形单元（2-39）对于平面应力三角形单元，应变矩阵B是常数53 将式（将式（2-9）和（）和（2-26）代入上式，）代入上式，即得平面应力三角形单元刚度矩阵。写成分块形式，有即得平面应力三角形单元刚度矩阵。写成分块形式，有（2-40）（2-9）（2-26）将式（2-9）和（2-26）代入上式，即得平面应54式中子矩阵为式中子矩阵为22矩阵，有矩阵，有（2-41）（2 2）平面应变三角形单元）平面应变三角形单元 对于平面应变问题，须将上式中的对于平面应变问题，须将上式中的E换为换为 ，换为换为 ，于是有：，于是有：，组合见式（，组合见式（2-40）其中，其中，bi(j,m)、ci(j,m)是形函数式（是形函数式（2-16）中的系数。）中的系数。式中子矩阵为22矩阵，有（2-41）（2）平面应变三角形55（2-42），组合见式（，组合见式（2-40）（3）示例）示例 平面应力直角三角形单元刚度矩阵平面应力直角三角形单元刚度矩阵 图图2-8示出一平面应力直角三角形单元，直角边长分别示出一平面应力直角三角形单元，直角边长分别为为a、b，厚度为，厚度为h，弹性模量为，弹性模量为E，泊松比为，泊松比为，计算单元，计算单元刚度矩阵。刚度矩阵。（2-42），组合见式（2-40）（3）示例 56第一步：计算第一步：计算bi、ci和单元面积和单元面积A 图图2-8（1-17）ijmabxyxijmyijmbijmcijmia0b0j0b0am00-b-a表表2-1 单元节点坐标和单元节点坐标和bi、ci值（值（i、j、m）参数参数节点节点单元面积单元面积:A=ab/2 计算步骤计算步骤第一步：计算bi、ci和单元面积A 图2-8（1-17）i57第二步：求子矩阵第二步：求子矩阵(由式（由式（2-41），算得），算得)其他从略。其他从略。第二步：求子矩阵(由式（2-41），算得)其他从略。58第三步：形成第三步：形成k 将将kii等按式（等按式（2-40）组集成）组集成k。（2-40）第三步：形成k（2-40）59（2-43a）1 2 3 4 5 6123456式中上边和右边的式中上边和右边的i、j、m表示子矩阵对应的节点号。表示子矩阵对应的节点号。（2-43a）1 2 60当当a=b时，即等腰直角三角形单元，有时，即等腰直角三角形单元，有（2-43b）1 2 3 4 5 6123456当a=b时，即等腰直角三角形单元，有（2-43b）1 61 子程序框图子程序框图SUBROUTINE SME3(NJ,NE,NI,ND2,E,AMU,H,LO,X,Y,SM,NTYPE)IMPLICIT DOUBLE PRECISION(a-h,o-z)DIMENSION LO(NE,3),X(NJ),Y(NJ),SM(ND2,ND2),B(3),C(3),D(3,3),YB(3,6),S(3,6),组集弹性刚度矩阵组集弹性刚度矩阵D（书式（书式1-19）计算应变矩阵计算应变矩阵B（书式（书式1-17、24、25、26）计算应力矩阵计算应力矩阵S（书式（书式1-29）计算单元刚度矩阵计算单元刚度矩阵k（书式（书式1-42）RETURNEND 子程序框图SUBROUTINE SME3(NJ,NE62 子程序子程序见附页见附页 SUBROUTING SME3()子程序见附页 SUBROUTING SME3(633、单元刚度矩阵性质、单元刚度矩阵性质njnj单元刚度矩阵单元刚度矩阵k的详细内容为的详细内容为（i、j是行列号是行列号）：）：3、单元刚度矩阵性质njnj单元刚度矩阵k的详细内容为64（2-38）单元刚度矩阵具有以下的性质：单元刚度矩阵具有以下的性质：（1）单元刚度矩阵中每个元素有明确的物理意义）单元刚度矩阵中每个元素有明确的物理意义例如，例如，kij表示当单元位移中第表示当单元位移中第j个元素为个元素为1(j=1)其余元素为零时，其余元素为零时，引起的单元力中的第引起的单元力中的第i个节点力个节点力Fi。把平衡方程（把平衡方程（2-38）写开）写开主对角线上元素主对角线上元素kii(i=1,nj)恒为正值。恒为正值。（2-38）单元刚度矩阵具有以下的性质：（1）单65（2）k的每一行或每一列元素之和为零的每一行或每一列元素之和为零F1=0F2=0F3=0Fi=0Fj=0Fnj=0rst图图2-6xyrst11以上式中第以上式中第i行为例，行为例，当所有节点沿当所有节点沿x向或向或y向都产生单位位移时，单元作平动运向都产生单位位移时，单元作平动运动，无应变，也无应力，因而单元结点力为零（不含初应动，无应变，也无应力，因而单元结点力为零（不含初应力）。力）。所以有所以有即，即，k的每一行元素之的每一行元素之和为零。由于对称性，和为零。由于对称性，每一列元素之和也为零。每一列元素之和也为零。（2）k的每一行或每一列元素之和为零F1=0F2=066（3）k是对称矩阵是对称矩阵 由由k单元的表达式，可见，由此可知单元的表达式，可见，由此可知k具有对称性。具有对称性。njnj对于主对角线元素对称。对称表达式：对于主对角线元素对称。对称表达式：kij =kji（3）k是对称矩阵njnj对于主对角线元素对称。对称表67证明：证明：kij表示当单元位移中第表示当单元位移中第j个元素为个元素为1(j=1)其余元素为零其余元素为零时，引起的单元力中的第时，引起的单元力中的第i个节点力个节点力Fi kji表示当单元位移中第表示当单元位移中第i个元素为个元素为1(i=1)其余元素为零其余元素为零时，引起的单元力中的第时，引起的单元力中的第j个节点力个节点力Fj第第 i自由度自由度 第第 j自由度自由度位移位移 i=1 j=1力力Fi=kijFj=kji虚功虚功Fi i=kijFj j=kji由虚功原理，得由虚功原理，得 kij=kji证明：kij表示当单元位移中第j个元素为1(j=1)68（4）单元刚度矩阵是奇异矩阵）单元刚度矩阵是奇异矩阵（即即k的行列式为零）的行列式为零）单元刚度矩阵是在单元处于平衡状态的前提下得出的。单元刚度矩阵是在单元处于平衡状态的前提下得出的。单元作为分离体看待，作用在它上面的外力（单元力）必单元作为分离体看待，作用在它上面的外力（单元力）必定是平衡力系。然而，定是平衡力系。然而，研究单元平衡时没有引入约束研究单元平衡时没有引入约束。承。承受平衡力系作用的无约束单元，其受平衡力系作用的无约束单元，其变形是确定的，但位移变形是确定的，但位移不是确定的不是确定的。所以出现性质（。所以出现性质（3）中的）中的“平动问题平动问题”，即，即单元可以发生任意的刚体运动。从数学上讲，方程（单元可以发生任意的刚体运动。从数学上讲，方程（2-38）的解不是唯一的或不能确定的。由此，）的解不是唯一的或不能确定的。由此，单元刚度矩阵单元刚度矩阵一定是奇异的一定是奇异的。（5）单元刚度矩阵是常量矩阵）单元刚度矩阵是常量矩阵单元力和单元位移成线性关系是基于弹性理论的结果。单元力和单元位移成线性关系是基于弹性理论的结果。（4）单元刚度矩阵是奇异矩阵（即k的行列式为零）（5）692.7 2.7 等价节点力等价节点力从前面单元分析可以看出：单元平衡所用到的量均属于节从前面单元分析可以看出：单元平衡所用到的量均属于节点的量，如单元位移、单元力。载荷亦应如此，点的量，如单元位移、单元力。载荷亦应如此，必须将体必须将体积力、表面力转化到节点上去，成为等价节点力积力、表面力转化到节点上去，成为等价节点力（载荷）。（载荷）。在第在第2.5节中已经得到了公式（节中已经得到了公式（2-35）和（）和（2-36）。）。（2-35）（2-36）这里，这里，Fd 就是体积力、表面力和集中力之和的总等价节就是体积力、表面力和集中力之和的总等价节点力。点力。2.7 等价节点力从前面单元分析可以看出：单元平衡所用到的量70（2-44）把总等价节点力把总等价节点力 Fd 分解成体积力、表面力和集中力的分解成体积力、表面力和集中力的等价节点力之和，有等价节点力之和，有 FV 单元上体积力的等价节点力单元上体积力的等价节点力 FS 单元上表面力的等价节点力单元上表面力的等价节点力 pC 单元上节点上的集中力单元上节点上的集中力注意到式（注意到式（2-35），得体积力等价节点力计算公式：），得体积力等价节点力计算公式：表面力的等价节点力计算公式：表面力的等价节点力计算公式：（2-45）（2-46）1、体积力的等价节点力、体积力的等价节点力2、表面力的等价节点力、表面力的等价节点力（2-44）把总等价节点力 Fd 分解成体积力、表面力和713、等价节点力计算举例、等价节点力计算举例（1）单元自重）单元自重 图图2-9所示平面应力三角形单元，单元厚度为所示平面应力三角形单元，单元厚度为h。单元。单元单位体积自重为单位体积自重为，自重指向，自重指向y轴的负方向。轴的负方向。图图2-9ijmPvixPvjxPvmxPviyPvjyPvmyxy（2-45）计算式计算式（2-21）3、等价节点力计算举例（1）单元自重 图2-972注意到形函数的性质注意到形函数的性质4：（2-23）得自重荷载的等价节点力得自重荷载的等价节点力（2-22）（i,j,m）根据体积力和式（根据体积力和式（2-45）、（）、（2-21）、（）、（2-22），得），得注意到形函数的性质4：（2-23）得自重荷载的等价节点力（73（2-47）上式表明：自重载荷的等价节点力为单元重量的上式表明：自重载荷的等价节点力为单元重量的1/3。子程序子程序见见 SUBJPE.FOR（2-47）上式表明：自重载荷的等价节点力为单元重量的1/374（2）均布面力）均布面力ijm图图2-10 xyqs单元边界上作用了均匀的分布力，单元边界上作用了均匀的分布力，如图如图2-10所示，其集度为所示，其集度为 qs。（2-46）（2-21）根据式（根据式（2-46）、（）、（2-21）和（）和（2-22）计算式计算式（2）均布面力ijm图2-10 xyqs单元边界上作用了均匀的75注意到形函数性质注意到形函数性质4：（2-23）得得（2-48）（2-22）均匀分布力的等价节点力为均匀分布力的等价节点力为 注意到形函数性质4：（2-23）得（2-48）（2-22）76 式（式（2-48）表明：在）表明：在ij边上受均布面力的平面问题边上受均布面力的平面问题三角形单元，其等价节点力等于将均布面力合力之半三角形单元，其等价节点力等于将均布面力合力之半简单地简化到简单地简化到i、j节点上，方向与分布力方向相同。节点上，方向与分布力方向相同。m节点上为零。节点上为零。子程序子程序见见 SUBJPE.FOR 式（2-48）表明：在ij边上受均布面力的平77（3）线性分布面力）线性分布面力ijm图图2-11xysqs 表面力集度在表面力集度在i点为点为qsx,qsyT，而在而在j点为点为0。设坐标轴。设坐标轴s的原点取在的原点取在i点，沿点，沿ij为正向，为正向，。ij边上任一点的面力集度边上任一点的面力集度 qs 计算式计算式（3）线性分布面力ijm图2-11xysqs 78ijm图图2-12xysl在在ij边上有：边上有：将将 qs 和上式代入式（和上式代入式（2-46），有），有由形函数的性质由形函数的性质3：ijm图2-12xysl在ij边上有：将qs和上式代入式79（2-49）（2-49）80 式（式（2-49）表明：）表明：ij边受线性分布面力：边受线性分布面力：i点为点为qsx,qsyT，j点为点为0时，其等价节点力可将总载荷的时，其等价节点力可将总载荷的2/3分配给分配给i点，点，1/3分配给分配给j点，点，m点为零得出。点为零得出。子程序子程序见见 SUBJPE.FOR 式（2-49）表明：ij边受线性分布面力：i点为qs81 单单元元上上的的体体积积力力和和表表面面力力向向结结点点的的移移置置都都是是符符合合直直观观的的静静力力等等效效原原理理的的，并并与与工工程程中中的的简简单单的的处处理理方方法法相相一一致致。应应当当指指出出，这这种种移移置置方方法法是是线线性性位位移移模模式式三三结结点点三三角角形形单单元元的的必必然然结结果果。对对于于非非线线性性位位移移模模式式的的单单元元，上上述述这这种种简简单单的的载载荷荷移移置置方方法法一一般般是是不不成成立立的的，而而应应当当采采用用公公式式（2-35）进行计算。）进行计算。（2-35）单元上的体积力和表面力向结点的移置都是符合直82本章小结本章小结：单元分析的主要任务是：单元分析的主要任务是：一、组集单元刚度矩阵；一、组集单元刚度矩阵；二、组集单元等价节点力矩阵。二、组集单元等价节点力矩阵。本章小结：单元分析的主要任务是：83作业：作业：2.以题以题1为依据，证明单元刚度矩阵的任何一行为依据，证明单元刚度矩阵的任何一行 （或列）元素的总和为零。（或列）元素的总和为零。1.证明常应变三角形单元发生刚体位移时，单元证明常应变三角形单元发生刚体位移时，单元 2.中将不产生应力。中将不产生应力。作业：2.以题1为依据，证明单元刚度矩阵的任何一行证明常应84