悬臂梁MATLAB有限元算例注释

1、用有限元法对悬臂梁分析的算例算例:h = 10mm。如下图所示的悬臂梁，受均布载荷q = 1N/ mm2作用。E = 2. 1 x 105n/mm2,卩=0.3厚度现用有限元法分析其位移及应力。jt5X | = 555l、52、开梁可视为平面应力状态，先按图示尺寸划分为均匀的三角形网格，共有8X 10= 80个单元,个节点，坐标轴以及单元与节点的编号如图。将均布载荷分配到各相应节点上，把有约束的节点 53、54、55视作固定铰链，建立如图所示的离散化计算模型。程序计算框图J （续左）处理根部约束，修改【K】【Q】程序中的函数功能介绍及源代码1. LinearTriangleElementSti

2、ffness(E,NU,t,xi,yi,xj,yj,xm,ym) 该函数用于计算平面应力情况下弹性模量为E、泊松比为NU厚度为t、第一个节点坐标为(xi,yi)、第二个节点坐标为(xj,yj)、第三个节 点坐标为(xm, ym)时的线性三角形元的单元刚度矩阵.该函数返回6X 6的单位刚度矩阵k.2. LinearTriangleAssemble(K,k,i,j,m) 该函数将连接节点 i,j,m 的线性三角形元的单元刚度矩阵 k 集成到整体刚度矩阵 K。每集成一个单元，该函数都将返回2NX 2N的整体刚度矩阵K.3. LinearTriangleElementStresses(E,NU,t,x

3、i,yi,xj,yj,xm,ym,u)-该函数计算在平面应力情况下 弹性模量为E、泊松比为NU厚度为t、第一个节点坐标为(xi,yi) 第二个节点坐标为(xj,yj)、第三个节点坐标 为(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;%三角形单元面积,单元节点应该按逆时针排序,保证每个三角形单元的面积都为正值(也可作为一个小函数：

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, gammam=x

5、j-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-1,2

8、*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) = K(2*

9、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) = K(

10、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-yj;

11、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=LinearTriangl

12、eElementStiffness(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 创建一个稀疏矩阵 for i=1:49if rem(i,5)%模取余，bool 型变量，非零即为真j=i;K=LinearTriangleAssemble(K,stifflike5,j,j+5,j

13、+6);%节点编号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=kQ;

14、 %高斯消元法，比克莱姆法则在计算速度上有绝对的优势！ 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,u)&#

15、39; 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(j+1

16、) 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,:)=LinearTriangleElem

17、entStresses(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,xn,

18、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+10)

19、 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