有限元 位移约束条件的引入

k31u1 k32v1 k331015u2 k34v2 L L k3,2nvn 3k331015
u2 3
其他方程不变 为此我们就建立了新的方程
K F
§4-6 有限元分析的实施步骤 根据前面的讨论，现以三角形常应变单元为例来说明应
用有限元法求解弹性力学平面问题的具体步骤。
①力学模型的确定根据工程实际情况确定问题的力学模型， 并按一定比例绘制结构图、注明尺寸、载荷和约束情况等。
节点的多少及其分布的疏密程度（即单元的大小），一般 要根据所要求的计算精度等方面来综合考虑。从计算结果的精 度上讲，当然是单元越小越好，但计算所需要的时间也要大大 增加。另外，在微机上进行有限元分析时，还要考虑计算机的 容量。因此，在保证计算精度的前提下，应力求采用较少的单 元。为了减少单元，在划分单元时，对于应力变化梯度较大的 部位单元可小一些，而在应力变化比较平缓的区域可以划分得 粗一些。
M M M M 0 k2n1,2 0 k2n,2
0L L L L L 0
0L L L L L 1L L L L L
k2,2n1 0
0 L L L L L k4,2n1 MM M M M M M
MM M M M M M
MM M M M M M
MM M M M M M
0 L L L L L k2n1,2n1 0 L L L L L k2n,2n1
第四章 平面问题的有限元分析
§4-4 等效节点力的计算
计算等效节点力
用虚功原理确定等效节点力
若三角形单元上作用有集中力g、分布力q(力/面积)和 体积力p(力/体积)，则根据静力等效原理，节点力所做的虚 功等于三种力所做的虚功。
m
* eT Re f * T g f * T qtds f * T ptdxdy
1k11 1015 Fy1
k11
1015 k21`
3k33
1015
Fy2
k31
k41
Fx3 Fy3
M M
M
M
M M
Fxn
k2n1,1
Fyn
2n1 k2n,1
k12 k22 k32 k42 M M M M k2n1,2 k2n,2
k13 k23 k33 1015 k43 M M M M k2n1,3 k2n,3
还有一点值得注意的是，单 元各边的长度不要相差太大，以 免出现过大的计算误差或出现病 态矩阵。例如，图4-12所示的(a) 、(b)两种单元划分，虽然都是同 样的四个节点，但(a)的划分方式 显然要比(b)的方式好。
(a)
(b)
图4-12
三. 节点的编号
节点编号时，应该注意要尽量使同一单元的相邻节点的号 码差尽可能地小，以便最大限度地缩小刚度矩阵的带宽，节 省存储、提高计算效率。
§4-5 边界条件的处理和整体刚度矩阵的修正， 计算实例
整体刚度矩阵[K]是奇异阵，必须考虑边界约束条件，排除弹性体的刚
体位移。消除了整体刚度矩阵的奇异性之后，才能从方程组 K F
中求解节点位移。 一般情况下，所考虑问题的边界往往已有一定的位移约束条件，排除
了刚体运动的可能性。否则，应当适当指定某些节点的位移值，以避免出 现刚体运动。在引用这些边界条件以后，待求节点未知量的数目和方程的 数目便可相应地减少。
②将计算对象进行离散化，即弹性体划分为许多三角形单元 ，并对节点进行编号。确定全部节点的坐标值，对单元进行 编号，并列出各单元三个节点的节点号。
③ 计算载荷的等效节点力（要求的输入信息）。
④ 由各单元的常数bi 、ci 、bj 、cj 、bm 、cm 及行列式2 ，
计算单元刚度矩阵。
返回
⑤ 组集整体刚度矩阵，即形成总刚的非零子矩阵。 ⑥ 处理约束，消除刚体位移。
⒈ 保持方程组为2n×2n系统，仅对[K]和{R}进行修正。 例如，若指定节点i在方向y的位移为vi ，则令[K]中的元素k2i, 2i 为1，而第2i行和第2i列的其余元素都为零。{R}中的第2i个元 素则用位移vi 的已知值代入，{R}中的其它各行元素均减去已 知节点位移的指定值和原来[K]中该行的相应列元素的乘积。
k2n,2n 2n2n vn 2n1
根据约束情况若在第一点的水平位移为: u1= β1，在第 二节点的水平位移为: u2 = β3，把节点所对应刚度矩阵的 行和列第一行和第一列及第三行和第三列, 除主对角元改 成1，其余的元素都改成零，同时把左端的{F}载荷列阵中
对应的行改为己知位移值β1,β3 ，其余的行都减去节点位
模型。
返回
例如，图4-11(a)所示受纯弯曲的梁，其结构对于x、y轴都是几何对 称的，而所受的载荷则是对于x轴对称，对于x轴反对称。可知，梁 的应力和变形也将具有同样的对称特性，所以只需取1/4梁进行计算 即可。取分离体如图4-11(b)所示，对于其它部分结构对此分离体的 影响，可以作相应的处理，即对处于y轴对称面内各节点的x方向位 移都设置为零，而对于在x轴反对称面上的各节点的x方向位移也都 设置为零。这些条件就等价于在图4-11(b)中相应节点位置处施加约 束，图中o点y方向施加的约束是为了消除刚体位移。
移值与原来刚度矩阵该行的相应列元素的乘积。
1
Fy1 1k21 3k23
1 0 0 k22
3 Fy2 1k41 3k43
0 0 0 k42
Fx3 1k51 3k53
Fy3 1k61 3k63
M
M
M M
M
M
Fxn
k1 2n1,1
k3 2n1,3
Fyn 1k2n,1 3k2n,3 2n1
LLLLL LLLLL LLLLL LLLLL MMMMM
k1, 2 n 1 k2, 2 n 1 k3, 2 n 1 k4, 2 n 1
M
MMMMM M
MMMMM M
MMMMM M
L L L L L k2n1,2n1 L L L L L k2n,2n1
k1, 2 n k2,2n
u1
v1
k3, 2 n
y
yRBiblioteka RRox
o
x
R
R
(b) 返回
(a)
图4-11
二. 节点的选择及单元的划分
节点的布置是与单元的划分互相联系的。通常，集中载荷 的作用点、分布载荷强度的突变点，分布载荷与自由边界的 分界点、支承点等都应该取为节点。并且，当物体是由不同 的材料组成时，厚度不同或材料不同的部分，也应该划分为 不同的单元。
时，对取绝对值，即可得到正确的计算结果。
五. 边界条件的处理及整体刚度矩阵的修正 整体刚度矩阵的奇异性可以通过考虑边界约束条件来排
除弹性体的刚体位移，以达到求解的目的。
返回
一般情况下，求解的问题的边界往往已有一点的位移约束 条件，本身已排除了刚体运动的可能性。否则的话，就必须适 当指定某些节点的位移值，以避免出现刚体位移。这里介绍两 种比较简单的引入已知节点位移的方法，这两种方法都可保持 原[K]矩阵的稀疏、带状和对称等特性。
q
g
Q f Ne f * N * e
p
f * T * eT N T
j i
代入上式，得
* eT Re * eT N T g N T qtds N T ptdxdy
Re NT g N T qtds N T ptdxdy
由此可知
由集中力引起的等效节点力 Ge NT g
k1,2n1 k2,2n1 k3,2n1 k4,2n1
M
MMMMM M
MMMMM M
MMMMM M
L L L L L k2n1,2n1 L L L L L k2n,2n1
k1,2n
k2,2n
u1
v1
k3,2n k4,2n
uv22
M
M
M
M
M
M
M M
k2n
1, 2 n
un
Fx3
Fy3
M
M
M M
M M
Fxn
k2n
1,1
Fyn 2n1 k2n,1
k12 k22 k32 k42 M M M M k2n1,2 k2n,2
k13 k23 k33 k43 M M M M k2n1,3 k2n,3
LLLLL LLLLL LLLLL LLLLL MMMMM
Lm
px
Lm py
形成载荷列阵{F}
把各单元上的等效节点力{R}e根据单元的编号迭加到载荷 列阵{F}对应行中
Rix
e
Riy
Re
R jx R jy
Rmx
Rmy
F 2n1
Re F0
{R } e = {F} e +{Q} e +{P} e
{F0} 表示作用在各节点上的集中力
N N
j j
gx gy
Lj Lj
gx gy
0
N
m
g
x
Lm
g
x
Nm
Nm g y Lm g y
表面分布力的等效节点力
由于
q
qx qy
Ni Li
Niqx
Niqy
Liqx
Li q y
Qe
N T qtds t
N N
j qx jqy
ds
平面问题的半带宽为
B =2 (d+1)
若采取带宽压缩存储，则整体刚度矩阵的存储量N 最多为 N =2nB = 4n (d+1)
其中 d为相邻节点的最大差值，n为节点总数。
例如在图4-13中，(a)与(b)的单元划分相同，且节点总数都等 于14，但两者的节点编号方式却完全不同。(a)是按长边进行编 号，d =7，N =488；而(b)是按短边进行编号，d =2，N =168。 显然(b)的编号方式可比(a)的编号方式节省280个存储单元。