有限元作業(yè)：三角形單元求解

《有限元作業(yè)》年級2015級學院機電工程學院專業(yè)名稱班級學號學生2016年05月如下圖所示為一受集中力P作用的結(jié)構(gòu)，彈性模量E為常量，泊松比V=1/6,厚度為I=1。按平面應力問題計算，運用有限元方法，分別采用三角形及四邊形單元求解，求節(jié)點位移及單元應力（要求三角形單元數(shù)量不少于4個，四邊形單元不少于2個）圖(一)圖（二）三角形單元求解圖（三）四邊形單元求解如圖劃分三角形單元，工分成四個分別為④如圖分別進行編號1、2、3、4、5、6，并建立坐標系編程進行求解，得出結(jié)果，其中假設(shè)力P=2000N調(diào)用Triangle2D3Node_Stiffness函數(shù)，求出單元剛度矩陣k1=1.0e+06*7.2857-3.0000-2.14290.8571-5.14292.1429-3.00007.28572.1429-5.14290.8571-2.1429-2.14292.14292.142900-2.14290.8571-5.142905.1429-0.85710-5.14290.85710-0.85715.142902.1429-2.1429-2.1429002.1429k2=1.0e+06*5.14290-5.14290.85710-0.857102.14292.1429-2.1429-2.14290-5.14292.14297.2857-3.0000-2.14290.85710.8571-2.1429-3.00007.28572.1429-5.14290-2.1429-2.14292.14292.14290-0.857100.8571-5.142905.1429k3=1.0e+06*2.14290-2.1429-2.142902.142905.1429-0.8571-5.14290.85710-2.1429-0.85717.28573.0000-5.1429-2.1429-2.1429-5.14293.00007.2857-0.8571-2.142900.8571-5.1429-0.85715.142902.14290-2.1429-2.142902.1429k4=1.0e+06*2.14290-2.1429-2.142902.142905.1429-0.8571-5.14290.85710-2.1429-0.85717.28573.0000-5.1429-2.1429-2.1429-5.14293.00007.2857-0.8571-2.142900.8571-5.1429-0.85715.142902.14290-2.1429-2.142902.1429調(diào)用Triangle2D3Node_Assembly函數(shù)，求出總體剛度矩陣求出的節(jié)點位移U=0000-0.00040.00080.00050.00100.00070.0023-0.00070.0026調(diào)用Triangle2D3Node_Stress函數(shù)，求出應力,S1、S2、S3、中求出的分別為Sx,Sy,SxyS1=1.0e+03*-4.4086-0.73483.5914S2=1.0e+03*4.4086-0.64050.4086S3=1.0e+03*1.8907-1.06012.1093S4=1.0e+03*-1.89072.10931.8907二、（1）如圖劃分四邊形單元，工分成四個分別為（2）如圖分別進行編號1、2、3、4、5、6，并建立坐標系（3）編程進行求解，得出結(jié)果，其中假設(shè)力P=2000N調(diào)用Quad2D4Node_Stiffness函數(shù)，求出單元剛度矩陣調(diào)用Quad2D4Node_Assembly函數(shù)，求出求出總體剛度矩陣求出節(jié)點位移U=00000.00120.0017-0.00120.00170.00160.0049-0.00170.0052調(diào)用Quad2D4Node_Stress函數(shù)，求出單元應力中的的S1、S2、S3分別為Sx,Sy,Sxy應力分量S1=1.0e+03*0.0000-0.24782.0000S2=1.0e+07*0.68564.1135-1.7137程序附錄一、1、三角形單元總程序：E=1e7;NU=1/6;t=1;ID=1;%調(diào)用Triangle2D3Node_Stiffness函數(shù)，求出單元剛度矩陣k1=Triangle2D3Node_Stiffness(E,NU,t,0,1,0,0,1,1,ID)k2=Triangle2D3Node_Stiffness(E,NU,t,0,0,1,0,1,1,ID)k3=Triangle2D3Node_Stiffness(E,NU,t,1,1,1,0,2,0,ID)k4=Triangle2D3Node_Stiffness(E,NU,t,2,0,2,1,1,1,ID)%調(diào)用Triangle2D3Node_Assembly函數(shù)，求出總體剛度矩陣KK=zeros(12,12);KK=Triangle2D3Node_Assembly(KK,k1,1,2,3);KK=Triangle2D3Node_Assembly(KK,k2,2,4,3);KK=Triangle2D3Node_Assembly(KK,k3,3,4,5);KK=Triangle2D3Node_Assembly(KK,k4,5,6,3)%邊界條件的處理及剛度方程求解k=KK(5:12,5:12)p=[0;0;0;0;0;0;0;2000]u=k\p%支反力的計算U=[0;0;0;0;u]%為節(jié)點位移P=KK*U%調(diào)用Triangle2D3Node_Strain函數(shù)，求出應變SN1、SN2、SN3中求出的分別為SNx,SNy,SNxyu1=[U(1);U(2);U(3);U(4);U(5);U(6)];u2=[U(3);U(4);U(7);U(8);U(5);U(6)];u3=[U(5);U(6);U(7);U(8);U(9);U(10)];u4=[U(9);U(10);U(11);U(12);U(5);U(6)];SN1=Triangle2D3Node_Strain(0,1,0,0,1,1,u1)SN2=Triangle2D3Node_Strain(0,0,1,0,1,1,u2)SN3=Triangle2D3Node_Strain(1,1,1,0,2,0,u3)SN4=Triangle2D3Node_Strain(2,0,2,1,1,1,u4)%調(diào)用Triangle2D3Node_Stress函數(shù)，求出應力,S1、S2、S3、中求出的分別為Sx,Sy,Sxyu1=[U(1);U(2);U(3);U(4);U(5);U(6)];u2=[U(3);U(4);U(7);U(8);U(5);U(6)];u3=[U(5);U(6);U(7);U(8);U(9);U(10)];u4=[U(9);U(10);U(11);U(12);U(5);U(6)];S1=Triangle2D3Node_Stress(E,NU,0,1,0,0,1,1,u1,ID)S2=Triangle2D3Node_Stress(E,NU,0,0,1,0,1,1,u2,ID)S3=Triangle2D3Node_Stress(E,NU,1,1,1,0,2,0,u3,ID)S4=Triangle2D3Node_Stress(E,NU,2,0,2,1,1,1,u4,ID)2、求剛度矩陣程序functionk=Triangle2D3Node_Stiffness(E,NU,t,xi,yi,xj,yj,xm,ym,ID)%該函數(shù)計算單元的剛度矩陣%輸入彈性模量E，泊松比NU，厚度t%輸入三個節(jié)點i、j、m的坐標xi,yi,xj,yj,xm,ym%輸入平面問題性質(zhì)指示參數(shù)ID(1為平面應力，2為平面應變)%輸出單元剛度矩陣k(6X6)%---------------------------------------------------------------A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=[betai0betaj0betam0;0gammai0gammaj0gammam;gammaibetaigammajbetajgammambetam]/(2*A);ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];endk=t*A*B'*D*B;求整體剛度矩陣functionz=Triangle2D3Node_Assembly(KK,k,i,j,m)%該函數(shù)進行單元剛度矩陣的組裝%輸入單元剛度矩陣k%輸入單元的節(jié)點編號I、j、m%輸出整體剛度矩陣KK%---------------------------------------------------------------DOF(1)=2*i-1;DOF(2)=2*i;DOF(3)=2*j-1;DOF(4)=2*j;DOF(5)=2*m-1;DOF(6)=2*m;forn1=1:6forn2=1:6KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);endendz=KK;4、求應變程序functionstrain=Triangle2D3Node_Strain(xi,yi,xj,yj,xm,ym,u)%該函數(shù)計算單元的應變%輸入三個節(jié)點i、j、m的坐標xi,yi,xj,yj,xm,ym%輸入單元的位移列陣u(6X1)%輸出單元的應力strain(3X1)，由于它為常應變單元，則單元的應變分量為SNx,SNy,SNz%---------------------------------------------------------------A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=[betai0betaj0betam0;0gammai0gammaj0gammam;gammaibetaigammajbetajgammambetam]/(2*A);strain=B*u;5、求應力程序functionstress=Triangle2D3Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,u,ID)%該函數(shù)計算單元的應力%輸入彈性模量E，泊松比NU，厚度t%輸入三個節(jié)點i、j、m的坐標xi,yi,xj,yj,xm,ym%輸入平面問題性質(zhì)指示參數(shù)ID(1為平面應力，2為平面應變)，單元的位移列陣u(6X1)%輸出單元的應力stress(3X1)，由于它為常應力單元，則單元的應力分量為Sx,Sy,Sxy%---------------------------------------------------------------A=(xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))/2;betai=yj-ym;betaj=ym-yi;betam=yi-yj;gammai=xm-xj;gammaj=xi-xm;gammam=xj-xi;B=[betai0betaj0betam0;0gammai0gammaj0gammam;gammaibetaigammajbetajgammambetam]/(2*A);ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];endstress=D*B*u;二、1、四邊形單元總程序：E=1e7;NU=1/6;h=1;ID=1;%調(diào)用Quad2D4Node_Stiffness函數(shù)，求出單元剛度矩陣k1=Quad2D4Node_Stiffness(E,NU,h,0,1,0,0,1,0,1,1,ID)k2=Quad2D4Node_Stiffness(E,NU,h,1,0,2,0,2,1,1,1,ID)%調(diào)用Quad2D4Node_Assembly函數(shù)，求出求出總體剛度矩陣KK=zeros(12,12);KK=Quad2D4Node_Assembly(KK,k1,1,2,3,4);KK=Quad2D4Node_Assembly(KK,k2,3,5,6,4)%邊界條件的處理及剛度方程求解k=KK(5:12,5:12)p=[0;0;0;0;0;0;0;2000]u=k\p%支反力的計算U=[0;0;0;0;u]%為節(jié)點位移P=KK*U%調(diào)用Quad2D4Node_Stress函數(shù)，求出單元應力中的的S1、S2、S3分別為Sx,Sy,Sxy應力分量u1=[U(1);U(2);U(3);U(4);U(5);U(6);U(7);U(8)];u2=[U(5);U(6);U(9);U(10);U(11);U(12);U(7);(8)];S1=Quad2D4Node_Stress(E,NU,0,1,0,0,1,0,1,1,u1,ID)S2=Quad2D4Node_Stress(E,NU,1,0,2,0,2,1,1,1,u2,ID)2、求剛度矩陣程序functionk=Quad2D4Node_Stiffness(E,NU,h,xi,yi,xj,yj,xm,ym,xp,yp,ID)%該函數(shù)計算單元的剛度矩陣%輸入彈性模量E，泊松比NU，厚度h%輸入4個節(jié)點i、j、m、p的坐標xi,yi,xj,yj,xm,ym,xp,yp%輸入平面問題性質(zhì)指示參數(shù)ID(1為平面應力，2為平面應變)%輸出單元剛度矩陣k(8X8)%---------------------------------------------------------------symsst;a=(yi*(s-1)+yj*(-1-s)+ym*(1+s)+yp*(1-s))/4;b=(yi*(t-1)+yj*(1-t)+ym*(1+t)+yp*(-1-t))/4;c=(xi*(t-1)+xj*(1-t)+xm*(1+t)+xp*(-1-t))/4;d=(xi*(s-1)+xj*(-1-s)+xm*(1+s)+xp*(1-s))/4;B1=[a*(t-1)/4-b*(s-1)/40;0c*(s-1)/4-d*(t-1)/4;c*(s-1)/4-d*(t-1)/4a*(t-1)/4-b*(s-1)/4];B2=[a*(1-t)/4-b*(-1-s)/40;0c*(-1-s)/4-d*(1-t)/4;c*(-1-s)/4-d*(1-t)/4a*(1-t)/4-b*(-1-s)/4];B3=[a*(t+1)/4-b*(s+1)/40;0c*(s+1)/4-d*(t+1)/4;c*(s+1)/4-d*(t+1)/4a*(t+1)/4-b*(s+1)/4];B4=[a*(-1-t)/4-b*(1-s)/40;0c*(1-s)/4-d*(-1-t)/4;c*(1-s)/4-d*(-1-t)/4a*(-1-t)/4-b*(1-s)/4];Bfirst=[B1B2B3B4];Jfirst=[01-tt-ss-1;t-10s+1-s-t;s-t-s-10t+1;1-ss+t-t-10];J=[xixjxmxp]*Jfirst*[yi;yj;ym;yp]/8;B=Bfirst/J;ifID==1D=(E/(1-NU*NU))*[1NU0;NU10;00(1-NU)/2];elseifID==2D=(E/(1+NU)/(1-2*NU))*[1-NUNU0;NU1-NU0;00(1-2*NU)/2];endBD=J*transpose(B)*D*B;r=int(int(BD,t,-1,1),s,-1,1);z=h*r;k=double(z);3、求總體剛度矩陣程序functionz=Quad2D4Node_Assembly(KK,k,i,j,m,p)%該函數(shù)進行單元剛度矩陣的組裝%輸入單元剛度矩陣k，單元的節(jié)點編號i、j、m、p%輸出整體剛度矩陣KK%---------------------------------------------------------------DOF(1)=2*i-1;DOF(2)=2*i;DOF(3)=2*j-1;DOF(4)=2*j;DOF(5)=2*m-1;DOF(6)=2*m;DOF(7)=2*p-1;DOF(8)=2*p;forn1=1:8forn2=1:8KK(DOF(n1),DOF(n2))=KK(DOF(n1),DOF(n2))+k(n1,n2);endendz=KK;4、求應力程序functionstress=Quad2D4Node_Stress(E,NU,xi,yi,xj,yj,xm,ym,xp,yp,u,ID)%該函數(shù)計算單元的應力%輸入彈性模量E，泊松比NU，厚度h，%輸入4個節(jié)點i、j、m、p的坐標xi,yi,xj,yj,xm,ym,xp,yp,%輸入平面問題性質(zhì)指示參數(shù)ID(1為平面應力，2為平面應變)%輸入單元的位移列陣u(8X1)%輸出單元的應力stress(3X1)%由于它為常應力單元，則單元的應力分量為Sx,Sy,Sxy%-------------------------------------------------