基于MATLAB的有限元法分析平面应力应变问题刘刚Word文档格式.docx
《基于MATLAB的有限元法分析平面应力应变问题刘刚Word文档格式.docx》由会员分享,可在线阅读,更多相关《基于MATLAB的有限元法分析平面应力应变问题刘刚Word文档格式.docx(17页珍藏版)》请在冰豆网上搜索。
![基于MATLAB的有限元法分析平面应力应变问题刘刚Word文档格式.docx](https://file1.bdocx.com/fileroot1/2022-11/27/821dc090-2098-4c3f-a817-c333764df2a8/821dc090-2098-4c3f-a817-c333764df2a81.gif)
E=
210000000
k1=LinearTriangleElementStiffness(E,NU,t,0,0,0.5,0.25,0,0.25,1)
k1=
1.0e+006*
Columns1through5
2.019200-1.0096-2.0192
05.7692-0.865400.8654
0-0.86541.44230-1.4423
-1.0096000.50481.0096
-2.01920.8654-1.44231.00963.4615
1.0096-5.76920.8654-0.5048-1.8750
Column6
1.0096
-5.7692
0.8654
-0.5048
-1.8750
6.2740
NU=0.3
NU=
0.3000
t=0.025
t=
0.0250
k2=LinearTriangleElementStiffness(E,NU,t,0,0,0.5,0,0.5,0.25,1)
k2=
1.44230-1.44230.86540
00.50481.0096-0.5048-1.0096
-1.44231.00963.4615-1.8750-2.0192
0.8654-0.5048-1.87506.27401.0096
0-1.0096-2.01921.00962.0192
-0.865400.8654-5.76920
-0.8654
0
5.7692
3.集成整体刚度矩阵8*8零矩阵
K=
00000000
K=LinearTriangleAssemble(K,k1,1,3,4)
2.01920000
05.769200-0.8654
00000
0-0.8654001.4423
-1.00960000
-2.01920.865400-1.4423
1.0096-5.7692000.8654
Columns6through8
-1.0096-2.01921.0096
00.8654-5.7692
000
0-1.44230.8654
0.50481.0096-0.5048
1.00963.4615-1.8750
-0.5048-1.87506.2740
K=LinearTriangleAssemble(K,k1,1,2,3)
1.0e+007*
0.403800-0.1010-0.20190-0.20190.1010
01.1538-0.086500-0.57690.0865-0.5769
0-0.08650.14420-0.14420.086500
-0.1010000.05050.1010-0.050500
-0.20190-0.14420.10100.4904-0.1875-0.14420.0865
0-0.57690.0865-0.0505-0.18750.67790.1010-0.0505
-0.20190.086500-0.14420.10100.3462-0.1875
0.1010-0.5769000.0865-0.0505-0.18750.6274
4.引入边界条件.用上一步得到的整体刚度矩阵,可以得到该结构的方程组如下形式
本题的边界条件:
将边界条件带入,得到:
5.解方程
分解上述方程组,提取总体刚度矩阵K的第3-6行的第3-6列作为子矩阵
Matlab命令
k=K(3:
6,3:
6)
k=
3.4615-1.8750-2.01920.8654
-1.87506.27401.0096-5.7692
-2.01921.00963.46150
0.8654-5.769206.2740
f=[9.375;
0;
9.375;
0]
f=
9.3750
u=k\f
u=
1.0e-005*
0.7111
0.1115
0.6531
0.0045
现在可以清楚的看出,节点2的水平位移和垂直位移分别是0.7111m和0.1115m。
节点3的水平位移和垂直位移分别是0.6531m和0.0045m。
6.后处理
用matlab命令求出节点1和节点4的支反力以及每个单元的应力。
首先建立总体节点位移矢量U,
U=[0;
u;
U=
F=K*U
F=
-9.3750
-5.6295
0.0000
5.6295
由以上知,节点1的水平反力和垂直反力分别是9.375kn(指向左边)和5.6295kn(作用力方向向下),节点4的水平反力和垂直反力分别是9.375kn(指向左边)和5.6295kn(作用力方向向下).满足力平衡条件。
接着,建立单元节点位移矢量
,然后调用matlab命令LinearTriangleElementStresses计算单元应力sigma1和sigma2
u1=[U
(1);
U
(2);
U(5);
U(6);
U(7);
U(8)]
u1=
u2=[U
(1);
U(3);
U(4);
U(6)]
u2=
sigma1=LinearTriangleElementStresses(E,NU,0.025,0,0,0.5,0.25,0,0.25,1,u1)
sigma1=
1.0e+003*
3.0144
0.9043
0.0072
sigma2=LinearTriangleElementStresses(E,NU,0.025,0,0,0.5,0,0.5,0.25,1,u2)
sigma2=
2.9856
-0.0036
-0.0072
由以上可知,单元1的应力
,
。
单元2的应力是
。
显然,在x方向的应力(拉应力)接近于正确的值3MPa(拉应力)。
接着调用LinearTriangleElementStresses函数计算每个单元的主应力和主应力方向角。
s1=LinearTriangleElementPStresses(sigma1)
s1=
0.0002
s2=LinearTriangleElementPStresses(sigma2)
s2=
-0.0001
主应力方向角
例2.考虑如图3.1所示的由均匀分布载荷和集中载荷作用的薄平板结构。
将平板离散化成12个线性三角单元,如图4所示。
假定E=210GPa,v=0.3,t=0.025m,w=100kN/m和P=12.5kN。
1.离散化
2.写出单元刚度矩阵
E=201e6;
NU=0.3;
t=0.025;
k1=LinearTriangleElementStiffness(E,NU,t,0,0.5,0.125,0.375,0.25,0.5,1);
k2=LinearTriangleElementStiffness(E,NU,t,0,0.5,0,0.25,0.125,0.375,1);
k3=LinearTriangleElementStiffness(E,NU,t,0.125,0.375,0.25,0.25,0.25,0.5,1);
k4=LinearTriangleElementStiffness(E,NU,t,0.125,0.375,0,0.25,0.25,0.25,1);
k5=LinearTriangleElementStiffness(E,NU,t,0,0.25,0.125,0.125,0.25,0.25,1);
k6=LinearTriangleElementStiffness(E,NU,t,0,0.25,0,0,0.125,0.125,1);
k7=LinearTriangleElementStiffness(E,NU,t,0.25,0.25,0.125,0.125,0.25,0,1);
k8=LinearTriangleElementStiffness(E,NU,t,0.125,0.125,0,0,0.25,0,1);
k9=LinearTriangleElementStiffness(E,NU,t,025,0.25,0.25,0,0.375,0.125,1);
k10=LinearTriangleElementStiffness(E,NU,t,0.25,0.25,0.375,0.125,0.5,0.25,1);
k11=LinearTriangleElementStiffness(E,NU,t,0.25,0,0.5,0,0.375,0.125,1);
k12=LinearTriangleElementStiffness(E,NU,t,0.375,0.125,0.5,0,0.5,0.25,1)
1.8637-0.8973-0.96630.8283-0.89730.0690
-0.89731.86370.9663-2.7610-0.06900.8973
-0.96630.96631.93270-0.9663-0.9663
0.8283-2.761005.5220-0.8283-2.7610
-0.8973-0.0690-0.9663-0.82831.86370.8973
0.06900.8973-0.9663-2.76100.89731.8637
3.集成整体刚度矩阵:
K=zero(22,22);
K=LinearTriangleAssemble(K,k1,1,3,2);
K=LinearTriangleAssemble(K,k2,1,4,3);
K=LinearTriangleAssemble(K,k3,3,5,2);
K=LinearTriangleAssemble(K,k4,3,4,5);
K=LinearTriangleAssemble(K,k5,4,6,5);
K=LinearTriangleAssemble(K,k6,4,7,6);
K=LinearTriangleAssemble(K,k7,5,6,8);
K=LinearTriangleAssemble(K,k8,6,7,8);
K=LinearTriangleAssemble(K,k9,5,8,9);
K=LinearTriangleAssemble(K,k10,5,9,10);
K=LinearTriangleAssemble(K,k11,8,11,9);
K=LinearTriangleAssemble(K,k12,9,11,10)
运行得
1.0e+008*
Columns1through7
0.0389-0.0187-0.00940.0007-0.03890.01870.0094
-0.01870.0389-0.00070.00940.0187-0.03890.0007
-0.0094-0.00070.03890.0187-0.0389-0.01870
0.00070.00940.01870.0389-0.0187-0.03890
-0.03890.0187-0.0389-0.01870.15580-0.0389
0.0187-0.0389-0.0187-0.038900.1558-0.0187
0.00940.000700-0.0389-0.01870.0779
-0.0007-0.009400-0.0187-0.03890
000.0094-0.0007-0.03890.0187-0.0187
000.0007-0.00940.0187-0.03890
000000-0.0389
0000000.0187
0000000.0094
000000-0.0007
0000000
Columns8through14
-0.0007000000
-0.0094000000
00.00940.00070000
0-0.0007-0.00940000
-0.0187-0.03890.01870000
-0.03890.0187-0.03890000
0-0.01870-0.03890.01870.0094-0.0007
0.077900.01870.0187-0.03890.0007-0.0094
00.0972-0.0093-0.0389-0.018700
0.0187-0.00930.0972-0.0187-0.038900
0.0187-0.0389-0.01870.15580-0.0389-0.0187
-0.0389-0.0187-0.038900.1558-0.0187-0.0389
0.000700-0.0389-0.01870.03890.0187
-0.009400-0.0187-0.03890.01870.0389
0-0.00090.0095-0.03890.0187-0.00940.0007
00.0095-0.03840.0187-0.0389-0.00070.0094
00.0004-0.00020000
0-0.00020.00040000
0-0.0094-0.00070000
00.00070.00940000
0000000
Columns15through21
-0.00090.00950.0004-0.0002-0.00940.00070
0.0095-0.0384-0.00020.0004-0.00070.00940
-0.03890.018700000
0.0187-0.038900000
-0.0094-0.000700000
0.00070.009400000
-1.94080.00951.9994-0.037700-0.0094
0.0095-5.6533-0.03775.7119000.0007
1.9994-0.0377-1.92190.0379-0.0389-0.0187-0.0389
-0.03775.71190.0379-5.6344-0.0187-0.03890.0187
00-0.0389-0.01870.03890.01870.0094
00-0.0187-0.03890.01870.0389-0.0007
-0.00940.0007-0.03890.01870.0094-0.00070.0389
-0.00070.00940.0187-0.03890.0007-0.0094-0.0187
Column22
-0.0007
0.0094
0.0187
-0.0389
0.0007
-0.0094
-0.0187
0.03890.0007
0.0389
4.引入边界条件:
U1x=U1y=U4x=U4y=U7x=U7y=0
F2x=F2y=F3x=F3y=F6x=F6y=F8x=F8y=F9x=F9y=F10x=F10y=F11x=F11y=0
F5x=0,F5y=-12.5
5.解方程:
k=[K(3:
6),K(3:
6,9:
12),K(3:
6,15:
22);
K(9:
12,3:
6),K(9:
12,9:
12),K(9:
12,15:
K(15:
22,3:
6),K(15:
22,9:
12),K(15:
22,15:
22)];
1.0e+008*
Columns1through8
0.0389-0.0187-0.00940.0007-0.03890.01870.0094-0.0007
-0.01870.0389-0.00070.00940.0187-0.03890.0007-0.0094
-0.0094-0.00070.03890.0187-0.0389-0.018700
0.00070.00940.01870.0389-0.0187-0.038900
-0.03890.0187-0.0389-0.01870.15580-0.0389-0.0187
0.0187-0.0389-0.0187-0.038900.1558-0.0187-0.0389
0.00940.000700-0.0389-0.01870.07790
-0.0007-0.009400-0.0187-0.038900.0779
000.0094-0.0007-0.03890.0187-0.01870
000.0007-0.00940.0187-0.038900.0187
000000-0.03890.0187
0000000.0187-0.0389
0000000.00940.0007
000000-0.0007-0.0094
Columns9through16
0.00940.0007000000
-0.0007-0.0094000000
-0.03890.0187000000
0.0187-0.0389000000
-0.01870-0.03890.01870.0094-0.000700