平面桁架程序的计算原理及程序编制优质PPT.ppt

平面桁架程序的计算原理及程序编制优质PPT.ppt

《平面桁架程序的计算原理及程序编制优质PPT.ppt》由会员分享，可在线阅读，更多相关《平面桁架程序的计算原理及程序编制优质PPT.ppt（69页珍藏版）》请在冰点文库上搜索。

节点自由度排序和节点平衡方程相对应。

根据虎克定律，对于任意一根杆，有：

考虑单根杆件在位移下产生的内力。

式中为杆的伸长，为杆的长度，为杆的内力，称为单根杆的刚度（单元刚度阵）。

单元局部刚度阵,用RD来表示单元刚度阵，于是,胡克定律可表达如下,由单根杆件的变形几何关系可得,角是杆件轴线与方向的夹角，由正向逆时针向转至杆轴的角度为正。

进一步有内力和节点位移之间的关系,对于每根杆件（以两端节点编号A和B定出角）应用上述公式，有,应用上述各杆内力和位移关系后，便可建立以位移为未知数的节点平衡方程。

平面桁架结构,节点平衡方程为,注意：

将各杆的内力用位移表示的方程代入上式，有,未知数为节点位移,写成矩阵形式,，即正则方程。

其中,结构刚度矩阵,而外力向量为：

采用高斯消元法求得位移为,利用每根杆的内力-位移关系计算杆内力,5-3最小总势能原理的应用,总势能由两部分组成结构的弹性应变能外力由于结构变位所产生的势能,最小总势能原理与位移法都是以位移为未知数使变形状态预先满足连续条件。

现对上述例题采用最小总势能原理进行求解。

对于整个桁架应变能是所有杆件应变能的叠加，即,由于桁架结构的杆件只能承受拉压力，所以单根杆件的应变能为,把上述值代入应变能表达式，得到,由上式可见，公式中只有位移的二次项，也就是说是位移的一个二次齐次函数，或者说是一个位移的二次型。

（位移的正定二次型，应变能总是正的）,对于该平面桁架，有,现在计算外力势能。

外力产生势能的原因是当节点发生位移时，外力要作功。

所作功的负值便是它们具有势能的改变量，如果取未变形位置外力的势能为零，有,平面桁架结构,1,2,3,4,P,P,将和相加，得到总势能。

由于是位移的一次函数，总势能就成为位移的二次非齐次函数。

根据最小总势能原理，在所有可能的位移状态中，真正发生的位移状态使总势能为最小。

即函数对自变量的偏微商为零，即,式中,现对各位移变量分别取偏微商后，得,注：

值得指出的是刚度矩阵中的系数只与结构本身的几何形态和约束条件有关，而与外载无关。

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位移未知数的确定,对于平面桁架，每个节点有两个自由度，把第个节点的水平位移、垂直位移分别记为，这样结构共有2NU个位移。

位移的方向与坐标轴相同为正，以这些位移作为未知数，并排列成一个列向量，称为结构总位移向量。

第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，系数排列成行向量BT，上式可以写成如下形式,其中,把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）ENDMODULEREALB（4）,DP（120）,R（1200）,C（4,4）INTEGERII（4）,V（120）,NW,IESG,NU,总体数据结构的设计,4NW4IESG2NU01X;

Y11100014100HI;

HL;

EF12100231001310010PX;

PY10,节点总数：

NDISP=2*NU=4,平面桁架结构,1,3,4,2,输入描述平面桁架本体结构和性质的数据结构本体描述数据（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,*）NODEXYDO10I=1,NWREAD（1,*）NODP_X（I）,NODP_Y（I）WRITE（2,*）I,NODP_X（I）,NODP_Y（I）CONTINUEWRITE（2,*）ELEMENTHLHREFDO20NT=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,输入本体结构和性质的数据数据,4NW4IESG2NU01X;

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:

=1IESGDO注：

对所有杆件循环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）USEPDATAINTEGERNT,II（4）,I,I0,J0I=BAR_HL（NT）-1I0=I+I!

