有限元讲义弹性力学平面问题有限单元法(四结点四边形等参元,八结点曲线四边形等参元,问题补充)分析

2.6 四结点四边形单元
(The four-node quadrilateral element)
前面介绍了四结点的矩形单元
其位移函数:
xy y x U 4321αααα+++=
xy y x V
8765αααα+++=
为双线性函数，应力，应变在单元内呈线性变化，比常应力三角形单元精度高。

但它对边界要求严格。

本节介绍的四结点四边形等参
元，它不但具有较高的精度，而且其网格划分也不受边界的影响。

对任意四边形单元（图见下面）若仍直接采用前面矩形单元的位移函数，在边界上它便不再是线性的（因边界不与x,y 轴一致），这样会使得相邻两单元在公共边界上的位移可能会出现不连续现象（非协调元），而使收敛性受到影响。

可以验证，利用坐标变换就能解决这个问题，即可以通过坐标变换将整体坐标中的四边形（图a ）变换成在局部坐标系中与四边形方向无关的边长为2的正方形。

正方形四个结点i,j,m,p 按反时钟顺序对应四边形的四个结点i j m p 。

正方形的 1-=η 和 1=η 二条边界，分别对应四边形的i ，j 边界和p,m 边界；ξ=-1和ξ=+1分别对应四边形的i ，p 边界和j ，m 边界。

如果用二组直线等分四边形的四个边界线段，使四边形绘成一个非正交网格，那么该非正交网格在正方形上对应着一个等距离的规则网格（见图a, b ）。

当然, 局部坐标上的A 点与整体坐标的A 点对应。

一、四结点四边形等参单元的形函数及坐标变换
由于可以将整体坐标下的四边形单元变换成局部坐标下的正方形单元，对于这种正方形单元，自然仍取形函数为： ξηαηαξαα2321+++=U ξηαηαξαα8765+++=V
引入边界条件，即可得位移函数：
∑=ijmp
i i U N U
i ijmp
i V N V ∑==
写成矩阵形式：
{}
{}[]{}e
e p i p i e
d N d N N N N V U f =⎥⎦
⎤⎢⎣⎡=⎭⎬⎫⎩⎨⎧=0
00 式中形函数: ()()()ηηξξηξi i i N ++=
114
1
, ()p m j i ,,, 按照等参元的定义,我们将坐标变换式亦取为: p p m m j j i i i ijmp
i x N x N x N x N x N x +++==∑
p p m m j j i i i ijmp
i y N y N y N y N y N y +++==
∑ ()162-- 式中形函数N 与位移函数中的完全一致。

可以验证，利用坐标变换式(2-6-1)，可以把整体坐标系中的任意四边形单元（图a ）变换成在局部坐标系中与四边形对应的边长为２的正方形。

因此可以将上述位移函数和形函数用于任意四边形单元，并将形函数中的ξ,η理解为任意四边形单元的局部坐标。

这样由位移函数可以得到单元各点的位移。

在四条边界上分别有ξ＝±和η＝±１，故边界上的位移呈线性变化，位移的连续性可得到保证。

于是，我们可以理解为：任意四边形单元是从基本的正方形单元变换过来的实际单元。

因此又称正方形单元为母体单元，或基本单元。

例题：
为了加深理解，现考察实际单元为矩形单元的坐标变换，在2.4节中，我们定义局部坐标与整体坐标的关系是：
()01
x x a
-=
ξ ()01y y b -=η
式中(x 0 , y 0 )为局部坐标原点。

由上第一式()01
x x a
-=
ξ得： ()()j i i j x x x x x a x ++-=+=2
1
210ξξ
将其重新组合：
()()j i x x x ξξ++-=12
1
121
()()()()()()()()p m j i x x x x ηξηξηξηξ+-++++-++--=
114
1114111411141
对照2.4中的形函数表达式,便知:
p p m m j j i i x N x N x N x N x +++= 自然同理可得: p p m m j j i i y N y N y N y N y +++=
由此知,矩形单元可以看作是四结点四边形单元的特例，自然，它也是等参元。

