平面三角形单元有限元程序的设计说明

. .

P

9 m 9 m

一、题目

如图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行，并建立坐标系，则

刚度矩阵的集成

建立与总刚度矩阵等维数的空矩阵，已变单元刚度矩阵的集成。 由单元分析已知节点、单元的排布规律，继而通过循环计算求得每个单元对应的节点序号。

通过循环逐个计算：（1）每个单元对应2种单元刚度矩阵中的哪一种；

（2）该单元对应总刚度矩阵的那几行哪几列

（3）将该单元的单元刚度矩阵加入总刚度矩阵的对应行列

循环又分为3层循环：（1）最外层：逐行计算

（2）中间层：该行逐个计算

（3）最里层：区分为第 奇/偶 数个计算

单元刚度的集成：[

][][][][][]'

'''''215656665656266256561661e Z e e e Z e Z e e e e k k k K k k k k k k +⋯++=⇓

=⇒==⇒==⇒=⨯⨯⨯⨯⨯⨯

边界约束的处理：划0置1法 X Y P X Y P 节点编号 单元编号

适用：这种方法适用于边界节点位移分量为已知(含为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 %受约束边界点数

NFORCE=1 %结点荷载个数

YOUNG=2e11 %弹性模量

POISS=0.25 %泊松比

THICK=0.1 %厚度

LNODS=[1 2 3;2 4 5;2 5 3;3 5 6;

4 7 8;4 8 5;

5 8 9;5 9 6;

6 9 10;

7 11 12;7 12 8;

8 12 13;

8 13 9;9 13 14;9 14 10;10 14 15;

11 16 17;11 17 12; 12 17 18; 12 18 13;

13 18 19; 13 19 14;14 19 20;14 20 15;

15 20 21;16 22 23;16 23 17;17 23 24;

17 24 18;18 24 25;18 25 19;25 19 26;

19 26 20;20 26 27;20 27 21;21 27 28] %单元定义数组（单元结点号）

%相应为单元结点号（编码）、按逆时针顺序输入

COORD=[0 0;-0.75 1.5;0.75 1.5;-1.5 3;0 3;

1.5 3;-

2.25 4.5;-0.75 4.5;0.75 4.5;

2.25 4.5;-3 6;-1.5 6;0 6;1.5 6;3 6;

-3.75 7.5;-2.25 7.5; -0.75 7.5;0.75 7.5;

2.25 7.5;

3.75 7.5;-

4.5 9;-3 9;

-1.5 9;0 9;1.5 9;3 9;4.5 9] %结点坐标数组

%坐标：x，y 坐标（共NPOIN 组）

FORCE=[1 0 -15] %结点力数组（受力结点编号, x 方向,y 方向）FIXED=[22 1 1;28 1 1] %约束信息（约束点，x 约束，y 约束）

%有约束为1，无约束为0

%********************************************************** %生成单元刚度矩阵并组成总体刚度矩阵

ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0

%********************************************************** for i=1:NELEM