有限元讲义.doc

§ 1.4 协调、非协调、广义协调及分电检验
1、4、0 引
以有限元数值分析的技术实现为目的本门课程，不仅要求学生能够进行实际的工程运算；另一方面也需要对解的收敛及精确性有所了解，是能从细节计算到理论性质都有所把握，这样，才能做到全面深入有助于对解结果得理论分析，此为基本之目的。
1、4、1 协调、非协调介绍
位移法有限元以Ritz的结构最小有限元为基础，该原理在数学上是一个泛函极值（变分）问题，系统势能可以表为以下数学形式：
面力势能
内力势能
体力势能
==1/2 - - (1)

＝0 。
表述为：在所有满足内部连续性和运动学边界条件的位移中满足平衡方程的位移使系统势能取驻值。如果驻值是极小点的，则平分行是稳定阶。
又：对于精确于问题的位移函数，系统势能的变分可求得关于问题应满足的所有微分方程：平衡方程边界条件（几何关系及物理方程是自然满足的）
遗憾的是精确位移难得寻找，故一般采用泛函的极小化序列逼近方法。类似于傅立叶级数逼近函数那样，把无穷维空间用有限空间去逼近。在有限元当中，当元素尺寸趋近于0时（即节点数目或节点自由度数趋于时），最后的解答若能无限逼近准确解，那么这样的位移函数（或形状函数）就称为收敛的，因此从收敛性及算收敛速度方面提出几点对形状函数的要求：
①、函数本身及其导数应在元素上连续，并含有常数部分；
②、元素之间的位移协调，不仅节点处的位移应当协调，沿整个内边界上的位移也应当协调（或称相容 ）。
③、多项式的项数越多越好，因用高次比低次多项式收敛快。
④、含有刚体位移（平动包含常数项，转动包含线性项）。
协调之： 即满足①、②条件的形状函数的元素，当然能满足3) 4)条件协调 元的收敛率就更高。

①
协调元的性质：
②
①
1) 能够以单调趋势逼近于正确解。如曲线①.
③
2) 势能总是大于最小状态，故解得上界。
④
3) 近似刚度k偏大，即元素偏“硬”。
②
4) 近似的位移偏小，即求得位移的下界。
能够以单调趋势逼近于正确解。如曲线②
.
势能总是大于最小状态，故解得上界。给点数

近似刚度k偏大，即元素偏“硬”。
近似的位移偏小，即求得位移的下界
非协调元：在弹性力学中，如板弯曲，相邻元素不仅要求位移本身连续，而且要求位移的导数连续（板弯边界上的相容性）。而在工程上能够保证导数相容的形变往往难以找到，以致工程上只能采用违反相容原则的一些形状函数，由违反相容原则的形函所构成的元素称为非协调元。
非协调元性质：
①、不能以最小位能原理作为它的理论基础。
②、解的趋势可能收敛，可能不收敛，（取决于网格划分）。对于收敛的趋势也未必满足单调性。可能收敛曲线如图中②。
Q
Q

m
j
i
k
例： 三角板元： 节点位移系统：｛｝＝

位移形函：
W=CL+CLL+CLL+CL+CLL+CLL+CL+
CLL+CLL+CLLL
L =( a+ bx + dy )/2
a= xy - xy ; b= y- y ; d= - x+ x 。
代入节点位移参数：
可得：w=[ N … N ]{}
变换一下写法：w=[ H HH H H H H HH H ] {}
H = L(3 – 2L) –7LLL
H = L ( dL– dL ) + ( d – d )LLL
H = L ( bL – bL ) + ( b – b) LLL
i +2=k
i

j k
H = 27 LLL
下标按循环计算
上述元素只能在结构上做到位移导数连续，在边界上其他点处，位移的法向导数并不连续，这因为：由于法向导数是一个完全的二次多项式，在元素的每条边上，其变化规律位一条二次抛物线，需要三个点上法向导数的相等条件才能维一确定，故相邻两条曲线一般不全重合。
故所举三角板弯元为非协调元。
例 ②书P的矩形元，由于坐标的交叉双乘积（不完备），可发现不该是w 或其导数， 都是连续的，这样只要节点的这些参数相同，边界上的这些是没有问题的，但展开 N的项，可以发现xy项，或者说缺少了代表热率变形的一项，因此，作为形状函数，是不能保证向正确的解答收敛，因而是非协调元。改进方案之一，是在节点处增加节点参数，并采用完全的埃尔米特三次多项式。

