第3章 平面结构问题的有限单元法
  弹性力学研究的问题可以分为平面问题和空间问题两大类。而平面问题又可以分为平面应力问题和平面应变问题。用弹性力学经典解法解决实际问题的主要困难在于求解偏微分方程的复杂性，所求解满足弹性力学基本方程的位移函数、应变函数、应力函数的表达式要覆盖整个区域，还要满足边界条件，因此求解这样的函数形式是十分困难的。而有限元法则将原来连续的弹性体离散化，其中最简单的就是采用三角形单元对弹性体进行划分。相应对每个三角形单元选择最简单的线性函数为位移模式，单元中任一点的位移可以通过3个节点的位移进行插值运算，这样整个区域中无限多个未知位移量就可以用有限个节点来表示，从而避免了求解覆盖整个区域的位移函数的困难。平面问题的有限元法，不仅可用来解决实际问题，而且通过其相对简单的概念，可以详细了解用有限元法对一般弹性体进行应力分析的基本原理和方法步骤，了解有限元法的性能特点、使用中应注意的问题，从而为学习后续各章节打下基础。从教学观点看，通过学习平面问题有限元的基本概念和方法，就能容易地理解有限元法中所通用的理论框架与核心问题。学习并掌握了有限元基本理论及分析方法后，就可将其思想推广到各种结构形式，并应用于工程结构设计中。