NOTETHENUMBERINGFORMULADO10I=1,210II（I）=I0+II=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）!

（*STIF*）,对角元地址数组：

半带宽数组：

刚度矩阵（二维存储）：

刚度矩阵（一维存储）：

存储每一行的对角元素在一维存储数组中R中的位置，存放在数组VBV:

BV+NDISP对角元素在一维数组R中的地址为VBV+I，而这个元素的值为RVBV+I,对角元地址数组：

一维刚度数组：

根据V可以求出：

（1）第i行主元在R中的位置为VBV+i，数值为RVBV+i。

如第3行，BV=0，i=3，V3=6，其值R6=r33。

（2）第i行半带宽（不计对角元）为：

BDWi=VBV+i-VBV+i-1-1,（3）在方阵中处于i行j列的元素rij在R中的地址为：

ix=VBV+i-i+j当i=j但是要求ixVBV+i-1（i=j），否则该元素为零，在R中并没有存放。

（4）第i行第一个非零元素所在列号为i1=VBV+i-1-VBV+i+i+1式中当i=1时VBV定义为零。

（5）第i行对角元地址可用下面的递推公式计算：

Vi=Vi-1+BDWi+1式中BDWi为第i行的半带宽。

半带宽选大的自然语言表示如下：

未知数总数NDISP已知安排数组V的存区半带宽存区清零，准备选大对杆件循环调出I0J0（NT,II）为计算NT号杆对半带宽作准备选出半带宽,根据自然语言表示，半带宽选大的程序段是:

BV:

=1;

FORI:

=0TONDISPDOVBV+I:

=0注：

半带宽选大之前，对V数组清零FORNT:

=1TOIESGDOBEGINI0J0（NT,II）杆位移对号,FORI:

=1TOIBDO!

IB为数组C的下标界偶，即单元位移向量的元素数BEGINFORJ:

=1TOIDO对杆件的行循环和列循环IA:

=III;

JA:

=IIJ;

转换到总刚度矩阵的行、列号IA与JAIFIAJATHENBEGINK:

=IA;

IA:

=JA;

=K;

!

保证总刚度矩阵元素的行号IA大于列号JAEND;

IF（JA0）AND（IA=NDISP）THEN防止出界BEGINIFVBV+IAIA-JATHENVBV+IA:

=IA-JA;

选大END;

END;

半带宽已选出,在此基础上就可以编制DIAGADR（NDISP,V）形成总刚度矩阵的对角元地址过程。

功能：

在给定总刚度矩阵R存放的基址IV，半带宽数组VBV:

BV+NDISP的条件下，求出总刚度矩阵各行的对角元地址，并仍就存放在VBV+1VBV+NDISP中。

SUBROUTINEASSEMB2（IB,II,C,R,NDISP,V,PHASE）USEPDATAINTEGERIB,NDISP,PHASE,I,J,K,IA,JA,II（4）,V（120）REALC（4,4）,R（1200）DO100I=1,IBDO50J=1,IIA=II（I）JA=II（J）!

TAKEOUTDISPNUMBERIF（IA.LT.JA）THENK=IAIA=JAJA=KENDIFIF（JA.GT.0）.AND.（IA.LE.NDISP）THENIF（PHASE.EQ.2）THENK=V（IA）-IA+JA!

（*THEADESSINR*）R（K）=R（K）+C（I,J）!

（*ACCUMULATING*）ELSEIF（V（IA）.LT.（IA-JA）V（IA）=IA-JAENDIFENDIF50CONTINUE100CONTINUERETURNEND!

ASSEMB2,SUBROUTINEDIAGADR（NDISP,V）USEPDATAINTEGERV（120）,IA!

FROMSEMI-BANDWIDTHTODIAGONALADRESSOFGLOBALMATRIX!

IV-BASEADRESSOFGLOBALMATRIX!

VBV.BV+NDISPSEMI-BANDINBUTRETURNWITHDIAGONALADRESS!

