有限元法對懸臂梁分析--MATLAB

1、用有限元法對懸臂梁分析的算例算例：h = 10mm。如下圖所示的懸臂梁，受均布載荷q = 1N/mm2作用。E = 2. 1 x 105n/mm2,卩=0.3厚度現(xiàn)用有限元法分析其位移及應力。ly10X405X | = 555l、52、開梁可視為平面應力狀態(tài)，先按圖示尺寸劃分為均勻的三角形網格，共有8X 10= 80個單元,個節(jié)點，坐標軸以及單元與節(jié)點的編號如圖。將均布載荷分配到各相應節(jié)點上，把有約束的節(jié)點 53、54、55視作固定鉸鏈，建立如圖所示的離散化計算模型。程序計算框圖J （續(xù)左）處理根部約束，修改【K】【Q】程序中的函數功能介紹及源代碼1. LinearTriangleElemen

2、tStiffness(E,NU,t,xi,yi,xj,yj,xm,ym) 該函數用于計算平面應力情況下彈性模量為E、泊松比為NU厚度為t、第一個節(jié)點坐標為(xi,yi)、第二個節(jié)點坐標為(xj,yj)、第三個節(jié) 點坐標為(xm, ym)時的線性三角形元的單元剛度矩陣.該函數返回6X 6的單位剛度矩陣k.2. LinearTriangleAssemble(K,k,i,j,m) 該函數將連接節(jié)點 i,j,m 的線性三角形元的單元剛度矩陣 k 集成到整體剛度矩陣 K。每集成一個單元，該函數都將返回2NX 2N的整體剛度矩陣K.3. LinearTriangleElementStresses(E,NU

3、,t,xi,yi,xj,yj,xm,ym,u)-該函數計算在平面應力情況下 彈性模量為E、泊松比為NU厚度為t、第一個節(jié)點坐標為(xi,yi) 第二個節(jié)點坐標為(xj,yj)、第三個節(jié)點坐標 為(xm, ym)以及單元位移矢量為 u時的單元應力。該函數返回單元應力矢量。函數源代碼：function y = LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym)A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj)/2;%三角形單元面積,單元節(jié)點應該按逆時針排序,保證每個三角形單元的面積都為正值(也可作為一個小

4、函數： LinearTriangleElementArea ) betai = yj-ym;betaj = ym-yi;betam = yi-yj;gammai = xm-xj;gammaj = xi-xm;gammam = xj-xi;B = betai 0 betaj 0 betam 0 ;0 gammai 0 gammaj 0 gammam ;gammai betai gammaj betaj gammam betam/(2*A);%B為應變矩陣，其中 betai=yi-ym,betaj=ym-yi,betam=yi-yj.gammai=xm-xj, gammaj=xi-xm, gamm