《有限元法概论》（第二版）P 172 中，是这样解释等参元的基本概念和推导方法的： 图形变换
四结点正方形(母元) 图形变换 四结点四边形(等参元) (ξ-η平面内) ───→ (x,y 平面内)
进行图形变换的关键是进行图形结点坐标之间的变换：
正方形结点坐标 坐标变换 四边形结点坐标 (ξi ,ηi ) ────→ (x,y) i=i,j,m,p i=i,j,m,p
为了实现上述结点坐标之间的变换，可利用母元的形函数，得出(ξ,η)和(x,y)之间的坐标变换式。

图形变换具有如下性质:
1. 母元中的坐标线对应于等参元的直线；
2. 四结点正方形母元对应于四个结点可以任意布置的直边四边形等参元；
3. 变换式(2-6-1)能保证相邻等参元的边界位移彼此协调。

二、几何矩阵[B]
已知单元的应变与结点位移之间的关系是:
{}{}d N N N N x y y x
p i p
i ⎥⎦
⎤
⎢⎣⎡⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂=
0000 ε ()262-- 形函数矩阵[N]只是局部坐标ξ,η的显函数,为求形函数对整体坐标x,y 的偏导数，必须
用复合函数求导公式:
ξ
ξξ∂∂∂∂+
∂∂∂∂=∂∂y
y N x x N N i i i ηηη∂∂∂∂+
∂∂∂∂=∂∂y
y N x x N N i i i ()362-- 或写成： []⎪⎪
⎭⎪
⎪⎬⎫
⎪⎪⎩⎪⎪⎨⎧∂∂∂∂=⎪⎪⎭
⎪⎪⎬
⎫⎪⎪⎩⎪⎪⎨⎧∂∂∂∂y N x N J N N i i i i ηξ ()a 362--
式中: []⎥⎥⎥
⎦
⎤
⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂=ηηξξx y
x y
x J ()b 362-- 称为雅可比矩阵，而把它的行列式称为雅可比行列式。

把式(2-6-1)代入[J]得:
[]⎪⎪⎭
⎪⎪⎬⎫
⎪⎪⎩⎪⎪⎨⎧⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣
⎡∂∂∂∂∂∂∂∂=∑∑∑∑p p
m m
j j
i i
p m j i
p m j i
ijmp i i ijmp i
i ijmp i i ijmp i
i
y x y x y x y x N N N N N N N N y N x
N
y N x N J ηη
η
η
ξξξξηηξξ 将形函数 ()()ηηξξi i i N ++=1141 ()p m j i ,,, 代入，分别对ξ,η求偏导，即可得到四结点四边形等参元的雅可比矩阵：
[]⎥⎦
⎤
⎢⎣⎡++++=ξβξ
βηβη
βB A B A J 432141 ()562--
式中常数记为: ∑=
ijmp i i i x A ηξ ∑=ijmp i i i y B ηξ ∑=
ijmp i i x ξβ1 ∑=ijmp
i i y ξβ2
∑=
ijmp i i x ηβ3 ∑=ijmp
i i y ηβ4 该雅可比矩阵的逆:
[]()()⎥⎦
⎤⎢⎣⎡++-+-+=
-ηβξβηβξ
βA A B B J
J 1324141
雅可比行列式:
()()()()[]ξβηββηβηβA A A J ++-++=
324116
1
()()()[]ηββξββββββ3421324116
1
B A A B -+-+-=
可以证明，如果四结点四边形的四个内角都小于180°的话，雅可比行列式|J|大于零，
其逆阵[J]-1
是存在的。

换句话说，为了使上述等参元能保持较好的精度，整体坐标系下所划分的任意四边形单元必须是凸四边形，即任意内角都不能大于180°。

四边形也不能太歪斜，否则会影响其精度。

