平面三角形单元有限元程序设计.docx
《平面三角形单元有限元程序设计.docx》由会员分享，可在线阅读，更多相关《平面三角形单元有限元程序设计.docx（13页珍藏版）》请在冰点文库上搜索。
平面三角形单元有限元程序设计
一、题目
如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖向载荷。
已知：
P=150N/m，E=200GPa，
=0.25，t=0.1m，忽略自重。
试计算薄板的位移及应力分布。
要求：
1.编写有限元计算机程序，计算节点位移及单元应力。
（划分三角形单元，单元数不得少于30个）；
2.采用有限元软件分析该问题（有限元软件网格与程序设计网格必须一致），详细给出有限元软件每一步的操作过程，并将结果与程序计算结果进行对比（任选取三个点，对比位移值）；
3.提交程序编写过程的详细报告及计算机程序；
4.所有同学参加答辩，并演示有限元计算程序。
有限元法中三节点三角形分析结构的步骤如下：
1）整理原始数据，如材料性质、荷载条件、约束条件等，离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2）单元分析，建立单元刚度矩阵。
3）整体分析，建立总刚矩阵。
4）建立整体结构的等效节点荷载和总荷载矩阵
5）边界条件处理。
6）解方程，求出节点位移。
7）求出各单元的单元应力。
8）计算结果整理。
一、程序设计
网格划分
如图，将薄板如图划分为6行，并建立坐标系，则
X
Y
P
X
Y
P
节点编号
单元编号
刚度矩阵的集成
建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。
由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。
通过循环逐个计算：
（1）每个单元对应2种单元刚度矩阵中的哪一种；
（2）该单元对应总刚度矩阵的那几行哪几列
（3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列
循环又分为3层循环：
（1）最外层：
逐行计算
（2）中间层：
该行逐个计算
（3）最里层：
区分为第奇/偶数个计算
单元刚度的集成：
边界约束的处理：
划0置1法
适用：
这种方法适用于边界节点位移分量为已知（含为0）的各种约束。
做法：
（1） 将总刚矩阵〔K〕中相应于已知位移行主对角线元素置1，其他元素改为零；同
时将载荷列阵｛R｝中相应元素用已知位移置换。
◎ 这样，由该方程求得的此位移值一定等于已知量。
（2） 将〔K〕中已知位移相应的列的非主对角成元素也置0，以保持〔K〕的对称性。
◎ 当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为
保证方程成立，须在方程右端减去已知位移对该方程的贡献——已知位移和相应总刚元素的乘积。
◎ 若约束为零位移约束时，此步则可省去。
特点：
（1） 经以上处理同样可以消除刚性位移（约束足够的前提下），去掉未知约束反力。
（2） 但这种方法不改变方程阶数，利于存贮。
（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:
3
ASTIF（（a（j）*2-1）:
a（j）*2,（a（k）*2-1）:
a（k）*2）=ASTIF（（a（j）*2-1）:
a（j）*2,（a（k）*2-1）:
a（k）*2）+ESTIF（j*2-1:
j*2,k*2-1:
k*2）;
%根据节点编号对应关系将单元刚度分块叠加到总刚
%度矩阵中
end
end
end
%**********************************************************
%将约束信息加入总体刚度矩阵（对角元素改一法）
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
end
%**********************************************************
%生成单元刚度矩阵并组成总体刚度矩阵
%**********************************************************
ifFIXED（i,3）==1
ASTIF（:
FIXED（i,1）*2）=0;%一列为零
ASTIF（FIXED（i,1）*2,:
）=0;%一行为零
ASTIF（FIXED（i,1）*2,FIXED（i,1）*2）=1;%对角元素为1
end
end
%**********************************************************
%生成荷载向量
ASLOD（1:
2*NPION）=0;%总体荷载向量置零
fori=1:
NFORCE
ASLOD（（FORCE（i,1）*2-1）:
FORCE（i,1）*2）=FORCE（i,2:
3）;
end
%**********************************************************
%求解内力
ASDISP=ASTIF\ASLOD'%计算节点位移向量
ELEDISP（1:
6）=0;%当前单元节点位移向量
fori=1:
NELEM
forj=1:
3
ELEDISP（j*2-1:
j*2）=ASDISP（LNODS（i,j）*2-1:
LNODS（i,j）*2）;
%取出当前单元的节点位移向量
end
i
STRESS=D*B1（:
:
i）*ELEDISP'%求内力
end
（程序计算结果和有限元软件得出的结果稍有偏差，可能是程序某些地方数据输入时出了问题，还在寻找具体原因）
2、有限元软件分析
建模
设置材料参数
边界约束
添加载荷
网格划分
升级会员 