§1、4、2 非协调元的排先检检验
协调元虽可以保证总位能从上往下地正确结果单调收敛，但往往过于复杂，使用麻烦。
在工程上往往使用形式简单的非协调元，自然，最小位能原理对此不再适用，那么在什么条件下，这类元素才能导致向正确解收敛呢？ Irons 提出了一个称作“拼片实验”（patch test） 的检验方法。实践表明，这种检验方法是有效的，但“拼片实验”的理论证明尚不清楚。拼片试验内容为“假设由若干元素拼成的一个任意拼片处于等应力状态，这时，其位移函数w( x , y ) 一般可用一m阶完全多项式函数 P( x , y ) 表示，（入在薄板问题中，m=2） ，而且，在这一拼片的边界上，也设置了符合等应变状态的位移边界条件。然后，将需要检验的某种元素按此条件进行计算。如果最后得到的有限元法解答能和P( x , y ) 一致，那么，称这种元素能够通过拼片试验，而通过拼片试验的元素将给出收敛的结果。注：P( x , y ) 至少应能代表各种等应变状态。
如： 拼片： 9节点参数三角变元
*拼片test实 际上成为非协调元的收敛准则

（1） ②

注释：
①、两种拼片均有满足各种等应变状态的边界条件。
②、用非协调元解上述问题。
③、对比①，不一致②。
④、一致①不一致②。
⑤、说明①通过，③的网格划分不通过，但误差不大，尚可用。
⑥、对于更不规则的网格，其误差可能更大，大到不宜采用。

§ 1、4、3 广义协调元简介
1、 概念： 广义协调元利用对变分原理的改造，使边界位移在位能条件降低，从而获得比非协调元有更好性质的单元，即对于不规则网格能通过分片检验；提高计算精度。
2、 原理： 用能量和加权线量法联合改造。
①、势能原理： ＝ － （b）
设单元边界真实变量为 ( x , y ),则在单元边界上每一点若均能满足，
　　　ｕ－= 0 ( c )
则单元协调，位移跨越单元连续，若不能保证（c）式成立，则为非协调元，其能量泛函需要做如下修正：
②　　　＝－ｕ－) d s ( d )
( d ) 式增加了单元不协调位移的能量贡献项。为Lagrange 乘子，其物理意义为单元边界的函数。( d )式可通过带约束的广义变分原理求解。
如果用加权残量法，放松单元向位移协调条件，即
　　　　　　　　　　　　　　　　　　　　　　　（ｅ）　　
式中T为权函数，以式（c）为约束条件建立单元位移形函的矩阵，可使单元向位移协调条件在某种积分意义下得到满足。
③、（e）式实际上是一个统一的形式，如取T权函数为函数，即得加权残方程为：　｜＝0　
④、全T＝n为边界应力，则权残方程为：

表示常量应力，等价于Irons的分片检验，是非协调远的收敛准则。
⑤、具体广义协调元的构造方法可参见：龙驭球著“新型有限元之引论”

　　　　　　　　§1.5　曲线坐标系及等参元概念
1、 5、0 引
　在工程有限元分析中，有几方面的原因需要需要采用高精度单元：
①、提高计算精度：
②、适应不规则的几何边界条件；
③、用较粗的单元数据准备，获得较高的计算效益；为获得高计算精度采用通常意义下的单元，本身具有某种缺陷
①、增加节点参数提高位移参数的阶次，使计算复杂，并增加求解规模
②、只适用于规则外形的结构，复杂边界难以准确描述。
因此，目前多采用等参元在二维问题中得到了极大的成功，以后，在三维问题和板壳计算中，也取得了很好的效果。
举例说明对受拉带孔板应力集中的计算结果（元与等参元对比）

单元类型

m
单元的相对边长
应力集中系数
K ＝ 3.02

误差

6
8
10
0.066
0.16
0.008
2.567
2.806
2.860
15
7.1
5.36
等参元
6
0.066

2.869
5.06

m等分
0w
等参元的优越性：
w
①、精度高（有线性，也有二次等）；
②、匹配各种复杂的直线及曲线边界；
③、构造单元方法统一，推广其它类型应用不困难。

