0 0 0 K 45 K 55
0 0 K 36 K 46 K 56 K 66
0 0 0 0 0 K 67 K 77
0 0 0 0 K 58 0 K 78 K88
A(1) A
A( 3) A( 2)
A(5) A(4)
0 0 0 0 0 0 0 0 0 0 0 0 0 0 K 0 0
j列
顶线下 Kmj,j
（2-2）
第 j列
Ki,j
第i 行
Kj,j 图5-12
1)列的列高之和加1，即
MAXA(j) = h1 +…+ hj-2 - hj-1 + 1
=（h1 +…+ hj-2 + 1）+ hj-1 = MAXA（j - 1）+ hj-1
因为永远有
MAXA(1)= 1， MAXA(2)= 2
故计算主元地址的公式可写为
MAXA(j+1)= MAXA(j)+ hj 式中， j = 2，3，…，N； （2-1）
平均法整理单元应力
两单元平均法：把两个相邻单元中的常应力加 以平均，用来表示公共边界中点处的应力。 绕结点平均法：把环绕某一结点的各单元常应 力加以平均，用以表示该结点的应力。在内结 点效果较好，而在边界结点可能很差，一般改 为应由内结点的应力外推计算出来。 （2）网格的细分 通过网格的细分，使每个单元的面积缩小，那 么尽管每个单元是应变、常应力单元，仍可较 好地反映结构中的应力变化，使得到的解答收 敛于问题的精确解。
2
x
○ ○
1 2 3 4
连接数组： 1号单元： LM=[0, 0, 1, 0, 0, 2] 2号单元： LM=[0, 0, 2, 3, 4, 5] 3号单元： LM=[3, 4, 5, 6, 7, 8] 4单元： LM=[0, 0, 1, 6, 7, 8]
0 0 3 6 ID 0 0 4 7 1 2 5 8
剪应力 xy计算结果 ( MPa)
考察点y(m) 有限元结果 弹性力学结果 误差 -1.25 -0.75 -0.25 0.25 16.2 31.2 37.2 33.7 10.9 26.7 34.6 34.6 5.3 4.5 2.6 -0.9 0.75 20.7 26.7 -6.0 1.25 3.6 10.9 -7.3
10
k12 k 22 k 32
k13 u1 1 k11 1010 k 23 u 2 F2 u k 33 F 3 3
10
k11 10 u1 k12 u2 k13u3 1 k11 10
A(9) A(8) A(7) A(6)
A(11) A(10)
A(15) A(14) A(13) A(12)
A(17) A(16)
A( 21) A( 20) A(19) A(18)
MAXA 1 2 4 6 10 12 16
18 22
变带宽存贮：按列存贮方式。从左到右，逐列存 放；对每一列，先存主对角线元素，然后由下而上顺 序存放，直到顶线下第一个元素为止。为避免混淆， 我们把存贮[K]的一维数组称为[A]。 实现变带宽存贮的关键问题是：总刚中元素 Kij在 一维数组 A中的地址是什么？为此，需要知道主元Kii 在A中的位置和相应列高hi。 主元位置：采用一个一维数组 MAXA 存主元在 A中 位置。 MAXA =[1,2,4,6,10,12,16,18，22]。 列高hj：第j行的左带宽。
从第j列的主对角线元素起到该列上方第一个非零 元素为止，所含元素的个数称为第j列的列高，记为hj ； 如果把第j列上方第1个非零元素的行号记为mj，则第j 列的列高为 hj = j - mj + 1 其实，hj就是第j行的左带宽，因而必有 UBW= max(hj)
j=1,2, …,N
利用节点位移信息数组 ID （去约束后节点位移自 由度编码），可容易地确定刚度矩阵 [K] 任何一列的 列高。
第2章 有限元程序设计方法
2.1 程序基本框图
开始
1、输入基本数据（结构描述）： 输入基本数据 （1）控制数据：如结点总数、单 计算单元刚度矩阵 元总数、约束条件总数等； （2）结点数据：如结点编号、结 形成总体刚度矩阵 点坐标、约束条件等； 形成结点荷载向量 （3）单元数据：如单元编号、单 元结点序号、单元的材料特性、 引入约束条件 几何特性等； 求解方程组，输出结点位移 （4）载荷数据：包括集中载荷、 分布载荷等。 计算单元应力，输出结果
y为挤压应力，属次要应 力。材料力学不考虑 , FEM与弹性力学误差较大。
2.2 提高计算精度的方法
（1） 计算结果的整理 计算结果包括位移和应力两个方面。在位移方 面，一般无须进行整理工作。应力结果则需要 整理。通常认为计算出的应力是三角形单元形 心处的应力。而相邻单元之间的应力存在突变， 甚至正、负符号都不相同。为了由计算结果推 算出结构内某一点的接接实际的应力，必须通 过某种平均计算。通常可采用两单元平均法或 绕结点平均法。
对 称
0 0 0 0 K 58 0 K 78 K88
顶线以上零元素无须存贮，仅顶线以下元素。
一维数组[A]存贮刚度矩阵[K]
K11 K
K12 K 22
0 K 23 K 33
K14 0 K 34 K 44
技巧要求较高，程序较长。
方阵形式的刚度矩阵[K]
UBW=4
顶线
K11 K
K12 K 22
0 K 23 K 33
K14 0 K 34 K 44
0 0 0 K 45 K 55
0 0 K 36 K 46 K 56 K 66
0 0 0 0 0 K 67 K 77
引入约束条件 求解方程组，输出结点位移 计算单元应力，输出结果 结束
[K]的存储；约束引入；求解
总刚存贮
全矩阵存贮法：不利于节省计算机的存贮空间， 很少采用。K[i,j] 对称三角存贮法：存贮上三角或下三角元素。 半带宽存贮法 ：存贮上三角形（或下三角形） 半带宽以内的元素 。 一维压缩存贮法 ：半带宽存贮中仍包含了许多 零元素。存贮每一行的第一个非零元素到主对 角线元素。
例：求图示框架结构h7=?。 利用ID数组得各单元的连接数组LM（假定小号为i） （1）ID数组
1 节点号：
2
3
4
Y
1 1 0 0 1 1 0 0 0 0 0 0
按列，遇1变0，遇0加1。
7
6
8
4
③ 3
4
④ 1 1
○
5 3
②
① 2
○ ○ ○
0 0 3 6 ID 0 0 4 7 1 2 5 8
Y
7
6
8 4 ③
4
3
②
5 3
④
1 1
○ ○ ○
①
2
○ ○ ○
2
a) 如果
ID（i, j ) = 0 则表明j号节点第i个自由度受有约束。 b) 如果
0 0 3 6 ID 0 0 4 7 1 2 5 8
ID( i, j ) ≠ 0 则j号节点第i个自由度不受约束。并且， j号节点第i 个位移分量在非约束节点位移列向量 f中的序号 就是: ID( i, j )
hj——刚度矩阵[K]第j列的列高。
一维数组A的总长度(S) ，即刚度矩阵K按变带宽存贮 的总存贮量 S = MAXA(N+1)- MAXA(1)
Ki,j在一维数组[A]中的地址 记Ki,j在一维数组A中的地址为AIJ。则由下图可知， A I J = MAXA[ J ] + J – I 其中，I = mj，mj+1，…，J。
UBW (d 1) ndf , ndf为一个结点的自由度数
结点编号：欲使最大半带宽UBW最小，必须注 意结点编号方法，使直接联系的相邻节点的最大点 号差最小。
例：计算下图半带宽。
结点数N=91，总刚[K]中的元素总数为： 82（91×2）×（91 ×2 ）=33124 最大半带宽UBW=(7+1) ×2=16，半带宽存储矩阵元素总数为 182 ×16=2912，约方阵元素的8.8%。
u1 1
5、线性方程组求解
求解方法常用：GAUSS消元法，QR分解法等。
其程序在一些专著中列出（例如见：徐士良编。
FORTRAN常用算法程序集。清华大学出版社，
1992）。在此不作详细介绍，其方法参阅[数值分
析]有关书籍。
6、单元应力
节点位移求单元应力。首先整体节点位移变换成单元 节点位移，然后再用物理方程求单元应力。
（1）半带宽存贮法
行号
UBW
行号
UBW
1 → 0 0 0 0 0 0 IR → 0 0 N→
列号 1
JC
1 → IR → N→ 1
（2） 变带宽存贮（一维压缩存贮）
等带宽存贮虽然已经节省了不少内存，但认真 研究半带宽内的元素，还有相当数量的零元素。在 平衡方程求解过程中，有些零元素只增加运算工作 量而对计算结果不产生影响。如果这些零元素不存、 不算，更能节省内存和运算时间，采用变带宽存贮 可以实现（也称一维数组存贮） 。变带宽存贮编程
{ } [S ]{ }
例1：对角受压的正方形薄板，载 荷沿厚度均匀分布，为2N/m。由 于对称性，取1/4部分作为计算对