力学WORD教案CHAPTER3弹性力学平面问题的有限元法

1、3 弹性力学平面问题的有限元法本章包括以下的内容：3.1 弹性力学平面问题的基本方程 3.2 单元位移函数 3.3 单元载荷移置 3.4 单元刚度矩阵 3.5 单元刚度矩阵的性质与物理意义 3.6 整体分析 3.7 约束条件的处理 3.8 整体刚度矩阵的特点与储备方法 3.9 方程组解法3.1 弹性力学平面问题的基本方程弹性力学是讨论弹性体在约束和外载荷作用下应力和变形分布规律的一门学科；在弹性 力学中针对微小的单元体建立基本方程,把复杂外形弹性体的受力和变形分析问题归结为偏 微分方程组的边值问题；弹性力学的基本方程包括平稳方程、几何方程、物理方程；弹性力学的基本假定如下：1）完全弹性, 2）

2、连续, 3）匀称, 4）各向同性, 5）小变形；3.1.1 基本变量 弹性力学中的基本变量为体力、面力、应力、位移、应变,各自的定义如下；体力体力是分布在物体体积内的力,例如重力和惯性力；面力 面力是分布在物体表面上的力,例如接触压力、流体压力；应力物体受到约束和外力作用,其内部将产生内力；物体内某一点的内力就是应力；图 3.1 如图 3.1 假想用通过物体内任意一点 p 的一个截面 mn 将物理分为、两部分；将部分撇开, 依据力的平稳原就,部分将在截面 mn上作用肯定的内力；在 mn截面上取包含p 点的微小面积 A ,作用于 A 面积上的内力为 Q ；令 A无限减小而趋于 p 点时,Q 的极

3、限 S就是物体在 p 点的应力；lim A 0 QA S应力 S 在其作用截面上的法向重量称为正应力,用 表示；在作用截面上的切向重量称为剪应力,用 表示；明显, 点 p 在不同截面上的应力是不同的；为分析点 p 的应力状态, 即通过 p 点的各个截面上的应力的大小和方向,在 p 点取出的一个平行六面体,六面体的各楞边平行于坐标轴；图 3.2 将每个上的应力分解为一个正应力和两个剪应力,分别与三个坐标轴平行；用六面体表面的应力重量来表示 p 点的应力状态；应力重量的下标商定如下：第一个下标表示应力的作用面,其次个下标表示应力的作用方向；xy,第一个下标x 表示剪应力作用在垂直于X 轴的面上,其

4、次个下标y 表示剪应力指向 Y 轴方向；正应力由于作用表面与作用方向垂直,用一个下标；x表示正应力作用于垂直于X 轴的面上,指向X 轴方向；应力重量的方向定义如下：假如某截面上的外法线是沿坐标轴的正方向,为正；假如某截面上的外法线是沿坐标轴的负方向,为正；这个截面上的应力重量以沿坐标轴正方向 这个截面上的应力重量以沿坐标轴负方向剪应力互等：xyyx,yzzy,zxxzx、y、z、xy、yz、zx物体内任意一点的应力状态可以用六个独立的应力重量来表示；位移位移就是位置的移动；物体内任意一点的位移,用位移在x,y,z 坐标轴上的投影u、v、w 表示；应变物体的外形转变可以归结为长度和角度的转变；各

5、线段的单位长度的伸缩,称为正应变,用 表示；zx 表示；两个垂直线段之间的直角的转变,用弧度表示,称为剪应变,用物体内任意一点的变形,可以用x、y、z、xy、yz、六个应变重量表示；3.1.2 平面应力和平面应变问题 弹性体在满意肯定条件时,其变形和应力的分布规律可以用在某一平面内的变形和应力的分布规律来代替,这类问题称为平面问题；平面问题分为平面应力问题和平面应变问题；1）平面应力问题设有很薄的等厚薄板,只在板边上受到平行于板面并且不沿厚度变化的面力,体力也平行于板面且不沿厚度变化；图 3.3 设板的厚度为t,在板面上：t0,zyzt0zzt0,zxz222由于平板很薄,外力不沿厚度变化,因

