7.4-单元刚度矩阵组装及整体分析

7.4单元刚度矩阵组装及整体分析

7.4,1单刚组装形成总刚依据全结构的平衡方程（为=8）可知，总体刚度矩阵是由单元刚度矩阵集合而成的.假如一个结构的计算模型分成片个单元，那么总体刚度矩阵可由各个单元的刚度矩阵组装而成，即[K]是由每个单元的刚度矩阵的每个系数按其脚标编号“对号入座”叠加而成的.这种叠加要统一的结构（总体）坐标系.将总体坐标轴分别用对某单元有求在同一总体坐标系下进行.假如各单元的刚度矩阵是在单元局部坐标下建立的，就必需要把它们转换到式中，（为和『分别是局部坐标系和总体坐标系下的单元结点位移向量；[T]为坐标转换阵,仅与两个坐标系的夹角有关，这样就有[Kr^Tfl^lT]是该单元在总体坐标系下的单元刚度矩阵.以后如不特别强调，总体坐标系下的各种物理参数均不加顶上的横杠.下面就通过简洁的例子来说明如何形成总体刚度矩阵.设有一个简洁的平面结构，选取6个结点，划分为4个单元.单元及结点编号如图3-27所示.每个结点有两个自由度.总体刚度矩阵的组装过程可分为下面几步:这样把位移到列阵中凡是斜支座的结点位移矢量都用局部坐标表示了.将式3T9中第1行左右两边前乘以［几］人a…扁，【即『….♦♦.一【郴••门…［皿班…㈤4一•3-20♦♦■.［刘……44由上式可见凡是边界点的斜支座，在刚度方程中对应于斜支座的位移和载荷向量均可干脆斜支座的局部坐标值，总刚度距阵中的相应行列需作相应的变换.上式的系数矩阵照旧是对称的，而且此方程中7结点位沿X，’轴表示，这样，限定＞’方向的位移就很便利了.实际计算中，并不须要建立结构总的位移方程组后再进坐标变换.而可以在形成单元刚度矩阵和结点载荷之后，就对斜支座点进行坐标变换，把变换后的单元刚度矩阵和结点载荷叠加入总刚度矩阵和总载荷的相应位置，最终叠加形成的也就是方程组3-20,即须要处理的结点，应当在单元计算中完成坐标变换后再叠加，当结构有不同的斜边约束时，都可以这样处理，只不过对不同边上的结点，应按不同的方向余弦矩阵变换就是了.总刚度平衡方程的求解应用有限元法，最终都是归结为解总体刚度平衡方程，它事实上是以总体刚度矩阵为系数矩阵的大型线性代数方程组,通过对结构施加位移边界条件，消退了结构的刚体位移，从而消退总体刚度矩阵的奇异性，解这个线性代数方程组可求出结位移我们已知，总体刚度矩阵具有大型、对称、稀疏、带状分布、正定、主元占优势的特点，稀疏表示刚度矩阵含有大量的零元素，带状表示非零元素集中在主对角线两侧.求解方程组应抓住上述特点，才能提高效率.首先，要为总刚度矩阵选择适当的存储方式，常用的有1整体存储总刚.总刚的全部元素以二维数组形式放在计算机内存中，存储效率最低，适用于小型问题的分析.2等带宽二维存储总刚.总刚的下三角或上三角的带内元素取最大半带宽以二维数组形式存放在计算机内存中.其行数同整体总刚，列数等于最大半带宽.3一维变带宽存储总刚.将总刚下三角实际半带宽内元素逐行存放在一个一维数组内.其次，依据总刚度矩阵的存储方式，选择适当的求解刚度方程的方法线性代数方程组的解法分干脆解法如高斯消去法、三角分解法和迭代解法如高斯-赛德尔迭代、超松弛迭代，本书主要探讨几种常用的直解法.1高斯消去法适合整体存储总刚.由于须要集合完整总刚，内存利用和计算效率都比较低.但高斯消去法原理和程序简洁，作为初学便于理解.2对称消元法利用刚度矩阵对称，有每次消元的子阵均对称的性质，对高斯消去法稍加改进形成.这样就是只需组装总刚的上三角或下三角部分3带消元法将对称消元法进一步改造，使之适合总刚的等带宽二维存储.

（4）因子化法（三角分解）又称Cholesky分解，适合一维变带宽存储总刚.这上方法储效率高，计算速度快，应用较为普遍.此外，还有一种方法，叫做波前法.波前法事实上也是一种改进的高斯消去法.它建立一个称为“波前”的空间，各单元刚度系数依次进入波前.一旦与某自由度有关的全部单元的刚度系数全部装入，便可将相应的变量消去.经过消元的方程的系数随即退出波前，存放在计算机的外存中.这样就可腾出空间装入新的刚度系数.所以，波前法不须要生成完整的总刚，而是边组装边消元，“成熟”一个消去一个.消元完成后，全部系数都已存储在计算机的外存或缓冲区中.回代时将各方程的系数按“先出后入”的依次调入内存求解.由此可见，这种方法是利用计算机充裕的外存资源，以多耗取机时来缓解内存不足的冲突,以便适应较大规模的问题.随着计算机技术的发展，内存资源不断扩大，对具有稀疏、带状性质的有限元刚度方程，这种以时间换取空间的方法得不偿失.另一方面，波前法的阐述和程序设计比较困难，且对多种单元并存的结构运用不便.所以，本书不拟介绍波前法.本书第九章将详细探讨适合整体存储总刚的高斯消去法和适合一维变带宽存储的因子化法以及有关的程序设计问题，以下仅列出这两种方法的梗概.