5、am=xj-xi.D = (E/(1-NU*NU)*1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2;%助彈性矩陣，分為平面應力問題和平面應變問題對于平面應力問題 D = (E/(1-NU*NU)*1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2;對于平面應變問題 E1=E/(1-NU*NU)， NU1=NU(/ 1-NU)y = t*A*B'*D*B;%單元剛度矩陣function y = LinearTriangleAssemble(K,k,i,j,m)K(2*i-1,2*i-1) = K(2*i-1,2*i-1) + k(1,1); K(2*i-1,2*i

6、) = K(2*i-1,2*i) + k(1,2);K(2*i-1,2*j-1) = K(2*i-1,2*j-1) + k(1,3); K(2*i-1,2*j) = K(2*i-1,2*j) + k(1,4);K(2*i-1,2*m-1) = K(2*i-1,2*m-1) + k(1,5); K(2*i-1,2*m) = K(2*i-1,2*m) + k(1,6);K(2*i,2*i-1) = K(2*i,2*i-1) + k(2,1);K(2*i,2*i) = K(2*i,2*i) + k(2,2);K(2*i,2*j-1) = K(2*i,2*j-1) + k(2,3);K(2*i,2*j

7、) = K(2*i,2*j) + k(2,4);K(2*i,2*m-1) = K(2*i,2*m-1) + k(2,5); K(2*i,2*m) = K(2*i,2*m) + k(2,6);K(2*j-1,2*i-1) = K(2*j-1,2*i-1) + k(3,1); K(2*j-1,2*i) = K(2*j-1,2*i) + k(3,2);K(2*j-1,2*j-1) = K(2*j-1,2*j-1) + k(3,3); K(2*j-1,2*j) = K(2*j-1,2*j) + k(3,4);K(2*j-1,2*m-1) = K(2*j-1,2*m-1) + k(3,5); K(2*j

8、-1,2*m) = K(2*j-1,2*m) + k(3,6);K(2*j,2*i-1) = K(2*j,2*i-1) + k(4,1);K(2*j,2*i) = K(2*j,2*i) + k(4,2);K(2*j,2*j-1) = K(2*j,2*j-1) + k(4,3);K(2*j,2*j) = K(2*j,2*j) + k(4,4);K(2*j,2*m-1) = K(2*j,2*m-1) + k(4,5); K(2*j,2*m) = K(2*j,2*m) + k(4,6);K(2*m-1,2*i-1) = K(2*m-1,2*i-1) + k(5,1); K(2*m-1,2*i) =

9、K(2*m-1,2*i) + k(5,2);K(2*m-1,2*j-1) = K(2*m-1,2*j-1) + k(5,3); K(2*m-1,2*j) = K(2*m-1,2*j) + k(5,4);K(2*m-1,2*m-1) = K(2*m-1,2*m-1) + k(5,5); K(2*m-1,2*m) = K(2*m-1,2*m) + k(5,6);K(2*m,2*i-1) = K(2*m,2*i-1) + k(6,1); K(2*m,2*i) = K(2*m,2*i) + k(6,2);K(2*m,2*j-1) = K(2*m,2*j-1) + k(6,3); K(2*m,2*j)

10、= K(2*m,2*j) + k(6,4);K(2*m,2*m-1) = K(2*m,2*m-1) + k(6,5); K(2*m,2*m) = K(2*m,2*m) + k(6,6);y = K; %對號入座，如前所述，每集成一次都將返回2NX 2N的整體剛度矩陣 K.此題為110 X 110function y = LinearTriangleElementStresses(E,NU,t,xi,yi,xj,yj,xm,ym,u)A = (xi*(yj-ym) + xj*(ym-yi) + xm*(yi-yj)/2;betai = yj-ym;betaj = ym-yi;betam = yi

11、-yj;gammai = xm-xj;gammaj = xi-xm;gammam = xj-xi;B = betai 0 betaj 0 betam 0 ;0 gammai 0 gammaj 0 gammam ;gammai betai gammaj betaj gammam betam/(2*A);D = (E/(1-NU*NU)*1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2;%平面應力和平面應變問題兩種情況y = D*B*u; %單元應力計算主程序源代碼%選取 2 個基本單%外部荷載，此處不E=21e7;NU=0.3;t=0.01; stifflike5=LinearTri

12、angleElementStiffness(E,NU,t,0.4,0.08,0.36,0.08,0.36,0.06,1) 元，調用M文件 stifflike1=LinearTriangleElementStiffness(E,NU,t,0.4,0.08,0.36,0.06,0.4,0.06,1) K=sparse(110,110); %creat a xishu matrix for total stiff 創(chuàng)建一個稀疏矩陣 for i=1:49if rem(i,5)%模取余，bool 型變量，非零即為真j=i;K=LinearTriangleAssemble(K,stifflike5,j,j

13、+5,j+6);%節(jié)點編號K=LinearTriangleAssemble(K,stifflike1,j,j+6,j+1);enden d%將每個單元剛度矩陣集成到總剛中K=full(K);%轉化稀疏矩陣 k=K(1:100,1:100);k=K,zeros(100,10);zeros(10,100),eye(10); k=sparse(k);%利用邊界條件簡化基本方程Q=sparse(2:10:92,1,-200,-400,-400,-400,-400,-400,-400,-400,-400,-400,110,1); 包括約束條件，通過形函數確定，是不是可以理解為梁的兩端為中間的一半呢？ d

14、=kQ; %高斯消元法，比克萊姆法則在計算速度上有絕對的優(yōu)勢！ x=0:0.04:0.4;plot(x,d(106:-10:6)%基本繪圖命令grid %帶網格 y=zeros(80,3); q=0;for i=1:49switch rem(i,5)case 1j=2*i;u=d(j-1) d(j) d(j+11) d(j+12) d(j+1) d(j+2);u=u'xl=0.4;yl=0.08;xm=0.36;ym=0.06;xn=0.4;yn=0.06; y(i+q,:)=LinearTriangleElementStresses(E,NU,t,xl,yl,xm,ym,xn,yn,

15、u)' xl=xl-0.04;xm=xm-0.04;xn=xn-0.04;case 2j=2*i;u=d(j-1) d(j) d(j+11) d(j+12) d(j+1) d(j+2);u=u'xl=0.4;yl=0.06;xm=0.36;ym=0.04;xn=0.4;yn=0.04; y(i+q,:)=LinearTriangleElementStresses(E,NU,t,xl,yl,xm,ym,xn,yn,u)' xl=xl-0.04;xm=xm-0.04;xn=xn-0.04;case 3j=2*i;u=d(j-1) d(j) d(j+11) d(j+12) d

16、(j+1) d(j+2);u=u'xl=0.4;yl=0.04;xm=0.36;ym=0.02;xn=0.4;yn=0.02; y(i+q,:)=LinearTriangleElementStresses(E,NU,t,xl,yl,xm,ym,xn,yn,u)' xl=xl-0.04;xm=xm-0.04;xn=xn-0.04;case 4j=2*i;u=d(j-1) d(j) d(j+11) d(j+12) d(j+1) d(j+2);u=u'xl=0.4;yl=0.02;xm=0.36;ym=0;xn=0.4;yn=0;y(i+q,:)=LinearTriangle

17、ElementStresses(E,NU,t,xl,yl,xm,ym,xn,yn,u)'xl=xl-0.04;xm=xm-0.04;xn=xn-0.04;otherwiseq=q+3;endendq=4;for i=1:49switch rem(i,5)case 1j=2*i;u=d(j-1) d(j) d(j+9) d(j+10) d(j+11) d(j+12);u=u'xl=0.4;yl=0.08;xm=0.36;ym=0.08;xn=0.36;yn=0.06;y(i+q,:)=LinearTriangleElementStresses(E,NU,t,xl,yl,xm,ym

18、,xn,yn,u)'xl=xl-0.04;xm=xm-0.04;xn=xn-0.04;case 2j=2*i;u=d(j-1) d(j) d(j+9) d(j+10) d(j+11) d(j+12);u=u'xl=0.4;yl=0.06;xm=0.36;ym=0.06;xn=0.36;yn=0.04; y(i+q,:)=LinearTriangleElementStresses(E,NU,t,xl,yl,xm,ym,xn,yn,u)'xl=xl-0.04;xm=xm-0.04;xn=xn-0.04;case 3j=2*i;u=d(j-1) d(j) d(j+9) d(j

19、+10) d(j+11) d(j+12);u=u'xl=0.4;yl=0.04;xm=0.36;ym=0.04;xn=0.36;yn=0.02;y(i+q,:)=LinearTriangleElementStresses(E,NU,t,xl,yl,xm,ym,xn,yn,u)'xl=xl-0.04;xm=xm-0.04;xn=xn-0.04;case 4j=2*i;u=d(j-1) d(j) d(j+9) d(j+10) d(j+11) d(j+12);u=u'xl=0.4;yl=0.02;xm=0.36;ym=0.02;xn=0.36;yn=0;y(i+q,:)=LinearTriangleElementStresses(E,NU,t,xl,yl,xm,ym,xn,yn,u)'xl=xl-0.04;xm=xm-0.04;xn=xn-0.04;otherwis