平面桁架程序的计算原理及程序编制

结构有限元程序设计基本原理——平面桁架程序的计算原理及程序编制5-01矩阵位移法5-02矩阵位移法算例5-03最小总势能原理的应用5-04矩阵位移法的求解步骤5-05结构计算简图的数据结构5-06位移未知数的确定5-07单根杆件的分析5-08结构总刚度矩阵的形成5-09杆件内力的计算5-10能量原理和矩阵位移法平面桁架程序的计算原理5-1矩阵位移法桁架是由离散杆件组成的构架结构，杆件的端点借助于无摩擦的铰连接起来。桁架主要靠各杆中的轴向拉力或压力来传递作用在桁架节点上的荷载，杆件的任何弯曲均忽略不计。用有限元分析桁架时，桁架中的每根杆件都是一个单元，称为杆单元。它是一维单元。不同分类：（1）平面桁架/空间桁架（2）静定桁架/超静定桁架求解方法力法：以力未知数，必须预先满足平衡条件，然后通过连续条件求解未知力；超静定基的选取。位移法：以位移为未知数，各杆件的变形由相连接的节点位移确定（变形协调条件），通过各个节点的平衡方程求出未知位移，再由位移计算出各杆件的内力；各节点的平衡方程也可由最小总势能原理推导。以平面桁架结构分析的程序设计为例，介绍结构分析和程序设计的方法。5-2矩阵位移法算例如图所示平面桁架，各根杆的截面积F相等，材料的弹性模量E相同，在两个单位力P=1的外荷载作用下，用位移法计算节点位移和各杆内力。平面桁架结构1234PP节点总数：NW=4可动节点数：NU=2位移未知数总数：NDISP=2*NU=4注意：节点自由度排序和节点平衡方程相对应。根据虎克定律，对于任意一根杆，有：考虑单根杆件在位移下产生的内力。式中为杆的伸长，为杆的长度，为杆的内力，称为单根杆的刚度(单元刚度阵)。单元局部刚度阵用RD来表示单元刚度阵，于是胡克定律可表达如下由单根杆件的变形几何关系可得

角是杆件轴线与方向的夹角，由正向逆时针向转至杆轴的角度为正。进一步有内力和节点位移之间的关系对于每根杆件（以两端节点编号A和B定出角）应用上述公式，有应用上述各杆内力和位移关系后，便可建立以位移为未知数的节点平衡方程。平面桁架结构节点平衡方程为

注意：节点自由度排序和节点平衡方程相对应。将各杆的内力用位移表示的方程代入上式，有未知数为节点位移写成矩阵形式，即正则方程。其中结构刚度矩阵而外力向量为：采用高斯消元法求得位移为利用每根杆的内力-位移关系计算杆内力5-3最小总势能原理的应用总势能由两部分组成结构的弹性应变能外力由于结构变位所产生的势能最小总势能原理与位移法都是以位移为未知数使变形状态预先满足连续条件。现对上述例题采用最小总势能原理进行求解。对于整个桁架应变能是所有杆件应变能的叠加，即由于桁架结构的杆件只能承受拉压力，所以单根杆件的应变能为把上述值代入应变能表达式，得到由上式可见，公式中只有位移的二次项，也就是说是位移的一个二次齐次函数，或者说是一个位移的二次型。（位移的正定二次型，应变能总是正的）对于该平面桁架，有杆号12345Laaasqrt(2)asqrt(2)a△Lv1u1-u2v2(u1+v1)/sqrt(2)(v2-u2)/sqrt(2)现在计算外力势能。外力产生势能的原因是当节点发生位移时，外力要作功。所作功的负值便是它们具有势能的改变量，如果取未变形位置外力的势能为零，有平面桁架结构1234PP将和相加，得到总势能。由于是位移的一次函数，总势能就成为位移的二次非齐次函数。根据最小总势能原理，在所有可能的位移状态中，真正发生的位移状态使总势能为最小。即函数对自变量的偏微商为零，即式中现对各位移变量分别取偏微商后，得注：值得指出的是刚度矩阵中的系数只与结构本身的几何形态和约束条件有关，而与外载无关。5-4矩阵位移法的求解步骤结构计算简图->（节点、单元编号，建立一个统一的坐标系等）分析节点位移的力学特性->确定位移未知数。建立每根杆件两端位移和内力的关系->（单根杆件的刚度矩阵）根据每根杆件上的上述关系建立结构可动节点的平衡方程->（结构总刚度矩阵）求解平衡方程，得到节点位移。根据求得的位移，利用每根杆件位移与内力关系计算各杆的内力。5-5结构计算简图的数据结构完整而确切描述一个平面桁架结构的数据有三个方面：（1）结构本体描述数据（NW,IESG,NU,X,Y,HL,HR）（2）性质数据（EF）（3）荷载数据（PX,PY）NW为节点总数IESG杆件总数NU可动节点总数X,Y节点坐标HL,HR每根杆件两端节点编号EF性质数据PX,PY外载荷数据5-6位移未知数的确定对于平面桁架，每个节点有两个自由度，把第个节点的水平位移、垂直位移分别记为，，这样结构共有2·NU个位移。位移的方向与坐标轴相同为正，以这些位移作为未知数，并排列成一个列向量，称为结构总位移向量。第i个节点的位移在总位移向量中占i0+1和i0+2位置，而i0=2*i-2。5-7单根杆件的分析根据虎克定律，对于任意一根杆，有：式中为杆的伸长，为杆的长度，为杆的内力，称为单根杆的刚度（单元刚度阵）。一、杆件在局部坐标系中的刚度矩阵