4.刚度计算
①、公式；　〔K〕＝
　　　　　　Ｋ=
Ｋ=
Ｋ=
②、Guass消元法消去中间节点，得聚缩刚阵；

③、程序系统不再介绍。

　　　　　　　§2　总体刚阵的形成及程序设计
从有限元的计算过程可知，形成单刚、集合总刚、置边界条件、求等效载荷列阵后，才能得到位移，应变及应力，课程正是按此工作进行，最终达到我们的目的。

　　 §2.1总刚的特征及其组装时所需要效能的问题
1 总刚特性及投放过程
总体结构刚阵的形成是各单元刚阵按节点序排列的集合矩阵（对号入座）。
①、对称性：依据 Betti互易原理（讲：在结构体的两点1、2作用单位位移位移，（１）先1点作用,再2点作用　
( 2 ) 先2点作用,再1点作用

②、带状性：（刚阵元素较集中在对角线的位置）一个节点自由度对应总刚的一行（一个节点不可能和结构的所有节点相连这一行总是和相关的节点刚阵元素相关，这也是总刚的组织进程反映出的特性。
③、稀疏性：　不相关的单元（指节点没有任何联系的单元）没有共同联系的刚性元素，故在对号累加中始终为0。
1
例：　　 节点（ 1 ，4）分属两个单元，
3
2
①
故刚阵中无K系数，同理
③
（1，6）、（2，6）、（3，4）、（4，6）

4
④
②

6
5

4 )投放过程：
总刚按节点排序（不与单元节点序相关）
⑥
⑤
④
③
②
①

⑥
⑤
④
③
②
①

考虑的问题：
总的按节点排序：
① ② ③ ④ ⑤ ⑥
④
⑥
⑤
③
②
①
2、投放。
2、 考虑的问题： ①坐标转换 ；②主从节点转换；③有效方式及其效率，依据总刚的特点，寻找节约的有序方法（由于对称可仅有一半）存贮效率措施使用存贮单元的多少，用总刚规模替标度量。
总刚规模：存贮总刚单元的总数及其存贮的形式。
*存贮方式 1、三角存放 （一维存） 总刚规模：KK=

一维存时，要有一个投放算法，使用时，要有一个调用算法。

n

m
方式：可一维存，可二维存（不好）切齐即可2维存
*存贮方式2、 等带宽存贮：
KK = m x n ；m 为最大带宽。 ×××××
按节点编号的最大差计算；由于带状性质 ××××
m<
*全希疏存放：
仅存贮各单元的非0元素及分解出的非0元素。需要记每一行的首尾特征，行元素的列位置等，则要有复杂的管理机构。管理机构也消耗单元，效率以远小于其它存方式消耗单元总数为恰当。
④求解效率：
一般不计编程的效率（程序编的优劣不计），只计存贮效率和求解效率。
求解效率当然以机时计算 ，效率成本（即存贮＋求解的消耗＝成本）。求解效率与存贮效率要综合评价，但存贮效率一般放在第一位考虑。
设：执行的时间为T，存贮量为S，则一个存贮与求解的成本估价是T, S的复杂函数，即：
C( S, T ) = T * P( S ) (P 是大于1的多项式)。
重要的存贮效率，故上述函数反映了节省存贮资源的重要性，因为减少存贮要求后，也许不再需要动用外存，从而节省了数据传输的费用。
T实际上也是S的函数，存入更多的元素，也就需要更多的计算，因此，该意义上讲也要省存贮。
3、 影响总刚规模的因素：
除存贮方式的因素之外，另一个影响总刚规模的因素是节点编号的顺序（按带状存贮时，也按全稀疏存贮（编号仅是一个因素），可能影响非0元的分解时的增加策略）
① 节点编号对带宽的影响：
变带宽
KK=17
总纲形式
3
2
1
带宽的主要影响因素是节点单元间的编号，差越大，则带宽宽，若差越小，各非0元集中在对角线附近，便会使带宽减少，这一点构成了我们过行节点编号的基本原则，举例说明：
最大带宽
m = 6
KK=36

1 ×
2 × ×
3 × × ×
4 × × × ×
5 × × × × ×
6 × × × × × ×
③
①

④
②
4
5
6

m = 4
KK=24
17
1 ×
2 × ×
3 × × ×
4 × × × ×
5 × × × × ×
6 　 　 × 　 × ×

③
①
3
2
1

④

②

6
5
4

