悬臂梁MATLAB有限元算例注释.docVIP

用有限元法对悬臂梁分析的算例 算例： 如下图所示的悬臂梁，受均布载荷q＝1N／mm2作用。E＝2．1×105N／mm2, μ＝0.3厚度h＝10mm。现用有限元法分析其位移及应力。 梁可视为平面应力状态，先按图示尺寸划分为均匀的三角形网格，共有8×10＝80个单元，5×ll＝55个节点，坐标轴以及单元与节点的编号如图。将均布载荷分配到各相应节点上，把有约束的节点5l、52、53、54、55视作固定铰链，建立如图所示的离散化计算模型。 程序计算框图： 将各单元刚阵按整体编号集成到整体刚阵K=0计算具有代表性的单元刚阵输入材料参数开 始 （续左） 将各单元刚阵按整体编号集成到整体刚阵 K=0 计算具有代表性的单元刚阵 输入材料参数 开 始 处理根部约束，修改【K 处理根部约束，修改【K】【Q】 求解[K][ 求解[K][δ]＝[Q] 整理[δ] 并画图计算单元应力，并输出 整理[δ] 并画图 计算单元应力，并输出 结束 （接右） 程序中的函数功能介绍及源代码 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*