注意⊿L=uB-uA，并且由平衡关系得到可把FA、FB排成列向量{N}，uA，uB，排成列向量{u}，系数阵排成矩阵RD，即设杆端分别得到平行于xoy坐标轴的位移uA，vA，uB，vB；则杆件的伸长量为二、位移的坐标转换把uA，vA，uB，vB

写成列向量{V}，系数排列成行向量{B}T，上式可以写成如下形式其中把⊿L代入上面的公式，得到此时，已用全局坐标系中的位移表达出杆的内力N。由于最终得到的平衡方程都是相对于全局坐标系建立的，上面所计算出的内力与杆轴方向一致。假设在杆端作用平行于杆轴的力FA、FB，则由平衡方程FA和FB在x轴、y轴方向的分量分别为：把、、、排列成列向量，则有由此令得到全局坐标系下的单元刚度矩阵5-8结构总刚度矩阵的形成假设平衡方程形式为展开上式可以写成式中nn=2*NU；P1，P2，…，Pnn为加在节点上的外力。5-9杆件内力的计算对任意一根杆件，如果杆的两端位移已知。这些位移在总位移向量中i0+1，i0+2，j0+1，j0+2的位置上，表示如下进而用下面公式计算结构内力平面桁架计算程序的编制矩阵位移法求解一般步骤：（1）结构计算简图。（节点、单元编号，建立一个统一的坐标系等）（2）分析节点位移的力学特性，确定位移未知数。（3）建立每根杆件两端位移和内力的关系。（单根杆件的刚度矩阵）（4）根据每根杆件上的上述关系建立结构可动节点的平衡方程。（结构总刚度矩阵）（5）求解平衡方程，得到节点位移。（6）根据求得的位移，利用每根杆件位移与内力关系计算各杆的内力。程序流程平面桁架计算程序的总功能表示{依据平面桁架及计算模型的有关数据，进行结构静力计算，并输出计算结果}根据平面桁架的特性进一步展开为：{总体数据结构的设计}{输入描述平面桁架本体结构和性质的数据}{输出结构图形}{计算各杆件的单元刚度矩阵}{形成结构总刚度矩阵}{总刚度矩阵三角化}{输入荷载数据，形成右端项总外力向量}{回代求解，求出总位移向量}{打印位移}{计算各杆内力}{输出各杆内力}根据平面桁架的特性进一步展开为：{总体数据结构的设计}{输入描述平面桁架本体结构和性质的数据}{输出结构图形}{计算各杆件的单元刚度矩阵}{形成结构总刚度矩阵}{总刚度矩阵三角化}{输入荷载数据，形成右端项总外力向量}{回代求解，求出总位移向量}{打印位移}{计算各杆内力}{输出各杆内力}{总体数据结构的设计}平面桁架程序的总体数据结构设计，包括全局量（标识符）的存储安排与说明。一个平面桁架结构的数据有三个方面：结构本体描述数据（NW,IESG,NU,X,Y,HL,HR）性质数据（EF）荷载数据（PX,PY）NW为节点总数IESG杆件总数NU可动节点总数X,Y节点坐标HL,HR每根杆件两端节点编号EF性质数据PX,PY外载荷数据MODULEPDATAREAL,PUBLIC::BAR_HL(100),BAR_HR(100),&BAR_TENSION(100),BAR_EF(100)REAL,PUBLIC::NODP_X(80),NODP_Y(80)ENDMODULE

