平面桁架的有限元法

平面桁架的有限元法一、建立有限元模型平面桁架是现成的有限元模型，只有节点力二、单元分析位移假设，几何方程，物理方程，单元刚度矩阵三、总刚度矩阵、总载荷列阵、引入支撑每个节点建立平衡方程集成技术（叠加法、对号入座）四、求解、结果分析，后处理一、建立有限元模型平面桁架是现成的有限元模型，只有节点力桁架原型=有限元模型yxP二、单元分析平面二力杆单元注意：一个单元的节点数目（2）；一个节点的自由度数目（2）；有整体坐标系和局部坐标系之分xoyxyuiiujjijbvv二、单元分析（先在局部坐标系中进行）问题：当节点位移一定，相应的节点力为多大？节点力与节点位移之间的比例常数为刚度系数目标：建立单元刚度矩阵。方法：结构力学的位移法加虚功原理，皆用节点位移表示。xoyxyuiiujjijbvv1、位移假设：设各点的位移，即位移场为xaavxaau4321421,,aaa其中为待定系数，由节点坐标和节点位移确定。将位移写成矩阵形式：4321100001aaaaxxvu}]{[}{aHfsxoyxyuiiujjijbvvvu在节点上满足：baavbaauavauijii432131写成矩阵形式：}]{[}{aAbexaavxaau4321bxxji,0xoyxyuiiujjijbvvvu432110000101000001aaaabbvuvujjii解线性代数方程组，得}]{[}{aAbe}{][}{1ebAa}]{[}{aHfs代入得141444212}{][][}{ebsAHf144212}{][}{efNf432110000101000001aaaabbvuvujjii4321100001aaaaxxvu1]][[)]([bsAHxN}]{[}{efNf节点位移与单元内位移的关系其中为形函数。用Machematica推导形函数：}{e在编程中只关心节点位移的系数矩阵xoyxyuiiujjijbvvvuxxHs100001][2、应变（从位移求应变）只有轴向应变3、应力（从应变求应力，本构方程）}]{[0}{eLvuBxuyvxu0}{xuEEivuivuxaavxaau43214、内力（从应力求内力）xuAEANiuyvxuAEAENSiie000}{二力杆只有轴向力0yv用Machematica推导应力矩阵：(EA的两个字母连写作参数名)}]{[}{eD用内力表示的物理方程yvxuEAEAN000EAEADe00][}]{[0}{eLBxuyvxu0}{NSe}]{][[}{eLeeBDS}]{[}{efNf4、单元刚度矩阵（用虚功原理）给节点一组虚位移，单元内即产生相应的虚应变，实际节点力（轴向力、横向力、弯矩）在虚位移上所做的功等于实际应力（横向剪应力可忽略不计）在虚应变上产生的应变能。用Machematica推导单刚矩阵[Ke]：beLeTLTeleTeTexδBDBlSR0}d]{][[][}~{d}{}~{}{}~{}~{e}~]{[}~{eLB}]{][[}{eLeeBDS00}{jieNNR144414}{][}{eeeδKRbLeTLexBDBK042222444d][][][][用Machematica推导，所得结果以矩阵形式显示三、总刚度矩阵、总载荷列阵、引入支撑1、坐标变换前面的推导是在局部坐标系中进行的，而总刚度矩阵、总载荷列阵、引入支撑等要在整体坐标中进行，因此，存在坐标变换问题。xoyxycossinsincosiiiiiiwuwwuuiiiiwuwucossinsincos}]{[}{iit三、总刚度矩阵、总载荷列阵、引入支撑1、坐标变换}]{[}{iitjijitt00}]{[}{eeTttT00][TTTTtt][][][][11}]{[}{eeRTR}]{[}{eekR}]{[}{eekR}]{[}]{[eekRT}]{][[}]{][[eeTkkTTTkTk]][][[][xoyxyuiiujjijbvvvucossinsincos][t三、总刚度矩阵、总载荷列阵、引入支撑2、在总结构和整体坐标中讨论问题，所需要的信息：（1）节点号；（2）单元长度；（2）单元倾角(单元起点到终点的有向线段与整体横轴之夹角）；（4）节点力所在的节点号和作用方向；（5）支撑所在的节点号和作用方向。3、总节点位移向量的排列规律按节点号和节点自由度的顺序（从小到大）排列。三、总刚度矩阵、总载荷列阵、引入支撑建立节点号数组（用Mathematica)：jm=Table[{0,0},{i,ne}];“jm数组名，ne总单元数;”jm={…，{i节点号，j节点号}，…};建立单元长度数组gc=Table[{0},{i,ne}];“gc数组名，ne总单元数;”gc={L1,L2,…Lne};建立单元方向角数组gj=Table[{0},{i,ne}];“gj数组名，ne总单元数;”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(jm[[e,k]]-1)+kk;e=单元号rreessTd})]{([对每个单元进行节点循环和自由度循环，从单元地址到结构地址搬家是单元到结构的映射(map)3、总节点位移向量的排列规律单元节点位移号：rr=2(k-1)+kk；k=1,2(单元节点)；kk=1,2(节点自由度)总节点位移号：ss=2(jm[[e,k]]-1)+kk；例如，已解得结构的总节点位移向量{d}求各单元的应力。用Mathematica：deg=Table[0,{i,4}];“整体坐标系中的单元节点位移向量”；del=Table[0,{i,4}];“局部坐标系中的单元节点位移向量”；del=se=Table[{0,0},{i,ne}];“单元内力数组ne行”;For[e=1,e=ne,e++,For[k=1,k=2,k++,For[kk=1,kk=2,kk++,rr=2(k-1)+kk；ss=2(jm[[e,k]]-1)+kk；deg[[rr]]=d[[ss]]]];del=Transpose[T[gj[[e]]]].deg;resultant=De.BL.del;se[[e,1]]=e;se[[e,2]]=resultant[[j,1]]];}]{][[}{eLeeBDS}{e},{Ne4、总节点载荷向量的排列规律按节点号和节点自由度的顺序（从小到大）排列。设总节点力向量为{Rz}单元节点位移号：rr=2(k-1)+kk；k=1,2(单元节点)；kk=1,2(节点自由度)；总节点位移号：ss=2(jm[[e,k]]-1)+kk；rreessssRTRzRz})]{([4、总节点载荷向量的排列规律单元节点力号：rr=2(k-1)+kk；k=1,2(单元节点)；kk=1,2(节点自由度)；总节点力号：ss=2(jm[[e,k]]-1)+kk；用Mathematica编程：nj=ne+1总节点数Pk=Table[0,{i,4}];“单元节点力向量”；Rz=Table[0,{i,2nj}];“总节点力向量2nj行”;For[e=1,e=ne,e++,Do[Pk[[i]]=….,{i,4}];“生成单元节点力”；For[k=1,k=2,k++,For[kk=1,kk=2,kk++,rr=2(k-1)+kk；ss=2(jm[[e,k]]-1)+kk；Rz[[ss]]=Rz[[ss]+Pk[[rr]];“叠加法”；]]];5、总刚度矩阵的排列规律（1）总刚度矩阵中的行和列均按节点号和节点自由度的顺序（从小到大）排列；（2）在总刚度矩阵中，与同一节点相连的单元的刚度要叠加。对一个单元而言：yjxjyixijjiiRRRRvuvukkkkkkkkkkkkkkkk44434241343332312423222114131211ijijij5、总刚度矩阵的排列规律单刚行号：r=2(i-1)+ii；i=1,2(单元节点)；ii=1,2(节点自由度)；单刚列号：s=2(j-1)+jj；j=1,2(单元节点)；jj=1,2(节点自由度)；总刚行号：rr=2(jm[[e,i]]-1)+ii；总刚列号：ss=2(jm[[e,j]]-1)+jj；jm=Table[{0,0},{i,ne}];“jm数组名，ne总单元数;”jm={…….，{i节点号，j节点号}，……};（2）在总刚度矩阵中，与同一节点相连的单元的刚度要叠加。用Mathematica编程：Kz=Table[0,{i,2nj},{j,2nj}];“开总刚度矩阵,nj总节点数”；For[e=1,e=ne,e++,“生成单刚,变坐标系”；For[i=1,i=2,i++,For[ii=1,ii=2,ii++,r=2(i-1)+ii;rr=2(jm[[e,i+1]]-1)+ii；For[j=1,j=2,j++,For[jj=1,jj=2,jj++,s=2(j-1)+jj;ss=2(jm[[e,j+1]]-1)+jj；“叠加法”；]]]]];;][TransposeTkeTek单刚行号：r=2(i-1)+ii；i=1,2(单元节点)；ii=1,2(节点自由度)；单刚列号：s=2(j-1)+jj；j=1,2(单元节点)；jj=1,2(节点自由度)；总刚行号：rr=2(jm[[e,i]]-1)+ii；总刚列号：ss=32(jm[[e,j]]-1)+jj；;]],[[]],[[]],[[srekssrrKzssrrKz6、引入支撑(1)位移为零的约束前处理时输入信息：支撑数目nz、在整体坐标系中的自由度号zc={…,2(jm[[e,k]]-1)+kk,…};“k=1,2;kk=1,2;”如图nz=6;zc={3,4,5,6,7,8};yxPFor[i=1,i=nz,i++,z=zc[[i]];Rz[[z]]=0;For[j=1,j=2nj,j++,Kz[[z,j]]=0;For[k=1,k=2nj,k++,Kz[[k,z]]=0]]；Kz[[z,z]]=1]6、引入支撑(1)位移为零的约束在Mathematica上实现如图:nz=6;zc={3,4,5,6,7,8};(2)给定位移的约束如图，在节点1给一向右的水平位移u1。已知位移数目nf=1;总自由度号和位移