6、此在整块板上有,z0,zx0,zy0 x、y、xy未知；剩下平行于XY 平面的三个应力重量2）平面应变问题设有很长的柱形体,支承情形不沿长度变化,在柱面上受到平行于横截面而且不沿长度变化的面力,体力也如此分布；图 3.4 以柱体的任一横截面为XY 平面,任一纵线为Z 轴；假定该柱体为无限长,就任一截面都可以看作对称面；由对称性,zx0,zy0,w0z0；y、xy,物体在 Z 方向处于自平稳状由于没有Z 方向的位移, Z 方向的应变未知量为平行于XY 平面的三个应力重量x、态；3.1.3 平稳方程 弹性力学中, 在物体中取出一个微小单元体建立平稳方程；平稳方程代表了力的平稳关系,建立了应力重量和

7、体力重量之间的关系；对于平面问题,在物体内的任意一点有,xyxX0（3-1）xyyxyY0yx3.1.4 几何方程 由几何方程可以得到位移和变形之间的关系；对于平面问题,在物体内的任意一点有,xyxuuv（3-2）x vyyyx刚体位移由位移 u=0,v=0 可以得到应变重量为零,反过来,应变重量为零就位移重量不为零；应变重量为零时的位移称为刚体位移；刚体位移代表了物体在平面内的移动和转动；由uuv00 x v0yyx可以得到刚体位移为以下形式,uu0y0可得,0可得,vv 0 x由u0,vxyuf1y,vf2x将f1, f2代入uvyxdf1ydf2xdydx积分后得到,f 1y u 0yf

8、2x v 0 x得到位移重量,uuu00 ,y,00时,物体内任意一点都沿x 方向移动相同的距离,可见u0代vv 0 x当0v0表物体在 x 方向上的刚体平移；当u00 ,v 00 ,0时,物体内任意一点都沿y 方向移动相同的距离,可见v 代表 0物体在 y 方向上的刚体平移；当u00 ,v 00 ,0时,可以假定0 ,此时的物体内任意一点P（x,y）的位移重量为,uy,vxP 点位移与 y 轴的夹角为 ,tgyytgxxP 点合成位移为,u2v2y2x2x2y2rr 为 P 点到原点的距离,可见 代表物体绕z 轴的刚体转动；3.1.5 物理方程弹性力学平面问题的物理方程由广义虎克定律得到；1

9、）平面应力问题的物理方程x1y（3-3）Exy1xE2 1yxyxyE平面应力问题有,z0Exyz2）平面应变问题的物理方程x12x1y（3-4）Ey12y1xExy2 1xyE平面应变问题有,z0 xyz在平面应力问题的物理方程中,将E 替换为1E2、替换为1,可以得到平面应变问题的物理方程；在平面应变问题的物理方程中,将E 替换为E 12、替换为 121,可以得到平面应力问题的物理方程；图 3.5 求解弹性力学平面问题,可以归结为在任意外形的平面区域 内已知掌握方程、在位移边界 S 上约束已知、 在应力边界 S 上受力条件已知的边值问题；然后以应力重量为基本未知量求解,或以位移作为基本未知

10、量求解；假如以位移作为未知量求解,求出位移后, 由几何方程可以运算出应变重量,得到物体的变形情形； 再由物理方程运算出应力重量,得到物体的内力分布,就完成了对弹性力学平面问题的分析；3.2 单元位移函数依据有限元法的基本思路,将弹性体离散成有限个单元体的组合,以结点的位移作为未知量；弹性体内实际的位移分布可以用单元内的位移分布函数来分块近似地表示；在单元内的位移变化可以假定一个函数来表示,这个函数称为单元位移函数、或单元位移模式；对于弹性力学平面问题,单元位移函数可以用多项式表示,ua 1a2xa3ya4x2a5xya6y2.（3-5）由单元形vb 1b2xb 3yb4x2b 5xyb 6y2

