您的当前位置:首页正文

悬臂梁MATLAB有限元算例注释

2023-07-31 来源:易榕旅网
用有限元法对悬臂梁分析的算例

算例:

如下图所示的悬臂梁,受均布载荷q=1N/mm2作用。E=2.1×105N/mm2, μ=0.3厚度h=10mm。现用有限元法分析其位移及应力。

梁可视为平面应力状态,先按图示尺寸划分为均匀的三角形网格,共有8×10=80个单元,5×ll=55个节点,坐标轴以及单元与节点的编号如图。将均布载荷分配到各相应节点上,把有约束的节点5l、52、53、54、55视作固定铰链,建立如图所示的离散化计算模型。

程序计算框图: (续左)

开 始 处理根部约束,修改【K】【Q】

输入材料参数 求解[K][δ]=[Q] 整理[δ] 并画图 计算具有代表性

的单元刚阵 计算单元应力,并输出 K<=0 结束

将各单元刚阵按

整体编号集成到

整体刚阵 (接右)

1 / 6

程序中的函数功能介绍及源代码

1. LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym)――该函数用于计算平面应力情况下弹性模量为E、泊松比为NU、厚度为t、第一个节点坐标为(xi,yi)、第二个节点坐标为(xj,yj)、第三个节点坐标为(xm,ym)时的线性三角形元的单元刚度矩阵.该函数返回6×6的单位刚度矩阵k.

2. LinearTriangleAssemble(K,k,i,j,m)――该函数将连接节点i,j,m的线性三角形元的单元刚度矩阵k集成

到整体刚度矩阵K。每集成一个单元,该函数都将返回2N×2N的整体刚度矩阵K.

3. LinearTriangleElementStresses(E,NU,t,xi,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;%三角形单元面积,单元节点应该按逆时针排序,保证每个三角形单元的面积都为正值(也可作为一个小函数: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=xj-xi.

D = (E/(1-NU*NU))*[1 NU 0 ; NU 1 0 ; 0 0 (1-NU)/2]; %D为弹性矩阵,分为平面应力问题和平面应变问题

对于平面应力问题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) = 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) = 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*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*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(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;%对号入座,如前所述,每集成一次都将返回2N×2N的整体刚度矩阵K.此题为110×110

2 / 6

function 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;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;%单元应力计算

主程序源代码

E=21e7;NU=0.3;t=0.01;

stifflike5=LinearTriangleElementStiffness(E,NU,t,0.4,0.08,0.36,0.08,0.36,0.06,1)%选取2个基本单元,调用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:49

if rem(i,5)%模取余, bool型变量,非零即为真 j=i;

K=LinearTriangleAssemble(K,stifflike5,j,j+5,j+6);%节点编号 K=LinearTriangleAssemble(K,stifflike1,j,j+6,j+1); end

end%将每个单元刚度矩阵集成到总刚中

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=k\\Q;%高斯消元法,比克莱姆法则在计算速度上有绝对的优势! x=0:0.04:0.4;

plot(x,d(106:-10:6))%基本绘图命令 grid%带网格 y=zeros(80,3); q=0; for i=1:49 switch rem(i,5) case 1 j=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)'; xl=xl-0.04;xm=xm-0.04;xn=xn-0.04; case 2 j=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;

3 / 6

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 3 j=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.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 4 j=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,:)=LinearTriangleElementStresses(E,NU,t,xl,yl,xm,ym,xn,yn,u)'; xl=xl-0.04;xm=xm-0.04;xn=xn-0.04; otherwise q=q+3; end end q=4; for i=1:49 switch rem(i,5) case 1 j=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,yn,u)'; xl=xl-0.04;xm=xm-0.04;xn=xn-0.04; case 2 j=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 3 j=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.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 4 j=2*i;

u=[d(j-1) d(j) d(j+9) d(j+10) d(j+11) d(j+12)]; u=u';

4 / 6

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; otherwise q=q+3; end

end% y(i+q,:)这是实现什么的?没见过这种用法,算法上应该就是通过节点位移实现指定单元的内力,这部分本人看的也晕晕的,望高人指点☺ N=y(73:80,1)

结果图及数据输出

悬臂梁轴线挠度图:

0

-0.1

-0.15

-0.2

-0.25

-0.3

-0.35 00.050.10.150.20.250.3 x/m一单元的单元刚阵 1.0e+006 *

0.8077 0 0 -0.4038 -0.8077 0.4038 0 2.3077 -0.3462 0 0.3462 -2.3077 0 -0.3462 0.5769 0 -0.5769 0.3462 -0.4038 0 0 0.2019 0.4038 -0.2019

w/m-0.050.350.4

5 / 6

-0.8077 0.3462 -0.5769 0.4038 1.3846 -0.7500

0.4038 -2.3077 0.3462 -0.2019 -0.7500 2.5096 五单元的单元刚阵

1.0e+006 *

0.5769 0 -0.5769 0.3462 0 -0.3462 0 0.2019 0.4038 -0.2019 -0.4038 0 -0.5769 0.4038 1.3846 -0.7500 -0.8077 0.3462 0.3462 -0.2019 -0.7500 2.5096 0.4038 -2.3077 0 -0.4038 -0.8077 0.4038 0.8077 0 -0.3462 0 0.3462 -2.3077 0 2.3077

根部73-80各单元应力计算结果如下(n/m2

):

1.0e+007 *

2.1119 -0.0621 -2.2816 -4.8824 5.0479 2.4065 0.0352 -2.37536 / 6

因篇幅问题不能全部显示,请点此查看更多更全内容