REALB(4),DP(120),R(1200),C(4,4)INTEGERII(4),V(120),NW,IESG,NU{总体数据结构的设计}4‘NW4‘IESG2‘NU01‘X;Y11100014100‘HI;HL;EF12100231001310010‘PX;PY10节点总数：NW=4可动节点数：NU=2位移未知数总数：NDISP=2*NU=4平面桁架结构1342{输入描述平面桁架本体结构和性质的数据}结构本体描述数据（NW,IESG,NU,X,Y,HL,HR）性质数据（EF）READ(1,*)NWREAD(1,*)IESGREAD(1,*)NUNDISP=NU+NUWRITE(2,*)'NW=',NW,'IESG=',IESG,'NU=',NUWRITE(2,*)'NODEXY‘DO10I=1,NWREAD(1,*)NODP_X(I),NODP_Y(I)WRITE(2,*)I,NODP_X(I),NODP_Y(I)CONTINUEWRITE(2,*)'ELEMENTHLHREF‘DO20NT=1,IESGREAD(1,*)BAR_HL(NT),BAR_HR(NT),BAR_EF(NT)WRITE(2,*)BAR_HL(NT),BAR_HR(NT),BAR_EF(NT)20CONTINUE!LISTDATA,NOWTOORGANIZETHEGLOBALMATRIX{输入本体结构和性质的数据数据}4‘NW4‘IESG2‘NU01‘X;Y11100014100‘HI;HL;EF12100231001310010‘PX;PY10根据平面桁架的特性进一步展开为：{总体数据结构的设计}{输入描述平面桁架本体结构和性质的数据}{输出结构图形}{计算各杆件的单元刚度矩阵}{形成结构总刚度矩阵}{总刚度矩阵三角化}{输入荷载数据，形成右端项总外力向量}{回代求解，求出总位移向量}{打印位移}{计算各杆内力}{输出各杆内力}{计算各杆件的单元刚度矩阵}相应位移向量SUBROUTINESTIF(NT,L,C,B,RD)USEPDATAINTEGERNT,I,JREALL,RD,C(4,4),B(4),CA,SA,DX,DYI=BAR_HL(NT)J=BAR_HR(NT) !LEFTANDRIGHTNODENUMBERDX=NODP_X(J)-NODP_X(I)DY=NODP_Y(J)-NODP_Y(I)L=SQRT(DX*DX+DY*DY)CA=DX/L;SA=DY/L !(*DIRECTINCOSINE*)B(1)=-CAB(2)=-SAB(3)=CAB(4)=SA!(*THEDISPLACEMENT->DEFORMATIONTRANSFORMATIONMATRIX*)RD=BAR_EF(NT)/L!(*BARSTIFFNESS*)DO10I=1,4DO10J=1,4C(I,J)=RD*B(I)*B(J)CONTINUE!(*CISTHESTIFFNESSMATRIXINGLOBALCOORDINATE*)RETURNEND !(*STIF*){计算各杆件的单元刚度矩阵}根据平面桁架的特性进一步展开为：{总体数据结构的设计}{输入描述平面桁架本体结构和性质的数据}{输出结构图形}{计算各杆件的单元刚度矩阵}{形成结构总刚度矩阵}{总刚度矩阵三角化}{输入荷载数据，形成右端项总外力向量}{回代求解，求出总位移向量}{打印位移}{计算各杆内力}{输出各杆内力}{形成结构总刚度矩阵}{对总刚度矩阵R清零}{FORNT:=1~IESGDO

注：对所有杆件循环BEGIN调出STIF(NT,T,C,B,RD)；计算杆件在全局坐标系下单元刚度矩阵C}调出I0J0(NT,II)；计算对号入坐信息数组II}调出ASSEMB2(IB,C,R,NDISP,BV,V)；注：结构总刚度矩阵R已形成END;}