11、.多项式中包含的项数越多,就越接近实际的位移分布,越精确；详细取多项,式来确定；即以结点位移来确定位移函数中的待定系数；图 3.6 如图 3.6 所示的 3 结点三角形单元,结点I、J、M 的坐标分别为x iyi、xjyj、xmym,结点位移分别为iu、v 、u 、v 、u 、vm；六个节点位移只能确定六个多项式的系数,所以3 结点三角形单元的位移函数如下,ua 1a2xa3y（3-6 ）va 4a 5xa 6y将 3 个结点上的坐标和位移重量代入公式（位移重量表示出来；3-6 ）就可以将六个待定系数用结点坐标和将水平位移重量和结点坐标代入（3-6 ）中的第一式,u i a 1 a 2 x i

12、 a 3 y iu j a 1 a 2 x j a 3 y ju m a 1 a 2 x m a 3 y m写成矩阵形式,u i1xiyia 1（3-7）uj1xjyja2um1xmyma3令1xiy i（3-8）1xjyjT,就有1x myma 1u iT1uja2a3umT1T*TT2 A,A 为三角形单元的面积；T 的相伴矩阵为,令T*xjymx myjyjy maxmxjT（3-9）xmy ix iymymy ixix mT*x iyjxjyiyiyjxjx i（3-10）Ta ib ic ia iamjajbjcjb ibb mj就a1amb mc mamc icjc m（3-11）

13、1aiajuia2b ibbmuj2Aa3c icjc mum同样,将垂直位移重量与结点坐标代入公式（3-6）中的其次式,可得,a41a iajamvi（3-12）a5b ibjb mvj2A6jac iccmvm将（ 3-11）、（3-12）代回（ 3-6）整理后可得,u1a ib ixc iy uiajbjxcjy uja mb mxc my um2 A1va ib ixc iy viajbjxcjyvj amb mxcmy vm2A令Ni1a ib ixciy（下标 i,j, m 轮换）2Au i可得uNi0eNj0jNm0vi（3-13）ujv0Ni0uN0Nmvjumvm单元内的位移

14、记为fiuivvi单元的结点位移记为ujjvjmumvm单元内的位移函数可以简写成,fNeN i 称为外形函数；（3-14）把N 称为外形矩阵,挑选单元位移函数应满意以下条件：1）反映单元的刚体位移与常量应变；2）相邻单元在公共边界上的位移连续,即单元之间不能重叠,也不能脱离；由（ 3-6）可以将单元位移表示成以下的形式,ua1a2xa52a3ya52a3yva4a 6ya52a3xa52a 3x反映了刚体位移和常应变；单元位移函数是线性插值函数,因此单元边界上各点的位移可以由两个结点的位移完全 确定；两个单元的边界共用两个结点,所以边界上的位移连续；外形函数 N i 具有以下性质：1）在单元

15、结点上外形函数的值为 1 或为 0；2）在单元中的任意一点上,三个外形函数之和等于 1；用 T 来运算三角形面积时,要留意单元结点的排列次序,当三个结点 i,j,m 取逆时针次序时,iA1T0；当三个结点i, j,m 取顺时针次序时,A1T0；22例题：如图3.7 所示等腰三角形单元,求其外形矩阵N ；解： 由axjymxmyjbiyjymcixmxj在公式中轮换下标可以运算得aixjymxmyj000a0,b imyjymja0a,cixmxj000a00,bjymyi000,ajxmy ixiym0cjxixma0a00a2,byiy0aaamxiyjxjyiaacmxjxi0aa三角形积

16、为Aa22外形函数为Ni1 aijb ixc iy 10ax0 xy2Aa2aNj1abjxcjy100ay2Aa2aNm1 amb mx0cmy x1a2axay 1xy2Aa2aa外形矩阵为Nx0y1y10yaaaaxyx000aaaa三角形面积的运算公式可得,A11xiyi11a01a21xjyj10a2221xmym100假如把三个结点按顺时针方向排列,即i（ a,0）, j（0,0）,m（0,a）A11xiyi11a01a21xjyj1002221xmym10a3.3 单元载荷移置有限元法的求解对象是单元的组合体,因此作用在弹性体上的外力,需要移置到相应的结点上成为结点载荷；载荷移置

