安徽建筑工程学院计算结构力学.ppt

第五章结构刚度矩阵与荷载向量,计算结构力学,5-1概述,以图示框架结构为例，设有n个未知量:,相应的结点荷载向量为:,则结构刚度方程可写为:,K-称为结构刚度矩阵(或称总刚度矩阵),本章讨论K的形成及程序设计，以及P的形成。,5-2应用能量原理形成结构刚度矩阵,其中:,结构在外荷载作用下的总势能可以写成：,NE是单元数,单元结点位移,整体坐标系中的单元刚度矩阵,NF单元自由度数,N结构未知量总数,C0或1。指明单元结点位移向量是由结构结点位移向量中的哪几个分量所组成,令:,(5)式反映了结构的离散过程，实际上表明了结构离散化后的变形协调条件。,将(5)式代入(3)式：,结构总势能为：,于是：,K是由k经过C变换后装配而成。,由于势能驻值原理等价于平衡方程，故装配总刚的有限元集合过程遵循平衡条件。,CNFxN是单元定位向量的增广写法：一个单元对应一个Ce，且有：,讨论,例1：,解：给单元结点编号，并写出各单元的定位向量：,1,2,由(5)式可写出各单元的结点位移向量与结构结点位移向量的关系式：,由此可得到各单元的C,如对第单元，可写出:,即Cij的行号与MWe的行号一致，把的序号作为C的列号j，便可由MWe得到Ce。,5-3按单元定位向量装配结构刚度矩阵,MW处理了约束，以及主从关系，无效未知量等特殊结点信息，也是C矩阵的实用(增广)写法。,MW是按单元结点编号顺序由结点的结构未知量编号顺序所组成的向量(列阵)。,单元定位向量可方便地指出单元的各个未知量在结构总体未知量中的对应位置(总体序号)。由此也就可以确定单元刚度矩阵中的元素在结构刚阵中的位置。,例2：求图示连续梁的结构刚度矩阵。,注意以下写法：,得到结构总刚度矩阵为：,主系数与副系数:相关未知量：相关结点：相关单元：,结构刚度矩阵的组成规律专有名词,凡未知量i的相关结点所在单元称为末知量i的相关单元。,与未知量i在同一单元的未知量叫做未知量i的相关未知量。若i，j相关，则Kij0若i，j不相关，则Kij=0,未知量i的相关未知量所在结点称为未知量i的相关结点。,5-4形成结构刚度矩阵的直接刚度法,不是列向量乘行向量,也不是向量的内积(点积),ai与bi不进行任何运算。,2、由单元定位向量的指标并积形成下标矩阵,如果将某个单元的定位向量代入上式，由上节中的例题可明显看出H中的元素就表示这个单元的刚度系数在结构刚度矩阵K中的下标。,解：如图,各单元的定位向量为：MW=12TMW=23T,例3：求图示连续梁的单元刚度矩阵在结构刚度矩阵中的下标矩阵。,根据并积定义：,清楚地表明了各单刚系数在总刚中的位置，参考前节例题,例4：求图示刚架中第单元的刚度系数对结构刚度矩阵的贡献。,式中圆括号内的元素就是第单元刚度系数在K中的下标。式中含零的元素说明单刚中此元素经C夹乘后为零,参考(5-2-9)式,不须叠加,只有H中的元素与结构K中下标一致时才进行叠加。这样便可根据单元定位向量的并积作为结构刚度矩阵K的下标直接来装配结构刚度矩阵。这就是直接刚度法。,上式并积的进一步说明：单元结点位移的序号为123456(I,J)单元的定位向量为000102(L,K)则意味着：,表示叠加到结构刚阵中去,由此可看出：由MW的并积形成下标矩阵，完全确立了单刚k中的元素在总刚K中位置，从而由数学的角度说明了用MW装配K的过程。,上述过程的FORTRAN程序模块可写成：L=MW(I)K=MW(J)ZK(L,K)=ZK(L,K)+DK(I,J),3、形成K的程序设计框图,本章新的变量和数组:NAI：EA或EI分组数(截面特性分组数)；DK(I,J)：单元刚度矩阵，其中I,J(16):单刚的行列号；ZK(L,K)：结构刚度矩阵,其中L,K(1N):总刚的行列号。,程序设计框图,4、形成结构刚度K的源程序设计,CTHEPROGRAMOFKJEXAMDIMENSIONJW(3,20),JTX(4,20),JH(2,20),*MW(6),JMH(20)REAL*8CX(20),SY(20),SL(20),EA(5),*X(20),Y(20),XSA(20,7),ZK(50,50),El(5)WRITE(*,*)FINDINGTHEMWOFEI.EMENTSOPEN(1,FILE=KJE.DAT),数据文件名为KJE.DAT,READ(1,*)NE,NJ,NJT,NAIREAD(1,*)(JH(I,J),I=1,2),J=1,NE)READ(1,*)(JTX(I,J),I=1,4),J=1,NJT)READ(1,*)(JMH(I),I=1,NE)READ(l,*)(EA(I),I=1,NAI)READ(l,*)(EI(I),I=1,NAI)READ(1,*)(X(I),I=1,NJ)READ(l,*)(Y(l),I=1,NJ)CALLQJW(NJ,NJT,JTX,JW,N)DO10M=1,NECALLQMW(M,NE,NJ,JH,JW,MW)WRITE(*,100)M,(MW(I),I=1,6),10,CALLDCH(NE,NJ,JH,X,Y,SL,CX,SY)CALLQXS(NE,NAI,JH,JMH,SL,CX,SY,EA,*EI,XSA)CALLKJX1(NE,N,XSA,NJ,JH,JW,ZK)WRITE(*,*)THISISSTRUCTURAL*MATRAIXZKWRITE(*,50)(ZK(I,J),I=I,N),J=I,N)100FORMAT(1X,ELEMT,NO=I5,/6X,MW*=,6I5)50FORMAT(9F9.1)STOPEND,SUBROUTINEKJX1(NE,N,XSA,NJ,JH,JW,ZK)DIMENSIONJH(2,NE),JW(3,NJ),MW(6)REAL*8DK(66),ZK(50,50),XSA(NE,7),*XS(7),ZK(50,50)CALLZERO2(ZK,50,50)DO20M=1,NEDO25I=1,7XS(I)=XSA(M,I)25CONTINUE,CALLDKX(XS,DK)WRITE(*,*)THISISELEMENT*MATRAIXWRHE(*,15)M15FORMAT(1X,ELEMENTNO.=,I5)WRHE(*,11)(DK(I,J),l=l,6),J=l,6)11FORMAT(1X,6FI0.2)CALLQMW(M,NE,NJ,JH,JW,MW)DO30I=1,6L=MW(I)IF(L.LE.0)GOTO30,DO40J=1,6K=MW(J)IF(K.LE.0)GOTO40ZK(L,K)=ZK(L,K)+DK(I,J)40CONTINUE30CONTINUE20CONTLNUERETURNEND,本程序是在“形成结构单元定位向量”的程序基础上扩充而成；本程序是用滚雪球的方法，扩充部分用红色表示本程序的总刚KNxN采用满阵存储。,数据文件KJE.DAT共8句NE与NJ分别为单元总数与结点总数NJT是特殊结点数NAI为截面特性分组数JH(2,NE)为单元两端结点编号数组JMH(NE)为截面特性分组号数组JTX(4,NJT)为特殊结点约束信息数组EA(NAI)为各组单元的EAEI(NAI)为各组单元的ElX(NJ)为各结点的X坐标Y(NJ)为各结点的Y坐标,5、形成结构刚度矩阵K的算例,目的：通过程序运行的中间步骤来说明直接刚度法的运算过程，要注意程序运行过程中所打印的中间结果：如单元定位向量、单元刚度矩阵等，通过对单元循环而直接形成结构刚度矩阵。,例5：形成图示三层框架的结构刚度矩阵。已知E=0.26x108KN/m2,梁和柱为矩形截面,尺寸bh分别为0.3x0.5m和0.3x0.6m。,解:先进行结点编号和单元编号计算梁、柱截面的有关物理量,梁:A=0.3x0.5=0.15m2EA=390000KNI=0.3x0.53/12=0.003l25m4EI=8125KNm2,柱：A=0.3x0.6=0.18m2EA=468000KNI=0.3x0.63/12=0.0054m4EI=14040KNm2,输人数据文件KJE.DAT内容9,8,2,23,1,4,2,5,3,6,4,7,5,8,6,1,2,3,4,5,67,1,1,1,8,1,1,11,1,1,1,1,1,2,2,2468000,39000014040,81250,6,0,6,0,6,0,610,10,7,7,4,4,0,0,部分结果为：K11=71240K17=-1354K22=156451K18=38176,5-5K的特性,1）KNxN，N由JW数组确定；2）K=KT，K是对称数组；3）若未知量i,j不相关，则Kij=0；4）带状特性，即非零元素分布在主对角线附近；5）K为稀疏阵，非零元素一般只占10左右；6）考虑约束处理后,K正定；7）主对角线上Kij0。,5-6等效结点力,对于K=P上式中P与一一对应在实际工程中，有：单元荷载、广义荷载（支座移动、变温），故必须研究上述情况下这些单元的等效荷载，推导等效荷载时采用虚功的不变性原理，即单元荷载与等效结点力在任何虚位移上所作的虚功相等，这实际上就是静力等效原则。静力等效移置只对所讨论单元的应力分布有影响，而不影响结构整体的应力分布。,1、用形函数求等效结点力,如图单元荷载Q，单元等效结点力FE，设单元发生虚位移v*,各结点虚位移为*，由虚功等效，得：,代入上式，并取特例x=l/2，考虑到Q力方向向下，故将Q代换成Q，便有：,故有：,上式正好是固端反力，但符号相反，由此可推出FE的另一计算公式：,式中F0为梁的固端反力，当q(x)为分布荷载时，由(2)式可推出,公式(3)是杆单元的特例，固端反力可由结构力学求出，列在附表，并编入程序。公式(3)仅对两端刚性连接的杆件，对于铰结点，需重推导形函数。,2、用固端反力求等效结点力FE,附表固定端梁的固端反力F0,附表(续)固定端梁的固端反力F0,关于附表，应注意此时z轴向里，Mzi逆时针为正；广义荷载列在7、8二栏；FE是在局部系里计算的，在由FE求P时，仍需进行坐标变换，由于,所以变换公式为：,5-7按单元定位向量装配荷载向量,列入P中的荷载有：1.直接作用在结点上的结点荷载，可直接装配；2.单元等效结点荷载FE。,将FE按单元的MW装配P是很简单的，因为FE的分量个数与MW的元素个数相等，所以在每个FE中的分量对应地在MW中有一个未知量编号，因此可以正确地叠加到P中去，若MW中对应的元素是0或-1，则相应的FE中分量就不必叠加。,形成右端项P的程序设计：可在形成结构刚阵的主程序中直接调用形成结点荷载的子程序SUB.YDXNLX(1,1)=NPJNLX(2,1)=NPMCALLYDX(1,N,NE,NJ,NJZ,1,NPJ,NPM,JH,NLX,JW,JZH,SL,CX,SY,PJZ,PMZ,FLZ,P)其中，NPJ：受载结点数NPM：受载单元数NLX(2，NL)：荷载组信息NL：荷载组数，每组有两个信息，即：受载结点数NPJ和受载单元数NPM,*,子程序SUB.YDX见教材P120。子程序中有两个输入语句，可加在数据文KJE.DAT的后面，作为第9句、第10句。第9句见P120，输入结点荷载信息：结点荷载信息数组PJZ(NPJ,2)，其中NPJ为受载结点数，第1列即PJZ(I,1)需输入两个信息JD.x式中，JD：荷载作用的结点号x：荷载作用的方向，在平面刚架情况下x=1