二、_弹性力学有限元法基本原理(一)

❖ 解决办法：将上述多项式系数广义坐标代换为单元节点位移广 义坐标（插值法）。
❖ 将三个节点坐标分别代入上述位移多项式：
u 1 2 x 3 y
ui 1 2 xi 3 yi u j 1 2 x j 3 y j
um 1 2 xm 3 ym
❖ 解上述方程，用节点位移表达多项式系数：
3）推论： 对三节点单元，形函数在单元内部和边界上均为线性函数， 形函数在单元一条边上的值，只跟边端点位置有关，与第三 点位置无关；某节点的形函数在对边上 恒为0。
4）形函数的几何意义
图2-5 形函数的几何意义
根据三角形单元形函数的上述性质也可以推断，单元的 假定位移在内部和边界上都线性分布，边界上的位移只 跟两端节点位移有关，保证了整个求解区域上位移的连 续性。
上一节以受轴向力的弹性杆为例讨论了有限元位移法 的基本原理和步骤，揭示了有限元法的本质特征。
本节讨论将该方法推广到解决弹性力学平面问题。弹 性区域离散化采用3节点三角形单元（T3单元）。
1、结构离散化
连续求解区域划分为有限个三角形区域（单元）； 或者连续体离散为有限个小单元的结合体，单元之间在节点处连接。 问题的未知量转化为节点位移。
下面用节点位移未知量作为待定参数（广义坐标），得到 其标准有限元形式。
重新构造单元内位移试探函数
❖离散结构中，节点位移分量是问 题的基本未知量。
❖在每个单元内通过对节点位移插 值，分片建立位移试探函数：
❖
对单元1有：u
L
L
s
u1
s L u2
矩阵形式为： u N d
其中：
N
L
L
s
s L
cL2 6
d T
4 5
和
cL2 6
d T
7 8
❖ 系统总势能等于三个单元应变能之和减去三个单元外力功。
❖ 为了适应矩阵形式运算，将单元势能矩阵表达式中的 d用整体节点 位移向量 D u1 u2 u3 u4 T 代替，同时单元刚度矩阵扩展成总
体规模（4×4）。则各单元相加后系统总势能为：
p
1 2
DT
3L 0
E
2
2 x
Adx
3L
qudx
0
❖单元应变表达式：
由应变~位移关系： x u,x u,s
得到单元上应变表达式： x
d ds
N d
Bd
上式就是用节点位移表达单元应变的公式。
其中
B
1 L
1 L
称为单元应变矩阵
❖计算一个单元内应变能(单元1）：
x
d ds
N d
Bd
U
L 0
1 2
▪ 离散化
考虑图2-1（a）所示的一维直杆结构，杆的总长度为3L。载 荷与图1-4相同。杆分为三个子区域，称为单元。单元之间的 连接点称为节点。这一步骤称为离散化。节点位移是问题的 基本未知量。
图2-1
（a）一维直杆的分域（截面积A，弹性模量E，轴向受力） （b）杆的有限单元
▪ 分片假设单元上位移试探函数
载荷与上节例题相同： q cx
载荷在三个单元内局部坐标下的表达式分别为：
q cs ， q c(L s) ， q c(2L s)
❖按外力功的积分式，分别计算三个单元的外力功 第1单元外力功为：
L
quds
0
L uT qds
0
L d T
0
N T
csds
cL2 6
d T
1 2
❖同理，第2、3单元外力功分别为：
Ni 称为单元的插值基函数或形函数。每个节点对应一个形函数。
❖ 位移模式的矩阵形式为：
u Niui N ju j Nmum v Nivi N jv j Nmvm
ui
vi
u
u
v
Ni
0
0 Ni
Nj 0
0 Nj
Nm 0
0 Nm
u v
j j
um
vm
Ni N j Nm ae Байду номын сангаасae
单元应力根据平面问题的物理方程得到：
x y
D
DBae
Sae
xy
其中：
S DB D Bi
Bj
B m
Si
Sj
Sm
S 称为应力矩阵。将平面应力或平面应变问题的弹性矩阵代
入，就可以具体计算出应力矩阵。
至此，已完成在二维弹性体区域上构造位移试探函数，并做好 计算系统总势能的准备。
5、利用最小势能原理建立有限元求解方程
4、用节点位移表达单元应变和应力
单元位移模式确定后，很容易得到用节点位移表示单元应 变和应力的表达式。
应用弹性力学平面问题几何方程的矩阵形式，得到：
0
x y
xy
x
0
y
u v
Lu
LNae
L
Ni
Nj
Nm ae
y x
Bi B j Bm ae Bae
矩阵 B 称为单元应变矩阵，其子块为：
c jv j
cmvm )
❖
将求得的广义坐标 1
~
代回前面位移多项式，整理得到以
6
单元节点位移分量表示的位移模式：
u Niui N ju j Nmum v Nivi N jv j Nmvm
u 1 2 x 3 y v 4 5x 6 y
其中
Ni
1 2A
(ai
bi x
ci y)
（i，j，m）
AE L
1
1
0
0
1 1 0 0
0 0 0 0
0 0 0 0
0 AE 0
L 0 0
0 1 1 0
0 1 1 0
0 0 0 0
0 AE 0
L 0 0
0 0 0 0
0 0 1 1
0
0
1
1
D
DT
cL2 6
1 2 0
cL2 6
0 4 5
cL2 6
00 7
0 0 8
U 1 dT kd
是插值基函数矩阵，称为“形函数”矩阵。
d
uu12
称为单元1的节点位移列阵。
❖其它两个单元也有同样的插值位移试探函数：
单元2： u N d
, d
u2 u3
单元3：
u Nd
,
d
uu43
N 同上
❖每个单元中，位移试探函数是位置坐标的简单函数（线性），任意 一点函数值取决于单元节点位移，且在单元节点上试探函数值等于待 定节点位移。
Se NT TtdS PSe ——单元面力等效节点力列阵；
E
x
2
Ads
1 2
L T
0x
AE xds
1 2
d T
L 0
BT
AE Bdsd
或 U 1 dT kd
2
其中：
k
L
0
T
B
AEBds
AE L
1 1
1
1
——称为杆单元刚度矩阵。
❖由于单元的几何、物理参数相同，上述单元应变能表达式对于3 个单元相同，只是节点位移列阵的分量不同。
❖计算单元上外力的功：
❖单元上的位移试探函数又称为单元位移模式。
显然，在使用插值试探函数情况下，整个杆上，由各单元位移函
数拼接而成的试探位移场是连续的。只要我们记住 u1 0 ，得
到的就是可能位移场。这样的位移试探场以未知节点位移作为待 定参量（广义坐标）。
下面进一步实施里兹法求解。
计算离散系统的势能泛函
p
把上述假设的分段线性位移场代入势能泛函：
p
3L 0
E 2
u,2x
Adx
3L
qudx
0
分段计算上式积分，再应用驻值条件，可立即求出待定系 数，位移场就完全确定，进而可求出各单元应力。
上述是原理上的里兹法有限元形式求解过程，其关键 是采用分片多项式拟合全域上的可能位移场——分片 试探函数。
a e
ai aj
ui
vi
uj
vj
um
vm T
am
——称为单元节点位移列阵。
图2-4 三节点三角形单元
在单元的三角形区域构造简单位移试探函数。 ❖ 采用x，y的一次多项式：
u 1 2 x 3 y v 4 5x 6 y
1 ~ 6 是待定系数，称为广义坐标。
❖ 显然，用这种形式的单元位移模式构造整个求解区域上的许 可位移场将非常困难。
定义上面方程的系数行列式为D：
1 xi D 1 xj
1 xm
yi yj 2A ym
（A为三角形单元的面积。）
用公式法（克莱姆法则）解出广义坐标 1 ~ 3 得：
ui 1 2 xi 3 yi
u j 1 2 x j 3 y j
um 1 2 xm 3 ym
1 xi D 1 xj
1 xm
图2-3 二维区域离散化
2、单元位移模式及插值函数
从离散结构中取出一个典型的分片区域（单元），研究二维区域 上分片多项式假设位移场（单元位移模式）。
单元如图2-4所示。单元节点采用局部编号i，j，m（逆时针旋 转），每个节点两个自由度（位移分量）：
ai
uvii
（i,j,m）
❖ 则每个单元六个节点位移分量：
yi yj 2A ym
1
1 D
ui uj
xi xj
yi yj
1 2A
(
ai
ui
a juj
amum )
um xm ym
1
2
11 D
ui uj
1 um
yi yj
1 2A
(bi
ui
bju j
bmum )
ym
1
3
11 D
xi xj
ui
1
u j 2 A (ciui c ju j cmum )
1 xm um
p
e
a eT
1 BTDBtdxdyae
e 2
e
a eT
NTftdxdy
e
aeT NTTtdS
e
se
❖ 在上式中定义：
p
1 T Dtdxdy
2
uTftdxdy
uT Ttds
S
BT DBtdxdy Ke ——单元刚度矩阵；
e
NT
e
ftdx
dy
Pbe
——单元体力等效节点力列阵；
第二单元 弹性力学有限元法基本原理（一）
第一节 里兹法的有限元形式
由于需要在整个求解区域上假设试探函数，经典里兹法在 解决实际问题时，尤其是几何形状复杂的二、三维问题， 具有局限性。
解决上述问题的办法是在求解区域上分片假设试探函数。 下面以一维直杆的分析为例子，研究基于里兹法的有限元
位移法基本原理和求解过程。
uu23 u4
cL3 3AE
13 23 27
❖结合单元位移模式 u N d 就得到整体上近似位移场。
❖ 单元应力由公式 x E x EBd 得到。
❖ 本问题有限元解与精确解的比较如图2-2所示。 图2-2 受轴向力杆的精确结果和有限元结果
有限元解与经典里兹解对比
第二节 常应变三角形单元解平面问题
首先按如下形式假定分片位移场
u b1 b2x 0 x L
u b3 b4 x L x 2L
u b5 b6 x 2L x 3L
bi 是待定常数。
为了使得上述假设位移场是“许可位移”，上述多项式待 定系数必须满足一定约束关系（3个），显然，该问题的 独立参量（广义坐标）只有3个。
上式中， 系数 ai ,bi , ci (i, j, m) 分别是行列式D第i,j,m行第1，2，3列
元素的代数余子式，其值取决于节点坐标。
同理，对y方向位移函数作代换，可求得广义坐标 4 ~ 6 ：
4
1 2A
(ai
vi
a jv j
amvm )
5
1 2A
(bi
vi
bjv j
bmvm )
6
1 2 A (civi
这个做法正体现了有限元法的实质。
上面形式的分片位移试探函数有下列缺点： 1) 必须对它进行调整，使其满足连续条件和边界约束条件； 2) 多项式系数作为广义坐标缺乏明显的物理意义。
因此，上述不是通常意义上标准的有限元形式，仍然具有 局限性，如对于二维以上的问题使各单元之间分片多项式 保持连续性很难处理。
x
Bi
LN i
0
y
0
y
Ni 0
x
N
i
0 Ni
x
0
N
i
y
0
Ni
y
N
i
x
（i, j, m）
❖ 将形函数分别代入上式，最后求得应变矩阵如下：
B
1 2A
b0i ci
0 ci bi
bj 0 cj
0 cj bj
bm 0 cm
0
c
m
bm
该应变矩阵的元素取决于节点坐标，是常数矩阵。因此，单元 内应变、应力是常数。该单元称为常应变单元。
弹性力学平面问题的总势能泛函表达式如下：
p
1 T Dtdxdy
2
uTftdxdy
uT Ttds
S
其中：
t ——弹性体厚度；
f ——平面问题体积力向量；
T ——边界面力向量。
应用前面在整个求解区域上分片假设的位移场（位 移模式），在有限元离散模型上对上述总势能泛函 进行分片计算（对各单元区域积分）并求和：
式中：
N Ni N j Nm 称为单元的形函数矩阵
3、形函数的性质和几何意义
1）相应于某节点 i 的形函数在i节点上值为1，在j，m节点上值为0。
Ni ( x j ,
y
j
)
ij
1 0
当j i 当j i
(i, j,m)
2）单元中任一点各形函数之和等于1，即：
Ni N j Nm 1
提示：三个形函数只有两个独立。该性质保证各节点函数值相同时， 插值得到单元上任意一点函数值亦相同。
2
简写为：
p
1 2
DT
KD DT R
应用势能驻值条件：
p
D
0
得到有限元求解方程——系统平衡方程：
KD R
1 1 0 0 u1
1
即：
AE 1
L 0
0
2 1 0
1 2 1
0
1
1
uu32 u4
cL2 6
6 12 8
引入约束条件：u1 0
即划去第一个方程，解出其余三个方程得到：