17、要满意静力等效原就；静力等效是指原载荷与结点载荷在任意虚位移上做的虚功相等；单元的虚位移可以用结点的虚位移*e表示为,（3-15）f*N*eXiY i令结点载荷为ReXjYjXmY m1）集中力的移置如图 3.7 所示,在单元内任意一点作用集中力PP xP y图 3.8 由虚功相等可得,*eTR e*eTNTP （3-16）由于虚位移是任意的,就R eNTP 例题 1：在均质、等厚的三角形单元ijm 的任意一点p（x p,y p）上作用有集中载荷；XiNi0P xY i0NiXjNj0Yj0NjP yXmNm0Y m0Nmxp,yp2）体力的移置令单元所受的匀称分布体力为pxy由虚功相等可得,

18、*eTR ep *TNTp tdxdy（ 3-17）R eNTtdxdy3）分布面力的移置设在单元的边上分布有面力PX,YT,同样可以得到结点载荷,R es NTPtds（3-18）例题 2：设有均质、等厚的三角形单元ijm ,受到沿 y 方向的重力载荷qy 的作用；求均布体力移置到各结点的载荷；XiN i0ciyc1AY0NiXjjNj0j0tdxdyqyY0NXmNm0Y m0NmXi0 ,Xj0,Xm0YNiqytdxdyqytNidxdyNidxdy1aibixciydxdy2A1aiAb iAxcciAycA1a ib ixc2A2A3Y1qyAt3同理,Yj1qyAt,Y m1qy

19、At33x 方向按三角形分布的载荷,求例题 3：在均质、等厚的三角形单元ijm 的 ij 边上作用有沿移置后的结点载荷；XiNi0ij 边上,以局部坐标表示Y0NiXjNj0jqxtdxdyYjs0N0XmNm0Y m0Nm取局部坐标s,在 i 点 s=0,在 j 点 s=l,L 为 ij 边的长度；在的插值函数为,Ni1s,N js,Nm0L1qtLLL载荷为q xqsqts2s3LtdsXiL 1sqs0LL2 L2 3 L06XjLsqstdsqt3 sL1 3qtL0LL3 L203.4 单元刚度矩阵 依据单元的位移函数,uNi0Nj0jNm0uiviujv0Ni0N0Nmvjumvm

20、由几何方程可以得到单元的应变表达式,u0Bui（3-19）v iux vv1bi0bj0bmuj0ci0cj0cmy2Avjcib icjbjcmbm（3-20）umyxv m记为Be,B 矩阵称为几何矩阵；jmB 矩阵可以表示为分块矩阵的形式BB iBBi1b i00ci2Acib i由物理方程,可以得到单元的应力表达式,DDBe（3-21）D 称为弹性矩阵,对于平面应力问题,D 1E101020012定义SDB为应力矩阵；将应力矩阵分块表示为,SS iSjS mS iDB i2E21b ic i1c ib i（3-22）b iciA 122应用虚功原理可以建立单元结点位移与结点力的关系矩阵

21、,单元刚度矩阵；虚功原理： 在外力作用下处于平稳状态的弹性体,假如发生了虚位移,就全部外力在虚位移上做的虚功等于内应力在虚应变上做的虚功；单元的结点力记为*FeBUieV iUjVjUmV mT单元的虚应变为*单元的外力虚功为,*eTFe单元的内力虚功为,*Ttdxdy由虚功原理可得,*eTFe*Ttdxdy（3-23）*e（3-24）*eT*eTB TTBSeDBe*eTF*eTB TDBtdxdyeFeBTDBtdxdye定义KeBTDBtdxdy为单元刚度矩阵；在 3 结点等厚三角形单元中B 和D 的重量均为常量,就单元刚度矩阵可以表示为,KeBTDBtA（3-25）单元刚度矩阵表示为分

22、块矩阵：KeKiiKijKim（3-26）KjiKjjKjmKrsKmiKmjKmmBrTDBsKrsKrx,sxKrx,sy（3-27）Kry,sxKry,sy3.5 单元刚度矩阵的性质与物理意义（一）单元刚度矩阵的物理意义Kix假设单元的结点位移如下：e100000 T由FeKee,得到结点力如下：（3-28）UiKix,ixV iKiy,ixUjKjx,ixVjKjy,ixUmKmx ,ixi 的水平方向上需要施加的结点力；VmKmy, ix, ix表示 i 结点在水平方向产生单位位移时,在结点Kiy, ix表示 i 结点在水平方向产生单位位移时,在结点i 的垂直方向上需要施加的结点力；

