运用MATLAB对桁架单元进行有限元分析

《运用MATLAB对桁架单元进行有限元分析》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、17-3-16下午12:54E:...第1页，共7页np=3;ne=2;nload=1;nb=4;nu=0;%np为节点数，ne为单元数，nload为外力数，nb为约束数，nu为泊松比np2=np*2;np3=np2-nb;%np2为不受约束时自由度是节点的两倍，np3为不受约束的节点自由度个数K=sym(zeros(np2,np2));%定义受整体刚度空数组F=sym(zeros(np2,1));%定义节点外力空数组KK=sym(zeros(np3,np3));%预置自由度刚度空数组FF=sym(zeros(np3,1));%预置自由度外力空数组symshAPEL%定

2、义未知正常数为变量xy=[0,h;sqrt(3)*h,0;0,0];%节点横纵坐标数组AA=[A;A];%单元杆件面积load=[2,2,-P];%荷载数组bound=[1,1,0;1,2,0;3,1,0;3,2,0];%边界条件数组IJ=[1,2;3,2];%各杆单元节点编号17-3-16下午12:54E:...第2页，共7页DW=zeros(1,4);forie=1:neip=IJ(ie,1);jp=IJ(ie,2);DW(1)=ip*2-1;DW(2)=ip*2;DW(3)=jp*2-1;DW(4)=jp*2;%给单元节点横纵方向编号x1=xy(ip,1);x2=

3、xy(jp,1);y1=xy(ip,2);y2=xy(jp,2);%杆单元节点x，y坐标L=sqrt((x2-x1)^2+(y2-y1)^2);%杆单元长度c=(x2-x1)/L;s=(y2-y1)/L;%c为余弦，s为正弦T=[c,s,0,0;-s,c,0,0;0,0,c,s;0,0,-s,c];%坐标转换矩阵（将局部坐标转换为整体坐标）A1=AA(ie);ke=[E*A1/L,0,-E*A1/L,0;0,0,0,0;-E*A1/L,0,E*A1/L,0;0,0,0,0];k=T.'*ke*T;%将局部坐标下单元刚度转换为整体坐标下单元刚度(转置后面加.可以去掉结果中

4、的虚数)fori=1:4i1=DW(i);forj=1:4j1=DW(j);K(i1,j1)=K(i1,j1)+k(i,j);%集成整体刚度矩阵17-3-16下午12:54E:...第3页，共7页endendendfori=1:nloadi1=(load(i,1)-1)*2+load(i,2);F(i1,1)=F(i1,1)+load(i,3);%将荷载按节点方向代码累加，计入外力荷载endFuu=sym(zeros(np2,1));%节点位移NR=zeros(np2,1);fori=1:nbi1=(bound(i,1)-1)*2+bound(i,2);NR(i1)=-

5、i1;%将有约束的节点的横纵节点编号挑出来变为负数uu(i1)=bound(i,3);%有约束的四个方向位移为0endj=0;fori=1:np2i1=NR(i);ifi1==0j=j+1;NR(i)=j;end17-3-16下午12:54E:...第4页，共7页end%此循环目的是为了将没有约束的节点以及方向挑出来fori=1:np2i1=NR(i);ifi1>0FF(i1)=F(i);end%将没有约束的节点方向的外力储存下来forj=1:np2j1=NR(j);if(i1>0&j1>0)KK(i1,j1)=K(i,j);end%将没有约束的节点的刚度储存下来end

6、endKKFFfori=1:np2i1=NR(i);forj=1:np2j1=NR(j);if(i1>0&j1<0)jj1=abs(j1);FF(i1)=FF(i1)-K(i,jj1)*uu(jj1);end17-3-16下午12:54E:...第5页，共7页%由于有约束方向的位移为0，所以无位移的地方FF[i1]不变endendu=sym(zeros(np2,1));u=KKFF%求解线性方程fori=1:np2i1=NR(i);ifi1>0uu(i)=u(i1);end%将求解的无约束方向的位移放进总的位移数组相应的位置中去enduudisp('位移')forip

7、=1:npstr=[ip;uu(ip*2-1);uu(ip*2)];disp(str)%输出节点位移endforie=1:neip=IJ(ie,1);jp=IJ(ie,2);DW(1)=ip*2-1;DW(2)=ip*2;DW(3)=jp*2-1;DW(4)=jp*2;x1=xy(ip,1);x2=xy(jp,1);y1=xy(ip,2);y2=xy(jp,2);17-3-16下午12:54E:...第6页，共7页L=sqrt((x2-x1)^2+(y2-y1)^2);c=(x2-x1)/L;s=(y2-y1)/L;T=[c,s,0,0;