有限元课程总结
一三节点三角形单元
1位移函数
移函数写成矩阵形式为：
确定六个待定系数
a4
v玉＞
矩阵形式如下：
J“= TV, 0 Nj 0 N m bJ _ 0 TV, 0
Nj 0
2单元刚度矩阵的计算
1）单元应变和节点位移的关系
由几何方程可以得到单元的应变表达式,
5 6 > = ----------------- b .
「2A '
7
厂 f
Mg
Y — A
”——，
Y Cd
As _
u i
匕?
宀=[N]{5丫
V7
u i
8x
dv
du dv
----- 1 ----
dy dx
J_
2A
C
C
i
bj 0 0
Cj C J b J
u j
V J
2）单元应力与单元节点位移的关系
[KJ = [B r ]T [D][B s ]
b r b s + —
c r c s t s 2 * s
“也+与仏
(T = i,jjn;s = i,jjn)
3) 单元刚度矩阵
卩心][K“] [K]J [K }i ] [K 〃]
[心][K mj ]
3载荷移置
1）集中力的移置
图3
由虚功相等可得，
（㈤丁附=（Q YJW ｛P ｝
由于虚位移是任意的，则皿｝"=["卩｛鬥
2）体力的移置
[S M D I B .] =
E
2A(1-Z /2) Mi Ci %
2 z
如图3所示,
令单元所受的均匀分布体力为｛〃｝=
Et
4(1 —“2)A
地C$ + [DfB i
% [K 加 [K
如
6
由虚功相等可得，
（{J*r）r{7?r =^}>f[N]r{p}tdxdy
{R}e =\\[N]r{p}tdxdy
3）分布面力的移置
设在单元的边上分布有面力{可二[片了r,同样可以得到结点载荷,
{R}e=\[N]T{P}tds
4.引入约束条件，修改刚度方程并求解
1）乘大数法处理边界条件
图3?4所示的结构的约束和载荷情况，如图3?7所示。结点1、4上有水平
方向的位移约束，结点4、6上有垂直方向的约束，结点3上作用有集中力（'，
匕）。
整体刚度矩阵[K]求出后，结构上的结点力可以表示为:
{F} = [K]{5}
根据力的平衡，结点上的结点力与结点载荷或约束反力平衡。用{?}表示结
点载荷和支杆反力，则可以得到结点的平衡方程：
[K]0}={P}
（3.4）
这样构成的结点平衡方程组，在右边向量{P}中存在未知量，因此在求解平衡
方程之前，要根据结点的位移约束情况修改方程（3-4）o先考虑结点n有水平方向位移约朿，与n结点水平方向对应的平衡方程为：
+ ^2w-1.2V l + …+ ?几_[.2幵-1冷+笛”_1.2必+??co
根据支承情况，方程(3-5)应该换成下面的方程:^=
° (3-6)
对比公式(3-5)和(3-6),在式(3-4)中应该做如下修正：
在［K ］矩阵中，第2n ?l 行的对角线元素?s 改为1,该行中全部非对角 线元素改为0；在｛P ｝中，第2n ?l 个元素改为0。为了保持［K ］矩阵的对称性，将 第2ml 列的全部非对角元素也改为0。
同理，如果结点n 在垂肓方向有位移约束，则(3-4)中的第2n 个方程修改 为，
=0
在［K ］矩阵中，第2n 行的对角线元素改为1,该行中全部非对角线元素改为
对图3-4所示结构的整体刚度在修改后可以得到以下的形式,
如果结点n 处存在一个已知非零的水平方向位移知，这时的约束条件为,
0；在｛P ｝中，第2n 个元素改为0。 为了保持［K ］矩阵的对称性，将第2n 列的全 部非对角元素也改为
0。 _
0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
u {
0 卩】
V
2
>
0 0 0 0 0
An-1
1 0 0 0 0
P” 0
10 0 0 0
* * * * * * *
* * *
*00 *00 *00
1 0 0 0 0 10 0 0 * * * 0 Er
0 T 0 0
(3-7)
在［K ］矩阵中，第2n ?l 行的对角线元素並12灯乘上一个大数A,向量{P} 中的对应换成人笛“一…—心，其余的系数保持不变。
方程改为，
^2n-\,\U \ + ^2n-\,2V \ + ???+ ^^2tt-\,2n-\U n + ^2n-\.2n V
n +???匕 ^^2n-\,2n-\ ( 3?8 )
A 的取值要足够大，例如取1010c 只有这样，方程(3-8)才能与方程(3-7) 等价。
二四面体单元
如图1所示的四面体单元，单元结点的编码为i,j,m,n o 每个结点的位移具有三 个分量u, v,w o 这样单元结点的位移列阵可表示成：
1T
单元的位移模式采用线性多项式
u = cc x + cc^x + oc 3y +
q =冬 + oc h x + cr 7 y + ^z 8 n
VP = 6Z 9 + 0()乂 + C] 1 y +
式中，为待定系数，由单元结点的位移和坐标决定。将四个结点的坐标(xi, yi, zi)> (xj,yj,zj)、(xm, ym, zm)> (xn, yn, zn)和结点位移
(ui, vi, wi)> (uj, vj, wj)>
m
(um, vm, wm)> (un, vn, wn)代入(2)式可得12个联立方程，解方程组便可求出。将这十二个系数回代到(2)式，则得到由结点位移和形函数表示的
单元内任一点的位移表达式:
二
NJ Njl N m I N n l]{3}e
= [N]{疔
式中，［I ］为三阶单位阵，［N ］为形函数矩阵。上式即为单元结点位移和单元任
意 点位移之间的关系。 1单元应变和应力
知道单元内任意一点位移后，可利用几何方程确定单元内该点的应变。
其中
单元的应力列阵：
{b} =［?］{￡}=辺旬伪=［s ］伪=区—Sj S m - s 十疗
式中：［S ］为四面体单元的应力矩阵，其分块形式为：
、
u Nj iv 0
Nj
0 0
0 0 i
0 N tt 0 0 N t 0 J
m
n
0 0 N m 0 0 N n
m
n
{易= [B]{S}e = B t -B j
B fn - B tl {5}
[讣
0 利 0 做 3c 0
bi
6V c z . 0 4 0 0 5 0
o 4
b t
0 d ： c i 0 b {
「
A ￡ ?
A&
6
A 】
[s z ]_ [Q]0] —
6A 3 V
4 6
A 》?
6
O
(i,j, m, n)
O
A^di
A Q C,
O
其中
4 — 1 —" - 1 —2 JLl A =
左(1 — Q ) —2(1 — “)
~ 36(1 -+- ")(1 一 2") 2单刚矩阵 对于四面体单元，利用虚功原理，采用类似平面问题的处理方法可以得到其单刚 矩阵
其中:[K]e 为单元刚度矩阵
山厂=JTfwr [Z)][ B\cbcclyclz { d>}
]6
{ d>}
写成分块形式为
g =川可 [D^B^cbcdydz =[J B]T
[Z?][J B]V
—k in
rin
式屮子矩阵[Krs]*下式计算
可以看出，单元刚度矩阵是由单元结点的坐标和单元材料的弹性常数所决定的, 是一
个常数矩阵。如果将空间弹性体划分为ne 个单元和n 个结点，再经过类似 于平而问题的组集过程，就可以得到弹性空间问题的平衡方程
{穴}=[尺]{￡}
—k Ji
kjj
k
mi
mm
k
mn 仏+&(g+d/J
+ A 2b r c s A\db + A 2h t .d s
+ &c/$ c r c s
+ AAd r d s +b r b s ) \d r c s + \c r d s
A
M + \d r d s + \d r c s
￡d$ +企如+g)
三平面四节点四边形单元
矩形单元也是一种常用的单元，它采用了比常应变三角形单元次数更高位的 模式，因而可以更好地反映弹性体中的位移状态和应力状态。
矩形单元1234如图3?1所示，其边长分别为2a 和2b,两边分别平行于x, y 轴。若取该矩形的四个角点为节点，因每个节点位移有两个分量，所以矩形单 元共有8个自由度。
A V4(V4)
图3-1 矩形单元1234
在局部坐标系中，节点i 的坐标是（i , i ）,其值分别为±1。取位移
模式
U =
+
6Z 3?7 + 心4百77
v = a 5 + c z
7 + a 疋 r/
由几何方程可以求得单元的应变
对于平血应力问题
阴(5)
"2 (6)
V2 (^2)
>
’
xy Q
—— <
-
乩一育勿一勿"
<
=
>
等 勿一昭旦
勿
lr 1_ 0
丄/莎
1 - b
[S,
E
Aab{l — /LC )
碍(1+ 77。)
M 彳(1+ %)
二^々77,(1 + ?)
4久(1 +氐) W (l + ￡o) 严碍(1 + %)
"4（
〃
2b
若将单元刚度矩阵写成分块形式
则其中的子矩阵可按下式进行计算
[k u] = JI[Q ]' [Q][6]山5
如果单元厚度t 是常量，则
同
样，对于平面应变问题，只要将上式屮的E 1- 即可。
四边形单元的节点位移与单元节点力之间的关系仍为
[可附十｝
四8节点六面体单元分析
一、形函数与坐标变换 1）形函数
Z, = — （1 T-厂匚厂）? （1 -4- 0三0）? （1 T-疋”疋）
k\ 1
k\2
k\4
E/ 1 ? 2
2）坐标变换
3）位移插值函数与几何矩阵
三、单元刚度矩阵与等效节点载荷向量
殆[v]?W〕
d
i
N\
N\
0...他0 0 0
他N2 0
d
i
d
i
d
x
6N\ 0
0 ON 2 0
0 dx
dx
dx
0 8N X 0
0 dN 2 0
....... 0 dN.
0 1
6
Sy
0 0
6N\ 0
dN 2 0
沁
dz
dz
dz
6N\
6N\ 0
dN 2 dN 2 0
6N 塔 0
労 dx dx
dx
dx
6N\ 6N\ 0
dN 2 8N 2
....... 0 沁
dz
Sy dz
Qy
dz
6N\ 0 6N\ dN 2 0 dN 2 ■ ■■■■■ ■
8z
dx
dz
dx
dz
8x
写成矩阵形式有:
单元刚度矩阵可以表示为
: [K e
]= JJJ[Bf [D] [B]?dv = JJJ[BY [D]
[B]?dxdydz
v e
v e
进一步写成数值积分形式为：
k]二￡￡￡恢苗训[功?恢护屛）]?卩（的,讣
/1側 上I J=1曰
单元体力载荷向量可以表示为
[dN t ]
< 、
dr
dx
> —D]
?v
ds _ ly 」v
f dN t
dN t dt J y
.dz 丿
底H川町?{/；}" = jjj [N] ■{F h}-\j\drdsdt
五 其他常用单元位移函数和自由度
平面三角形单
元
平面应力或应
变
u = a x
+ a 2^
+ a 3ri + a^rj
2°5+必+如77 +。向
单元名称及 适用情况
杆单元 桁架
单元图形
位移模式
平面梁单元 平面刚架
（荟^点）起始瑤
u = a x + a 2x
v = a 3 + a 4x + a 5x 2 + a 6x 3
u = a x +a 2x + a 3y
v = a 4 +a 5x + a h y
平面四边形单
元
平面应力或应
变
w = + a 2x + a 3y 2
+ a 4x + a 5xy
+ a 6y 2 + 如兀'
2 2 + a s x y + a 9xy + a iO y
3 +a n x 3y + a ]2xy 3
w = a }L } +a 2L 2 +?3L 3 +a 4￡2 + a 5 厶
3
厶 +a 6 L|L 2
+ ?(厶2厶;一厶厶;) + a K (L y L\ -厶厶;)
+ a () (L| L~2 — L 2L\)
u = a 】+ a 2e + a^ri + a 4^ v = a 5+ a b ￡ + tz 77 +
w = a g +a]o￡ + 4]〃 + Q]2：
u = a } + a 2c + 57 + +
5科+唧：+ 6? + a 亦
v = a 9+ + q [〃 + mg
+ 5印+ 5皿+ gQ +仏釦： w
= a l7+a ls ￡ + a lg ^ + a 2O C
+ a?何 + a 22i^ +幺 + 5品
有限元中，梁单元的节点有6个自由度，壳单元节点有5个自由度，而体单元 节点有3个自由度。
四面体单元
三维应力
六面体单元 三维应力
矩形板单元 薄板弯曲问题
三角形板单元 薄板弯曲问题
. 一、题目 如图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
有限元课程总结 一三节点三角形单元 1位移函数 移函数写成矩阵形式为： 确定六个待定系数 a4 v玉＞ 矩阵形式如下： J“= TV, 0 Nj 0 N m bJ _ 0 TV, 0 Nj 0 2单元刚度矩阵的计算 1）单元应变和节点位移的关系 由几何方程可以得到单元的应变表达式, 5 6 > = ----------------- b . 「2A ' 7 厂 f Mg Y — A ”——， Y Cd As _ u i 匕? 宀=[N]{5丫 V7 u i 8x dv du dv ----- 1 ---- dy dx J_ 2A C C i bj 0 0 Cj C J b J u j V J
2）单元应力与单元节点位移的关系 [KJ = [B r ]T [D][B s ] b r b s + — c r c s t s 2 * s “也+与仏 (T = i,jjn;s = i,jjn) 3) 单元刚度矩阵 卩心][K“] [K]J [K }i ] [K 〃] [心][K mj ] 3载荷移置 1）集中力的移置 图3 由虚功相等可得， （㈤丁附=（Q YJW ｛P ｝ 由于虚位移是任意的，则皿｝"=["卩｛鬥 2）体力的移置 [S M D I B .] = E 2A(1-Z /2) Mi Ci % 2 z 如图3所示, 令单元所受的均匀分布体力为｛〃｝= Et 4(1 —“2)A 地C$ + [DfB i % [K 加 [K 如 6
由虚功相等可得， （{J*r）r{7?r =^}>f[N]r{p}tdxdy {R}e =\\[N]r{p}tdxdy 3）分布面力的移置 设在单元的边上分布有面力{可二[片了r,同样可以得到结点载荷, {R}e=\[N]T{P}tds 4.引入约束条件，修改刚度方程并求解 1）乘大数法处理边界条件 图3?4所示的结构的约束和载荷情况，如图3?7所示。结点1、4上有水平 方向的位移约束，结点4、6上有垂直方向的约束，结点3上作用有集中力（'， 匕）。 整体刚度矩阵[K]求出后，结构上的结点力可以表示为: {F} = [K]{5} 根据力的平衡，结点上的结点力与结点载荷或约束反力平衡。用{?}表示结 点载荷和支杆反力，则可以得到结点的平衡方程： [K]0}={P} （3.4） 这样构成的结点平衡方程组，在右边向量{P}中存在未知量，因此在求解平衡
第一章 1.有限单元法求得的解为：[ ]3 A.精确解 B.解析解 C.近似解 D.整数解 2.弹性力学问题的基本解法有：[ ] ABD A.按位移求解 B.按应力求解 C.按单元刚度求解 D. 混合求解 E.按整体刚度求解 23.弹性力学问题的基本解法有：按位移求解，按应力求解和[ ]3 A. 按单元刚度求解 B. 按整体刚度求解 C. 混合求解 D.按平衡方程求解 24.弹性力学问题的基本解法有：按位移求解，混合求解和[]4 A. 按平衡方程求解 B. 按单元刚度求解 C. 按整体刚度求解 D. 按应力求解 25.弹性力学问题的基本解法有：按应力求解，混合求解和[ ]2 A. 按整体刚度求解 B. 按位移求解 C. 按单元刚度求解 D. 按平衡方程求解 3.用弹性力学经典解法解决实际问题的主要困难在于：[ ]4 A.对弹性体离散化的复杂性 B.刚度矩阵求解的困难性 C.受力分析的复杂性 D.求解偏微分方程的复杂性 4.用三角形单元的节点位移，可以表示单元中的：[ ]BDE A.弯矩 B.应变 C.扭矩 D.应力 E.结点力 26.用三角形单元的节点位移，可以表示单元中的应变，应力和[ ]3 A. 扭矩 B. 弯矩 C. 结点力 D.外力 27.用三角形单元的节点位移，可以表示单元中的应变，结点力和[ ]4 A.弯矩 B. 外力 C. 扭矩 D. 应力 28. 用三角形单元的节点位移，可以表示单元中的应力，结点力和[ ]4， A. 外力 B. 扭矩 C. 弯矩 D. 应变 5.将各个单元集合成离散化的结构模型进行整体分析，问题最后归结为求解[ ]。2 A.结点位移 B. 以结点位移为未知量的线性方程组 C.整体刚度矩阵 D.单元刚度矩阵 6.对于三角形三结点单元，其结点按照[]顺序进行排列。3 A.从左至右 B. 顺时针 C. 逆时针 D.以上均可 7.对于三角形三结点单元，每个结点位移在单元平面内有[ ]个分量 2 A.1 B.2 C.3 D.4 8.对于三角形三结点单元，共有[ ]个位移分量。4 A.3 B.4 C.5 D.6 9.形函数 N在结点i上的值等于[ ]。2 i A.0 B.1 C. -1 D.2 10.在单元中任意一点，三个形函数之和等于[ ]2 A.0 B.1 C.2 D.3 11.有了单元的位移模式，就可以应用[ ]求得单元的应变3 A.平衡微分方程 B.物理方程 C. 几何方程 D.积分方程 12.单元应力矩阵[S]与弹性矩阵[D]和单元应变矩阵[B]的关系是：[ ]C A. [S]= [D]+ [B] B. [S]= [D]—[B] C. [S]= [D] [B] D. [S]= [D]/ [B] 13.三角形三结点单元中，单元应力矩阵[S]是一个[ ]4 A.对称矩阵 B.零矩阵 C.非常数矩阵 D.常数矩阵 14.三角形三结点单元的应力分量为[ ]1 A.常量 B.变量 C.零 D.不确定
第三章平面问题的有限元法本章通过三角形常应变单元，介绍有限元法应用于弹性体应力分析的基本原理和方法：包括弹性体的离散化，单元特性的分析，刚度矩阵的建立，等效节点力的计算，解答的收敛性以及实施步骤和注意事项，同时对形函数的性质也作简要的叙述。 第一节三角形常应变单元 一、结构离散化 用有限元法分析弹性力学平面问题，第一步就是把原来的连续的弹性体离散化。 (a) (b) 图3.1 弹性体和有限元模型 将整个结构（平板）划分成有限个三角形。这样的三角形称为单元（三角形单元）。 三角形单元的顶点取为节点。3节点三角形单元用边界节点间的直线段来近似板的曲线边界。 这些三角形在其节点处相互连接，组成一个单元集合体，以代替原来的弹性体。 注：1. 全部节点和全部单元一般由1开始按自然顺序编号。2. 节点编码：
总码-----------用于整体分析，如1,2，…，按自然顺序编号局部码--------用于单元分析，i，j，m 要求按逆时针方向的顺序进行编码 每个单元的节点局部码i，j，m和节点总码有一一对应的关系 3. 单元间不能有重叠 4. 一个单元的任一顶点不许为另一单元任一边的内点 5. 所有作用在单元上的载荷，包括集中载荷、表面载荷和体积力，都按虚功等效的原则移置到节点上，成为等效节点载荷。
二、 位移模式 1. 单元节点位移列阵 i u 图 3.2 平面三角形单元 设单元e 的节点号码为i ，j ，m 。由弹性力学平面可知，单元内任意一点有两个位移分量u ，v ，记为 {}T f u v ????= 故每个节点也有两个位移分量，因此称节点自由度为2。3个节点得位移分量分别是 ,,,,,m m i i j j u v u v u v ，用列阵表示为 {} i e i i e j j j m m m u v u v u v δδδδ?? ???????????????????????????? ?????? == （3-1）
第八章 平面问题的有限元分析及三角形单元的应用 第一节 概述 分析弹性力学平面问题时，最简单的单元式由三个结点组成的三角形单元。当用以分析平面应力问题时，可将其视为三角板；当用以分析平面应变问题时，则可式为三棱柱。各单元在结点处为铰结。图8-1所示位移悬臂梁离散为三角形单元的组合体 以矩阵形式列出弹性力学平面问题的基本量和基本方程。 谈形体所受体力分量可表示为 [ ] T y x y x p p p p p =??? ? ????= (8-1) 所受面力分量可表示为 [ ] T y x y x p p p p p =??? ? ????= （8-2） 体内任一点应力分量可表示为 []T xy y x τδδδ= （8-3） 任一点的应变分量可表示为 []T xy y x γεεε= （8-4） 任一点的位移分量可表示为 []T v u =δ （8-5） 弹性力学平面问题的几何方程的矩阵表达式为 ?? ???? ???????? ??? ???+??????=????????????=x u y v y v x u xy y x εεεε （8-6） 平面应力问题的物理方程的矩阵表达式为 ? ???? ? ?????????? ????????? ?--= ????? ? ??????xy y x xy y x E γεεμμμ μτσσ210 0010112 （8-7） 或简写成为 εσD = （8-8） 式中
???? ?? ? ?????? ?--=210 0010112μμμ μ E D （8-9） 称为弹性矩阵。 平面应变问题的物理方程也可写成式（8-8），但须将式（8-9）中的E 换成 2 1μ -E ，μ换成 2 1μμ -，因此得出 ???? ?? ????????? ?? ?-----+-= )1(2210 00110 11)21)(1()1(2 2 μμμμμμ μμμE D （8-10） 平衡微分方程及边界条件也可以用矩阵表示，但弹性力学有限元位移法中，通常用虚功 方程代替平衡微分方程和应力边界条件。虚功方程的矩阵表达式为 ?????***=+tdxdy tds p f ptdxdy f T T σε （8-11） 式中：[ ] T v u f ** * =，表示虚位移； []T xy x x * ***=γεεε，表示与虚位移相对应的虚应变。 为了便于计算，作用于弹性体上的体力和面力替换为作用在结点上的集中力，即等效结 点荷载。设作用于各个结点上的外力分量用如下列阵来表示 []T n n V U V U V U F ?=2211 与这些结点外力分量相对应得结点虚位移分量列阵为 []T n n v u v u v u * ******?=2211δ 则外力在虚位移上做的虚功为 F v V u U v V u U v V u U T n n n n ** *****=++?++++δ22221111 如平面弹性体的厚度为t ，该虚功除以t ，即可得出单位厚度薄板上的外力虚功。于是，式（8-11）所示虚功方程可写成 ??**=tdxdy F T T σεδ （8-11） 虚功方程不仅仅应用于弹性力学，也可用于塑性力学。其应用条件是：只要变形体的全部外力和应力满足平衡方程；位移是微小的，并满足边界条件，位移与应变满足几何方程。
三角形常应变单元程序的编制与使用 有限元法是求解微分方程边值问题的一种通用数值方法，该方法是一种基于变分法（或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。 有限元分析的基本步骤可归纳为三大步：结构离散、单元分析和整体分析。对于平面问题，结构离散常用的网格形状有三角形、矩形、任意四边形，以三个顶点为节点的三角形单元是最简单的平面单元，它较矩形或四边形对曲边边界有更好的适应性，而矩形或四边形单元较三节点三角 形有更高的计算精度。 Matlab语言是进行矩阵运算的强大工具，因 此，用Matlab语言编写有限元中平面问题的程序 有优越性。本章将详细介绍如何利用Matlab语言 编制三角形常应变单元的计算程序，程序流程图见 图1。 有限元法中三节点三角形分析结构的步骤如 下： 1）整理原始数据，如材料性质、荷载条件、约 束条件等，离散结构并进行单元编码、结点 编码、结点位移编码、选取坐标系。 2）单元分析，建立单元刚度矩阵。 3）整体分析，建立总刚矩阵。 4）建立整体结构的等效节点荷载和总荷载矩 阵 5）边界条件处理。 6）解方程，求出节点位移。 7）求出各单元的单元应力。 8）计算结果整理。计算结果整理包括位移和应 力两个方面；位移计算结果一般不需要特别 的处理，利用计算出的节点位移分量，就可 画出结构任意方向的位移云图；而应力解的 误差表现在单元内部不满足平衡方程，单元与单元边界处应力一般不连续，在边界上应力解一般与力的边界条件不相符合。图1 程序流程图
1.1 程序说明 %******************************************************************* % 三角形常应变单元求解结构主程序 %******************************************************************* ●功能：运用有限元法中三角形常应变单元解平面问题的计算主程序。 ●基本思想：单元结点按右手法则顺序编号。 ●荷载类型：可计算结点荷载。 ●说明：主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算 法的全过程，即实现程序流程图的程序表达。 %----------------------------------------------------------------------------------------------------- 1 程序准备 format short e %设定输出类型 clear all %清除所有已定义变量 clc %清屏 ●说明： format short e －设定计算过程中显示在屏幕上的数字类型为短格式、科学计数法； clear all －清除所有已定义变量，目的是在本程序的运行过程中，不会发生变量名相同等可能使计算出错的情况； clc －清屏，使屏幕在本程序运行开始时 %----------------------------------------------------------------------------------------------------- 2 全局变量定义 global NNODE NPION NELEM NVFIX NFORCE COORD LNODS YOUNG POISS THICK global FORCE FIXED global BMATX DMATX SMATX AREA global ASTIF ASLOD ASDISP global FP1 ●说明： NNODE—单元结点数，NPION—总结点数，NELEM—单元数，NVFIX—受约束边界点数，NFORCE—结点力数，COORD—结构结点坐标数组，LNODS —单元定义数组，YOUNG—弹性模量，POISS—泊松比，THICK—厚度
【ANSYS 算例】4.7.1(3) 基于3节点三角形单元的矩形薄板分析 如图4-20所示为一矩形薄平板，在右端部受集中力100 000N F =作用，材料常数为：弹性模量7110Pa E =?、泊松比1/3μ=，板的厚度为0.1m t =，在ANSYS 平台上，按平面应力问题完成相应的力学分析。 (a) 问题描述 (a) 有限元分析模型 图4–20 右端部受集中力作用的平面问题(高深梁) 解答 在ANSYS 平台上，完成的分析如下。 1. 基于图形界面的交互式操作(step by step) (1) 进入ANSYS(设定工作目录和工作文件) 程序 → ANSYS Interactive →Working directory (设置工作目录) →Initial jobname (设置工作文件名): 2D3Node →Run → OK (2) 设置计算类型 ANSYS Main Menu : Preferences… → Structural → OK (3) 选择单元类型 ANSYS Main Menu : Preprocessor →Element Type →Add/Edit/Delete… →Add… →Solid ：Quad 4node 42 →OK (返回到Element Types 窗口) → Options… →K3: Plane Strs w/thk(带厚度的平面应力问题) →OK →Close (4) 定义材料参数 ANSYS Main Menu : Preprocessor →Material Props →Material Models →Structural →Linear →Elastic → Isotropic: EX:1.0e7 (弹性模量)，PRXY: 0.33333333 (泊松比) → OK → 鼠标点击该窗口右上角的“ ”来关闭该窗口 (5) 定义实常数以确定平面问题的厚度 ANSYS Main Menu: Preprocessor →Real Constant s… →Add/Edit/Delete →Add →Type 1→ OK →Real Constant Set No: 1 (第1号实常数), THK: 0.1 (平面问题的厚度) →OK →Close (6) 生成单元模型 生成4个节点 ANSYS Main Menu: Preprocessor →Modeling → Create → Nodes → On Working Plane →输入节点1的x,y,z 坐标(2,1,0)，回车→输入节点2的x,y,z 坐标(2,0,0)，回车→输入节点3的x,y,z 坐标(0,1,0)，回车→输入节点4的x,y,z 坐标(0,0,0)，回车→OK 定义单元属性 ANSYS Main Menu: Preprocessor →Modeling → Create → Elements → Elem Attributes →Element type number:1 →Material number:1→Real constant set number:1 →OK 生成单元 ANSYS Main Menu: Preprocessor →Modeling → Create → Elements → User Numbered → Thru Nodes →Number to assign to element:1→Pick nodes:2,3,4→OK →Number to assign to element:2→Pick nodes:3,2,1→OK (7) 模型施加约束和外载 左边两个节点施加X,Y 方向的位移约束 ANSYS Main Menu: Solution → Define Loads → Apply →Structural → Displacement → On
单元节点和积分节点的联系和区别 有限元方法的实质： 通过变分原理极小值转化为矩阵的极小值（(变分原理)→ (最小势能原理) (虚功原理) 变分原理：把一个物理学问题（或其他学科的问题）用变分法化为求泛函极值（或驻值）的问题，就称为该物理问题 （或其他学科的问题）的变分原理。如果建立了一个新的变分原理，它解除了原有的某问题变分原理的某些约束条件，就称为该问题的广义变分原理；如果解除了所有的约束条件，就称为无条件广义变分原理，或称为完全的广义变分原理。1964年，钱伟长教授明确提出了引进拉格朗日乘子（Lagrange multiplier ）把有约束条件的变分原理化为较少（或没有）约束条件的变分原理的方法。 最小势能原理：最小势能原理就是说当一个体系的势能最小时，系统会处于稳定平衡状态。举个例子来说，一个小球在曲面上运动，当到达曲面的最低点位置时，系统就会趋向于稳定平衡。势能最小原理与虚功原理本质上是一致的。宇宙万物，如果其势能未达到“最小”（局部概念），它总要设法变化到其“相对”最小的势能位置。举个例子：一个物体置于高山上，它相对于地面来说有正的势能（非最小），因而它总有向地面运动的“能力”（向地面“跃迁”）（其力学本质是其处于一种不稳平衡状态）。因此，它试图（也只有）向下运动，才能保证其达到一个相对平稳的状态。λ 最小势能原理是势能驻值原理在线弹性范围里的特殊情况。对于一般性问题：真实位移状态使结构的势能取驻值（一阶变分为零），在线弹性问题中取最小值。形象的说，当你在一百米高的钢丝绳上走的时候你总是希望尽早回到地上，但其实只要你不动你也是平衡的，因为驻值也可以是极大值（此时称为随遇平衡）。而当你在一百米高的大楼里的办公室里时，你并不害怕，因为周围的物体的势能均不比你小，此时驻值取的是极小值而不是最小值。在有限元的理论中，最小势能原理是在所有满足给定边界条件的位移时，满足平衡微分方程的位移使得势能取得最小值。公式如下： 在江见鲸的有限元就很清晰的介绍了有限元的原理，再参考汪新老师的ppt ，理解就十分透彻。以下是我在网上搜索到的关于有限元的书，很值得一看。另外cook 的书也很好，是陈贡发老师介绍的，很值得一看。 1《The Finite Element Method 》O.C.Zienkiewicz ，R.Taylor 著，第五版，三卷本，有中文译本 《有限元法》(英)监凯维奇著（第四版中译本1985年出版，上下册，尹泽勇等译，权威著作，有限元研究者必读；第五版译名改成了《有限单元法》，曾攀译，2008年出版~~） 2.《Nonlinear Finite Element for Continua and Structures 》T .Belytschko 等著，有中文译本 《连续体和结构的非线性有限元》庄茁译，清华大学出版社，固体力学非线性有限元的集大成之作~~ ：外力对系统作的功势能：W U 21???? ??--+??? ??=+=∏???dS dV dV W U S T V T V T σf u f u εσ
北方工业大学 高等有限元课程总结 姓名:韩双鹏学号: 2011303310105 专业班级:结构研-11 系(部、院):建筑工程学院 2012 年5 月25 日
高等有限元学习总结——六节点三角形等参数单元 1 概述 从弹性力学基本方程到有限元原理再到最新进展，经过本课程的学习，比较系统的掌握了有限元相关内容，更学习到了一种方法、一些生活中的哲理。首先从大方向掌握所学内容，避免迷失在局部造成一叶遮目不见泰山之悲剧，比如弹性力学原理从大方向说就是三类方程，以及其在各类问题中的应用；其次了解了科研的相关过程及创新之处，从已知的东西到无知的领域，正如老师所说，能成功地把某一领域的东西搬到相关领域，这就是一大创造，比如有限元中将梁弯曲的理论研究厚板弯曲问题，由有限元标准单元到等参元的研究等；再有，我们生活中的常识、学习中的某些东西值得我们细细品味，也许这就是平时所说的小事反应大道理，老师的理论：“很多想法都是错误的”“很好想到的方法也许很难走通”“有缺陷的东西才更体现出美”“平衡的理论，吃点亏也许是福”等等，受益匪浅。不再一一赘述，本文将取其中的一个知识点，总结六节点三角形等参单元的相关内容。 我们知道，无论三节点或者六节点三角形单元还是四节点或者八节点矩形单元，它们形状简单、规则但计算精度低，且对于复杂边界的适应性差，难以很好的拟合曲边边界，解决这一问题的通用方法是细分边界，以直代曲，利用更多的简单单元去拟合边界复杂的区域。但这样处理仍存在折线代替曲线所带来的误差，且这种误差不能通过提高单元位移函数的精度来补偿。那么能否构造出单元形状任意、边界适应性好、计算精度高的曲边单元，以便在给定的精度下用较少数目的单元去解决实际问题？这就是有限元中一类重要的单元——等参数单元。本文将总结等参数单元的基本概念，并以六节点三角形单元为例讲述等参元实现过程中的三种变换，以及该等参元的收敛性等问题。 2 等参数单元及实现过程 2.1 等参数单元概念 由于实际问题的复杂性，通常需要使用一些形状不规整和形状复杂的单元来离散边界形状复杂的原问题。如下图所示（a）中为常见的几何形状不规整的实际单元，称为实际单元，也称为参数单元。（b）中为对应的形状规整的单元，称为标准单元。对于形状复杂的实际单元的单元分析，若仍采用前面介绍的方法进行，则在单元位移函数的建立和单元刚度矩阵计算方面会遇到许多困难。由此可考虑利用前面介绍过的形状规整的标准单元的单元分析来研究实际单元，几何形状的不同可认为是坐标变换的结果。
有限元八种三维单元介绍 有限元三维体单元常见单元有四面体4、10节点单元、六面体8、20、27节点单元、三棱柱6、15节点单元。我们在2000年新问世的四面体20节点单元。下面分别介绍如下： 1 四面体4节点单元（常应变单元、一次单元），见图一。 单元内部的位移插值函数为一次多项式，即只含常数项和Z Y X ,,四项。应变是位移的偏导数，故在单元内部，应力和应变为常数，位移和应力收敛速度都很慢，是非常落后的单元。 图一 四面体4节点单元（常应变单元） 2 四面体10节点单元（二次单元），见图二。 用体积坐标定义的单元：单元内位移插值函数为二次完全多项式，即含常数项和Z Y X ,,,YZ XZ XY Z Y X ,,,,,222十项，在单元内部，应力和应变为一次完全多项式，位移收敛速度很快，但应力收敛速度仍较慢。由于整体加密使用的节点数太多，而局部加密生成的单元奇异，刚度阵病态，故应力集中问题中很难得到精度较高的解，在不考虑应力集中、疲劳寿命的问题中，由于该单元使用节点较少、几何适应性强，被人们经常使用。 用直角坐标定义的单元：由六面体20节点单元通过节点重合退化得到。这种单元误差较大，无法求节点应力,只能求出 GAUSS 积分点的应力值，不推荐使用。 3 四面体20节点单元（三次单元），见图三。 用体积坐标定义的单元，单元内位移插值函数为完全三次多项式，即含常数项和Z Y X ,,, YZ XZ XY Z Y X ,,,,,222，XYZ Y Z X Z Z Y X Y Z X Y X Z Y X ,,,,,,,,,222222333二十项， 在单元内部，应力和应变为完全二次多项式，位移和应力收敛速度都很快，精度最高、几何适应性强，在应力集中、疲劳寿命分析问题中使用是非常有用和令人放心的单元。 4 三棱柱6节点单元（一次单元），见图四。 与四面体4 节点单元类似。
【MATLAB 算例】4.7.1(2) 基于3节点三角形单元的矩形薄板分析(T riangle2D3Node) 如图4-20所示为一矩形薄平板，在右端部受集中力100 000F N =作用，材料常数为：弹性模量7110E Pa =?，泊松比13 μ=，板的厚度0.1t m =。基于MA TLAB 平台求解该结构的节点位移、支反力以及单元应力。 图4-20 解答：对该问题进行有限元分析的过程如下。 （1）结构的离散化与编号 将结构离散为二个3节点三角形单元，单元编号及节点编号如图4-20(b)所示。 （2）计算各单元的刚度矩阵(以国际标准单位) 首先在MA TLAB 环境下，输入弹性模量E 、泊松比NU 、薄板厚度t 和平面应力问题性质指示参数ID ，然后针对单元1和单元2，分别两次调用函数Triangle2D3Node_Stiffness ，就可以得到单元的刚度矩阵k1(6×6)和k2(6×6)。 >> E=1e7; >> NU=1/3; >> t=0.1; >> ID=1; >> k1=Triangle2D3Node_Stiffness(E,NU,t,2,0,0,1,0,0,ID) k1 = 1.0e+006 * 0.2813 0 0 0.1875 -0.2813 -0.1875 0 0.0938 0.1875 0 -0.1875 -0.0938 0 0.1875 0.3750 0 -0.3750 -0.1875 0.1875 0 0 1.1250 -0.1875 -1.1250 -0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750 -0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188 >>k2=Triangle2D3Node_Stiffness(E,NU,t,0,1,2,0,2,1,ID) k2 = 1.0e+006 * 0.2813 0 0 0.1875 -0.2813 -0.1875 0 0.0938 0.1875 0 -0.1875 -0.0938 0 0.1875 0.3750 0 -0.3750 -0.1875 0.1875 0 0 1.1250 -0.1875 -1.1250 -0.2813 -0.1875 -0.3750 -0.1875 0.6563 0.3750 -0.1875 -0.0938 -0.1875 -1.1250 0.3750 1.2188 （3） 建立整体刚度方程 由于该结构共有4个节点，则总共的自由度数为8，因此，结构总的刚度矩阵为KK (8×8)，先对KK 清零，然后两次调用函数Triangle2D3Node_Assembly 进行刚度矩阵的组装。 >>KK = zeros(8,8); >>KK=Triangle2D3Node_Assembly(KK,k1,2,3,4); >>KK=Triangle2D3Node_Assembly(KK,k2,3,2,1) KK = 1.0e+006 * Columns 1 through 6 0.6563 0.3750 -0.3750 -0.1875 -0.2813 -0.1875 0.3750 1.2188 -0.1875 -1.1250 -0.1875 -0.0938 -0.3750 -0.1875 0.6563 0 0 0.3750 -0.1875 -1.1250 0 1.2188 0.3750 0 -0.2813 -0.1875 0 0.3750 0.6563 0 -0.1875 -0.0938 0.3750 0 0 1.2188
一：用三结点三角形平面单元计算平面结构的应力和位移。 1，设计说明书 计算简图，网格划分，单元及结点的编号如下图所示。由于结构对称，去四分之一结构分析。其中 E=2e10pa,mu=0.167,h=1m. 变量注释： Node ------- 节点定义 gElement ---- 单元定义
gMaterial --- 材料定义，包括弹性模量，泊松比和厚度 gBC1 -------- 约束条件 gNF --------- 集中力 gk------------总刚 gDelta-------结点位移 子程序注释： PlaneStructualModel ———定义有限元模型 SolveModel ———————求解有限元模型 DisplayResults ——————显示计算结果 k = StiffnessMatrix( ie )———计算单元刚度 AssembleStiffnessMatrix( ie, k )—形成总刚 es = ElementStress( ie )————计算单元应力 function exam1 % 输入参数：无 % 输出结果：节点位移和单元应力 PlaneStructualModel ; % 定义有限元模型 SolveModel ; % 求解有限元模型 DisplayResults ; % 显示计算结果 return ; function PlaneStructualModel % 定义平面结构的有限元模型 % 输入参数：无 % 说明： % 该函数定义平面结构的有限元模型数据： % gNode ------- 节点定义 % gElement ---- 单元定义 % gMaterial --- 材料定义，包括弹性模量，泊松比和厚度% gBC1 -------- 约束条件 % gNF --------- 集中力 global gNode gElement gMaterial gBC1 gNF % 节点坐标 % x y gNode = [0.0, 2.0 % 节点1 0.0, 1.0 % 节点2 1.0, 1.0 % 节点3 0.0, 0.0 % 节点4 1.0, 0.0 % 节点5 2.0, 0.0] ; % 节点6 % 单元定义 % 节点1 节点2 节点3 材料号 gElement = [3, 1, 2, 1 % 单元1
六节点三角形单元 单元节点位移向量： T e v u v u v u v u v u v u } {}{665544332211=δ 单元内任意一点的位移是单元节点位移的插值函数 已知u(x,y)和v(x,y)在六个顶点的函数值，可分别设u(x,y)和v(x,y)为具有六个待定系数的 插值多项式 2 652 432126524321y B xy B x B y B x B B v y A xy A x A y A x A A u +++++=+++++= 将已知条件代入，可解得其中的待定系数 u(x,y)和v(x,y)也可用Lagrenge 插值公式表示为： ∑=+++++==6 1665544332211k k k u N u N u N u N u N u N u N u ∑=+++++==61 665544332211k k k v N v N v N v N v N v N v N v 写成矩阵形式： ???????? ? ?????????? ?????????????????????????=??????66554433221165 4 3 2 1 654321v u v u v u v u v u v u 0 000000N N N N N N N N N N N N v u
缩写为： {}e ][δδN v u =? ?? ???= 单元内任意一点的应变 x v y u y v x u xy y x ??+ ??=??=??= γεε 矩阵形式： {}{}e N H v u H x y y x δε]][[][v u 0 0=??????=??? ???????? ?? ? ???????????? ???? = }[B]{}{ e δε=∴ []12 365 4 3 2 1 65 4 3 2 1 654321 0000000 0 0]][[][?=?? ??????????? ? ???????????????? ==B B B B B B N N N N N N N N N N N N x y y x N H B []?? ?? ??? ?????????????????=x N y N y N x N B k k k k k 0 0 k = 1,2,3...6 N k (k=1,2,3...6)为六节点三角形单元的形函数，N k 的表达式： 26524321),(y a xy a x a y a x a a y x N k k k k k k k +++++= k =1,2,3 (6) k i y x N y x N i i k k k k ≠== 0),(1 ),( y a x a a y N y a x a a x N k k k k k k k k 6535422 2 ++=??++=??∴ 单元内任意一点的应力：}]{][[}]{[}{e B D D δεσ==
《有限元作业》 年级2015级 学院机电工程学院 专业名称 班级学号 学生姓名 2016年05月
如下图所示为一受集中力P作用的结构，弹性模量E为常量，泊松比V=1/6,厚度为I=1。按平面应力问题计算，运用有限元方法，分别采用三角形及四边形单元求解，求节点位移及单元应力（要求三角形单元数量不少于4个，四边形单元不少于2个） 图(一) 图（二）三角形单元求解 图（三）四边形单元求解
（1）如图划分三角形单元，工分成四个分别为???④ （2）如图分别进行编号1、2、3、4、5、6，并建立坐标系（3）编程进行求解，得出结果，其中假设力P=2000N 调用Triangle2D3Node_Stiffness函数，求出单元刚度矩阵k1 = +06 * 0 0 0 0 0 0 0 0 k2 = +06 * 0 0 0 0 0 0 0 0 k3 = +06 * 0 0 0 0 0 0 0 0 k4 = +06 *
0 0 0 0 0 0 0 0 调用Triangle2D3Node_Assembly函数，求出总体刚度矩阵 求出的节点位移 U = 调用Triangle2D3Node_Stress函数，求出应力,S1、S2、S3、中求出的分别为Sx,Sy,Sxy S1 = +03 * S2 = +03 *
S3 = +03 * S4 = +03 * 二、 （1）如图划分四边形单元，工分成四个分别为?? （2）如图分别进行编号1、2、3、4、5、6，并建立坐标系（3）编程进行求解，得出结果，其中假设力P=2000N 调用 Quad2D4Node_Stiffness函数，求出单元刚度矩阵 调用Quad2D4Node_Assembly函数，求出求出总体刚度矩阵
八节点等参单元平面有限元 分析程序 土木工程学院 2011.2
目录 1.概述 (1) 2.编程思想 (2) 2.1.八节点矩形单元介绍 (2) 2.2.有限元分析的模块组织 (5) 3.程序变量及函数说明 (6) 3.1.主要变量说明： (6) 3.2.主要函数说明 (7) 4.程序流程图 (8) 5.程序应用与ANSYS分析的比较 (9) 5.1.问题说明 (9) 5.2.ANSYS分析结果 (10) 5.3.自编程序分析结果 (12) 5.4.结果比较分析 (12) 参考文献 (14) 附录程序源代码 (15)
《计算力学》课程大作业 1.概述 通常情况下的有限元分析过程是运用可视化分析软件(如ANSYS、SAP等)进行前处理和后处理，而中间的计算部分一般采用自己编制的程序来运算。具有较强数值计算和处理能力的Fortran语言是传统有限元计算的首选语言。随着有限元技术的逐步成熟,它被应用在越来越复杂的问题处理中，但在实际应用中也暴露出一些问题。有时网格离散化的区域较大,而又限于研究精度的要求，使得划分的网格数目极其庞大，结点数可多达数万个，从而造成计算中要运算的数据量巨大，程序运行的时间较长的弊端，这就延长了问题解决的时间，使得求解效率降低。因为运行周期长，不利于程序的调试，特别是对于要计算多种运行工况时的情况；同时大数据量处理对计算机的内存和CPU 提出了更高的要求，而在实际应用中，单靠计算机硬件水平的提高来解决问题的能力是有限的。因此，必须寻找新的编程语言。 随着有限元前后处理的不断发展和完善，以及大型工程分析软件对有限元接口的要求，有限元分析程序不应只满足解题功能，它还应满足软件工程所要求的结构化程序设计条件，能够对存储进行动态分配，以充分利用计算机资源，它还应很容易地与其它软件如CAD 的实体造型，优化设计等接口。现在可编写工程应用软件的计算机语言较多，其中C语言是一个较为优秀的语言，很容易满足现在有限元分析程序编程的要求。 C语言最初是为操作系统、编译器以及文字处理等编程而发明的。随着不断完善，它已应用到其它领域，包括工程应用软件的编程。近年来，C语言已经成为计算机领域最普及的一个编程语言，几乎世界上所有的计算机都装有C的编译器，从PC机到巨型机到超巨型的并行机，C与所有的硬件和操作系统联系在一起。用C 编写的程序，可移植性极好，几乎不用作多少修改，就可在任何一台装有ANSI、C编译器的计算机上运行。C既是高级语言，也是低级语言，也就是说，可用它作数值计算，也可用它对计算机存储进行操作。
P 9 m 9 m 一、题目 如图1所示，一个厚度均匀的三角形薄板，在顶点作用沿板厚方向均匀分布的竖向载荷。已知：P=150N/m，E=200GPa，=，t=，忽略自重。试计算薄板的位移及应力分布。 要求： 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）最里层：区分为第奇/偶数个计算 单元刚度的集成： [][] [][] [][] ' ' ' ' ' ' 2 1 56 56 6 6 56 56 2 6 6 2 56 56 1 6 6 1 e Z e e e Z e Z e e e e k k k K k k k k k k + ? + + = ? = ? = = ? = = ? = ? ? ? ? ? ? 边界约束的处理：划0置1法 适用：这种方法适用于边界节点位移分量为已知(含为0)的各种约束。 做法： （1）将总刚矩阵〔K〕中相应于已知位移行主对角线元素置1，其他元素改为零；同 时将载荷列阵｛R｝中相应元素用已知位移置换。 ◎这样，由该方程求得的此位移值一定等于已知量。 （2）将〔K〕中已知位移相应的列的非主对角成元素也置0，以保持〔K〕的对称性。 ◎当然，在已知位移分量不为零的情况下，这样做就改变了方程左端的数值，为 保证方程成立，须在方程右端减去已知位移对该方程的贡献——已知位移和相应总刚元素的乘积。◎若约束为零位移约束时，此步则可省去。 特点： （1）经以上处理同样可以消除刚性位移(约束足够的前提下)，去掉未知约束反力。 （2）但这种方法不改变方程阶数，利于存贮。