NDISP:

GLOBALVARIABLEDO10IA=1,NDISP10V（IA）=V（IA）+V（IA-1）+1RETURNEND!

DIAGADR,总刚度矩阵三角化修改的平方根法运算工作量小且可以处理满阵存储的系数矩阵，可以使用较小的存储空间。

每增加一个右端项时只须增加一次回代求解，而无须重新分解系数矩阵。

具体实施在形成结构总刚度矩阵之后。

NDISPTHEORDEROFMATRIXTOBELDLTED!

VBV.BV+NDISPTHEARRAYOFDIAGONALADREDD,THEMATRIXTOBE!

LDLTEDISINRVBV+1.VBV+NDISP;

T1.BDWTHEWORKINGARRAY;

BDWISTHEMAX-BANDOFRSUBROUTINELDLT1（NDISP,V,R,T）,根据平面桁架的特性进一步展开为：

总体数据结构的设计输入描述平面桁架本体结构和性质的数据输出结构图形计算各杆件的单元刚度矩阵形成结构总刚度矩阵总刚度矩阵三角化输入荷载数据，形成右端项总外力向量回代求解，求出总位移向量打印位移计算各杆内力输出各杆内力,输入荷载数据，形成右端项总外力向量WRITE（2,*）NODEDPXDPYDO70I=1,NUJ=I+I-2DP（J+1）=0.0READ（1,*）DP（J+1）,DP（J+2）WRITE（2,*）DP（J+1）,DP（J+2）CONTINUE!

THELOADONNODESHAVEBEENREADIN,4NW4IESG2NU01X;

PY10,回代求解，求出总位移向量调用SLVEQ1（NDISP,V,R,DP）计算结构的总位移向量。

输出节点位移：

WRITE（2,*）NODEXYUVDO80I=1,NUJ=I+I-2WRITE（2,*）I,NODP_X（I）,NODP_Y（I）,DP（J+1）,DP（J+2）80CONTINUE,根据平面桁架的特性进一步展开为：

总体数据结构的设计输入描述平面桁架本体结构和性质的数据输出结构图形计算各杆件的单元刚度矩阵形成结构总刚度矩阵总刚度矩阵三角化输入荷载数据，形成右端项总外力向量回代求解，求出总位移向量打印位移计算各杆内力输出各杆内力,计算杆件内力：

其中,杆件的伸长量为,表示为矩阵形式,杆单元内力为,计算杆件内力的流程表示如下：

已知总位移向量DP开设局部变量BEGIN调出STIF（NT,L,C,B,RD）;

调出I0J0（NT,II）;

计算杆件内力END;

计算第NT号杆件的内力过程是：

SUBROUTINEINTFORCE（NDISP,NT,DP,T）USEPDATAINTEGERNT,I,J,II（4）REALDP（120）,C（4,4）,B（4）,T,RD,U!

GIVENTTHEELEMENT_NO,RETURNTENSION.DPISGLOBALDISPCALLSTIF（NT,L,C,B,RD）!

GETRDANDBCALLI0J0（NT,II）!

GETELEMENTDISPNUMBERINGVECTORIIT=0.0!

SET0.0BEFOREACCUMULATEDO10I=1,4J=II（I）IF（J.GT.0）.AND.（J.LT.NDISP）THENU=DP（J）!

THERESPECTIVEDISPDEPENDSONISITFIXEDOFNOTELSEU=0.0ENDIFT=T+B（I）*U!

ACCUMULATETHEENLONGATIONOFBAR10CONTINUET=T*RDRETURNEND,计算和输出杆件内力：

WRITE（2,*）BARHLHREFLTENSIONDO90NT=1,IESGL=BAR_TENSION（NT）CALLINTFORCE（NDISP,NT,DP,BAR_TENSION（NT）WRITE（2,*）NT,BAR_HL（NT）,BAR_HR（NT）,BAR_EF（NT）,L,&

BAR_TENSION（NT）90CONTINUE,本章结束！