对一个形状复杂的弹性体，要想用某种函数来描述结构整体内部任一点的位移是不大可能的。但是如果把弹性体离散为有限个单元体，就可以利用其节点的位移来构造出单元的位移插值函数，即位移函数。平面问题中用来离散的单元类型有多种，最常用的单元是三角形单元和矩形单元。图3.1所示为三角形三节点单元和四边形四节点单元。图3.2所示两种平面问题分别采用了矩形单元和三角形单元划分网格。本章仅就平面三角形常应变单元展开详细介绍，至于其他类型单元在第4章统一介绍。
图3.1 三角形与四边形单元
(a) 三角形单元; (b) 四边形单元
图3.2 三角形与四边形单元划分网格
(a) 四边形单元; (b) 三角形单元
3.1 平面三角形常应变单元位移模式
图3.3为任一三角形单元。设其顶点为节点，节点编号按右手法则依次为i、j、m。每个节点在平面内存在沿x轴和y轴的位移u、v，单图3.3 三角形单元元的3个节点共有6项位移分量。单元内任一点(x, y)的位移u(x,y)、v(x,y)可假设为坐标x、y的某种函数，即可选用适当的位移函数表示，这种函数也称为位移模式。
一种简单的线性位移函数为汽车结构有限元分析第3章 平面结构问题的有限单元法u=α1+α2x+α3y
v=α4+α5x+α6y(3.1)式中，α1、α2, …、α6为6个待定常数，可以由单元的节点位移确定。
设节点i,j,m的坐标分别为(xi, yi)、(xj, yj)、(xm, ym)，其节点位移为(ui,vi)、(uj,vj)、(um,vm)，将它们代入上式得: ui=α1+α2xi+α3yi, vi=α4+α5xi+α6yi
uj=α1+α2xj+α3yj, vj=α4+α5xj+α6yj
um=α1+α2xm+α3ym, vm=α4+α5xm+α6ym(3.2)联立求解上述公式左边的3个方程，可以求出待定常数α1、α2、α3:α1=12Auixiyi
ujxjyj
umxmym,  α2=12A1uiyi
1ujyj
1umym,  α3=12A1xiui
1xjuj
1xmum(3.3)式中，A为三角形单元ijm的面积:A=121xiyi
1xjyj
1xmym  为使求得的面积A为正值，单元节点编号i,j,m的次序必须是逆时针转向。至于将哪个节点作为起始节点i，则没有关系。将式(3.3)代入式(3.1)的第一式，整理后得u=12A［(ai+bix+ciy)ui+(aj+bjx+cjy)uj
+(am+bmx+cmy)um］(3.4)同理可得v=12A［(ai+bix+ciy)vi+(aj+bjx+cjy)vj
+(am+bmx+cmy)vm］(3.5)式中ai=xjym-xmyj
bi=yj-ym
ci=-xj+xm (i=i,j,m)(3.6)下标i,j,m按顺序轮换。如令Ni=12A(ai+bix+ciy) (i=i,j,m)(3.7)就得到由节点位移表达单元内任一点位移的插值公式，即位移模式的另一形式:u=Niui+Njuj+Nmum
v=Nivi+Njvj+Nmvm (i=i,j,m)(3.8)函数Ni表示单元内部的位移分布形态，故Ni、Nj、Nm可称为单元的形状函数，简称为形函数。通常把式(3.8)写成矩阵形式:f=u
v=Ni0Nj0Nm0
0Ni0Nj0Nmδe=Nδe(3.9)式中,N称为形函数矩阵,其维数为2×6:N=Ni0Nj0Nm0
0Ni0Nj0Nm(3.10)3.2 单元应变和应力
有了单元的位移模式，就可以用几何方程求得单元的应变。将式(3.4)、式(3.5)代入式(2.24), 得到应变和节点位移的关系式: εx
εy
γxy=12Abi0bj0bm0
0ci0cj0cm
cibicjbjcmbmui
vi
uj
vj
um
vm(3.11)上式简写成ε=Bδe(3.12)式中，B为单元应变转换矩阵，它可以写成分块形式:B=［Bi Bj Bm］其中子矩阵为Bi=12Abi0
0ci
cibi (i=i,j,m)(3.13)由于矩阵B的元素都是常数，则由式(3.11)可知，单元内各点的应变也都是常数，这是采用了线性位移模式的必然结果。所以三节点三角形单元也称为常应变单元。
将式(3.12)代入物理方程(2.25)，得σ=DBδe=Sδe(3.14)这就是应力与节点位移的关系式。其中S为单元应力矩阵，也可写成分块形式: S=DB=［Si Sj Sm］对于平面应力问题，其子矩阵为Si=DBi=E2(1-μ2)Abiμci
μbici
(1-μ)ci/2(1-μ)bi/2 (i=i,j,m)(3.15)对于平面应变问题，只要把式(3.15)中的E换为E/(1-μ2)，把μ换为μ/(1-μ)，就可得到相应的子矩阵。
显然，S矩阵是常数矩阵，单元中的应力分量也是常量。这就造成在相邻单元的公共边上存在着应力突变现象。但是随着网格的细分，这种应力突变将会急剧减小，有限元的解答将收敛于正确结果。
3.3 单元平衡方程与单元刚度矩阵
整个结构处于平衡状态，所划分出的一个小单元体同样处于平衡状态，而结构的平衡条件可通过节点的平衡条件表示。有限元的任务就是要建立和求解整个弹性体的节点位移和节点力之间关系的平衡方程。为此首先要建立每一个单元的节点位移和节点力之间关系的平衡方程。单元平衡方程可以利用最小势能原理建立，也可以利用虚功原理求解。下面采用虚功方程来建立单元的平衡方程。
作用在每个单元上的载荷都移置到各节点之后，各单元所受到的力就只有通过节点传递的节点力。图3.4表示作用于单元e上的节点力Fe以及相应的应力分量σ，它们使单元处于平衡。在每个节点处，沿x和y方向有2个节点力分量，每个单元共6个节点力分量，组成单元节点力列阵为Fe=［Xi Yi Xj Yj Xm Ym］T(3.16)  设想单元节点发生了虚位移f*=［u* v*］T，相应的节点虚位移如图3.5所示，其6个分量组成单元节点虚位移列阵为
δ*=［u*i v*i u*j v*j u*m v*m］T(3.17) 图3.4 三角形单元节点力
图3.5 三角形单元节点虚位移
在单元内部引起的虚应变为ε*=［ε*x ε*y ε*z］T(3.18)根据虚功原理，外力虚功等于内力虚功。所以节点力在节点的虚位移上所做的虚功应等于单元内部应力在虚应变上所做的虚功。这就是单元保持平衡状态所必须满足的条件，即单元的平衡条件。δ*eTFe=ε*Tσtdxdy(3.19)式中，t为单元厚度。由式(3.12)和式(3.14)可知，ε*=Bδ*e, σ=DBδe，将它们代入式(3.19), 并因虚位移是任意的，可以在等号两边同时消去，得Fe=BTDBtdxdyδe(3.20)令ke=BTDBtdxdy(3.21)则Fe=keδe(3.22)这就建立了单元节点力与节点位移之间的关系，通常称为单元刚度方程，其实质代表的是单元平衡方程。式中的ke称为单元刚度矩阵。由于D中的元素和t对某一给定的单元来说都是常量，而且在选取线性位移函数的情况下，B中的元素也是常量，所以式(3.21)可写成:ke=BTDBtdxdy=BTDBtA(3.23) 式中,A=dxdy为该单元的面积。通常，单元刚度矩阵式(3.23)多用分块形式来推导并表达，以便于编制计算程序，ke由9个分块组成: ke=BTi
BTj
BTmD［Bi Bj Bm］tA=kiikijkim
kjikjjkjm
kmikmjkmme(3.24)其中每个子块是2×2的子矩阵:(krs)e=(Br)TDBstA=kxxrskxyrs
kyxrskyyrse (r,s=i, j, m)对于平面应力问题(krs)e=12Abr0cr
0crbr·E1-μ21对
μ1称
001-μ2·12Abs0
0cs
csbs·tA
=Et4(1-μ2)Abrbs+1-μ2crcsμbrcs+1-μ2crbs
μcrbs+1-μ2brcscrcs+1-μ2brbs
 (r,s=i, j, m)(3.25)对于平面应变问题，把上式中的E换为E1-μ2, μ换为μ1-μ即可。于是得(krs)e=E(1-μ)t4(1+μ)(1-2μ)Abrbs+1-2μ2(1-μ)crcsμ1-μbrcs+1-2μ2(1-μ)crbs
