基于MATLAB的有限元法分析平面应力应变问题刘刚Word文档格式.docx
《基于MATLAB的有限元法分析平面应力应变问题刘刚Word文档格式.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的有限元法分析平面应力应变问题刘刚Word文档格式.docx（17页珍藏版）》请在冰豆网上搜索。
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
升级会员 