梁单元有限元法分析关键词：梁单元有限元分析1.摘要：二维平面梁单元是梁单元中最简单的单元之一，这次作业旨在学习如何运用有限元分析法分析梁单元。

2.目的：运用MATLAB软件分析二维梁单元。

3.题目：设一方形的截面梁，截面每边长为5cm,长度为10m,在左端约束固定，在右端施以一个沿y方向的集中力ω=100N,求其挠度与转角。

3.建立有限元分析模型：（1）.结构离散化：单元的选择：由于为悬臂梁，且横向的长度远远小于轴向的长度，所以在这选择平面梁单元；单元的数量：将这个梁从中间划分为两个单元；建立坐标系，坐标系包括结构的整体坐标系与单元的局部坐标系；（2.）建立平面梁单元的位移模式：建立整体坐标系：建立一个有两个单元组成的模型，由于X方向的位移U1,U2,U3太小所以我们略去这三个自由度的变化；节点坐标码：单元编码：同样出1号单元，建立局部坐标系：4.具体的MATLAB求解过程与结果：>> clearx1=0;x2=sym('L');x=sym('x');j=0:3;v=x.^jv =[ 1, x, x^2, x^3]>> %计算形函数矩阵m=...[1 x1 x1^2 x1^30 1 2*x1 3*x1^21 x2 x2^2 x2^30 1 2*x2 3*x2^2]m =[ 1, 0, 0, 0][ 0, 1, 0, 0][ 1, L, L^2, L^3][ 0, 1, 2*L, 3*L^2]>> mm=inv(m)mm =[ 1, 0, 0, 0] [ 0, 1, 0, 0] [ -3/L^2, -2/L, 3/L^2, -1/L][ 2/L^3, 1/L^2, -2/L^3, 1/L^2]>> mm=inv(m);N =[ (2*x^3)/L^3 - (3*x^2)/L^2 + 1, x - (2*x^2)/L + x^3/L^2, (3*x^2)/L^2 - (2*x^3)/L^3, x^3/L^2 - x^2/L]>> %N=[1 x x^2 x^3]*（inv(m)）%由最小势能原理导出刚度矩阵keB=diff(N,2) %梁单元的单元应变矩阵是形函数矩阵的2介导数（由梁的应变能得出）B =[ (12*x)/L^3 - 6/L^2, (6*x)/L^2 - 4/L, 6/L^2 - (12*x)/L^3, (6*x)/L^2 - 2/L]>> k=transpose(B)*(B);ke=int(k,0,'L') %从0到L上积分k矩阵ke =[ 12/L^3, 6/L^2, -12/L^3, 6/L^2][ 6/L^2, 4/L, -6/L^2, 2/L][ -12/L^3, -6/L^2, 12/L^3, -6/L^2][ 6/L^2, 2/L, -6/L^2, 4/L]>> %Element1:E=4.0e11,I=bh^3/12=5.2e-7EI=4.0e11*5.2e-7EI =208000>> ke1=EI*subs(ke,'L',5)ke1 =19968 49920 -19968 4992049920 166400 -49920 83200-19968 -49920 19968 -4992049920 83200 -49920 166400>> %由上市我们就计算出了1号单元刚度矩阵ke,由于分成两个单元，所以L=10/2=5>> %同理，我们运用上述办法得到2号单元的刚度矩阵ke2>> clearx2=5;x=sym('x');j=0:3;v=x.^jv =[ 1, x, x^2, x^3]>> m=...[1 x2 x2^2 x2^30 1 2*x2 3*x2^21 x3 x3^2 x3^30 1 2*x3 3*x3^2]m =[ 1, 5, 25, 125][ 0, 1, 10, 75][ 1, L, L^2, L^3][ 0, 1, 2*L, 3*L^2]>> mm=inv(m);N=v*mmN =[ (2*x^3)/(L - 5)^3 + (30*L*x)/(L - 5)^3 - (x^2*(3*L + 15))/(L - 5)^3 + (L^2*(L - 15))/(L - 5)^3, x^3/(L - 5)^2 -(5*L^2)/(L - 5)^2 - (x^2*(2*L + 5))/(L - 5)^2 + (L*x*(L + 10))/(L - 5)^2, (75*L - 125)/(L - 5)^3 - (2*x^3)/(L - 5)^3 - (30*L*x)/(L - 5)^3 + (x^2*(3*L + 15))/(L - 5)^3, x^3/(L - 5)^2 - (25*L)/(L - 5)^2 + (x*(10*L + 25))/(L - 5)^2 - (x^2*(L + 10))/(L - 5)^2]>> B=diff(N,2)B =[ (12*x)/(L - 5)^3 - (2*(3*L + 15))/(L - 5)^3, (6*x)/(L - 5)^2 - (2*(2*L + 5))/(L - 5)^2, (2*(3*L + 15))/(L - 5)^3 -(12*x)/(L - 5)^3, (6*x)/(L - 5)^2 - (2*(L + 10))/(L - 5)^2]>> k=transpose(B)*(B);ke =[ 12/(L - 5)^3, 6/(L - 5)^2, -12/(L - 5)^3, 6/(L - 5)^2][ 6/(L - 5)^2, 4/(L - 5), -6/(L - 5)^2, 2/(L - 5)][ -12/(L - 5)^3, -6/(L - 5)^2, 12/(L - 5)^3, -6/(L - 5)^2][ 6/(L - 5)^2, 2/(L - 5), -6/(L - 5)^2, 4/(L - 5)]>> %Element1:E=4.0e11,I=bh^3/12=5.2e-7EI=4.0e11*5.2e-7EI =208000>> ke2=EI*subs(ke,'L',10)ke2 =19968 49920 -19968 4992049920 166400 -49920 83200-19968 -49920 19968 -4992049920 83200 -49920 166400>> %由此我们也得到了2号单元的刚度矩阵ke2>> %由于ke1,ke2都是在各自的局部坐标下得到的，所以我们必须把他们向整体坐标系做变换>> %局部坐标系想整体坐标系的转换>> T=eye(4,4) %定义坐标变换矩阵T =1 0 0 00 1 0 00 0 1 00 0 0 1>> %由于局部坐标系与整体坐标系的的夹角为零度，所以得到的T矩阵是一个4行4列的单位阵>> ke1=ke2ke1 =19968 49920 -19968 4992049920 166400 -49920 83200-19968 -49920 19968 -4992049920 83200 -49920 166400>> %由于运算问题，这里必须再次，定义ke1,而我们得到的ke2恰好等于之前的ke1 >> ke1=T*ke1*T';>> ke2=T*ke2*T';>> %系统分析F=[K]u>> %首先我们要在这里对整体刚度矩阵组集：直接法>> G1=...[1 0 0 0 0 00 1 0 0 0 00 0 1 0 0 00 0 0 1 0 0];>> G2=...[0 0 1 0 0 00 0 0 1 0 00 0 0 0 1 00 0 0 0 0 1];>> K1=G1'*ke1*G1K1 =19968 49920 -19968 49920 0 049920 166400 -49920 83200 0 0-19968 -49920 19968 -49920 0 049920 83200 -49920 166400 0 00 0 0 0 0 00 0 0 0 0 0 >> K2=G2'*ke2*G2K2 =0 0 0 0 0 00 0 0 0 0 00 0 19968 49920 -19968 499200 0 49920 166400 -49920 832000 0 -19968 -49920 19968 -499200 0 49920 83200 -49920 166400 >> K=K1+K2K =19968 49920 -19968 49920 0 049920 166400 -49920 83200 0 0 -19968 -49920 39936 0 -19968 4992049920 83200 0 332800 -49920 832000 0 -19968 -49920 19968 -499200 0 49920 83200 -49920 166400>> %引入约束条件>> %v1=0,xta1=0相当于>> K(1,:)=0;K(:,1)=0;>> KK =0 0 0 0 0 00 0 0 0 0 00 0 39936 0 -19968 499200 0 0 332800 -49920 832000 0 -19968 -49920 19968 -499200 0 49920 83200 -49920 166400 >>F=[0 0 0 0 -100 0]' %节点外载荷F =-100>>%求解系统方程，得到所有节点的位移>>%排除V1,与Xta1的影响>> KX=K(3:6,3:6)KX =39936 0 -19968 499200 332800 -49920 83200-19968 -49920 19968 -4992049920 83200 -49920 166400>> FX=F(3:6,1)FX =-100>> u=inv(KX)*FXu =-0.0501-0.0180-0.1603-0.0240>>%上述得到了2，3节点的挠曲与转角其中中间点（2）的挠曲与转角位：-0.1603 -0.0240右端点（3）的挠曲与转角位：-0.0501 -0.01805.参考文献：1.弹性力学与及有限元法基础教程韩清凯孙伟编著东北大学出版社2009.062.百度文库魏磊学号：200818932011.04~05。