μ1-μcrbs+1-2μ2(1-μ)brcscrcs+1-2μ2(1-μ)brbs
 (r,s=i, j, m)(3.26)  单元刚度矩阵具有以下性质。
(1) 单元刚度矩阵中每个元素有明确的物理意义。其物理意义是单位节点位移分量所引起的节点力。例如，kmn是表示当单元第n个自由度产生单位位移而其他自由度固定时，在第m个自由度产生的节点力。
(2) ke是对称矩阵。其元素之间有如下关系: krs=ksr，这个特性是由弹性力学中功的互等定理所决定的。
(3) ke是奇异矩阵。ke在无约束的条件下，单元可作刚体运动。根据行列式性质，可知|k|值也为零。
由上可见，单元刚度矩阵ke反映了单元节点力与节点位移之间的关系。ke中的每一个元素表明该单元的各节点分别沿坐标轴方向发生单位位移时(其余节点位移分量为零)引起的各节点的节点力分量。ke中诸元素的值越大就表示单元的刚度(即抵抗单元变形的能力)越大。单元刚度矩阵ke取决于: ①单元的位移函数(因为B是由位移函数求偏导所得);  ②单元的几何参数(即单元的形状、大小和方位);  ③单元的材料性质(就是指材料的弹性常数E和μ)，但ke不随坐标的平行移动而改变。
3.4 单元等效节点力
作用在弹性体上的载荷可分为集中力、体积力(自重、惯性力及离心力等)和面力三种。在用有限元法求解时，需要把单元所承受的载荷全部移置到节点上，形成等效节点载荷。这种移置是按静力等效原则进行的。所谓静力等效原则是指原载荷与移置后的等效节点载荷，在弹性体产生任何虚位移的过程中，所做的虚功相等。在一定的位移函数下，这种移置的结果是唯一的，而且总能符合通常对刚体而言的静力等效原则，即原载荷与等效节点载荷在任一坐标轴上投影之和应相等，对任一坐标轴的力矩之和也应相等。
设单元e上任一点m(x, y)作用一集中力R, R在坐标轴x、y方向的分量分别为Rx、Ry，即:R=［Rx Ry］T。假设集中力R向单元的三个节点移置而得到单元等效节点载荷，其单元节点载荷列阵Pe可表示为Pe=Pi
Pj
Pme=PexiPeyiPexjPeyjPexmPeymT现在，假设该单元发生一个任意的虚位移，则m点相应的虚位移为f=［u v］T=Nδ*e而该单元各节点处相应的虚位移为δe=［uiviujvjumvm］T根据虚功原理有(δe)TPe=f*TR=(Nδe)TR=(δe)TNTR由于δe是任意的，故(δe)T也是任意的，因此上式两边与它相乘的矩阵应相等。于是得:Pe=NTR(3.27)或Pe=Ni0Nj0Nm0
0Ni0Nj0NmTRx
Ry〗
=［NiRxNiRyNjRxNjRyNmRxNmRy］T(3.28)这就是单元内作用一集中力向节点移置的公式。利用此公式，可以类似求得单元内分布的体力或单元边界上的面力向节点移置的公式。
设上述单元受有体力作用，单元内单位体积内的体力为P，在坐标轴x、y方向的分量为Px、Py，即P=［Px Py］T。可将微体积tdxdy上的体力Ptdxdy看作是一集中载荷，利用式(3.27)并积分得到体力向节点移置的公式: (PV)e=NTPtdxdy(3.29)式中积分是对三角形单元的面积进行的。
如果上述单元的某一边界s上受有分布的面力q，其分量为qx、qy，即q=［qxqy］T。可将边界上微面积tds上的面力qtds看作集中载荷，利用式(3.27)并积分即可得到面力向节点移置的公式:(Ps)e=∫NTqtds(3.30)式中积分是对受有分布面力的边界线s进行的。
总之，前述式(3.27)、式(3.29)和式(3.30)分别表示作用在单元体内的集中力、体力和单元边界上的分布面力向节点移置成等效节点载荷的积分公式。无论位移函数是线性函数或非线性函数，均可适用。当然，选用的位移函数不同，所得的结果就不一样。所以，上述三个积分公式也称为载荷移置的一般公式。
以体力为例，设单元受有总体积力W=［WxWy］T，当单元划分得较小时，单位体积内的体力Pe可认为在单元内是均布的:Pe=［INiINjINm］TdxdytP(3.31)而形函数Ni、Nj、Nm与微面积dxdy乘积的积分为Nidxdy=12A(ai+bix+ciy)dxdy
=12Aaidxdy+bixdxdy+ciydxdy由三角形形心计算公式可得形心的坐标xc、yc为xc=xdxdydxdy=13(xi+xj+xm)
yc=ydxdydxdy=13(yi+yj+ym)由dxdy=A，得xdxdy=A3(xi+xj+xm)
ydxdy=A3(yi+yj+ym)将式(3.6)及式(3.7)代入以上诸式，积分得: Nidxdy=12A(xjym-xmyj)A+(yj-ym)A3(xi+xj+xm)
 +(xm-xj)A3(yi+yj+ym)
