梁单元有限元法分析
关键词：梁单元 有限元 分析
1. 摘要：二维平面梁单元是梁单元中最简单的单元之一，这次作业旨在学习如何运用有限元分析法分析梁单元。
2. 目的：运用MATLAB 软件分析二维梁单元。
3. 题目：设一方形的截面梁，截面每边长为5cm, 长度为10m, 在左端约束固定，在右端施以一个沿y 方向的集中力ω=100N,求其挠度与转角。
3. 建立有限元分析模型： （1）. 结构离散化：
单元的选择：由于为悬臂梁，且横向的长度远远小于轴向的长度，所以在这选择平面梁单元； 单元的数量：将这个梁从中间划分为两个单元；
建立坐标系，坐标系包括结构的整体坐标系与单元的局部坐标系； （2. ）建立平面梁单元的位移模式： 建立整体坐标系：
建立一个有两个单元组成的模型，由于X 方向的位移U1,U2,U3太小所以我们略去这三个自由度的变化； 节点坐标码：
单元编码：
同样出1号单元，建立局部坐标系：
4. 具体的MATLAB 求解过程与结果： >> clear x1=0; x2=sym('L'); x=sym('x'); j=0:3; v=x.^j v =
[ 1, x, x^2, x^3]
>> %计算形函数矩阵 m=...
[1 x1 x1^2 x1^3 0 1 2*x1 3*x1^2 1 x2 x2^2 x2^3 0 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)） %由最小势能原理导出刚度矩阵ke
B=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-7 EI=4.0e11*5.2e-7 EI =
208000
>> ke1=EI*subs(ke,'L',5) ke1 =
19968 49920 -19968 49920 49920 166400 -49920 83200 -19968 -49920 19968 -49920 49920 83200 -49920 166400
>> %由上市我们就计算出了1号单元刚度矩阵ke, 由于分成两个单元，所以L=10/2=5 >> %同理，我们运用上述办法得到2号单元的刚度矩阵ke2 >> clear x2=5;
x=sym('x'); j=0:3; v=x.^j v =
[ 1, x, x^2, x^3] >> m=... [1 x2 x2^2 x2^3 0 1 2*x2 3*x2^2 1 x3 x3^2 x3^3 0 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*mm N =
[ (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-7 EI=4.0e11*5.2e-7 EI =
208000
>> ke2=EI*subs(ke,'L',10) ke2 =
19968 49920 -19968 49920 49920 166400 -49920 83200 -19968 -49920 19968 -49920 49920 83200 -49920 166400
>> %由此我们也得到了2号单元的刚度矩阵ke2
>> %由于ke1,ke2都是在各自的局部坐标下得到的，所以我们必须把他们向整体坐标系做变换 >> %局部坐标系想整体坐标系的转换 >> T=eye(4,4) %定义坐标变换矩阵 T =
1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1
>> %由于局部坐标系与整体坐标系的的夹角为零度，所以得到的T 矩阵是一个4行4列的单位阵 >> ke1=ke2 ke1 =
19968 49920 -19968 49920 49920 166400 -49920 83200
-19968 -49920 19968 -49920 49920 83200 -49920 166400
>> %由于运算问题，这里必须再次，定义ke1, 而我们得到的ke2恰好等于之前的ke1 >> ke1=T*ke1*T'; >> ke2=T*ke2*T'; >> %系统分析F=[K]u
>> %首先我们要在这里对整体刚度矩阵组集：直接法 >> G1=... [1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0]; >> G2=... [0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1]; >> K1=G1'*ke1*G1 K1 =
19968 49920 49920 166400 -19968 -49920 49920 83200 0 0 0 0
>> K2=G2'*ke2*G2 K2 =
0 0 0 0 0 0 0 0 0 0 0 0
>> K=K1+K2 K =
-19968 49920 -49920 83200 19968 -49920 -49920 166400 0 0 0 0 0 0 0 0 19968 49920 49920 166400 -19968 -49920 49920 83200 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -19968 49920 -49920 83200 19968 -49920 -49920 166400
19968 49920 -19968 49920 0 0 49920 166400 -49920 83200 0 0 -19968 -49920 39936 0 -19968 49920 49920 83200 0 332800 -49920 83200 0 0 -19968 -49920 19968 -49920 0 0 49920 83200 -49920 166400
>> %引入约束条件 >> %v1=0,xta1=0相当于 >> K(1,:)=0;K(:,1)=0; >> K K =
0 0 0 0 0 0 0 0 39936 0 0 0 0 0 -19968 0 0 49920
>>F=[0 0 0 0 -100 0]' %节点外载荷 F = 0 0 0 0 -100 0
>>%求解系统方程，得到所有节点的位移 >>%排除V1, 与Xta1的影响 >> KX=K(3:6,3:6) KX =
39936 0 -19968 0 332800 -49920 -19968 -49920 19968 49920 83200 -49920
>> FX=F(3:6,1)
0 0 0 0 0 -19968 332800 -49920 -49920 19968 83200 -49920 49920 83200 -49920 166400 0 0 49920 83200 -49920 166400
FX = 0 0 -100 0 >> u=inv(KX)*FX u =
-0.0501 -0.0180 -0.1603 -0.0240
>>%上述得到了2，3节点的挠曲与转角
其中中间点（2）的挠曲与转角位：-0.1603 -0.0240 右端点（3）的挠曲与转角位：-0.0501 -0.0180 5. 参考文献：
1. 弹性力学与及有限元法基础教程 韩清凯 孙伟 编著 2. 百度文库
东北大学出版社 魏磊
学号：20081893 2011.04~05 2009.06
梁单元有限元法分析
关键词：梁单元 有限元 分析
1. 摘要：二维平面梁单元是梁单元中最简单的单元之一，这次作业旨在学习如何运用有限元分析法分析梁单元。
2. 目的：运用MATLAB 软件分析二维梁单元。
3. 题目：设一方形的截面梁，截面每边长为5cm, 长度为10m, 在左端约束固定，在右端施以一个沿y 方向的集中力ω=100N,求其挠度与转角。
3. 建立有限元分析模型： （1）. 结构离散化：
单元的选择：由于为悬臂梁，且横向的长度远远小于轴向的长度，所以在这选择平面梁单元； 单元的数量：将这个梁从中间划分为两个单元；
建立坐标系，坐标系包括结构的整体坐标系与单元的局部坐标系； （2. ）建立平面梁单元的位移模式： 建立整体坐标系：
建立一个有两个单元组成的模型，由于X 方向的位移U1,U2,U3太小所以我们略去这三个自由度的变化； 节点坐标码：
单元编码：
同样出1号单元，建立局部坐标系：
4. 具体的MATLAB 求解过程与结果： >> clear x1=0; x2=sym('L'); x=sym('x'); j=0:3; v=x.^j v =
[ 1, x, x^2, x^3]
>> %计算形函数矩阵 m=...
[1 x1 x1^2 x1^3 0 1 2*x1 3*x1^2 1 x2 x2^2 x2^3 0 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)） %由最小势能原理导出刚度矩阵ke
B=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-7 EI=4.0e11*5.2e-7 EI =
208000
>> ke1=EI*subs(ke,'L',5) ke1 =
19968 49920 -19968 49920 49920 166400 -49920 83200 -19968 -49920 19968 -49920 49920 83200 -49920 166400
>> %由上市我们就计算出了1号单元刚度矩阵ke, 由于分成两个单元，所以L=10/2=5 >> %同理，我们运用上述办法得到2号单元的刚度矩阵ke2 >> clear x2=5;
x=sym('x'); j=0:3; v=x.^j v =
[ 1, x, x^2, x^3] >> m=... [1 x2 x2^2 x2^3 0 1 2*x2 3*x2^2 1 x3 x3^2 x3^3 0 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*mm N =
[ (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-7 EI=4.0e11*5.2e-7 EI =
208000
>> ke2=EI*subs(ke,'L',10) ke2 =
19968 49920 -19968 49920 49920 166400 -49920 83200 -19968 -49920 19968 -49920 49920 83200 -49920 166400
>> %由此我们也得到了2号单元的刚度矩阵ke2
>> %由于ke1,ke2都是在各自的局部坐标下得到的，所以我们必须把他们向整体坐标系做变换 >> %局部坐标系想整体坐标系的转换 >> T=eye(4,4) %定义坐标变换矩阵 T =
1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1
>> %由于局部坐标系与整体坐标系的的夹角为零度，所以得到的T 矩阵是一个4行4列的单位阵 >> ke1=ke2 ke1 =
19968 49920 -19968 49920 49920 166400 -49920 83200
-19968 -49920 19968 -49920 49920 83200 -49920 166400
>> %由于运算问题，这里必须再次，定义ke1, 而我们得到的ke2恰好等于之前的ke1 >> ke1=T*ke1*T'; >> ke2=T*ke2*T'; >> %系统分析F=[K]u
>> %首先我们要在这里对整体刚度矩阵组集：直接法 >> G1=... [1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0]; >> G2=... [0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1]; >> K1=G1'*ke1*G1 K1 =
19968 49920 49920 166400 -19968 -49920 49920 83200 0 0 0 0
>> K2=G2'*ke2*G2 K2 =
0 0 0 0 0 0 0 0 0 0 0 0
>> K=K1+K2 K =
-19968 49920 -49920 83200 19968 -49920 -49920 166400 0 0 0 0 0 0 0 0 19968 49920 49920 166400 -19968 -49920 49920 83200 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -19968 49920 -49920 83200 19968 -49920 -49920 166400
19968 49920 -19968 49920 0 0 49920 166400 -49920 83200 0 0 -19968 -49920 39936 0 -19968 49920 49920 83200 0 332800 -49920 83200 0 0 -19968 -49920 19968 -49920 0 0 49920 83200 -49920 166400
>> %引入约束条件 >> %v1=0,xta1=0相当于 >> K(1,:)=0;K(:,1)=0; >> K K =
0 0 0 0 0 0 0 0 39936 0 0 0 0 0 -19968 0 0 49920
>>F=[0 0 0 0 -100 0]' %节点外载荷 F = 0 0 0 0 -100 0
>>%求解系统方程，得到所有节点的位移 >>%排除V1, 与Xta1的影响 >> KX=K(3:6,3:6) KX =
39936 0 -19968 0 332800 -49920 -19968 -49920 19968 49920 83200 -49920
>> FX=F(3:6,1)
0 0 0 0 0 -19968 332800 -49920 -49920 19968 83200 -49920 49920 83200 -49920 166400 0 0 49920 83200 -49920 166400
FX = 0 0 -100 0 >> u=inv(KX)*FX u =
-0.0501 -0.0180 -0.1603 -0.0240
>>%上述得到了2，3节点的挠曲与转角
其中中间点（2）的挠曲与转角位：-0.1603 -0.0240 右端点（3）的挠曲与转角位：-0.0501 -0.0180 5. 参考文献：
1. 弹性力学与及有限元法基础教程 韩清凯 孙伟 编著 2. 百度文库
东北大学出版社 魏磊
学号：20081893 2011.04~05 2009.06