1、高斯消去法高斯循序消去法的一般公式:对于n阶线性代数方程［KIO=（回，需进行力一1次消元.接受循序消去时，第川次消元以一101内吁】）人apf〃・・=R俨F_¥VK伊*3-21w=l.2•…D想・・J=a+L+2j•彳・次消元后的m行元素作为主元行，人必为主元，对第行元素I＞用）的消元公式为式中LW L…等的上角码,表示该元素是经过第m次消元后得到的结果.同样，可以m把经过m次消元后的系数矩阵和载荷阵分别记为［勺由及”式3—21表时第田次消元是在经m-l次消元的基础上进行的.消元过程中，主元及被消元素的位置可见图3-30a.图中阴影部分已完成消元过程的元素，主元行以下的矩阵为待消部分,在进行第m次时，1-m行元素的消元过程已经完成，其中的元素就是消元最终得到的上三角阵中的元素.m行发下的元素消元过程尚未结束，连同01行元素在内构成一个待消的方阵.消元共需进行n-l次.消元完成后，即可回代求解.我们把消元最终结果记为11=为上三角阵，回代公式可写作il\…3-22i=%12•

3.

2.D一.为一回代过程自后向前进行.当回代求解“时，~已经解得.回代示意图见图3-30b,阴影部分为已求得解答的部分.

2.三角分解法总体刚度平衡方程囚⑹=9中,[K]是对称、正定矩阵，因而可做如下分解[K]^[L][D][L其中3-23]r图7-30高斯消去法豆4Ai Aiw=J J昨⑼是单位上三角矩阵，【匕】二【乙】【01代入整本结构平衡方程记（

①二（力，则［〃/）=（乃即由AiMl・■向下回代.由其中第一个方程解得41,再由其次个方程+A22y2=巴解得L,依此类推可求得{Y}.又由向上回代，可得a=匕，由Q=」得,依此类推可求得⑶.由上述过程可见，三角分解法求解线性代数方程组的关键是对系数矩阵进行三角分解.求解内力由平衡方程组［均@=（力解出位移｛力后，从中分别出各单元的结点位移｛,再通过方程（3-3）.（3-4）和（3-6）等计算各单元的应变（/,、应力（6,和结点力等内力456图7-271按单元局部编号依次形成单元刚度矩阵.图7-27中所示的单元

③，结点的局部编号依次为讣.形成的单元刚度矩阵以子矩阵的形式给出是32将单元结点的局部编号GJ，*换成总体编号，相应的把单元刚度矩阵中的子矩阵的下标也换成总体编号.对下图3-27所示单元

③的刚度矩阵转换成总体编号后为◎73上⑵3将转换后的单元刚度矩阵的各子矩阵，投放到总体刚度矩阵的对应位置上.单元

③的各子矩阵投放后状况如下:

（4）将全部的单元都执行上述的1,2,3步,便可得到总体刚度矩阵，如式（3-9）.其中右上3-9角的上标S表示第S单元所累加上的子矩阵.

（5）从式（3-9）可看出，总体刚度矩阵中的子矩阵AB是单元刚度矩阵的子矩阵转换成总体编号后具有相同的下标J,J的那些子矩阵的累加.总体刚度矩阵第』行的非零子矩阵是由与结点）相联系的那些单元的子矩阵向这行投放所构成的.742结点平衡方程我们首先用结构力学方法建立结点平衡方程.连续介质用有限元法离散以后，取出其中随意一个结点L从环绕［点各单元移置而来的结点载荷为｛印=工｛々｝‘rZ式中，表示对环绕结点7的全部单元求和，环绕结点7的各单元施加于结点了的结点力为£（用，.因此，结点的平衡方程可表示为£（用=（用f（3-10）以［K］代入平衡方程，得到以结点位移表示的结点,的平衡方程，对于每个结点，都可列出平衡方程,于是得到整个结构的平衡方程组如下【灯⑷=（丹式中，［K］为整体刚度矩阵，｛6为全部结点位移组成的向量，为全部结点载荷组成的向量.当然，假如各点的载荷向量也是在单元局部坐标下建立的，在合成以前，也应把它们转换到统一的结构（总体）坐标系下，即⑻J⑷