23、挑选不同的单元结点位移,可以得到单元刚度矩阵中每个元素的物理含义：Krx,sx表示 s 结点在水平方向产生单位位移时,在结点r 的水平方向上需要施加的结点力；Kry,sx表示 s 结点在水平方向产生单位位移时,在结点r 的垂直方向上需要施加的结点力；Krx,sy表示 s 结点在垂直方向产生单位位移时,在结点r 的水平方向上需要施加的结点力；Kry,sy表示 s 结点在垂直方向产生单位位移时,在结点r 的垂直方向上需要施加的结点力；因此单元刚度矩阵中每个元素都可以懂得为刚度系数,加的力；（二）单元刚度矩阵的性质1）对称性 利用分块矩阵的性质证明如下：KrsB rTDBsKsrBsTDB r即在结

24、点产生单位位移时需要施KsrTBsTDBrTB rTDTBsB rTDBsKrs即KeKeT2）奇特性即单元刚度矩阵的行列式为零,Ke0；将定单元产生了x 方向的刚体移动,e101010 T,此时对应的单元结点力为零；011由行列式的性质可000Ke1000100可以得到,在单元刚度矩阵中1,3,5 列中对应行的系数相加为零,知,Ke0；01 T,可以得同样假如假定单元产生了y 方向上的刚体位移e010到,在单元刚度矩阵中2,4,6 列中对应行的系数相加为零；3.6 整体分析得到了单元刚度矩阵后,要将单元组成一个整体结构,依据结点载荷平稳的原就进行分析,即整体分析；整体分析包括以下 4 个步骤

25、：1）建立整体刚度矩阵,2）依据支承条件修改整体刚度矩阵,3）解方程组,求出结点的位移,4）依据结点位移,求出单元的应变和应力；在这里把结点位移作为基本未知量求解；成；如何得到整体刚度矩阵？基本方法是刚度集成法,即整体刚度矩阵是单元刚度矩阵的集图 3.9 如图 3.9 所示,一个划分为 6 个结点、 4 个单元的结构；得到了每个单元的单元刚度矩阵后,要集成为整体刚度矩阵；3.6.1 刚度集成法的物理意义由单元刚度矩阵的物理意义可知,起的单元结点力；单元刚度矩阵的系数是由单元结点产生单位位移时引在如图 3.9 所示的结构中,使结点 3 产生单位位移时,在单元（1）中的结点 2 上引起结点力；由于

26、结点 2、3 同时属于单元（1）、（3）,在单元（ 2）中的结点 2 上同样也引起结点力,因此,在整体结构中当结点 3 产生位移时,结点 2 上的结点力应当是单元（1）、（2）在结点 2 上的结点力的叠加；刚体集成法即结构中的结点力是相关单元结点力的叠加,整体刚度矩阵的系数是相关单元的单元刚度矩阵系数的集成；结点 3 在整体刚度矩阵的对应系数,应当是单元 （1）、（3）、（4）中对应系数的集成；3.6.2 刚度矩阵集成的规章1）将单元刚度矩阵中的每个分块放到在整体刚度矩阵中的对应位置上,得到单元的扩大刚度矩阵；单元刚度矩阵系数取决于单元结点的局部编号次序,必需知道单元结点的局部编号与该结点在整

27、体结构中的总体编号之间的关系,才能得到单元刚度矩阵中的每个分块在整体刚度矩阵中的位置； 将单元刚度矩阵中的每个分块按总体编码次序重新排列后,可以得到单元的扩大矩阵；假定单元结点的局部编号与整体的对应关系如下：单元编号单元结点局部编号单元结点整体编号6 1 i 3 1 j 1 1 m 2 2 i 5 2 j 2 2 m 4 3 i 5 3 j 3 3 m 2 4 i 3 4 j 5 4 m 6 单元（ 2）的单元扩大矩阵 K2的分块矩阵形式如下,只列出非零的分块：局部 编号1jm 11i整体1 2 3 4 5 编号Kjj2Kjm2Kji21j1 m 12 1i3 Kmj2Kmm2Kmi24 5

