平面桁架的有限元法

1、平面桁架的有限元法一、建立有限元模型 平面桁架是现成的有限元模型， 只有节点力二、单元分析 位移假设，几何方程，物理方程， 单元刚度矩阵三、总刚度矩阵、总载荷列阵、引入支撑 每个节点建立平衡方程 集成技术（叠加法、对号入座）四、求解、结果分析，后处理一、建立有限元模型平面桁架是现成的有限元模型， 只有节点力桁架原型=有限元模型yxP二、单元分析 平面二力杆单元注意：一个单元的节点数目（ 2）； 一个节点的自由度数目（2）； 有整体坐标系和局部坐标系之分xoyxyuiiujjijbvv二、单元分析（先在局部坐标系中进行） 问题：当节点位移一定，相应的节点力为多大？节点力与节点位移之间的比例常数为

2、刚度系数目标：建立单元刚度矩阵。方法：结构力学的位移法加虚功原理， 皆用节点位移表示。xoyxyuiiujjijbvv1、位移假设： 设各点的位移，即位移场为xaavxaau4321421,aaa其中为待定系数，由节点坐标和节点位移确定。将位移写成矩阵形式：4321100001aaaaxxvuaHfsxoyxyuiiujjijbvvvu在节点上满足：baavbaauavauijii432131写成矩阵形式：aAbexaavxaau4321bxxji , 0 xoyxyuiiujjijbvvvu432110000101000001aaaabbvuvujjii解线性代数方程组，得aAbe1ebAa