!确定这根杆件两端点（HL端和HR端）位移在总位移向量中的编号SUBROUTINEI0J0(NT,II)USEPDATA

INTEGERNT,II(4),I,I0,J0I=BAR_HL(NT)-1I0=I+I!NOTETHENUMBERINGFORMULADO10I=1,210

II(I)=I0+I

I=BAR_HR(NT)-1J0=I+IDO20I=3,4II(I)=J0+I-2RETURNEND!(*I0J0*)

ui在总位移向量中是第2*（i-1）+1号uj在总位移向量中是第2*（j-1）+1号SUBROUTINESTIF(NT,L,C,B,RD)USEPDATAINTEGERNT,I,JREALL,RD,C(4,4),B(4),CA,SA,DX,DYI=BAR_HL(NT)J=BAR_HR(NT) !LEFTANDRIGHTNODENUMBERDX=NODP_X(J)-NODP_X(I)DY=NODP_Y(J)-NODP_Y(I)L=SQRT(DX*DX+DY*DY)CA=DX/L;SA=DY/L !(*DIRECTINCOSINE*)B(1)=-CAB(2)=-SAB(3)=CAB(4)=SA!(*THEDISPLACEMENT->DEFORMATIONTRANSFORMATIONMATRIX*)RD=BAR_EF(NT)/L!(*BARSTIFFNESS*)DO10I=1,4DO10J=1,4C(I,J)=RD*B(I)*B(J)CONTINUE!(*CISTHESTIFFNESSMATRIXINGLOBALCOORDINATE*)RETURNEND !(*STIF*)13581011182225330112106327对角元地址数组：半带宽数组：3.2

-1.04.3

0.82.0

4.00.18.8

4.14.0

2.1

-0.81.20.30002.8

1.5003.1

0.703.7

-0.200-1.20-2.14.17.2刚度矩阵（二维存储）：3.24.30.82.04.00.18.84.14.02.11.20.300002.81.5003.10.703.70004.17.2刚度矩阵（一维存储）：存储每一行的对角元素在一维存储数组中R中的位置，存放在数组

V[BV:BV+NDISP]对角元素在一维数组R中的地址为V[BV+I]，而这个元素的值为

