四、弹性力学平面问题
的有限元分析及程序
四、弹性力学平面问题
的有限元分析及程序
引 言
常应变三角形单元
矩形双线性单元
平面问题程序 (一 )
平面等参数单元
平面问题程序 (二 )
Wilson 非协调元
杆系问题以结点作为分割单元的“结点”是很自
然的，但对于平面问题，待分析物体是连续的，并
不存在实际结点。要将物体“拆”成单元，必须用
一些假想的线或面作人为地分割。实际计算时，可
将连续体分成多种形状单元，为讨论简单，现暂时
规定只用一种单元来分割。
平面问题有限单元法可用的单元很多，作为初学，
先介绍两种最简单的单元：三角形和矩形。然后再
介绍高级些的单元“等参数单元”。
将物体进行分割时，必须保证相邻单元具有公共
边界。假定相邻单元仅在一些点（顶点或顶点加边
中点）相连接。这些点即为“结点”。
4.1 引 言
4.2 常应变三角形单元
4.2.1 面积坐标
三角形单元中任一点 P可用直
角坐标 (x,y) 表示。
P
21
3
y
x
如图所示连 P1,P2,P3，则
可得三个小三角形。它们和大三
角形 ?123的面积比，记作 Li（ =
?Pjk/ ?123），称为面积坐标。
三个面积坐标显然 L1+ L2 + L3= 1，只有两个是独
立的。三角形中任一点 P的 位臵可用面积坐标 L1,L2
确定。
当 P点在 1时 L2 = L3= 0,L1= 1。余类推。可见面积
坐标具有“形函数”的性质。
4.2 常应变三角形单元
4.2.2 位移模式
由于面积坐标有形函数性质，
因此根据试凑法可得
P
21
3
y
x
形函数 = Ni=Li = 面积坐标
1) 面积坐标和直角坐标关系
如果结点 i 位移为 ui,vi，则
单元位移模式（位移场）为
u=? Niui ； v=? Nivi
?? 2
1
1
1
1232
33
22
11
???
yx
yx
yx
33
221
1
1
1
2
1
yx
yx
yx
L
?
?
11
332
1
1
1
2
1
yx
yx
yx
?
22
113
4.2 常应变三角形单元
2) 矩阵表达
P
21
3
y
x
)(2 1 iiii cxbaL ??? ?
kjkji xyyxa ??
kji yyb ?? kji xxc ???
321,,k,j,i ?
? ? ? ? ? ? ? ?? ?TT3T2T1 dddd e ?
? ? ? ?iii vud ?T
? ? ? ? ? ? ? ?? ?321 NNNN ? ? ? ? ? ? ?? ?
edNvud ??
T
? ? ?
?
?
?
?
?
?
i
i
i L
L
N
0
0
321,,i ?
4.2 常应变三角形单元
4.2.3 单元列式
1) 微分算子矩阵
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
xy
yx
A
0
0
2) 应变、应力矩阵
? ? ? ? ? ? ? ?? ?edBdA ?? T?
? ?D弹性矩阵：
? ? ? ?? ? ? ?? ?edSD ?? ??
? ? ? ? ? ? ? ?? ?321 BBBB ?
? ?
?
?
?
?
?
?
?
?
?
?
?
ii
i
i
i
bc
c
b
B 0
0
2
1
?
? ? ? ?? ? ? ?
?
?
?
?
?
?
?
?
?
?
??
ii
i
i
ii
bc
c
b
DBDS 0
0
2
1
?
平面应力问题
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
ii
ii
ii
'
i
bc
cc
bb
DS
2
1
2
1 ??
?
?
式中
)-(12 2??
E'D ?
平面应变时
???? ??? 1-1 2 ;EE
4.2 常应变三角形单元
由此可见，单元应变、应力都是常量。
当所分析的问题具有初应变时，单元的弹性应变
为 [e]=[? ]-[? 0]，应力为 [?]=[D][e]。
3) 单元应变能
? ? ? ???
?
?? dA
2
1 T tU
? ? ? ? ? ?? ?ee dBSdtU TT
2
??
将上述应变、应力代入
4) 单元外力势能
? ? ? ?
? ? ? ? ? ? ? ?d l )(
dA
TT
T
?
?
??
??
ijl
ee
bf
dtdF
dFtP
?
?
第一项体积力、第二项结
点力、第三项表面力的外
力势。
? ? ? ?
? ? ? ? ? ? ? ?e
l
e
bf
dNtF
NFtP
ij
d l ) ](
dA[
TT
T
?
?
??
??
?
?
代入位移后，经整理可得
4.2 常应变三角形单元
5) 令总势能一阶变分等于零，推导单元刚度方程
当有初应变时推导结果如何？
? ? ? ? ? ?BStk e T??
? ? ? ? ? ?jieij BStk T??
321,,j,i ?
? ? ? ? ? ?
? ? ? ?d s )(
dA
T
T
?
?
?
?
ijl
be
Nt
FNtP
?
?
6) 单元刚度矩阵、等效荷载矩阵当有初应变时结果如何？
具体显式表达式见教材
P。 47 式（ 3,2-39）
? ? ? ?? ? ? ? ? ?
? ? ? ? ? ? 0ds )(
dA
T
TT
???
??
?
?
ijl
be
NF
FNdBSt
?
???
?
? ? ? ?? ? ? ? ? ? ? ?
? ? ? ?ds )(
dA
T
TT
?
?
?
??
ijl
be
Nt
FNtFdBSt
?
?
?
? ? ? ?? ? ? ? ? ? ? ?
? ? ? ?? ? ? ? ? ?ds )dA (
dA
T
0
T
TT
??
?
??
??
ijl
be
NtDBt
FNtFdBSt
??
?
?
?
4.2 常应变三角形单元
7) 关于等效结点荷载
等效结点荷载可用公式积分计算，但由于形函数
的图形是一平面（边界处为一直线），因此可证明也
可按杠杆原理通过静力等效来求。
如 P.48 图 3-4所示。
4.2.4 解答的收敛性准则
1) 位移模式（也称位移函数）必须包含刚体位移。
2) 位移模式必须包含常应变位移。
3) 位移模式必须保证单元间位移协调。
1),2) 对平面问题也即要求具有常数项和坐标一
次项，这称作“完备性准则”。
3) 称作“协调性准则”。既完备又协调的单元一
定是收敛的。但不等于说非协调单元一定不收敛。
4.3 矩形双线性单元
三角形单元划分灵活，能较好拟合边界复杂（如
曲线边界）物体计算。但是，单元应变、应力是常量，
对一般问题精度较低，要提高精度就的增加结点、增
加未知量，为此讨论其他单元。其一为本节单元。
1) 自然坐标
2a
2b
2
2
图示矩形单元，设 ?=x/a，
?=y/b，则转换成正则单元。
2) 形函数
2
2
2
3 4
1
?
?
由形函数的性质“本点 1，它点零”，利用试凑法
可设,N1=a（ 1- ?)(1- ?）它满足“它点零条件”。
再令本点为 1，可得 a =1/4，代回可的形函数 N1。
同理可得,Ni=1/4(1+?0)(1+?0) (i=1,2,3,4)。
式中 ?0 = ??i ； ?0=??i 。
请大家验证 Ni是否满足形函数性质。
4.3 矩形双线性单元
3) 位移模式
u=?Niui ； v=?Nivi。
或以矩阵表示为
2
2 2
3 4
1
?
?
可以用势能原理，也可以用虚位移原理。一经建立
单元位移模式后，剩下的工作和杆系、三角形单元类
似，因此这里从略。
? ? ? ?? ?edNd ?
? ? ? ? ? ?? ?41 NNN ????
? ? ?
?
?
?
?
?
?
i
i
i N
N
N
0
0
[d]=[u,v]T
单元结点
位移矩阵
4) 关于单元列式
5) 关于单元特性结果请看 P,53 式 (3,6-13~15)。
6) 关于计算结果的整理
里兹法已经知道：位移结果比应力、内力结果精
度高。位移达到满意结果，有几何方程求应变，再
由物理方程求应力，结果精度较差。上述三角形单
元常应力，矩形单元应力线性变化，许多工程问题
的应力是复杂的。为更好标征性，需要对计算结果
进行整理。常用处理方法有两种。
4.3 矩形双线性单元
6-1) 绕结点平均法
以交于同一结点各单元此结点处某应力分量的代
数平均值，作为此结点该实际应力的近似值。
对于边界处的结点，由内结点结果的外插得到。
6-2) 两单元平均法
三角形单元时,以两相邻单元应力平均值作为边
中点的应力近似值。 矩形单元时,以两相邻单元公
共边两端结点四个应力的平均值作为边中点的应力
近似值。对于边界处的结点，同样由内结点结果的
外插得到。
1) 程序功能
本程序可用三角形或
矩形单元计算平面应力
问题。当计算平面应变
问题时需要自行转换弹
性常数。
本程序为了减少计算
数据的准备，对规则问
题具有做、单元结点编
码等自动生成功能。
本程序荷载生成功能
较弱，请自行修改。
4.4 平面问题程序 (一 )
本程序可以用来计算
如墙梁、剪力墙（可以
带孔洞）等结构。
2）程序数据文件说明
2-1）基本数据
结点位移数，单元结点
数，结点总数，最大半
带宽，总约束位移码数，
单元总数，点的坐标数，
规则标志，问题标志，
弹性常数及厚度。
2-2) 结点坐标
如果不规则,按结点号
顺序读入全部结点的坐
标值。
如果规则无孔：孔标志,
X方向单元数,Y方向单
元数,X方向单元长度，
Y方向单元长度。
如果规则有孔：
控制结点数，生成结点
类数。
结点号,X,Y
4.4 平面问题程序 (一 )
起点号，终点号，生成
的点数，相邻点号差值，
“相邻两点间距”。
2-3）读入结点荷载值
有荷载的结点数
结点号,X方向荷载值,
Y方向荷载值。
2-4）单元的整体结点码
如果规则无孔,不需要
输入
如果规则有孔：
待修改
如果不规则：
按单元类型读入单元
整体位移码。
2-5) 读入全部零位移约
束的位移码。
如果问题类型不等于零
2-6）读入第二种材料的
弹性常数，厚度，首单
元号，终单元号，循环
步长。
如果问题类型等于零
没有第六组
4.4 平面问题程序 (一 )
2-7）一算例数据
2,4,441,46,22,400,2,0,0,6.93e4
,0.3,2.0e-2
1
20,20,0.05,0.1
21
421,0.0,0.05
422,0.0,0.1
423,0.0,0.1
424,0.0,0.1
425,0.0,0.1
426,0.0,0.1
427,0.0,0.1
428,0.0,0.1
429,0.0,0.1
430,0.0,0.1
431,0.0,0.1
432,0.0,0.1 结
433,0.0,0.1 点
434,0.0,0.1 荷
435,0.0,0.1 载
436,0.0,0.1 数
437,0.0,0.1 据
438,0.0,0.1
439,0.0,0.1
440,0.0,0.1
441,0.0,0.05
运行程序查看计算结果
1,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42
4.5 平面等参数单元
三角形和矩形单元是最简单的单元形式，前已提
及，由于位移模式是线性和双线性的，精度较低。对
于曲线边界问题，还有以直边代替曲边的离散误差。
为此，介绍本节等参数单元族。
首先以四结点等参元为例进行介绍，然后再介绍
其他等参元。
2
2
2
3 4
1
?
?
x
y ?
?
1
2
3 4
母单元
4.5.1 单元描述
为克服矩形单元不能拟合曲线边界，用图示任意
四边形单元，但在直角坐标下要描绘属于单元的点比
较困难。
自然坐标下的图示单元（母单元）形状都是规则图
形，考虑到母单元形函数性质，并设 i 点的坐标为：
(xi,yi)，则由 x=?Nixi ； y=?Niyi 可以描绘子单元形
状。 为什麽？
子单元（等参元）
4.5 平面等参数单元
x
y ?
?
1
2
3 4
母单元
利用上述转换公式可看出：
1) 母单元正交坐标线映射后成为图
示子单元斜角坐标线。仍为直线。
2) 子单元有两套坐标系：整体 x,y坐
标，和局部 ?,?坐标。
3) ?,?坐标又是母单元正交坐标。
4) 根据子单元结点坐标情况，规则
母单元可映射出任意四边形单元。
5) 相邻单元映射后仍然连续。
子单元（等参元）
2
2
2
3 4
1
?
?
4.5.2 单元位移模式
设 i 点的位移为,(ui,vi)，则单
元位移场为 u=?Niui； v=?Nivi。
和矩形单元一样，
可用矩阵表示为,? ? ? ?? ?edNd ?
4.5 平面等参数单元
x
y ?
?
1
2
3 4
母单元
4.5.3 坐标系间的转换关系
2
2
2
3 4
1
?
?
y
y
x
x
y
y
x
x
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
???
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
yyy
xxx
由复合函数求导数的规则可得用矩阵表示则为：
?
?
?
??
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
??
?
?
?
?
?
?
y
x
yx
yx
??
??
?
?
?
?
?
??
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
??
?
?
?
?
?
?
?
?
??
??
yy
xx
y
x
引入记号 称作雅可比矩阵其逆矩阵为：
4.5.4 其他等参元母单元形函数
1) 八结点等参元（矩形族）
5
6
7
8
子单元（等参元）
6
7
8
5
用试凑法 N1可设为：
N1=a(1-?)(1-?)(1+?+?)；它能自
动满足它点为零。本点为 1得：
a=-1/4
同理可得角结点 Ni通式为：
i -1/4(1+?0)(1+? 0)(1-? 0 -? 0)
(i=1,2,3,4).式中 ?0 = ??I;? 0=??i,
用试凑法 N5可设为：
N5=a(1-?2)(1-?) ；它能自动满足
它点为零。本点为 1得,a=1/2
同理可得边中点 Ni通式为：
N5,7=1/2(1-?2)(1+?0)
N6,8=1/2(1-?2)(1+?0)
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
yx
yx
J ? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
yy
xxJ
??
??
1? ?
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
??
??
xx
yy
J
y
xx
J
11
映射后
子单元可以
是曲边单元
4.5 平面等参数单元
2) 六结点 (三角形族 ) 单元如图所示
3 1
2
1
1
4
56 L
1=?
L2=?
母单元
子单元
3 1
2
4
56
?
?
x
y
设 N1=a?(?-1/2)；它满足它点为零
条件，为使本点为 1，可得 a =2。同
理可得
Ni=Li(2Li-1) (i=1,2,3)
设 N4=aL1L3；它满足它点为零条
件，为使本点为 1，可得 a =4。由此
可得
N4=4L1L3
N5=4L2L1
N6=4L3L2
映射后
子单元可以
是曲边单元
4.5 平面等参数单元
3) 常用的两族单元
矩形：四、八、十二结点等参元。
三角形：六、十结点等参元。
问题：试建立十、十二结点母单元形函数。
各种单元分别包含坐标几次完整多项式？
能否给出建立两族形函数的一般方法？
4.5.5 单元描述和位移模式
1) 单元描述 2) 位移模式
? ? ? ? ? ?? ?TTT1 ne xxx ????
? ? ? ? ? ?? ?edNvud ?? T
? ? ? ? ? ?? ?nNNN ???? 1
? ? ?
?
?
?
?
?
?
i
i
i N
N
N
0
0
? ? ? ? ? ?? ?TTT1 ne ddd ????
? ? ? ? ? ?? ?exNyxx ?? T
结点数
4.5 平面等参数单元
4.5.6 等参元单元分析
记微分符号,
?? ?
??
?
?? D;
xD x
记
iyiyixix NDD;NDD ??
1) 应变矩阵
? ? ? ? ? ? ? ?? ?edBdA ?? T?
? ? ? ? ? ?? ?nBBB ???? 1
? ?
?
?
?
?
?
?
?
?
?
?
?
i
x
i
y
i
y
i
x
i
DD
D
D
B 0
0
?? ??
i
ii
i
ii yNy;xNx
? ? ? ? ?
?
?
?
?
?
?
?
? ?? xDxD
yDyD
JJ
??
??11
? ?
)
(1
i
i
i
x
NyDD
NyDDJD
??
??
?
? ?
? ?
)
(-1
i
i
i
y
NxDD
NxDDJD
??
??
?
? ?
4.5 平面等参数单元
2) 应力矩阵
? ?D弹性矩阵：
? ? ? ?? ? ? ?? ?edSD ?? ??
3) 单元应变能
? ? ? ???
A
tU dA
2
1 T ??
4) 单元外力势能
? ? ? ?
? ? ? ? ? ? ? ?e
l
e
A
bf
dNtF
NFtP
ij
d s ) ](
dA[
TT
T
?
?
??
??
?
5) dA,ds 的计算
为用虚位移原理推导
单元刚度方程，必须解决
dA,ds 的计算。
母单元规则微元体 d?d?
映射后变成图示（曲边）
四边形。
x
y ?
?
1
2
3 4
5
6
7
8
4.5 平面等参数单元
如图示，此微面积为
x
y ?
?
1
2
3 4
5
6
7
8
? ?
? ? ??
????
dd
dddddA
Jd e t
J
?
???
坐标的积分上下限均为 -1,1。
沿边线的积分（ ?=1为例）
?
??
?
??
??
d][dds 1 / 2
2
1
2
1
11
??
?? ?
??
?
??? yx
一般情况见 P.72 式 (3,6-23)。
? ? ?
?
?
?
?
?
?
xDxD
yDxD
J
??
??
??
i
ii xNx
??
i
ii yNy
4.5 平面等参数单元
有了上述结果，经虚位移原理或势能原理即可推
得式 (3,6-24)~(3,6-26)单元刚度和等效荷载结果。
4.5.7 数值积分
三角形和矩形单元可以写出刚度显式表达式，但
对于等参元，由于两套坐标的转换，导致刚度、荷载
的被积表达式十分复杂，一般不可能积出显式结果。
只能用数值积分由程序来得到。
目前常用的是高斯积分（矩形族）和哈默尔积分
（三角形族）。它们的积分点位臵、加权系数等见表
3-1,3-2（ P.74~76)。
其积分公式见式（ 3,6-40）、（ 3,6-41）。
4.5 平面等参数单元
4.5.8 作等参元分析时应注意的问题
等参元分析中要用 det[J]-1,可见雅可比行列式等于
零将导致刚度矩阵等无法积分，使分析失效。因此要
避免以下可能使 det[J]=0的情况：
1) 子单元边界不能过于扭曲。
2) 矩形子单元不能退化成三角形。
3) 子单元角顶处单元边线切线角角不能等于 1800。
上述情况如 P.81 图 3-39 示意。
此外,子单元边界上结点应尽可能是或接近等分点，
避免产生奇异单元。
可能情况下应采用直边子单元，这样可使雅可比
矩阵简单，提高计算效率。
4.5 平面等参数单元
4.5.9 离散化时应注意的问题
除对等参元应注意上述问题外，任何有限元分析
都还应注意以下几点：
1) 相互邻接的单元大小应尽可能均匀。
2) 单元最大尺寸与最小尺寸之比应尽可能接近一，最
多不应大于二。
3) 应合理编码，使单元结点间的整体编号差值最小。
4) 应尽可能使各界点的单元数目相同,如 P.82 图 3,6-42
左图示意,
4.6 平面问题程序 (二 )
应用本程序时,数据文件的准备见 P.217。
运行程序查看计算结果
4.7 Wilson非协调单元
至今为止所介绍的单元都是能保证收敛的协调单
元。单元计算结果的精度，取决于位移模式中坐标完
全多项式的次数。为改进精度 Wilson提出非协调的单
元，简单介绍如下。
在矩形双线性单元基础上，增加一位移修正项：
? ? ? ?? ? ? ?? ? ee aNdNd ??
? ? ?
?
?
?
?
?
?
21
21
00
00
NN
NNN
? ? ? ?vvuue aaaaa 2121?
21 1 ???N
22 1 ???N
显然在四个结点处修正项等于零，因此它只影响单元
内部位移，可见 [a]e是单元内部位移参数。所求得的
结点位移将仍直接是实际问题的近似值。
4.7 Wilson非协调单元
根据这一位移模式，可得应变矩阵为
? ? ? ?? ? ? ?? ? ee aNdNd ???
经虚位移原理做单元列式，可得
? ? ? ?? ? ? ?? ? ee aBdB ??? ?
? ? ? ?
? ? ? ?
? ?
? ?
? ? ? ?
? ? ??
?
?
?
??
?
?
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
a
e
ee
e
e
a
e
P
PF
a
d
kk
kk
T
1
1
双线性
修正项
相关项
相关项
结点力
双线性
等效
修正项
等效
各项的具体计算公式见下页。
4.7 Wilson非协调单元
? ? ? ? ? ?? ???
Ve DDBk dV
T
双线性
修正项
等效
? ? ? ? ? ?? ???
V BDBk dV
TT
1
相关项
? ? ? ? ? ?? ???
Va
BDBk dVT
修正项
? ? ? ? ? ? ? ? ? ??? ??
?
?
SV b
a
e NFNP dSdV
TT
? ? ? ? ? ? ? ? ? ??? ??
?
?
SV be
NFNP dSdV TT
双线性
等效
4.7 Wilson非协调单元
由于 [a]e是单元内部位移参数，与单元集合体（整
体）的其他单元无关，因此在集成之前可做如下消去
? ? ? ?
? ? ? ?
? ?
? ?
? ? ? ?
? ? ??
?
?
?
??
?
?
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
a
e
eee
a
e
P
PF
a
d
kk
kk
T
1
1
? ? ? ? ? ? ? ? ? ? )( T11 eaeae dkPka ?? ?
? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ?aeaee
eae
PkkPF
dkkkk
1
1
T
1
1
1 )(
?
?
??
??
单元刚度
单元等效荷载
4.7 Wilson非协调单元
由此可见做了上述消去（也称静力凝聚）后，单元刚
度方程形式上和双线性单元完全一样。
问题：试增加 PSTE 程序的单元库，使其能用非
协调元进行分析。
能否在四边形等参元基础上来作非协调元？
协调元是肯定收敛的，非协调元收敛性不能保证。
下一章将介绍如何检查非协调元的收敛性。这里只给
出结论：
1) 当单元网格为矩形或平行四边形时，分析将试收敛
的。
2) 任意四边形时，以 ?=?=0点的 [J]作为单元的各点
的 [J]，这样计算的修正单元刚度也能使分析收敛。
的有限元分析及程序
四、弹性力学平面问题
的有限元分析及程序
引 言
常应变三角形单元
矩形双线性单元
平面问题程序 (一 )
平面等参数单元
平面问题程序 (二 )
Wilson 非协调元
杆系问题以结点作为分割单元的“结点”是很自
然的，但对于平面问题，待分析物体是连续的，并
不存在实际结点。要将物体“拆”成单元，必须用
一些假想的线或面作人为地分割。实际计算时，可
将连续体分成多种形状单元，为讨论简单，现暂时
规定只用一种单元来分割。
平面问题有限单元法可用的单元很多，作为初学，
先介绍两种最简单的单元：三角形和矩形。然后再
介绍高级些的单元“等参数单元”。
将物体进行分割时，必须保证相邻单元具有公共
边界。假定相邻单元仅在一些点（顶点或顶点加边
中点）相连接。这些点即为“结点”。
4.1 引 言
4.2 常应变三角形单元
4.2.1 面积坐标
三角形单元中任一点 P可用直
角坐标 (x,y) 表示。
P
21
3
y
x
如图所示连 P1,P2,P3，则
可得三个小三角形。它们和大三
角形 ?123的面积比，记作 Li（ =
?Pjk/ ?123），称为面积坐标。
三个面积坐标显然 L1+ L2 + L3= 1，只有两个是独
立的。三角形中任一点 P的 位臵可用面积坐标 L1,L2
确定。
当 P点在 1时 L2 = L3= 0,L1= 1。余类推。可见面积
坐标具有“形函数”的性质。
4.2 常应变三角形单元
4.2.2 位移模式
由于面积坐标有形函数性质，
因此根据试凑法可得
P
21
3
y
x
形函数 = Ni=Li = 面积坐标
1) 面积坐标和直角坐标关系
如果结点 i 位移为 ui,vi，则
单元位移模式（位移场）为
u=? Niui ； v=? Nivi
?? 2
1
1
1
1232
33
22
11
???
yx
yx
yx
33
221
1
1
1
2
1
yx
yx
yx
L
?
?
11
332
1
1
1
2
1
yx
yx
yx
?
22
113
4.2 常应变三角形单元
2) 矩阵表达
P
21
3
y
x
)(2 1 iiii cxbaL ??? ?
kjkji xyyxa ??
kji yyb ?? kji xxc ???
321,,k,j,i ?
? ? ? ? ? ? ? ?? ?TT3T2T1 dddd e ?
? ? ? ?iii vud ?T
? ? ? ? ? ? ? ?? ?321 NNNN ? ? ? ? ? ? ?? ?
edNvud ??
T
? ? ?
?
?
?
?
?
?
i
i
i L
L
N
0
0
321,,i ?
4.2 常应变三角形单元
4.2.3 单元列式
1) 微分算子矩阵
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
xy
yx
A
0
0
2) 应变、应力矩阵
? ? ? ? ? ? ? ?? ?edBdA ?? T?
? ?D弹性矩阵：
? ? ? ?? ? ? ?? ?edSD ?? ??
? ? ? ? ? ? ? ?? ?321 BBBB ?
? ?
?
?
?
?
?
?
?
?
?
?
?
ii
i
i
i
bc
c
b
B 0
0
2
1
?
? ? ? ?? ? ? ?
?
?
?
?
?
?
?
?
?
?
??
ii
i
i
ii
bc
c
b
DBDS 0
0
2
1
?
平面应力问题
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
ii
ii
ii
'
i
bc
cc
bb
DS
2
1
2
1 ??
?
?
式中
)-(12 2??
E'D ?
平面应变时
???? ??? 1-1 2 ;EE
4.2 常应变三角形单元
由此可见，单元应变、应力都是常量。
当所分析的问题具有初应变时，单元的弹性应变
为 [e]=[? ]-[? 0]，应力为 [?]=[D][e]。
3) 单元应变能
? ? ? ???
?
?? dA
2
1 T tU
? ? ? ? ? ?? ?ee dBSdtU TT
2
??
将上述应变、应力代入
4) 单元外力势能
? ? ? ?
? ? ? ? ? ? ? ?d l )(
dA
TT
T
?
?
??
??
ijl
ee
bf
dtdF
dFtP
?
?
第一项体积力、第二项结
点力、第三项表面力的外
力势。
? ? ? ?
? ? ? ? ? ? ? ?e
l
e
bf
dNtF
NFtP
ij
d l ) ](
dA[
TT
T
?
?
??
??
?
?
代入位移后，经整理可得
4.2 常应变三角形单元
5) 令总势能一阶变分等于零，推导单元刚度方程
当有初应变时推导结果如何？
? ? ? ? ? ?BStk e T??
? ? ? ? ? ?jieij BStk T??
321,,j,i ?
? ? ? ? ? ?
? ? ? ?d s )(
dA
T
T
?
?
?
?
ijl
be
Nt
FNtP
?
?
6) 单元刚度矩阵、等效荷载矩阵当有初应变时结果如何？
具体显式表达式见教材
P。 47 式（ 3,2-39）
? ? ? ?? ? ? ? ? ?
? ? ? ? ? ? 0ds )(
dA
T
TT
???
??
?
?
ijl
be
NF
FNdBSt
?
???
?
? ? ? ?? ? ? ? ? ? ? ?
? ? ? ?ds )(
dA
T
TT
?
?
?
??
ijl
be
Nt
FNtFdBSt
?
?
?
? ? ? ?? ? ? ? ? ? ? ?
? ? ? ?? ? ? ? ? ?ds )dA (
dA
T
0
T
TT
??
?
??
??
ijl
be
NtDBt
FNtFdBSt
??
?
?
?
4.2 常应变三角形单元
7) 关于等效结点荷载
等效结点荷载可用公式积分计算，但由于形函数
的图形是一平面（边界处为一直线），因此可证明也
可按杠杆原理通过静力等效来求。
如 P.48 图 3-4所示。
4.2.4 解答的收敛性准则
1) 位移模式（也称位移函数）必须包含刚体位移。
2) 位移模式必须包含常应变位移。
3) 位移模式必须保证单元间位移协调。
1),2) 对平面问题也即要求具有常数项和坐标一
次项，这称作“完备性准则”。
3) 称作“协调性准则”。既完备又协调的单元一
定是收敛的。但不等于说非协调单元一定不收敛。
4.3 矩形双线性单元
三角形单元划分灵活，能较好拟合边界复杂（如
曲线边界）物体计算。但是，单元应变、应力是常量，
对一般问题精度较低，要提高精度就的增加结点、增
加未知量，为此讨论其他单元。其一为本节单元。
1) 自然坐标
2a
2b
2
2
图示矩形单元，设 ?=x/a，
?=y/b，则转换成正则单元。
2) 形函数
2
2
2
3 4
1
?
?
由形函数的性质“本点 1，它点零”，利用试凑法
可设,N1=a（ 1- ?)(1- ?）它满足“它点零条件”。
再令本点为 1，可得 a =1/4，代回可的形函数 N1。
同理可得,Ni=1/4(1+?0)(1+?0) (i=1,2,3,4)。
式中 ?0 = ??i ； ?0=??i 。
请大家验证 Ni是否满足形函数性质。
4.3 矩形双线性单元
3) 位移模式
u=?Niui ； v=?Nivi。
或以矩阵表示为
2
2 2
3 4
1
?
?
可以用势能原理，也可以用虚位移原理。一经建立
单元位移模式后，剩下的工作和杆系、三角形单元类
似，因此这里从略。
? ? ? ?? ?edNd ?
? ? ? ? ? ?? ?41 NNN ????
? ? ?
?
?
?
?
?
?
i
i
i N
N
N
0
0
[d]=[u,v]T
单元结点
位移矩阵
4) 关于单元列式
5) 关于单元特性结果请看 P,53 式 (3,6-13~15)。
6) 关于计算结果的整理
里兹法已经知道：位移结果比应力、内力结果精
度高。位移达到满意结果，有几何方程求应变，再
由物理方程求应力，结果精度较差。上述三角形单
元常应力，矩形单元应力线性变化，许多工程问题
的应力是复杂的。为更好标征性，需要对计算结果
进行整理。常用处理方法有两种。
4.3 矩形双线性单元
6-1) 绕结点平均法
以交于同一结点各单元此结点处某应力分量的代
数平均值，作为此结点该实际应力的近似值。
对于边界处的结点，由内结点结果的外插得到。
6-2) 两单元平均法
三角形单元时,以两相邻单元应力平均值作为边
中点的应力近似值。 矩形单元时,以两相邻单元公
共边两端结点四个应力的平均值作为边中点的应力
近似值。对于边界处的结点，同样由内结点结果的
外插得到。
1) 程序功能
本程序可用三角形或
矩形单元计算平面应力
问题。当计算平面应变
问题时需要自行转换弹
性常数。
本程序为了减少计算
数据的准备，对规则问
题具有做、单元结点编
码等自动生成功能。
本程序荷载生成功能
较弱，请自行修改。
4.4 平面问题程序 (一 )
本程序可以用来计算
如墙梁、剪力墙（可以
带孔洞）等结构。
2）程序数据文件说明
2-1）基本数据
结点位移数，单元结点
数，结点总数，最大半
带宽，总约束位移码数，
单元总数，点的坐标数，
规则标志，问题标志，
弹性常数及厚度。
2-2) 结点坐标
如果不规则,按结点号
顺序读入全部结点的坐
标值。
如果规则无孔：孔标志,
X方向单元数,Y方向单
元数,X方向单元长度，
Y方向单元长度。
如果规则有孔：
控制结点数，生成结点
类数。
结点号,X,Y
4.4 平面问题程序 (一 )
起点号，终点号，生成
的点数，相邻点号差值，
“相邻两点间距”。
2-3）读入结点荷载值
有荷载的结点数
结点号,X方向荷载值,
Y方向荷载值。
2-4）单元的整体结点码
如果规则无孔,不需要
输入
如果规则有孔：
待修改
如果不规则：
按单元类型读入单元
整体位移码。
2-5) 读入全部零位移约
束的位移码。
如果问题类型不等于零
2-6）读入第二种材料的
弹性常数，厚度，首单
元号，终单元号，循环
步长。
如果问题类型等于零
没有第六组
4.4 平面问题程序 (一 )
2-7）一算例数据
2,4,441,46,22,400,2,0,0,6.93e4
,0.3,2.0e-2
1
20,20,0.05,0.1
21
421,0.0,0.05
422,0.0,0.1
423,0.0,0.1
424,0.0,0.1
425,0.0,0.1
426,0.0,0.1
427,0.0,0.1
428,0.0,0.1
429,0.0,0.1
430,0.0,0.1
431,0.0,0.1
432,0.0,0.1 结
433,0.0,0.1 点
434,0.0,0.1 荷
435,0.0,0.1 载
436,0.0,0.1 数
437,0.0,0.1 据
438,0.0,0.1
439,0.0,0.1
440,0.0,0.1
441,0.0,0.05
运行程序查看计算结果
1,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42
4.5 平面等参数单元
三角形和矩形单元是最简单的单元形式，前已提
及，由于位移模式是线性和双线性的，精度较低。对
于曲线边界问题，还有以直边代替曲边的离散误差。
为此，介绍本节等参数单元族。
首先以四结点等参元为例进行介绍，然后再介绍
其他等参元。
2
2
2
3 4
1
?
?
x
y ?
?
1
2
3 4
母单元
4.5.1 单元描述
为克服矩形单元不能拟合曲线边界，用图示任意
四边形单元，但在直角坐标下要描绘属于单元的点比
较困难。
自然坐标下的图示单元（母单元）形状都是规则图
形，考虑到母单元形函数性质，并设 i 点的坐标为：
(xi,yi)，则由 x=?Nixi ； y=?Niyi 可以描绘子单元形
状。 为什麽？
子单元（等参元）
4.5 平面等参数单元
x
y ?
?
1
2
3 4
母单元
利用上述转换公式可看出：
1) 母单元正交坐标线映射后成为图
示子单元斜角坐标线。仍为直线。
2) 子单元有两套坐标系：整体 x,y坐
标，和局部 ?,?坐标。
3) ?,?坐标又是母单元正交坐标。
4) 根据子单元结点坐标情况，规则
母单元可映射出任意四边形单元。
5) 相邻单元映射后仍然连续。
子单元（等参元）
2
2
2
3 4
1
?
?
4.5.2 单元位移模式
设 i 点的位移为,(ui,vi)，则单
元位移场为 u=?Niui； v=?Nivi。
和矩形单元一样，
可用矩阵表示为,? ? ? ?? ?edNd ?
4.5 平面等参数单元
x
y ?
?
1
2
3 4
母单元
4.5.3 坐标系间的转换关系
2
2
2
3 4
1
?
?
y
y
x
x
y
y
x
x
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
???
???
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
yyy
xxx
由复合函数求导数的规则可得用矩阵表示则为：
?
?
?
??
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
??
?
?
?
?
?
?
y
x
yx
yx
??
??
?
?
?
?
?
??
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
??
?
?
?
?
?
?
?
?
??
??
yy
xx
y
x
引入记号 称作雅可比矩阵其逆矩阵为：
4.5.4 其他等参元母单元形函数
1) 八结点等参元（矩形族）
5
6
7
8
子单元（等参元）
6
7
8
5
用试凑法 N1可设为：
N1=a(1-?)(1-?)(1+?+?)；它能自
动满足它点为零。本点为 1得：
a=-1/4
同理可得角结点 Ni通式为：
i -1/4(1+?0)(1+? 0)(1-? 0 -? 0)
(i=1,2,3,4).式中 ?0 = ??I;? 0=??i,
用试凑法 N5可设为：
N5=a(1-?2)(1-?) ；它能自动满足
它点为零。本点为 1得,a=1/2
同理可得边中点 Ni通式为：
N5,7=1/2(1-?2)(1+?0)
N6,8=1/2(1-?2)(1+?0)
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
yx
yx
J ? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
yy
xxJ
??
??
1? ?
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
??
??
??
??
xx
yy
J
y
xx
J
11
映射后
子单元可以
是曲边单元
4.5 平面等参数单元
2) 六结点 (三角形族 ) 单元如图所示
3 1
2
1
1
4
56 L
1=?
L2=?
母单元
子单元
3 1
2
4
56
?
?
x
y
设 N1=a?(?-1/2)；它满足它点为零
条件，为使本点为 1，可得 a =2。同
理可得
Ni=Li(2Li-1) (i=1,2,3)
设 N4=aL1L3；它满足它点为零条
件，为使本点为 1，可得 a =4。由此
可得
N4=4L1L3
N5=4L2L1
N6=4L3L2
映射后
子单元可以
是曲边单元
4.5 平面等参数单元
3) 常用的两族单元
矩形：四、八、十二结点等参元。
三角形：六、十结点等参元。
问题：试建立十、十二结点母单元形函数。
各种单元分别包含坐标几次完整多项式？
能否给出建立两族形函数的一般方法？
4.5.5 单元描述和位移模式
1) 单元描述 2) 位移模式
? ? ? ? ? ?? ?TTT1 ne xxx ????
? ? ? ? ? ?? ?edNvud ?? T
? ? ? ? ? ?? ?nNNN ???? 1
? ? ?
?
?
?
?
?
?
i
i
i N
N
N
0
0
? ? ? ? ? ?? ?TTT1 ne ddd ????
? ? ? ? ? ?? ?exNyxx ?? T
结点数
4.5 平面等参数单元
4.5.6 等参元单元分析
记微分符号,
?? ?
??
?
?? D;
xD x
记
iyiyixix NDD;NDD ??
1) 应变矩阵
? ? ? ? ? ? ? ?? ?edBdA ?? T?
? ? ? ? ? ?? ?nBBB ???? 1
? ?
?
?
?
?
?
?
?
?
?
?
?
i
x
i
y
i
y
i
x
i
DD
D
D
B 0
0
?? ??
i
ii
i
ii yNy;xNx
? ? ? ? ?
?
?
?
?
?
?
?
? ?? xDxD
yDyD
JJ
??
??11
? ?
)
(1
i
i
i
x
NyDD
NyDDJD
??
??
?
? ?
? ?
)
(-1
i
i
i
y
NxDD
NxDDJD
??
??
?
? ?
4.5 平面等参数单元
2) 应力矩阵
? ?D弹性矩阵：
? ? ? ?? ? ? ?? ?edSD ?? ??
3) 单元应变能
? ? ? ???
A
tU dA
2
1 T ??
4) 单元外力势能
? ? ? ?
? ? ? ? ? ? ? ?e
l
e
A
bf
dNtF
NFtP
ij
d s ) ](
dA[
TT
T
?
?
??
??
?
5) dA,ds 的计算
为用虚位移原理推导
单元刚度方程，必须解决
dA,ds 的计算。
母单元规则微元体 d?d?
映射后变成图示（曲边）
四边形。
x
y ?
?
1
2
3 4
5
6
7
8
4.5 平面等参数单元
如图示，此微面积为
x
y ?
?
1
2
3 4
5
6
7
8
? ?
? ? ??
????
dd
dddddA
Jd e t
J
?
???
坐标的积分上下限均为 -1,1。
沿边线的积分（ ?=1为例）
?
??
?
??
??
d][dds 1 / 2
2
1
2
1
11
??
?? ?
??
?
??? yx
一般情况见 P.72 式 (3,6-23)。
? ? ?
?
?
?
?
?
?
xDxD
yDxD
J
??
??
??
i
ii xNx
??
i
ii yNy
4.5 平面等参数单元
有了上述结果，经虚位移原理或势能原理即可推
得式 (3,6-24)~(3,6-26)单元刚度和等效荷载结果。
4.5.7 数值积分
三角形和矩形单元可以写出刚度显式表达式，但
对于等参元，由于两套坐标的转换，导致刚度、荷载
的被积表达式十分复杂，一般不可能积出显式结果。
只能用数值积分由程序来得到。
目前常用的是高斯积分（矩形族）和哈默尔积分
（三角形族）。它们的积分点位臵、加权系数等见表
3-1,3-2（ P.74~76)。
其积分公式见式（ 3,6-40）、（ 3,6-41）。
4.5 平面等参数单元
4.5.8 作等参元分析时应注意的问题
等参元分析中要用 det[J]-1,可见雅可比行列式等于
零将导致刚度矩阵等无法积分，使分析失效。因此要
避免以下可能使 det[J]=0的情况：
1) 子单元边界不能过于扭曲。
2) 矩形子单元不能退化成三角形。
3) 子单元角顶处单元边线切线角角不能等于 1800。
上述情况如 P.81 图 3-39 示意。
此外,子单元边界上结点应尽可能是或接近等分点，
避免产生奇异单元。
可能情况下应采用直边子单元，这样可使雅可比
矩阵简单，提高计算效率。
4.5 平面等参数单元
4.5.9 离散化时应注意的问题
除对等参元应注意上述问题外，任何有限元分析
都还应注意以下几点：
1) 相互邻接的单元大小应尽可能均匀。
2) 单元最大尺寸与最小尺寸之比应尽可能接近一，最
多不应大于二。
3) 应合理编码，使单元结点间的整体编号差值最小。
4) 应尽可能使各界点的单元数目相同,如 P.82 图 3,6-42
左图示意,
4.6 平面问题程序 (二 )
应用本程序时,数据文件的准备见 P.217。
运行程序查看计算结果
4.7 Wilson非协调单元
至今为止所介绍的单元都是能保证收敛的协调单
元。单元计算结果的精度，取决于位移模式中坐标完
全多项式的次数。为改进精度 Wilson提出非协调的单
元，简单介绍如下。
在矩形双线性单元基础上，增加一位移修正项：
? ? ? ?? ? ? ?? ? ee aNdNd ??
? ? ?
?
?
?
?
?
?
21
21
00
00
NN
NNN
? ? ? ?vvuue aaaaa 2121?
21 1 ???N
22 1 ???N
显然在四个结点处修正项等于零，因此它只影响单元
内部位移，可见 [a]e是单元内部位移参数。所求得的
结点位移将仍直接是实际问题的近似值。
4.7 Wilson非协调单元
根据这一位移模式，可得应变矩阵为
? ? ? ?? ? ? ?? ? ee aNdNd ???
经虚位移原理做单元列式，可得
? ? ? ?? ? ? ?? ? ee aBdB ??? ?
? ? ? ?
? ? ? ?
? ?
? ?
? ? ? ?
? ? ??
?
?
?
??
?
?
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
a
e
ee
e
e
a
e
P
PF
a
d
kk
kk
T
1
1
双线性
修正项
相关项
相关项
结点力
双线性
等效
修正项
等效
各项的具体计算公式见下页。
4.7 Wilson非协调单元
? ? ? ? ? ?? ???
Ve DDBk dV
T
双线性
修正项
等效
? ? ? ? ? ?? ???
V BDBk dV
TT
1
相关项
? ? ? ? ? ?? ???
Va
BDBk dVT
修正项
? ? ? ? ? ? ? ? ? ??? ??
?
?
SV b
a
e NFNP dSdV
TT
? ? ? ? ? ? ? ? ? ??? ??
?
?
SV be
NFNP dSdV TT
双线性
等效
4.7 Wilson非协调单元
由于 [a]e是单元内部位移参数，与单元集合体（整
体）的其他单元无关，因此在集成之前可做如下消去
? ? ? ?
? ? ? ?
? ?
? ?
? ? ? ?
? ? ??
?
?
?
??
?
?
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
a
e
eee
a
e
P
PF
a
d
kk
kk
T
1
1
? ? ? ? ? ? ? ? ? ? )( T11 eaeae dkPka ?? ?
? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ?aeaee
eae
PkkPF
dkkkk
1
1
T
1
1
1 )(
?
?
??
??
单元刚度
单元等效荷载
4.7 Wilson非协调单元
由此可见做了上述消去（也称静力凝聚）后，单元刚
度方程形式上和双线性单元完全一样。
问题：试增加 PSTE 程序的单元库，使其能用非
协调元进行分析。
能否在四边形等参元基础上来作非协调元？
协调元是肯定收敛的，非协调元收敛性不能保证。
下一章将介绍如何检查非协调元的收敛性。这里只给
出结论：
1) 当单元网格为矩形或平行四边形时，分析将试收敛
的。
2) 任意四边形时，以 ?=?=0点的 [J]作为单元的各点
的 [J]，这样计算的修正单元刚度也能使分析收敛。
