《有限元基础教程》_【MATLAB算例】基于节点四边形单元的矩形薄板分析(QuadDNode)
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
基于4节点四边形单元的矩形薄板分析(Quad2D4Node)
如图4-21所示的一个薄平板,在右端部受集中力F 作用,其中的参数为:
75110Pa,=1/3,=0.1m,=110N E t F μ=⨯⨯。基于MA TLAB 平台,按平面应力问题计算各
个节点位移、支座反力以及单元的应力。
(a) 问题描述 (b) 有限元分析模型
图4-21 右端部受集中力作用的薄平板
解答:对该问题进行有限元分析的过程如下。 (1)结构的离散化与编号
将结构离散为二个4节点矩形单元,单元编号及节点编号如图4-21(b)所示,连接关系见表4-6,节点的几何坐标见表4-7,载荷F 按静力等效原则向节点1,2移置。
表4-6 结构的单元连接关系
单元号 节点号 1 2
3 5 6
4 1 3 4 2
表4-7 节点的坐标
节点
节点坐标/m x y
1 2 3 4 5 6
2 1 2 0 1 1 1 0 0 1 0 0
节点位移列阵
[]1
1223344556
6T
u v u v u v u v u v u v =q (4-186)
节点外载列阵
000000000022
T
F F
⎡⎤
=--
⎢⎥⎣
⎦
F
(4-187)
约束的支反力列阵
556
60
0000000T
x y x y R R R R ⎡⎤=⎣⎦R (4-188)
总的节点载荷列阵
556
600000022
T
x y x y F F
R R R R ⎡
⎤
=+==-
-
⎢⎥⎣
⎦
P F R R (4-189)
其中,(
)55,x y R R 和()
66,x y R R 分别为节点5和节点6的两个方向的支反力。 (2)计算各单元的刚度矩阵(以国际标准单位)
首先在MATLAB 环境下,输入弹性模量E 、泊松比NU 、薄板厚度h 和平面应力问题性质指示参数ID ,然后针对单元1和单元2,分别两次调用函数Quad2D4Node_Stiffness ,就可以得到单元的刚度矩阵k1(8×8)和k2(8×8)。
>> E=1e7; >> NU=1/3; >> t=0.1; >> ID=1;
>> k1= Quad2D4Node_Stiffness(E,NU,t, 1,1, 0,1,0,0,1,0, ID); >>k2=Quad2D4Node_Stiffness(E,NU,t, 2,1,1,1,1,0,2,0, ID);
(3) 建立整体刚度方程
由于该结构共有6个节点,则总共的自由度数为12,因此,结构总的刚度矩阵为KK (12×12),先对KK 清零,然后两次调用函数Quad2D4Node_Assembly 进行刚度矩阵的组装。
>>KK = zeros(12,12);
>>KK= Quad2D4Node_Assembly(KK,k1, 3,5,6,4); >>KK= Quad2D4Node_Assembly(KK,k2, 1,3,4,2);
(4) 边界条件的处理及刚度方程求解
由图4-21(b)可以看出,节点5和节点6的两个方向的位移将为零,即
55660,0,0,0u v u v ====,因此,将针对节点1,2,3和4的位移进行求解。节点1,2,
3和4的位移将对应KK 矩阵中的前8行和前8列,则需从KK (12×12)中提出,置给k ,然后生成对应的载荷列阵p ,再采用高斯消去法进行求解。注意:MATLAB 中的反斜线符号“\”就是采用高斯消去法。
>>k=KK(1:8,1:8);
>>p=[0;-50000;0;-50000;0;0;0;0]; >>u=k\p
u = 0.8000 -2.5333 -0.8000 -2.5333 0.6000 -0.8667 -0.6000 -0.8667 [将列排成行]
由此可以看出,所求得的结果为:
11220.800 0, 2.533 3,0.800 0, 2.533 3,u v u v ==-=-=- 33440.600 0,0.866 7,0.600 0,0.866 7u v u v ==-=-=-。
(5)支反力的计算
在得到整个结构的节点位移后,由原整体刚度方程就可以计算出对应的支反力。先将
上面得到的位移结果与位移边界条件的节点位移进行组合(注意位置关系),可以得到整体的
位移列阵U(12×1),再代回原整体刚度方程,计算出所有的节点力P(12×1),按式(4-189)的对应关系就可以找到对应的支反力。
>>U=[u;0;0;0;0]; >>P=KK*U P = 1.0e+005 *
-0.0000 -0.5000 0.0000 -0.5000 0.0000 0.0000 [将列排成行] 0.0000 -0.0000 -2.0000 0.5000 2.0000 0.5000 [将列排成行]
由式(4-189)的对应关系,可以得到对应的支反力为
5566200 000,50 000,200 000,50 000x y x y R R R R =-===。
(6)各单元的应力计算
先从整体位移列阵U(12×1)中提取出单元的位移列阵,然后,调用计算单元应力的函数Quad2D4Node_Stress ,就可以得到各个单元的应力分量。
>>u1=[ U(5);U(6) ;U(9);U(10);U(11);U(12);U(7);U(8)]
u1 = 0.6000 -0.8667 0 0 0 0 -0.6000 -0.8667 [将列排成行]
>>stress1=Quad2D4Node_Stress(E,NU, 1,1,0,1,0,0,1,0, u1,ID) stress1 = 1.0e+005 *
0.0000 0.0000 -10.0000 [将列排成行]
>>u2=[ U(1);U(2) ;U(5);U(6); U(7);U(8);U(3);U(4)]
u2 = 0.8000 -2.5333 0.6000 -0.8667 [将列排成行] -0.6000 -0.8667 -0.8000 -2.5333 [将列排成行]
>>stress2=Quad2D4Node_Stress(E,NU, 2,1,1,1,1,0,2,0, u2,ID) stress2 = 1.0e+005 *
0.0000 0.0000 -10.0000 [将列排成行]
可以看出:计算得到的单元1的应力分量为0Pa,x σ=0Pa,y σ= 1 000 000Pa xy τ=-;单元2的应力分量为0Pa,x σ= 0Pa,y σ= 1 000 000Pa xy τ=-。