（2）式中，（户是总体坐标系下的结点载荷向量，为坐标转换阵.位移边界条件在有限元法对结构进行整体分析时，建立了整体刚度矩阵［K］,也得到了结构的刚度平衡方程，即=（丹.结构刚度方程的求解相当于总刚［K］求逆的过程.但是，从数学上看，未经处理的总刚是对称、半正定的奇异矩阵，它的行列式值为零，不能立刻求逆.从物理意义看，在进行整体分析时，结构是处于自由状态，在结点载荷的作用下，结构可以产生随意的刚体位移.所以，在已知结点载荷的条件下，仍不能通过平衡方程惟一地解出结点位移（

①.为了使问题可解，必需对结构加以足够的位移约束，也就是应用位移边界条件.首先要通过施加适当的约束，消退结构的钢体位移，再依据问题要求设定其他已知位移.所以，处理位移边界条件在有限元分析步骤中特别重要.万一n万一方约束的种类包括使某些自由度上位移为零，一，或给定其位移值，”一“•还有给定支承刚度等，本书涉及前两种.处理约束的方法，常用的有删行删列法、分块法、置大数法和置“1”法等，下面分别予以介绍.

1、删行删列法若结构的某些结点位移值为零时（即与刚性支座连接点的位移），则可将总体刚度矩阵中相应的行列、删行删列划掉，然后将矩阵压缩即可求解.这种方法的优点是道理简洁.假如删去的行列很多，则总体刚度矩阵的阶数可大大缩小.通常用人工计算时常接受该方法.若用计算机算题，在程序编制上必带来麻烦，因为刚度矩阵压缩以后，刚度矩阵中各元素的下标必全变更.因而一般计算机算题不太接受.

2.分块法为了理解这个方法，我们把方程［打⑷二（丹分块如下

（3）G R其中，%假设是给定的结点位移；为是无约束的（自由）结点位移.因而日是已知的结点力；玛是未知的结点力.方程（3-11）可以写为即囱-网鬲⑶⑵体/4+［0；=白-313其中，不是奇异的，因而可以解方程3-12得出因4=2⑸⑶14一旦知道了心，就可以由方程3-13求得未知结点力玛.在全部给定的结点自由度都等于零的特别状况下，我们可以删除对应于0的各行和各列即删行删列法，故可把方程简写为［鸟］居4=3T

53.置“1”法由于全部给定的结点位移通常都不能在位移向量5的起先或终了，故分块法的编号方法是很麻烦的.因此，为了引入给定的边界条件%,可以接受下述等价的方法.可以把方程3-12和3-13合在一起写为【[]口-【％】£]41]5]_pl工0]mJkJ-l J3-16在实际计算中，方程3-16所示的过程可以在不重新排列所述方程的状况下用下述分块的方法为进行.步骤1假如把/给定为/,则载荷向量P可以修改为；Pi=g_KQ,j=1,2,…，为卸为结点自由度总数.步骤

（2）除对角线元素以外，使［K］中对应于的行和列为零，而对角线元素为1,即4==0,i=・・・L2,E5=1步骤

（3）在载荷向量中引入规定的$/值，即P对全部规定的结点位移7均应反复运用上述过程（步骤

（1）到

（3））,应当指出，由于这个过程保持了方程的对称性，因此，［K］可以按带状存储，而且几乎不会增加编制程序的工作量.

4.置大数法置大数法的思路是在总体刚度矩阵中，把指定位移所对应的行和列的对角元素乘上一个很大的数，如1°”，此行其他元素保持不变，同时把该行对应的载荷项也相应地用来代替，这里4为指定位移，于是原平衡方程组变为除第7行外，其他各行仍保持原来的平衡特性，而第1个方程式绽开为心+…+io%4+…二］15降4,由于上式中的10乂区比其他项的系数大得多，求和后可略去其小量，则上式变为”W©I勺a=IO M即4=

4.这样就用近似方程组代替原方程组，得到近似满足边界条件的解.当指定位移为零时，只要将对角元素乘上一个大数，而相应的载荷项经证明可以不置零.删行删列法适用于指定零位移点，而置大数法适用于给定位移（包括零位移）.

5.斜支座的处理对于简洁的约束状况（如限定某些结点位移为零或取得给定数值），可以用前述置大数法处理.有的结构在直角坐标系内建立了位移方程组，但在某个斜边上受有法向约束.如图3-28所示正方形固支板，受均布横向载荷，对此，可利用对称性而只计算其1/8,如图中ABC部分，其中AC为固支边，按对称性,8=0AB边上有％，但在BC边上应限定绕BC的转用等于零.为处理此类斜边上的约束，须对斜边上的结点做坐标变换.若结构的总体坐标系为寸户»,为斜支座的局部坐标系（见图3-29）.对于边界结点L须限定〉方向位移M=°,为此，将边界结点了的位移及载荷都变换到局部坐标x〉’轴系.设X轴与斜支座的X’轴夹角为°,逆时针为正，y则依据其次单中坐标转换关系有cos a sin acos asin a［冏-sinacos a=.或写成-sm acosa低）图7-28（即=W3-17与位移关系相同有3-18将上两式带入结构刚度方程有…中/…/%Ka••••••••勺/Q…QW…3-19。