28、Kij2Kim2Kii26 2）将全部单元的扩大矩阵相加得到整体刚度矩阵；KK 1K2K3 K4整体刚度矩阵如下所示：整体1 jj 1 2 3 4 jm25 6 im4编号KKjm 1Kji 1 1 KKji22 Kmj 1 Kmm 1Kmi 1 3 Kij1 +Kjj2 +Kmj3 Kmm2+Kmi3 K+Kmm3 Kii 1 Kji3 Kim1 Kjm3+Kjj3+Kij44 +Kii4 Kmj2Kmi25 Kij2Kij3 Kim2Kii2Kjm46 +Kim3+Kji4 +Kii3Kmm4+Kjj4Kmi4Kmj43.7 约束条件的处理图 3.9 所示的结构的约束和载荷情形,如图3.1

29、0 所示；结点1、4 上有水平方向的位移约束,结点4、6 上有垂直方向的约束,结点3 上作用有集中力（Px,Py）；图 3.10 整体刚度矩阵 K 求出后,结构上的结点力可以表示为：FK用 P 表示结点载荷和支依据力的平稳, 结点上的结点力与结点载荷或约束反力平稳；杆反力,就可以得到结点的平稳方程：KP（3-29）这样构成的结点平稳方程组,在右边向量 P 中存在未知量, 因此在求解平稳方程之前,要依据结点的位移约束情形修改方程3-29；先考虑结点n 有水平方向位移约束,与n 结点水平方向对应的平稳方程为：K2 n,11 u1K2n,12v 1.K2n1,2n1unK2n,12nvn.P 2n1

30、（3-30）依据支承情形,方程（3-30）应当换成下面的方程：un0（3-31）对比公式（ 3-30）和（ 3-31）,在式（ 3-29）中应当做如下修正：在K 矩阵中, 第 2n-1 行的对角线元素 K 2 n ,1 2 n 1 改为 1,该行中全部非对角线元素改为0；在 P 中,第 2n-1 个元素改为 0；为了保持 K 矩阵的对称性,将第 2n-1 列的全部非对角元素也改为 0；同理,假如结点 n 在垂直方向有位移约束,就（3-29）中的第 2n 个方程修改为,v n 0在K 矩阵中,第 2n 行的对角线元素改为 1,该行中全部非对角线元素改为 0；在 P中,第 2n 个元素改为00；为

31、了保持 K 矩阵的对称性, 将第 2n 列的全部非对角元素也改为0；00u 1P 100v 1P 200u2P 300v 2P 400000000unP 2n1（3-32）0000010000000010000v nP 2n00000000对图 3.9 所示结构的整体刚度在修改后可以得到以下的形式,100000000000Et（3-33）*00*0*00*0*00*0对*00*0*00*01000002称10000*0*0*01假如结点n 处存在一个已知非零的水平方向位移u ,这时的约束条件为,*n成unu*（3-34）K2n1,2n1乘上一个大数A ,向量 P 中的对应换n在K 矩阵中,第2

32、n-1 行的对角线元素AK2n,12n1u*,其余的系数保持不变；n方程改为,K2n,11 u 1K2n,12v 1.AK2n,12n1u nK2n1,2nv n.AK2n,12n1 u*（3-35）nA 的取值要足够大,例如取1010；只有这样,方程（3-35）才能与方程（3-34）等价；假如结点n 处存在一个已知非零的垂直方向位移* v ,这时的约束条件为,vnv*；n也可以采纳同样的方法修改整体刚度矩阵；3.8 整体刚度矩阵的特点与储备方法用有限元方法分析复杂工程问题时,结点的数目比较多,整体刚度矩阵的阶数通常也是很高的； 那么,是否在进行运算时要储存整体刚度矩阵的全部元素？能否依据整体