3、aHfs代入得141444212ebsAHf144212efNf432110000101000001aaaabbvuvujjii4321100001aaaaxxvu1)(bsAHxNefNf节点位移与单元内位移的关系其中为形函数。用Machematica推导形函数：e在编程中只关心节点位移的系数矩阵xoyxyuiiujjijbvvvuxxHs1000012、应变（从位移求应变） 只有轴向应变3、应力（从应变求应力，本构方程）0eLvuBxuyvxu0 xuEEivuivuxaavxaau43214、内力（从应力求内力）xuAEANiuyvxuAEAENSiie000二力杆只有轴向力0yv用Ma

4、chematica推导应力矩阵：(EA的两个字母连写作参数名)eD用内力表示的物理方程yvxuEAEAN000EAEADe000eLBxuyvxu0NSeeLeeBDSefNf4、单元刚度矩阵（用虚功原理） 给节点一组虚位移，单元内即产生相应的虚应变，实际节点力（轴向力、横向力、弯矩）在虚位移上所做的功等于实际应力（横向剪应力可忽略不计）在虚应变上产生的应变能。用Machematica推导单刚矩阵Ke：beLeTLTeleTeTexBDBlSR0ddeeLBeLeeBDS00jieNNR144414eeeKRbLeTLexBDBK042222444d用Machematica推导，所得结果以矩阵

5、形式显示三、总刚度矩阵、总载荷列阵、引入支撑1、坐标变换 前面的推导是在局部坐标系中进行的，而总刚度矩阵、总载荷列阵、引入支撑等要在整体坐标中进行，因此，存在坐标变换问题。xoyxycossinsincosiiiiiiwuwwuuiiiiwuwucossinsincosiit三、总刚度矩阵、总载荷列阵、引入支撑1、坐标变换iitjijitt00eeTttT00TTTTtt11eeRTReekReekReekRTeeTkkTTTkTkxoyxyuiiujjijbvvvucossinsincost三、总刚度矩阵、总载荷列阵、引入支撑2、在总结构和整体坐标中讨论问题，所需要的信息：（1）节点号；（2

6、）单元长度；（2）单元倾角(单元起点到终点的有向线段与整体横轴之夹角）；（4）节点力所在的节点号和作用方向；（5）支撑所在的节点号和作用方向。3、总节点位移向量的排列规律 按节点号和节点自由度的顺序（从小到大）排列。三、总刚度矩阵、总载荷列阵、引入支撑建立节点号数组（用Mathematica)：jm=Table0, 0, i, ne; “jm数组名，ne 总单元数;”jm=，i节点号， j节点号，;建立单元长度数组gc=Table0, i, ne; “gc数组名，ne 总单元数;”gc=L1, L2, Lne;建立单元方向角数组gj=Table0, i, ne; “gj数组名，ne 总单元数;

7、”gj=af1, af2, afne;3、总节点位移向量的排列规律jm=1, 2，2, 3，i节点号， j节点号，;jjiievuvu单元节点位移排列规律：rr =2(k-1)+kkk =1, 2 (单元节点);kk =1, 2, (节点自由度)k=1k=2kk=1kk=2总节点位移向量d的排列规律(从小到大排序）：ss=2(jme, k-1)+kk;e=单元号rreessTd)(对每个单元进行节点循环和自由度循环，从单元地址到结构地址搬家是单元到结构的映射(map)3、总节点位移向量的排列规律单元节点位移号：rr =2(k-1)+kk；k =1, 2 (单元节点)； kk =1, 2(节点自

8、由度)总节点位移号：ss=2(jme, k-1)+kk；例如，已解得结构的总节点位移向量d求各单元的应力。 用Mathematica：deg=Table0, i, 4; “整体坐标系中的单元节点位移向量” ；del=Table0, i, 4; “局部坐标系中的单元节点位移向量” ；del=se=Table0, 0, i, ne; “单元内力数组 ne行” ;Fore=1, e=ne, e+, Fork=1, k=2, k+, Forkk=1, kk=2, kk+, rr =2(k-1)+kk； ss=2(jme, k-1)+kk；degrr=dss; del=TransposeTgje.deg

9、; resultant=De.BL.del; see, 1=e; see, 2=resultantj,1 ;eLeeBDSe,Ne4、总节点载荷向量的排列规律按节点号和节点自由度的顺序（从小到大）排列。设总节点力向量为Rz单元节点位移号：rr =2(k-1)+kk；k =1, 2(单元节点)； kk =1, 2 (节点自由度)；总节点位移号：ss=2(jme, k-1)+kk；rreessssRTRzRz)(4、总节点载荷向量的排列规律单元节点力号：rr =2(k-1)+kk；k =1, 2(单元节点)；kk =1, 2 (节点自由度)；总节点力号：ss=2(jme, k-1)+kk；用Mat

10、hematica编程： nj =ne+1总节点数Pk=Table0, i, 4; “单元节点力向量” ；Rz=Table0, i, 2nj; “总节点力向量2nj行” ;Fore=1, e=ne, e+, DoPki=., i, 4; “生成单元节点力” ； Fork=1, k=2, k+, Forkk=1, kk=2, kk+, rr =2(k-1)+kk； ss=2(jme, k-1)+kk； Rzss=Rzss+Pkrr; “叠加法” ； ;5、总刚度矩阵的排列规律（1）总刚度矩阵中的行和列均按节点号和节点自由度的顺序（从小到大）排列；（2）在总刚度矩阵中，与同一节点相连的单元的刚度要叠

11、加。对一个单元而言：yjxjyixijjiiRRRRvuvukkkkkkkkkkkkkkkk44434241343332312423222114131211 i jijij5、总刚度矩阵的排列规律单刚行号：r =2(i-1)+ii；i =1, 2(单元节点)； ii =1, 2(节点自由度)；单刚列号：s =2(j-1)+jj；j=1, 2(单元节点)； jj=1, 2(节点自由度)；总刚行号：rr=2(jme, i-1)+ii；总刚列号：ss=2(jme, j-1)+jj； jm=Table0, 0, i, ne; “jm数组名，ne总单元数;” jm=.，i节点号，j节点号，;（2）在总刚

12、度矩阵中，与同一节点相连的单元的刚度要叠加。用Mathematica编程：Kz=Table0, i, 2nj, j,2nj; “开总刚度矩阵, nj 总节点数” ；Fore=1, e=ne, e+, “生成单刚,变坐标系” ； Fori=1, i=2, i+, Forii=1, ii=2, ii+, r =2(i-1)+ii; rr=2(jme, i+1-1)+ii； Forj=1, j=2, j+, Forjj=1, jj=2, jj+, s =2(j-1)+jj; ss=2(jme, j+1-1)+jj； “叠加法” ； ; Transpose TkeTek单刚行号：r =2(i-1)+i

13、i；i =1, 2(单元节点)； ii =1, 2 (节点自由度)；单刚列号：s =2(j-1)+jj；j=1, 2(单元节点)； jj=1, 2 (节点自由度)；总刚行号：rr=2(jme, i-1)+ii；总刚列号：ss=32(jme, j-1)+jj；;,srekssrrKzssrrKz6、引入支撑(1) 位移为零的约束前处理时输入信息：支撑数目nz、在整体坐标系中的自由度号zc=, 2(jme, k-1)+kk, ; “k=1, 2; kk=1, 2;”如图 nz=6; zc=3, 4, 5, 6, 7, 8;yxPFori=1, i=nz, i+, z=zci; Rzz=0; For

14、j=1, j=2nj, j+, Kzz, j=0; Fork=1, k=2nj, k+, Kzk, z=0； Kzz,z=1 6、引入支撑(1) 位移为零的约束在Mathematica上实现如图: nz=6; zc=3, 4, 5, 6, 7, 8;(2) 给定位移的约束如图，在节点1给一向右的水平位移u1。 已知位移数目nf=1; 总自由度号和位移值 zf=1, u1;Fori=1, i=nf, i+, z=zfi, 1; zz= zfi, 2; Rzz= zz Kzz,z1015; Kzz,z= Kzz,z1015 yxP算例1：求各杆的轴力，EA不相同，1杆长L。yxP计算结果：各杆内力（解析解）算例2： 计算结果：各节点位移：xyP求：P作用点的位移各杆应力各杆内力：节点号，u/mm，