平面三角形单元有限元程序设计Word下载.docx
- 文档编号:20069447
- 上传时间:2023-01-16
- 格式:DOCX
- 页数:13
- 大小:226.32KB
平面三角形单元有限元程序设计Word下载.docx
《平面三角形单元有限元程序设计Word下载.docx》由会员分享,可在线阅读,更多相关《平面三角形单元有限元程序设计Word下载.docx(13页珍藏版)》请在冰豆网上搜索。
P
X
Y
刚度矩阵的集成
建立与总刚度矩阵等维数的空矩阵,已变单元刚度矩阵的集成。
由单元分析已知节点、单元的排布规律,继而通过循环计算求得每个单元对应的节点序号。
通过循环逐个计算:
(1)每个单元对应2种单元刚度矩阵中的哪一种;
(2)该单元对应总刚度矩阵的那几行哪几列
(3)将该单元的单元刚度矩阵加入总刚度矩阵的对应行列
循环又分为3层循环:
(1)最外层:
逐行计算
(2)中间层:
该行逐个计算
(3)最里层:
区分为第奇/偶数个计算
单元刚度的集成:
边界约束的处理:
划0置1法
适用:
这种方法适用于边界节点位移分量为已知(含为0)的各种约束。
做法:
(1)
将总刚矩阵〔K〕中相应于已知位移行主对角线元素置1,其他元素改为零;
同
时将载荷列阵{R}中相应元素用已知位移置换。
◎
这样,由该方程求得的此位移值一定等于已知量。
(2)
将〔K〕中已知位移相应的列的非主对角成元素也置0,以保持〔K〕的对称性。
当然,在已知位移分量不为零的情况下,这样做就改变了方程左端的数值,为
保证方程成立,须在方程右端减去已知位移对该方程的贡献——已知位移和相应总刚元素的乘积。
若约束为零位移约束时,此步则可省去。
特点:
经以上处理同样可以消除刚性位移(约束足够的前提下),去掉未知约束反力。
但这种方法不改变方程阶数,利于存贮。
(3)
不过,若是要求出约束反力,仍要重新计算各个划去的总刚元素。
程序如下:
变量说明
NNODE单元节点数
NPION总结点数
NELEM单元数
NVFIX受约束边界点数
FIXED约束信息数组
NFORCE节点力数
FORCE节点力数组
COORD结构节点坐标数组
LNODS单元定义数组
YOUNG弹性模量
POISS泊松比
THICK厚度
B单元应变矩阵(3*6)
D单元弹性矩阵(3*3)
S单元应力矩阵(3*6)
A单元面积
ESTIF单元刚度矩阵
ASTIF总体刚度矩阵
ASLOD总体荷载向量
ASDISP节点位移向量
ELEDISP单元节点位移向量
STRESS单元应力
%**********************************************************
%初始化
clear
formatshorte%设定输出类型
clear%清除内存变量
NELEM=36%单元个数(单元编码总数)
NPION=28%结点个数(结点编码总数)
NVFIX=2%受约束边界点数
NFORCE=1%结点荷载个数
YOUNG=2e11%弹性模量
POISS=0.25%泊松比
THICK=0.1%厚度
LNODS=[123;
245;
253;
356;
478;
485;
589;
596;
6910;
71112;
7128;
81213;
8139;
91314;
91410;
101415;
111617;
111712;
121718;
121813;
131819;
131914;
141920;
142015;
152021;
162223;
162317;
172324;
172418;
182425;
182519;
251926;
192620;
202627;
202721;
212728]%单元定义数组(单元结点号)
%相应为单元结点号(编码)、按逆时针顺序输入
COORD=[00;
-0.751.5;
0.751.5;
-1.53;
03;
1.53;
-2.254.5;
-0.754.5;
0.754.5;
2.254.5;
-36;
-1.56;
06;
1.56;
36;
-3.757.5;
-2.257.5;
-0.757.5;
0.757.5;
2.257.5;
3.757.5;
-4.59;
-39;
-1.59;
09;
1.59;
39;
4.59]%结点坐标数组
%坐标:
x,y坐标(共NPOIN组)
FORCE=[10-15]%结点力数组(受力结点编号,x方向,y方向)
FIXED=[2211;
2811]%约束信息(约束点,x约束,y约束)
%有约束为1,无约束为0
%生成单元刚度矩阵并组成总体刚度矩阵
ASTIF=zeros(2*NPION,2*NPION);
%生成特定大小总体刚度矩阵并置0
fori=1:
NELEM
%生成弹性矩阵D
D=[1POISS0;
POISS10;
00(1-POISS)/2]*YOUNG/(1-POISS^2)
%计算当前单元的面积
A=-det([1COORD(LNODS(i,1),1)COORD(LNODS(i,1),2);
1COORD(LNODS(i,2),1)COORD(LNODS(i,2),2);
1COORD(LNODS(i,3),1)COORD(LNODS(i,3),2)])/2
%生成应变矩阵B
forj=0:
2
b(j+1)=COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);
c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1);
end
B=[b
(1)0b
(2)0b(3)0;
0c
(1)0c
(2)0c(3);
c
(1)b
(1)c
(2)b
(2)c(3)b(3)]/(2*A);
B1(:
:
i)=B;
%**********************************************************
%求应力矩阵S=D*B
S=D*B;
ESTIF=B'
*S*THICK*A;
%求解单元刚度矩阵
a=LNODS(i,:
);
%临时向量,用来记录当前单元的节点编号
forj=1:
3
fork=1:
ASTIF((a(j)*2-1):
a(j)*2,(a(k)*2-1):
a(k)*2)=ASTIF((a(j)*2-1):
a(k)*2)+ESTIF(j*2-1:
j*2,k*2-1:
k*2);
%根据节点编号对应关系将单元刚度分块叠加到总刚
%度矩阵中
%将约束信息加入总体刚度矩阵(对角元素改一法)
fori=1:
NVFIX
ifFIXED(i,2)==1
ASTIF(:
(FIXED(i,1)*2-1))=0;
%一列为零
ASTIF((FIXED(i,1)*2-1),:
)=0;
%一行为零
ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1;
%对角元素为1
ifFIXED(i,3)==1
ASTIF(:
FIXED(i,1)*2)=0;
ASTIF(FIXED(i,1)*2,:
ASTIF(FIXED(i,1)*2,FIXED(i,1)*2)=1;
%生成荷载向量
ASLOD(1:
2*NPION)=0;
%总体荷载向量置零
NFORCE
ASLOD((FORCE(i,1)*2-1):
FORCE(i,1)*2)=FORCE(i,2:
3);
%求解内力
ASDISP=ASTIF\ASLOD'
%计算节点位移向量
ELEDISP(1:
6)=0;
%当前单元节点位移向量
ELEDISP(j*2-1:
j*2)=ASDISP(LNODS(i,j)*2-1:
LNODS(i,j)*2);
%取出当前单元的节点位移向量
i
STRESS=D*B1(:
:
i)*ELEDISP'
%求内力
(程序计算结果和有限元软件得出的结果稍有偏差,可能是程序某些地方数据输入时出了问题,还在寻找具体原因)
2、有限元软件分析
建模
设置材料参数
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 平面三角形 单元 有限元 程序设计
![提示](https://static.bdocx.com/images/bang_tan.gif)