利用雅可比的逆矩阵，即可求出整体坐标系下形函数的偏导数：
[][]()()),,,(11 1141p m j i i i i i i i i i i J N N J y N x N ⎭⎬⎫⎩⎨⎧++=⎪⎪⎭
⎪⎪⎬⎫⎪⎪⎩⎪
⎪⎨⎧∂∂∂∂=⎪⎪⎭⎪⎪
⎬⎫
⎪⎪⎩⎪⎪⎨⎧∂∂∂∂--ξξηηηξηξ ()662-- 求出全部偏导，即代回（2-6-2）右侧，即可得到几何矩阵[B], [B]是ηξ,的函数,即:
[][]⎥⎥
⎥⎥⎥⎥
⎥⎦
⎤⎢⎢⎢⎢
⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡∂∂∂∂∂∂∂∂
=x N y N x N y N x N y N y N y N y N x N x N x N N x y
y x B p p j j i
i p j i p
j i 0
000000
0 将(2-6-6)代入即可获得[B], [B]是ηξ,的函数。

三、单元刚度矩阵
获得[B]后,便可由单刚的一般表达式:
[][][][]⎰⎰=dxdy B D B t K T
求出四结点四边形的单元刚度矩阵。

在按上述公式作积分运算时,必须把面积元dxdy 变换成ηξd d ,图a 上的面积元abdc 的面积等于矢量→
ab 与矢量→
ac 的矢量积的模，
即微元
→
→
⨯=ac ab dA
沿ξ轴对应于d ξ的矢量增量是: →
→→
∂∂+∂∂=j d y i d x ab ξξ
ξξ 沿η轴对应于η的矢量增量是: →→→
∂∂+∂∂=j d y i d x ac ηη
ηη 式中
→
→
j
i , 是坐标x,y 的单位矢,注意到:
0=⨯=⨯→
→
→
→j j i i 1=⨯→
→j i 则有:
⎪⎭⎫ ⎝⎛∂∂+∂∂⨯⎪⎭⎫ ⎝⎛∂∂+∂∂=→→→→j d y
i d x j d y i d x dA ηηηηξξξξ
ηξηξηξηξd d J j i d d x y y x =⨯⎪⎪⎭
⎫
⎝⎛∂∂∂∂-∂∂∂∂=→→
因此刚度矩阵的积分式:
[][][]⎰⎰--=111
1 ηξd td J D B K T ()762--
在计算单元刚度矩阵[K]中元素时，由于被积函数中出现了雅可比行列式，使得它用解
析法很难求其积分，故常采用高斯数值积分法.
四、数值积分
1．一维数值积分
()⎰b
a d F ξξ
基本思想:构造一个多项目式()ξψ，使在i ξ(i=1,2……n)上有 ()()i i F ξξψ= ， 然后用近似函数()ξψ的积分
()⎰b
a
d ξξψ来近似原被积函数()ξF 的积分()⎰b
a d F ξξ。

i
ξ称为积分点或取样点，积分点i ξ的数值和位置决定了()ξψ近似()ξF 的程度，亦即决定数值积分的精度。

对于n 个积分点,按照积分点位置的不同选择，通常采用两种不同的数值积分方法，Newton-Cotes 积分和高斯积分方案。

二者方法基本相同，只是前者的积分点ξi 是等间距分布,而后者不是等间距分布。

高斯积分的积分点位置由下述方法确定：
① 定义n 次多项式()()()()n P ξξξξξξξ---= 21 ② 由下列条件确定n 个积分点位置
()⎰
=b
a
i d P 0ξξξ 11,0-=n i
由上二式可见，()ξP
有以下性质：
① 在积分点上()0=i P ξ； ② 多项式()ξP
与1210
... ,,,-n ξξξξ
在（a, b ）域内正交。

由此可见n 个积分点的位置
ξi 是在求积域（a, b ）内与
1210 ... ,,,-n ξξξξ正交的n 次多项式()ξP 构成方程
()⎰=b
a i
d P 0ξξξ的解。

于是,被积函数()ξF 可由2n-1次多项式()ξψ来近似。

()()()()∑∑=-=-+=n
i n i i i n i
P F l
11
1
ξξβξξξψ
()ξi l 是n-1阶Lagrange 插值函数 ()()()()()()()n n i i i l ξξξξξξξξξξξξξ------= 121
用
()⎰b a
d ξξψ近似()⎰b
a d F ξξ, 并考虑上面确定n 个积分点位置的约定条件得:
()()⎰∑=+=b
a n
i i i R F H d F 1
ξξξ 式中 ()⎰-=
b a n i
i d l H ξξ1
H i 称为积分的权系数(加权系数, 权)，可见加权系数H i 与被积函数()ξF 无关，只与积分点个数和位置有关。

为便于计算积分点位置ξi 和加权系数H i , 常将上式中的积分限规格化，即令: 1-=a 1=b
由此计算出的i ξ , H i 对应于原积分域(a, b)的关系为:
i b a b a ξ22--+ i H a
b 2
- 对于多重积分，按重积分规则，计算内层时，保持外层为常数，逐层计算即可。

有关数
值积分更详尽的资料可参阅其教科书。

2．等参元计算中数值积分阶次的选择
数值积分中，如何选择积分阶次将直接影响计算精度、工作量，选择不当则有可能使计算失败。

选择积分阶次的原则首先是要保证积分的精度能满足所求问题的要求。

如一维问题：
设：插值函数中的多项式阶数为p ， 微分算子中导数的最高阶为m ，则有限元得到的近似能量是2(p-m)次多项式。

若被积函数是2(p-m)次多项式，应选高斯积分的阶次为n=p-m+1。

三结点三角形单元（线性单元）刚度矩阵中被积函数是常数，故只需一个积分点。

双线性单元 p-m+1=1, 但插值函数将包含有二次项（ξ2 ,ξη,η2
）中的ξη项，所以要达到精确积分应采用2×2阶的高斯积分。

有的有限元书中列出了不同插值函数的高斯积分点数n 及相应积分点坐标i ξ和权系数H i 的值，可供编写有限元程序时参考。

五、四边形等参元的荷载等效变换
四边形等参元等效结点荷载的计算，仍然利用局部坐标体系。

1．集中力
设在单元上任意一点C 作用有一集中力{}⎭
⎬⎫
⎩⎨⎧=y x P P P ,根据荷载等效变换的一般公式，其等效结点荷载计算公式是： {}[]{}P N R T
C =
式中，
[]T
C N 是形函数[N]在集中力作用点C 上的值。

2．体积力
设在单元上作用有单位体积力
{}⎭
⎬⎫
⎩⎨⎧=y x W W W ，其等效结点荷载的计算公式是：
1
1
11
{}[]{}T R N W t J d --=⎰⎰
ξdη
3．分布力
设单元的某一边界上承受的分布表面力是
{}⎭
⎬⎫
⎩⎨⎧=y x q q q ，其等效结点荷载的计算公式是，
当分布力作用在单元ij, mp 边界上时：
ξξ
ξd y x t q N R T
⎰-∂∂+∂∂=1
12
2)()(}{][}{ 当分布力作用在单元ij, mp 边界上时：
ηη
ηd y
x t q N R T ⎰-∂∂+∂∂=1
1
22)()(
}{][}{
2.7 八结点曲线四边形等参元
在常规有限元程序的单元库内，目前工程上对平面问题最适用的是八结点等参元。

下面作一简要介绍。

单元及结点编号如图a所示，位移和力列向量仍采用与前面类似的排列方式。

为了构造其位移函数，按照前面等参元的概念，可先构造八结点正方形母体单元的位移函数（图b），并取为双二次插值位移函数：
2
8
2
7
2
6
5
2
4
3
2
1
ξη
α
η
ξ
α
η
α
ξη
α
ξ
α
η
α
ξ
α
α+
+
+
+
+
+
+
=
U
2
16
2
15
2
14
13
2
12
11
10
9
ξη
α
η
ξ
α
η
α
ξη
α
ξ
α
η
α
ξ
α
α+
+
+
+
+
+
+
=
V
代入边界条件，得到用结点位移{d}表示的位移场:
⎪
⎪
⎭
⎪⎪
⎬
⎫
=
=
∑
∑
=
=
8
1
8
1
i
i
i
i
i
i
V
N
V
U
N
U
()1
7
2-
-
式中的形函数:
()()()
()()⎪
⎭
⎪
⎬
⎫
=
+
+
-
-
=
=
-
+
+
+
=
8,7,6,5
1
1
2
1
4,3,2,1
1
1
1
4
1
2
2
2
2j
N
i
N
i
i
j
j
j
i
i
i
i
i
η
η
ξ
ξ
η
ξ
ξ
η
η
η
ξ
ξ
η
η
ξ
ξ
()2
7
2-
-
按照等参元的定义，实际单元映射成正方的母体单元的坐标变换式与式2-7-2等同,得:
⎪
⎪
⎭
⎪⎪
⎬
⎫
=
=
∑
∑
=
=
8
1
8
1
i
i
i
i
i
i
y
N
y
x
N
x
()3
7
2-
-
前面我们在推导单元刚度矩阵时，都是假定单元结点编号从左下角开始，逆时针编号。

从有限元理论上说，其编号应该是可以任意的。

下面给出一种不同结点编号的例子。

实际单元基本单元
2.8 几个问题的补充 一、变厚问题 各单元取不同t 。

二、不同材料问题
各单元取不同E,μ
三、平面应力与平面应变问题
在前面三角形单元的推导中，我们假定其为平面应力问题。

工程中，还有另一类情形─平面应变状态。

例如，在对坝体或遂洞等长柱体进行分析时，如果取xoy 坐标平面与其横截平面行,而
Z 轴与其长度方向一致（如图）
那么，由于所考察物体在Z 方向的尺寸很大，且又受
到平行于xoy 平面，且不沿长度方向变化的荷载作用，就
可认为各个横截面应处于同样的状况，即近似认为Z 方向
的位移分量W=0, (位移与Z 无关)
于是，由弹性力学知，在六个应力分量中也仅有三
个独立分量σx ,σy 和τxy, 而()y x z
σσμσ+= 不
独立。

并可得到平面应变问题的物理方程。

⎪⎪⎭⎫ ⎝⎛---=y x x E σμμσμε112 ， ⎪⎪⎭
⎫ ⎝⎛---=x y y E σμμσμε112 21112μμμτγ-⎪⎪⎭⎫ ⎝⎛-+=E
xy xy
比较平面应力问题：
()y x x E μσσε-=1， ()x y y E
μσσε-=1 ， )1(2μτγ-=E xy xy 得知只需把应力问题中的E 换成μ-1E ，μ换成μμ
-1即得应变问题。

所以在这类问题的程序设计中，通常可以同时求解应力和应变问题，只需设置一开关变量便可以实现。

四、各向异性材料
在弹性矩阵[D]中反应。

如正交各向异性时的弹性矩阵[D]为：
[]
⎥
⎥
⎥
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎢
⎢
⎢
⎣
⎡
-
-
-
-
=
xy
y
x
y
y
x
y
x
y
x
y
x
y
x
x
G
E
E
E
E
D
1
1
1
1
μ
μ
μ
μ
μ
μ
μ
μ
μ
μ
（见凯维奇《有限元法》，中P102～104，英 P98～109 ）
五、设置不同类型的单元
在工程中，同一构件经常要用到一种以上材料。

如利用角、槽钢等在开洞板内作加强筋用。

另一种情况是钢筋混凝土构件的全过程分析中，纵向受拉筋的单元划分，如图。

混凝土被划成若干平面三角形单元，而将纵向受力钢筋（或箍筋）当作线单元（图中红线所示），也可把钢筋等效成与混凝土叠合的三角形单元，但此时均为两种不同材料。

有限元分析时，可认为结构是由若干（面）单元（三角形或矩形）和线（杆）单元（二力杆）共同组成，划分网格时，若碰到线单元，便将结点取在线单元上，使面元和杆元使用共同的结点。

引进杆元后，并不增加结构的自由度（未知数），只是
装配总刚和计算过程中多了杆元单刚。

如本来是三个三角形
的结点，可能还有两个杆元汇交如同一个结点（见上图）。

杆元（见右图）单刚的一般形式可取为：
[]
⎥
⎥
⎥
⎥
⎦
⎤
⎢
⎢
⎢
⎢
⎣
⎡
-
-
-
-
=
α
α
α
α
α
α
α
α
α
α
α
α
α
α
2
2
2
2
2
2
sin
sin
cos
cos
sin
sin
cos
sin
sin
cos
cos
sin
cos
cos
symmetric
l
EA
K e
l
应力矩阵
[][]α
α
α
αsin
cos
sin
cos-
-
=l E
S
六、温度应力问题
通常是将各单元由于温度改变所产生的应力：
[][]{}
⎰⎰tdxdy
D
B T
ε当作外力，化成等效结点荷载加到右端项中求解。

式中：{}[]T
T
T0 0
α
α
ε=。