基于matlab的有限元法分析平面应力应变问题刘刚.docx
- 文档编号:29898392
- 上传时间:2023-08-03
- 格式:DOCX
- 页数:20
- 大小:161.25KB
基于matlab的有限元法分析平面应力应变问题刘刚.docx
《基于matlab的有限元法分析平面应力应变问题刘刚.docx》由会员分享,可在线阅读,更多相关《基于matlab的有限元法分析平面应力应变问题刘刚.docx(20页珍藏版)》请在冰豆网上搜索。
基于matlab的有限元法分析平面应力应变问题刘刚
姓名:
刘刚学号:
15
平面应力应变分析有限元法
Abstruct:
本文通过对平面应力/应变问题的简要理论阐述,使读者对要分析的问题有大致的印象,然后结合两个实例,通过MATLAB软件的计算,将有限元分析平面应力/应变问题的过程形象的展示给读者,让人一目了然,快速了解有限元解决这类问题的方法和步骤!
一.基本理论
有限元法的基本思路和基本原则以结构力学中的位移法为基础,把复杂的结构或连续体看成有限个单元的组合,各单元彼此在节点出连接而组成整体。
把连续体分成有限个单元和节点,称为离散化。
先对单元进行特性分析,然后根据节点处的平衡和协调条件建立方程,综合后做整体分析。
这样一分一合,先离散再综合的过程,就是把复杂结构或连续体的计算问题转化简单单元分析与综合问题。
因此,一般的有限揭发包括三个主要步骤:
离散化单元分析整体分析。
二.用到的函数
1.LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p)
(KkIf)
(ku)
(kuA)
(ENUt)
三.实例
例1.考虑如图所示的受均布载荷作用的薄平板结构。
将平板离散化成两个线性三角元,假定E=200GPa,v=,t=0.025m,w=3000kN/m.
1.离散化
2.写出单元刚度矩阵
通过matlab的LinearTriangleElementStiffness函数,得到两个单元刚度矩阵
和
,每个矩阵都是6
6的。
>>E=210e6
E=
0
>>k1=LinearTriangleElementStiffness(E,NU,t,0,0,,,0,,1)
k1=
+006*
Columns1through5
00
00
00
00
Column6
>>NU=
NU=
>>t=
t=
>>k2=LinearTriangleElementStiffness(E,NU,t,0,0,,0,,,1)
k2=
+006*
Columns1through5
00
0
0
00
Column6
0
0
3.集成整体刚度矩阵8*8零矩阵
K=
00000000
00000000
00000000
00000000
00000000
00000000
00000000
00000000
>>K=LinearTriangleAssemble(K,k1,1,3,4)
K=
+006*
Columns1through5
0000
000
00000
00000
000
0000
00
00
Columns6through8
0
000
000
0
>>K=LinearTriangleAssemble(K,k1,1,2,3)
K=
+007*
000
000
0000
0000
0
0
00
00
4.引入边界条件.用上一步得到的整体刚度矩阵,可以得到该结构的方程组如下形式
本题的边界条件:
将边界条件带入,得到:
5.解方程
分解上述方程组,提取总体刚度矩阵K的第3-6行的第3-6列作为子矩阵
Matlab命令
>>k=K(3:
6,3:
6)
k=
+006*
0
0
>>f=[;0;;0]
f=
0
0
>>u=k\f
u=
*
现在可以清楚的看出,节点2的水平位移和垂直位移分别是0.7111m和0.1115m。
节点3的水平位移和垂直位移分别是0.6531m和0.0045m。
6.后处理
用matlab命令求出节点1和节点4的支反力以及每个单元的应力。
首先建立总体节点位移矢量U,
U=[0;0;u;0;0]
U=
*
0
0
0
0
>>F=K*U
F=
由以上知,节点1的水平反力和垂直反力分别是(指向左边)和(作用力方向向下),节点4的水平反力和垂直反力分别是(指向左边)和(作用力方向向下).满足力平衡条件。
接着,建立单元节点位移矢量
,然后调用matlab命令LinearTriangleElementStresses计算单元应力sigma1和sigma2
>>u1=[U
(1);U
(2);U(5);U(6);U(7);U(8)]
u1=
*
0
0
0
0
>>u2=[U
(1);U
(2);U(3);U(4);U(5);U(6)]
u2=
*
0
0
>>sigma1=LinearTriangleElementStresses(E,NU,,0,0,,,0,,1,u1)
sigma1=
+003*
>>sigma2=LinearTriangleElementStresses(E,NU,,0,0,,0,,,1,u2)
sigma2=
+003*
由以上可知,单元1的应力
,
,
。
单元2的应力是
。
显然,在x方向的应力(拉应力)接近于正确的值3MPa(拉应力)。
接着调用LinearTriangleElementStresses函数计算每个单元的主应力和主应力方向角。
>>s1=LinearTriangleElementPStresses(sigma1)
s1=
+003*
>>s2=LinearTriangleElementPStresses(sigma2)
s2=
+003*
,
主应力方向角
,
,
例2.考虑如图所示的由均匀分布载荷和集中载荷作用的薄平板结构。
将平板离散化成12个线性三角单元,如图4所示。
假定E=210GPa,v=,t=0.025m,w=100kN/m和P=。
1.离散化
2.写出单元刚度矩阵
>>E=201e6;
>>NU=;
>>t=;
>>k1=LinearTriangleElementStiffness(E,NU,t,0,,,,,,1);
>>k2=LinearTriangleElementStiffness(E,NU,t,0,,0,,,,1);
>>k3=LinearTriangleElementStiffness(E,NU,t,,,,,,,1);
>>k4=LinearTriangleElementStiffness(E,NU,t,,,0,,,,1);
>>k5=LinearTriangleElementStiffness(E,NU,t,0,,,,,,1);
>>k6=LinearTriangleElementStiffness(E,NU,t,0,,0,0,,,1);
>>k7=LinearTriangleElementStiffness(E,NU,t,,,,,,0,1);
>>k8=LinearTriangleElementStiffness(E,NU,t,,,0,0,,0,1);
>>k9=LinearTriangleElementStiffness(E,NU,t,025,,,0,,,1);
>>k10=LinearTriangleElementStiffness(E,NU,t,,,,,,,1);
>>k11=LinearTriangleElementStiffness(E,NU,t,,0,,0,,,1);
>>k12=LinearTriangleElementStiffness(E,NU,t,,,,0,,,1)
k1=
+006*
0
0
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)
运行得
+008*
Columns1through7
0
0
0
0
00
000
00
000
000000
000000
000000
000000
0000000
0000000
0000000
0000000
0000000
0000000
0000000
0000000
Columns8through14
000000
000000
00000
00000
0000
0000
00
0
000
00
0
0
00
00
0
0
00000
00000
00000
00000
0000000
0000000
Columns15through21
0000000
0000000
0000000
0000000
0000000
0000000
0000000
0000000
0
0
00000
00000
00000
00000
00
00
00
00
Column22
0
0
0
0
0
0
0
0
0
0
0
0
0
0
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=
5.解方程:
>>k=[K(3:
6,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:
22);K(15:
22,3:
6),K(15:
22,9:
12),K(15:
22,15:
22)];
K=
+008*
Columns1through8
00
00
0
0
000
000
000
000
000000
000000
000000
000000
00000000
00000000
00000000
00000000
00000000
00000000
00000000
00000000
Columns9through16
00000000
00000000
000000
000000
000000
000000
000
000
00
00
0
0
00
00
0000
0000
000000
000000
000000
000000
Columns17through22
000000
000000
000000
000000
000000
000000
000000
000000
00
00
000000
000000
000000
000000
00
00
>>f=[0;0;0;0;0;;0;0;0;0;0;0;0;0;0;0];
>>u=k\f
u=
*[]T
6.后处理
>>U=[0;0;u(1:
4);0;0;u(5:
8);0;0;u(9:
16)];
>>F=K*U
F=
[00]T
三.小结
通过这次练习,对有限元结合MATLAB软件解决平面应力/应变问题的方法和过程逐渐清晰,特别是在应用MATLAB软件的过程中学到了很多东西。
同时,通过此次课程设计,我对以前学过的有限元理论课程知识也有了进一步了解。
一开始还没有头绪,耗了半天也没正处东西来,后来吧课本拿来仔细看模仿学习,逐渐对软件有所熟悉,对解题思路更加清晰,我深刻的认识到自己动手自己找方法的重要性。
以后要加强这方面的能力。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 基于 matlab 有限元 分析 平面 应力 应变 问题