结构有限元程序设计基本原理——平面桁架程序的计算原理及程序编制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=41122uvWuv注意:节点自由度排序和节点平衡方程相对应。根据虎克定律,对于任意一根杆,有:EFNLL考虑单根杆件在位移下产生的内力。式中为杆的伸长,为杆的长度,为杆的内力,称为单根杆的刚度(单元刚度阵)。LLNEFL单元局部刚度阵用RD来表示单元刚度阵,于是EFRDL胡克定律可表达如下NRDL由单根杆件的变形几何关系可得cossin()cos()sinBABALuvuuvv角是杆件轴线与方向的夹角,由正向逆时针向转至杆轴的角度为正。xxAAuvBBuvABxy进一步有内力和节点位移之间的关系()cos()sinBABAEFNuuvvL对于每根杆件(以两端节点编号A和B定出角)应用上述公式,有11vNEFa122uuNEFa23vNEFa1142uvNEFa2152vuNEFa应用上述各杆内力和位移关系后,便可建立以位移为未知数的节点平衡方程。平面桁架结构节点平衡方程为024cos450NNP014cos450NN025cos450NNP035cos450NN注意:节点自由度排序和节点平衡方程相对应。将各杆的内力用位移表示的方程代入上式,有11211102222EFuvuPa11111002222EFuvPa12211012222EFuuvPa22110012222EFuvPa未知数为节点位移写成矩阵形式RWP,即正则方程。其中11121314212223243132333441424344111102222111002222/111012222110012222rrrrrrrrREFarrrrrrrr结构刚度矩阵而外力向量为:123400PPPPPPP采用高斯消元法求得位移为13.824/uaEF11.000/vaEF21uu21vv利用每根杆的内力-位移关系计算杆内力11N20N31N41.412N51.412N5-3最小总势能原理的应用•总势能由两部分组成结构的弹性应变能外力由于结构变位所产生的势能UV最小总势能原理与位移法都是以位移为未知数使变形状态预先满足连续条件。现对上述例题采用最小总势能原理进行求解。对于整个桁架应变能是所有杆件应变能的叠加,即5211()()/2kkkkUEFLL由于桁架结构的杆件只能承受拉压力,所以单根杆件的应变能为2()/2EFULL把上述值代入应变能表达式,得到2222211221222()()/22()/222EFUvuuvuvvua由上式可见,公式中只有位移的二次项,也就是说是位移的一个二次齐次函数,或者说是一个位移的二次型。(位移的正定二次型,应变能总是正的)对于该平面桁架,有杆号12345Laaasqrt(2)asqrt(2)a△Lv1u1-u2v2(u1+v1)/sqrt(2)(v2-u2)/sqrt(2)现在计算外力势能。外力产生势能的原因是当节点发生位移时,外力要作功。所作功的负值便是它们具有势能的改变量,如果取未变形位置外力的势能为零,有P12VPuPu平面桁架结构1234PP将和相加,得到总势能。由于是位移的一次函数,总势能就成为位移的二次非齐次函数。UVV根据最小总势能原理,在所有可能的位移状态中,真正发生的位移状态使总势能为最小。即函数对自变量的偏微商为零,即()0(1,2,3,4)iUViw11213242wuwvWwuwv式中现对各位移变量分别取偏微商后,得11211102222EFuvuPa11111002222EFuvPa12211012222EFuuvPa22110012222EFuvPa注:值得指出的是刚度矩阵中的系数只与结构本身的几何形态和约束条件有关,而与外载无关。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个位移。位移的方向与坐标轴相同为正,以这些位移作为未知数,并排列成一个列向量,称为结构总位移向量。iiuiv1112232iinunuNDISPuwvwuwvWuvuvw第i个节点的位移在总位移向量中占i0+1和i0+2位置,而i0=2*i-2。5-7单根杆件的分析根据虎克定律,对于任意一根杆,有:EFNLL式中为杆的伸长,为杆的长度,为杆的内力,称为单根杆的刚度(单元刚度阵)。LLN/EFL一、杆件在局部坐标系中的刚度矩阵注意⊿L=uB-uA,并且由平衡关系得到()ABABEFFFNuuL可把FA、FB排成列向量{N},uA,uB,排成列向量{u},系数阵排成矩阵RD,即AABBFuNuFuEFEFLLRDEFEFLLRDuN设杆端分别得到平行于xoy坐标轴的位移uA,vA,uB,vB;则杆件的伸长量为cossin()cos()sinBABALuvuuvv二、位移的坐标转换把uA,vA,uB,vB写成列向量{V},系数排列成行向量{B}T,上式可以写成如下形式TLBV其中cossincossinBAABBuvVuv把⊿L代入上面的公式,得到cossincossinAABBuvEFNRDBVuLv此时,已用全局坐标系中的位移表达出杆的内力N。由于最终得到的平衡方程都是相对于全局坐标系建立的,上面所计算出的内力与杆轴方向一致。假设在杆端作用平行于杆轴的力FA、FB,则由平衡方程ABFFNFA和FB在x轴、y轴方向的分量分别为:coscosAxAFFNsinsinAyAFFNcoscosBxBFFNsinsinByBFFN把、、、排列成列向量,则有AxFAyFBxFByFPcossincossinAxAyBxByFFPNFF由此TTPBNBRDBV令TCBRDB得到PCV全局坐标系下的单元刚度矩阵5-8结构总刚度矩阵的形成假设平衡方程形式为RWP展开上式可以写成11112211211222221122nnnnnnnnnnnnnnnnnnnnrwrwrwPrwrwrwPrwrwrwP式中nn=2*NU;P1,P2,…,Pnn为加在节点上的外力。5-9杆件内力的计算对任意一根杆件,如果杆的两端位移已知。这些位移在总位移向量中i0+1,i0+2,j0+1,j0+2的位置上,表示如下01020102iijjwwVww进而用下面公式计算结构内力NRDBV平面桁架计算程序的编制矩阵位移法求解一般步骤:(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),BA