文档介绍：一、题目
如图1所示,一个厚度均匀的三角形薄板,在顶点作用沿板厚方向均匀分布的竖向载荷。已知:P=150N/m,E=200GPa,=0、25,t=0、1m,忽略自重。试计算薄板的位移及应力分布。
要求:
编写有限元计算机程序,计算节点位移及单元应力。(划分三角形单元,单元数不得少于30个);
采用有限元软件分析该问题(有限元软件网格与程序设计网格必须一致),详细给出有限元软件每一步的操作过程,并将结果与程序计算结果进行对比(任选取三个点,对比位移值);
提交程序编写过程的详细报告及计算机程序;
所有同学参加答辩,并演示有限元计算程序。
有限元法中三节点三角形分析结构的步骤如下:
1)整理原始数据,如材料性质、荷载条件、约束条件等,离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2)单元分析,建立单元刚度矩阵。
3)整体分析,建立总刚矩阵。
4)建立整体结构的等效节点荷载与总荷载矩阵
5)边界条件处理。
6)解方程,求出节点位移。
7)求出各单元的单元应力。
8)计算结果整理。
一、程序设计
网格划分
如图,将薄板如图划分为6行,并建立坐标系,则
P
P
X
X
Y
Y
单元编号
节点编号
刚度矩阵的集成
建立与总刚度矩阵等维数的空矩阵,已变单元刚度矩阵的集成。
由单元分析已知节点、单元的排布规律,继而通过循环计算求得每个单元对应的节点序号。
通过循环逐个计算:(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
format short e %设定输出类型
clear %清除内存变量
NELEM=36 %单元个数(单元编码总数)
NPION=28 %结点个数(结点编码总数)
NVFIX=2