33、刚度矩阵 的特点提高运算效率？整体刚度矩阵具有以下几个显著的特点：对称性,稀疏性,非零系数带形分布；1）对称性由单元刚度矩阵的对称性和整体刚度矩阵的集成规章,利用对称性,只储存整体矩阵上三角部分的系数即可；2）稀疏性可知整体刚度矩阵必为对称矩阵；单元刚度矩阵的多数元素为零,非零元素的个数只占较小的部分；如图 3.11 所示的结构,结点 2 只和通过单元联接的 1、3、4、5 结点相关, 结点 5 只和通过单元联接的 2、3、4、6、8、9 结点相关；由单元刚度矩阵的物理意义和整体刚度矩阵的形成方式可知,相关结点 2、3、4、6、8、9 及结点 5 本身产生位移时,才使结点5 产生结点力,其余结

34、点产生位移时不在该结点处引起结点力；在用分块形式表示的整体矩阵中,与相关结点对应的分块矩阵具有非零的元素,其它位置上的分块矩阵的元素为零,如图 3.12 所示；图 3.11 图 3.12 整体刚度矩阵的分块矩阵示意3）非零元素带形分布在图 3.12 中,明显可以看出,整体刚度矩阵的非零元素分布在以对角线为中心的带形区域内,这种矩阵称为带形矩阵；在包括对角线元素的半个带形区域内,每行具有的元素个数叫做半带宽,用 hbd 表示；d 相邻结点的编码的最大 差值 1 2图 3.11 所示结构的相邻结点编码的最大差值为 4,所以半带宽为 10；二维等带宽储备设整体刚度矩阵 K 为一个 n n 的矩阵,最

35、大半带宽为 d；利用带形矩阵的特点和对称性,只需要储存以 d 为固定带宽的上半带的元素,称为二维等带宽储备；进行储备时,把整体刚度矩阵 K 每行中的上半带元素取出,储存在另一个矩阵 K* 的对应行中,得到一个n d 矩阵 K * ；把元素在 K 矩阵中的行、列编码记为r、s,在矩阵 K*中的行、列编码记为r*、s*,对应关系如下：r*=r s *=s-r+1 图 3.13a 图 3.13b 如图 3.13a所示的最大半带宽为d 的整体刚度矩阵K ,采纳二维等带宽储备后得到如图 3.13b所示的矩阵 K*；用新的方法储备后,K 矩阵中的对角线元素储存在新矩阵中的第 1 列中, K 矩阵中的r 行

36、元素仍旧储存在新矩阵的r 行中, K 矩阵中的s 列元素就依据新的列编码储存在新矩阵的不同列中；采纳二维等带宽储备,需要储存的元素数量与 K 矩阵中的总元素数量之比为 d ；所存n储的元素数量取决于最大半带宽 d 的值, d 的值就由单元结点的编码方式打算；虽然在采纳二维等带宽储备时,仍旧会储存一些零元素,但是采纳这种方法时元素寻址很便利；图 3.14a 图 3.14b 对于同样的有限元单元网格,依据图 3.14a的结点编码,最大的半带宽为 14；依据图3.14b 的结点编码,最大的半带宽为 18；依据图 3.11 的结点编码,最大的半带宽为 10；3.9 线性方程组解法由于有限元分析需要使用

37、较多的单元,线性方程组的阶数很高,有限元求解的效率很大程度上取决于线性方程组的解法；率是关键；利用矩阵的对称、 稀疏、 带状分布等特点提高方程求解效线性方程组的解法分为两大类：直接解法,迭代解法；直接解法以高斯消去法为基础,包括高斯消去法、等带宽高斯消去法、三角分解法,以及适用于大型方程组求解的分块算法和波前法等；迭代算法有高斯-赛德尔迭代、超放松迭代和共轭梯度法等；在方程组的阶数不是特殊高时,通常采纳直接解法；当方程组的阶数过高时,为防止舍 入误差和消元时有效数缺失等对运算精度的影响,可以挑选迭代方法；这里给出了用Fortran 语言编写的等带宽高斯消去法的代码,其中 NROW 为矩阵行的数目, NHBW 为最大半带宽；SUBROUTINE SOLVERBNROW,NHBW,STIFF,DISPL C Band elimination method C . : C : SOLVE EQUA TIONS KS*H=Q