matlab编程平面梁问题有限元分析

1、2019.4.28平面梁问题有限元分析作业题作业题：平面刚架结构如图所示）由矩形厳面梁组咸，其中水平梁和垂直累各一根，各梁的材料和戴面尺寸均相同，弹性模旨7额陥，載面面积J4xh-（4.5xir1）KpxlO-1）xlOJiiiJp截面惯性矩43町用*,茸它几何尺寸（单1i;in）在图中给岀©试求裁面Z的水平位移和I转角.PlOJNh/-Im.理论解图IF面刚架红构图前处理程序clear;clcXY=1,0,02,2,03,2,-2;%XY为N行3列，N为节点总数ELB=1,1,2,12,2,3,1;%ELB为MB行4列，MB为单元总数b=4.5e-2;h=2e-1;%截面宽b和高h

2、A1=b*h;IZ1=3e-5;%惯性矩EAIZ=200e9,A1,IZ1;%弹性模量与截面积和惯性矩ELPQ=1,1,-2e4,22,1,-1e4,2;%ELPQ第一列为受载荷单元号，第二列为长度C（见下表）,三列为载荷大小，四列为载荷类型SU=1,02,03,0;%SU第一列为已知节点位移对应的方程号，第2列为节点位移数值，约束即等于0子程序1functionPO,ii,jj=dxjd(ELB,XY,ELPQ1)PO=zeros(6,1);k=ELPQ1(1);ii=ELB(k,2);jj=ELB(k,3);dxy=XY(ii,2),XY(ii,3)XY(jj,2),XY(jj,3);DY

3、=dxy(2,2)-dxy(1,2);DX=dxy(2,1)-dxy(1,1);L=sqrt(DXT+DYA2);C=ELPQ1(2);Q=ELPQ1(3);type=ELPQ1(4);switchtypecase1P0(2)=P0(2)+0.5*Q*C*(2-2*S2/LA2+CA3/LA3);P0(3)=PO(3)+Q*CA2*(6-8*C/L+3*CA2/LA2)/12;PO(5)=PO(5)+Q*C-PO(2);PO(6)=PO(6)-Q*CA3*(4-3*C/L)/12/L;case2D=L-C;PO(2)=PO(2)+Q*(L+2*C)*DA2/LA3;PO(3)=PO(3)+Q*

4、C*DA2/LA2;PO(5)=PO(5)+Q-PO(2);PO(6)=PO(6)+Q*D*CA2/LA2;case3D=L-C;PO(1)=PO(1)+Q*D/L;PO(4)=PO(4)+Q*C/L;case4PO(2)=PO(2)+7*Q*L/20;PO(3)=PO(3)+Q*LA2/20;PO(5)=PO(5)+3*Q*L/20;PO(6)=PO(6)+Q*LA2/30;end子程序2functionKE,T=ketb(dxy,E,A,IZ)DY=dxy(2,2)-dxy(1,2);DX=dxy(2,1)-dxy(1,1);L=sqrt(DXA2+DYA2);S=DY/L;C=DX/L;

5、a1=IZ/L;a2=a1/L;a3=E/L;KE=a3*A,0,0,-A,0,00,12*a2,6*a1,0,-12*a2,6*a10,6*a1,4*IZ,0,-6*a1,2*IZ-A,0,0,A,0,00,-12*a2,-6*a1,0,12*a2,-6*a10,6*a1,2*IZ,0,-6*a1,4*IZ;t=C,S,0;-S,C,0;0,0,1;t1=zeros(3,3);T=t,t1;t1,t;end子程序functionKZ=kztb(XY,ELB,EAIZ)N,M=size(XY);KZ=zeros(9,9);MB,m=size(ELB);fork=1:2ii=ELB(k,2);j

6、j=ELB(k,3);LTB=ELB(k,4);dxy=XY(ii,2),XY(ii,3)XY(jj,2),XY(jj,3);E=EAIZ(LTB,1);A=EAIZ(LTB,2);IZ=EAIZ(LTB,3);KE,T=ketb(dxy,E,A,IZ);CN=3*ii-2,3*ii-1,3*ii,3*jj-2,3*jj-1,3*jj;KE=(T')*KE*T;fori=1:6forj=1:6KZ(CN(i),CN(j)=KZ(CN(i),CN(j)+KE(i,j);endendend子程序functionU=weiyi(KZ,P,SU)LR,m=size(SU);fork=1:LRi

7、=SU(k,1);P(i)=1e10*SU(k,2);KZ(i,i)=1e10*KZ(i,i);endU=KZP;end子程序functionP=ydx(XY,ELB,EAIZ,ELPQ)N,m=size(XY)P=zeros(3*N,1)LPQ,m=size(ELPQ)forj=1:LPQELPQ1=ELPQ(j,:)k=ELPQ(j,1)PO,ii,jj=dxjd(ELB,XY,ELPQ1)LTB=ELB(k,4);dxy=XY(ii,2),XY(ii,3)XY(jj,2),XY(jj,3);E=EAIZ(LTB,1);A=EAIZ(LTB,2);IZ=EAIZ(LTB,3);KE,T=k

8、etb(dxy,E,A,IZ);PO=(T')*PO;CN=3*ii-2,3*ii-1,3*ii,3*jj-2,3*jj-1,3*jjfori=1:6P(CN(i)=P(CN(i)+PO(i);endend后处理程序clc;N,M=size(XY);DU=zeros(N,4);fori=1:NDU(i,1)=i;DU(i,2)=U(3*i-2);DU(i,3)=U(3*i-1);DU(i,4)=U(3*i);endMB,M=size(ELB);fori=1:MBP1(i,1)=i;end'节点号水平位移竖直位移转角'DU主程序>>GJINPUT;>>KZ=kztb(XY,ELB,EAIZ);>>P=ydx(XY,ELB,EAIZ,ELPQ);>>U=weiyi(KZ,P,SU);>