=16［xi(yj-ym)+xj(ym-yi)+xm(yi-yj)］
=13A同理可得: Njdxdy=13A,  Nmdxdy=13A则式(3.31)的计算结果为Pe=13［III］TAtP因为AtP=W=［WxWy］T，所以Pe=［PexiPeyiPexjPeyjPexmPeym］T
=13Wx13Wy13Wx13Wy13Wx13WyT这就是说，单元上的总体积力平均分配到三个节点上即得等效节点载荷。
如果单元所受的总体积力W是自重，则其等效节点载荷列阵为Pe=-W3［010101］T(3.32)其中的负号表示等效节点载荷分量的方向与坐标的y轴方向相反。
实际上，只要位移函数是线性的，对于作用在边界上的分布力或集中力可更简单地直接按静力学原则把载荷移置到相邻两节点上(而不是移置到三个节点上)，其结果与按虚功等效原则移置所得完全一致。例如某单元边界ij边上作用了三角形分布面力，移置后的单元等效节点载荷列阵为Pe=［PeiPejPem］T=12qjtl13230T(3.33)式中，l表示三角形单元边界ij的长度，即把总载荷12qjtl的三分之一移置到节点i，三分之二移置到节点j，而等效节点载荷的方向与原载荷相同。
如果单元边界线上的面力为均布载荷，则总载荷平均分配到边界两端的节点上。如果在边界线上作用有一集中载荷，则可按集中载荷作用点至边界两端节点之间的距离成反比分配至两端节点上。一般在结构离散时，多在集中载荷作用处安置一个节点，这样就可免去移置载荷的手续。这里，必须强调指出，上述非节点载荷的简单移置法只在位移函数为线性的条件下才适用，对于非线性位移函数，则载荷的移置必须按前述普遍公式(3.27)、公式(3.29)、公式(3.30)进行计算。现行有限元通用程序中体力、面力都可直接施加，程序会自动处理，但为了施加集中载荷，一般需在建模时预留节点位置，以便加载。
3.5 整体平衡方程与总刚度矩阵
单元平衡方程建立了单元节点力与节点位移之间的关系。对每个单元都可以建立单元平衡方程，将所有单元平衡方程集合在一起，就得到结构的整体平衡方程:Kδ=R(3.34)式中,δ是整个结构上节点位移分量的列阵，设节点个数为n，每个节点有2个位移分量，δ=［ui vi］T是二阶向量，按照从小到大的节点编号列阵为 δ=［δT1 δT2  δT3 … δTn］T(3.35)式中,R是整个结构上节点力分量的列阵，每个节点有2个节点力分量，Ri=［Rxi Ryi］T，则R=［RT1 RT2 RT3 … RTn］T(3.36)K为结构的整体刚度矩阵，一般称为总刚度矩阵，其维数为2n×2n，可写成分块形式: K11K12……K1n
K21K22……K2n