1 ×
2 × ×
3 × × ×
4 　 × × ×
5 　 　 × × ×
6 　 　 　 × × ×
5
3
1

15
m = 3
KK=18
③
①

④
②

2
4
6

由此可见，取定了简化模型之后，对节点如何编号是一个很重要的问题，可直接越小到存贮规模的大小和机时的多少，但如何把节点的编号的顺序便的较优，是一个很复杂的问题，这将在以后专门研究这里给出一个编号的原则（工程准则），即是“优先沿结构较宽的方向按坐标进行编号。
② 节点单元编号的自动形成
　　所谓节点单元编号就是给XE,ME赋值，对于一般结构的分析问题，都有很多的节点单元若仅用人力在图纸上划分单元和计算节点坐标，则是很麻烦的，也是易于出错的，因此，这部分工作尽可能地去让计算机完成。
自动编序（又称自动划分单元），当然也需要把结构的简化模型取定，然后给某个编序原则，（当然应尽量靠近编排优序的原则，当然也可以依靠算法去优化），这一原则也需要设计者事先在草图上示意性标处，以便检验这原则的实施正确性。（也可以用图显方式对照草图检验），最后即可编制程序。
举一例：
y
x
自动说明自动编序实施方法：
1
5
9
13
17

①
⑦
(13)
(19)

③
⑨
(15)
(21)
②
⑧
(14)
(20)
2
6
10
14
18

JA(JB=3)
④
⑩
(16)
(22)

3
7
11
15
19

⑤
(11)
(17)
(23)

⑥
(12)
(18)
(24)

4
8
12
16
20

JB (分4等)

SUBROU XEME( JA , JB , A , B , XE , ME , NP) NP=(JA + 1) * (JB + 1)
DIMENDION XE( NP , 1 ) ME( NP , 3 )
D 1 N = 1 NP
XE( N , 1 ) = B * ( JB - ( N -1 )/( JA + 1 )
XE( N ,2 ) = -A * ( JA + 1 ) * (FLOAT ( N - 1 )/( JA + 1) - ( N -1 )/( JA + 1 ))
D 2 I = 1 JA
D 2 J = 1 JB
KI = 2 * ( I + ( J - 1 ) * JA ) 产生偶数单元编号KI = 2 、8、14、20第一行
ME ( KI , 1 ) = I + ( J - 1 ) * ( JA + 1 ) + 1
ME ( KI , 2 ) = ME( KI 1 ) - 1
ME( KI , 3 ) = ME( KI , 1 ) + JA + 1
ME( KI - 1 , 1 ) = ME ( KI , 1 )
ME( KI - 1 , 2 ) = ME ( KI , 3 ) - 1
ME( KI - 1 , 3 ) = ME ( KI , 3 )
RETURN
END
§2.2 总刚压缩存贮的实现
由上述分析可知，变带宽存贮 KK 最小，显然最省机时，但由于带宽是变的，故只能用一维数组依次存放总刚阵中各行的第1 个非0元素到对角线的各元素，令该数组称为AKM(kk).
要使用AKM(kk)，还必须另外一个整形数组来记忆原总刚阵中〔K〕中每一行（对应一个自由度） 的非0 元素到对角线的个数，该数组计为KDKM(NID),NID 为总自由度数。
在实用中，我们并不是真正去记忆每行的非0 元素数，而是记忆每行非0元素的累加数，也即每行最后一个元素在一维数组中的位置，这样给使用带来了很大的方便。所以在机器程序中记忆位置的编号，经常采用后一种方法。按这种方法，例2的KDKM（I）的元素应为1，3，5，9，13，17（＝KK）。而第3个例题则应为1，3，6，9，12，15（＝KK）。显然这时KDKM( NID ) = KK (最后一个数)。
1． 计算原理： P，①，②，③ 点。
解释：行快是以每一个节点所对应的自由度数为行数，列号从1到最后一列，2维是 Row = 2 ; 3 维是 Row = 3 .
A、 只有与节点之构成的相关点（与节点相关的那些单元的编号），才在总刚
阵中的第i 行块上提供非0子块。
B、 第i行块上非0子块的分布情况与节点 i的相关节点编号的差值有关。带宽对于下 阵，则是相关节点最小编号列对角线子块。对上阵存贮，则是相关节点中最大编号列对角线子块（从右到左）。
2、框图：书P 。