R[V[BV+I]]1358101119232634对角元地址数组：3.24.30.82.04.00.18.84.14.02.11.20.300002.81.5003.10.703.70004.17.2一维刚度数组：根据V可以求出：（1）第i行主元在R中的位置为V[BV+i]，数值为R[V[BV+i]]。如第3行，BV=0，i=3，V[3]=6，其值R[6]=r33。（2）第i行半带宽（不计对角元）为：BDW[i]=V[BV+i]-V[BV+i-1]-1（3）在方阵中处于i行j列的元素rij在R中的地址为：ix=V[BV+i]-i+j当i>=j但是要求ix>V[BV+i-1]（i>=j），否则该元素为零，在R中并没有存放。（4）第i行第一个非零元素所在列号为i1=V[BV+i-1]-V[BV+i]+i+1式中当i=1时V[BV]定义为零。（5）第i行对角元地址可用下面的递推公式计算：V[i]=V[i-1]+BDW[i]+1式中BDW[i]为第i行的半带宽。半带宽选大的自然语言表示如下：{未知数总数NDISP已知}{安排数组V的存区}{半带宽存区清零，准备选大}{对杆件循环}{调出I0J0（NT,II）}{为计算NT号杆对半带宽作准备}{选出半带宽}根据自然语言表示，半带宽选大的程序段是:BV:=1;FORI:=0TONDISPDOV[BV+I]:=0注：半带宽选大之前，对V数组清零FORNT:=1TOIESGDOBEGINI0J0(NT,II)杆位移对号FORI:=1TOIBDO!IB为数组C的下标界偶，即单元位移向量的元素数BEGINFORJ:=1TOIDO对杆件的行循环和列循环IA:=II[I];JA:=II[J];转换到总刚度矩阵的行、列号IA与JAIFIA<JATHENBEGINK:=IA;IA:=JA;JA:=K;!保证总刚度矩阵元素的行号IA大于列号JAEND;IF(JA>0)AND(IA<=NDISP)THEN防止出界BEGINIFV[BV+IA]<IA-JATHENV[BV+IA]:=IA-JA;!选大END;END;END;{半带宽已选出}在此基础上就可以编制DIAGADR(NDISP,V)形成总刚度矩阵的对角元地址过程。功能：在给定总刚度矩阵R存放的基址IV，半带宽数组V[BV:BV+NDISP]的条件下，求出总刚度矩阵各行的对角元地址，并仍就存放在V[BV+1]~V[BV+NDISP]中。 SUBROUTINEASSEMB2(IB,II,C,R,NDISP,V,PHASE) USEPDATA INTEGERIB,NDISP,PHASE,I,J,K,IA,JA,II(4),V(120) REALC(4,4),R(1200) DO100I=1,IB DO50J=1,I IA=II(I) JA=II(J) !TAKEOUTDISPNUMBER IF(IA.LT.JA)THEN K=IA IA=JA JA=K ENDIF IF((JA.GT.0).AND.(IA.LE.NDISP))THEN IF(PHASE.EQ.2)THEN K=V(IA)-IA+JA!(*THEADESSINR*) R(K)=R(K)+C(I,J)!(*ACCUMULATING*) ELSE IF(V(IA).LT.(IA-JA))V(IA)=IA-JA ENDIF ENDIF50CONTINUE100CONTINUE RETURN END!ASSEMB2 SUBROUTINEDIAGADR(NDISP,V)USEPDATAINTEGERV(120),IA!FROMSEMI-BANDWIDTHTODIAGONALADRESSOFGLOBALMATRIX!IV--BASEADRESSOFGLOBALMATRIX!V[BV..BV+NDISP]SEMI-BANDINBUTRETURNWITHDIAGONALADRESS!NDISP:GLOBALVARIABLEDO10IA=1,NDISP10V(IA)=V(IA)+V(IA-1)+1RETURNEND!DIAGADR{总刚度矩阵三角化}修改的平方根法运算工作量小且可以处理满阵存储的系数矩阵，可以使用较小的存储空间。每增加一个右端项时只须增加一次回代求解，而无须重新分解系数矩阵。具体实施在形成结构总刚度矩阵之后。!NDISP—THEORDEROFMATRIXTOBELDLTED!V[BV..BV+NDISP]—THEARRAYOFDIAGONALADREDD,THEMATRIXTOBE!LDLTEDISINR[V[BV]+1..V[BV+NDISP]];!T[1..BDW]—THEWORKINGARRAY;BDWISTHEMAX-BANDOFRSUBROUTINELDLT1(NDISP,V,R,T)根据平面桁架的特性进一步展开为：{总体数据结构的设计}{输入描述平面桁架本体结构和性质的数据}{输出结构图形}{计算各杆件的单元刚度矩阵}{形成结构总刚度矩阵}{总刚度矩阵三角化}{输入荷载数据，形成右端项总外力向量}{回代求解，求出总位移向量}{打印位移}{计算各杆内力}{输出各杆内力}{输入荷载数据，形成右端项总外力向量}WRITE(2,*)"NODEDPXDPY"DO70I=1,NU J=I+I-2 DP(J+1)=0.0 READ(1,*)DP(J+1),DP(J+2) WRITE(2,*)DP(J+1),DP(J+2)CONTINUE!THELOADONNODESHAVEBEENREADIN4‘NW4‘IESG2‘NU01‘X;Y11100014100‘HI;HL;EF12100231001310010‘PX;PY10{回代求解，求出总位移向量}调用SLVEQ1(NDISP,V,R,DP)计算结构的总位移向量。

输出节点位移：WRITE(2,*) 'NODEXY