Kn1Kn2……Knnδ1
δ2


δn=R1
R2


Rn(3.37)子矩阵Kij是2×2矩阵。组集的原则是将每个单元的子矩阵(Krs)2×2置于整体刚度矩阵的第r行、第s列上。根据各单元节点的局部编号(i,j,m)与整体编号(1,2,3,…,n)的对应关系，把各Ke子块的下标换成整体编号，“对号入座”地把各Ke的子块放入K中相应位置，凡是下标相同的子块，应放入K中同一位置叠加，如此形成总刚度矩阵。
总刚度矩阵具有以下一些性质。
(1) 总刚度矩阵是对称矩阵。单元刚度矩阵是对称矩阵，由它按对称方式组集的总刚度矩阵必然也是对称矩阵。
(2) 总刚度矩阵呈稀疏带状分布。因为任一节点只与绕它的相连单元发生联系，所以K中的每一行含有大量的零元素，而非零的元素往往分布在对角线主元素的附近。
(3) 总刚度矩阵是奇异矩阵。因为弹性体在外力作用下处于平衡，在平面问题中应满足3个平衡方程，反映在K中就是存在3个线性相关的行列式，因而是奇异的。只有在排除刚体位移后，K才是正定的。
根据总刚度矩阵的对称、稀疏带状性质，按照一定规则进行节点编号，可以节省大量计算机存储量和运算工作量。
有限元法的核心是: 它不是寻求整个域上连续的解析解，而是去寻找在每个子域(单元)上，近似满足基本方程的分片插值函数解。它将物体人为地划分成一个个单元，在单元分析的基础上，再集合单元，对整体结构进行综合分析。
3.6 边界条件处理
由于总刚度矩阵是奇异矩阵，不存在逆阵，为了求解整体平衡方程，必须消除刚体位移，引入位移约束条件。从力学角度来看，进行边界条件约束处理，消除结构由于外力作用而产生刚体位移，使刚度矩阵从半正定阵成为正定阵，使解具有唯一性。如果物体本身有位移约束，可根据实际情况直接引入；如果物体仅有应力边界条件而无位移约束条件，则需人为的加上消除刚体位移的约束条件。例如对平面问题，就要引入3个不在同一点的x、y方向的位移约束。有限元法中的边界条件也是假定在节点上受到约束。位移约束可分为零位移约束和非零位移约束两类，相当于在节点上有支承或节点位移为已知非零量。有限元法中通常采用两种方法，即划行划列法和乘大数法来引入位移约束条件。前者适用于简单问题的手算练习，后者适合于实际问题的计算机处理。
由于K矩阵中与零位移对应的行和列上的元素，在求解其余位移时不起作用，因而可将它们划去，这使得修正后的总刚度矩阵降低了阶数，行列式不等于零，从而可以解出其余节点位移，这就是划行划列法。乘大数法是将K矩阵中与指定位移对应的主对角元素乘上一个非常大的数，K矩阵中其他元素不变，同时将载荷列阵R中对应元素换为节点位移指定值与扩大了的主对角元素的乘积。用以上方法进行位移约束条件处理后，总刚度方程得到修正，方程可解。
由上可见，进行约束处理是方便的，但要注意的是如何正确地确定这些约束。对于一般固定的结构物，只要根据实际结构的具体约束情况与变形情况，加上相应的约束条件(包括强迫位移)即可；对于承受自相平衡的力系作用，又无实际位移约束条件的结构物，在平面问题中，可人为地规定某一节点s的两个位移分量为零，再规定另一节点在绕s点转动的切线方向的位移分量为零即可。也就是说，所加约束的数量应正好使结构消除刚体运动，而约束所在位置则可以任意选定。第2章已指出，所加的约束位置不同，解出的节点位移值是不同的(相差一刚体位移值)，但结构内各点的相对位移值却仍然不变，因此并不影响应力值。