有限元单元法及程序设计

⎥ ⎥
2
kjj 2 ⎥⎦ 3
(2) 两端支承条件的引入
先不考虑约束条件,得到整体刚度矩阵后,将其主对角线元素 k ii 改为 1,第 i 行,第 j 列 其余元素改为 0,对应的载荷元素也改为 0.
(3) 非结点荷载的处理
利用等效结点荷载进行分析:
1 各结点(包括两端结点)加约束,阻止结点转动,其约束力矩分别为交于该结点的各相关单
0
0
0⎤⎡
Xi
e
⎤
⎢⎥ ⎢
⎥⎢ ⎥
⎢ Yi ⎥ ⎢− sinα cosα 0 0 ⎢⎥ ⎢
0 0⎥⎢ Yi ⎥ ⎥⎢ ⎥
⎢ ⎢
M
i
⎥ ⎥
⎢ =⎢
0
01 0
0
0⎥ ⎥
⎢ ⎢
M
i
⎥ ⎥
⎢ ⎢
X
j
⎥ ⎥
⎢0 ⎢
0
0
cosα
sin α
0⎥⎢ ⎥⎢
X
j
⎥ ⎥
⎢ ⎢
Yj
⎥ ⎥
⎢0 ⎢
0
0
− sinα
cosα
0
⎥ ⎥
单元①: K①= ⎢⎢kji kjj 0⎥⎥ 2
单元② :
K②= ⎢⎢0
kii
kij
⎥ ⎥
2
⎢⎣ 0 0 0⎥⎦ 3
⎢⎣0 kji kjj⎥⎦ 3
②将单元贡献矩阵想叠加,形成整体刚度矩阵
⎡kii
K=
K①+
K②=
⎢ ⎢
kji
1
⎢⎣ 0
123 k1
ij
k k 1
2
jj + ii
k2 ji
0 ⎤1
kji
2
⎥ ⎥ ⎥ ⎥
⎡ui
e
⎤
⎢⎥
⎢ vi ⎥ ⎢⎥
⎢θ ⎢
i
⎥ ⎥
⎢u ⎢
j
⎥ ⎥
⎢ ⎢
v
j
⎥ ⎥
⎢⎣θ
j
⎥ ⎦
l⎦
(1-17)
上式可简写成
e
ee
F =k δ
(1-18)
其中单元刚度矩阵为
⎡ EA
⎢ ⎢
l
⎢0
⎢
e
⎢ ⎢
0
k
=
⎢ ⎢−
EA
⎢l
⎢ ⎢
0
⎢
⎢0
⎣
0
12EI l3 6EI l2
0
12EI −
l3 6EI l2
4.根据国家标准（GB-1526-89）规定的程序流程图标准化符号及规定 ：
a) 图表示程序流程图的起点和终点； b) 图表示数据信息的输入和输出； c) 图表示数据进行系列运算之前要完成的数据预置； d) 图表示判断条件； e) 图表示各种处理功能，如数学运算方式等； f) 图表示流程的路径和指向。
件全为零。将该支承条件引入到整体刚度方程，得
⎡ P1 ⎤ ⎡ K11 K12 K13 K14 ⎤ ⎡ 0 ⎤
⎢⎥⎢
⎥⎢ ⎥
⎢P2 ⎥ ⎢K 21 K 22 K 23 K 24 ⎥ ⎢∆ 2 ⎥
⎢ ⎥=⎢
⎥⎢ ⎥
⎢ ⎢
P3
⎥ ⎥
⎢ ⎢
K
31
K 32
K 33
K
34
⎥ ⎥
⎢ ⎢
∆
3
⎥ ⎥
⎢⎣P4 ⎥⎦ ⎢⎣K 41 K 42 K 43 K 44 ⎥⎦ ⎢⎣ 0 ⎥⎦
[ ] ∆ = ∆1 ∆ 2 ∆3 ∆ 4 T [ ] = u1v1θ1 u2v2θ 2 u3v3θ 3 u4v4θ 4 T
Pix、Piy、Piθ分别代表作用在结点 i（i=1，2，3，4）上的水平力、竖向力和力偶。规定， 结点力 Pix、Piy 的正方向与整体坐标系 x、y 的正方向相同，Piθ以顺时针指向为正；结点
量。
为了导出整体坐标系中杆端力 Xi、Yi、Mi 和局部坐标系中 X i、Yi、Z i 之间的关系，将 Xi、 Yi 分别向 x、y 轴上投影，可得
X i = X i cosα + Yi sinα ⎫ ⎬
a)
Yi = − X i sinα + Yi cosα ⎭
式中，α表示由 x 轴到 x 轴之间的夹角，以顺时针为正。
Xi
=
EA l
(u
i
−uj)
=
EA l
u
i
−
EA l
u
j
⎫ ⎪⎪
X
j
=
−
EA l
(ui
−uj)
=
−
EA l
ui
+
EA l
u
j
⎬ ⎪ ⎪⎭
（a）
其次考虑杆端弯矩 M i、M j 与杆端剪力 Yi、Y j 与杆端转角θ i、θ j 和横向位移 vi、v j 的关系。
根据结构力学位移法的转角位移方程，并按照本节规定的符号和正负号，可得
l 0
0 EA l 0
0
⎣
0
12EI l3 6EI l2
0
12EI − l3
6EI l2
0
6EI l2 4EI l
0
6EI − l2 2EI
l
EA −
l 0
0 EA l 0
0
0
12EI −
l3 6EI − l2
0
12EI l3 6EI
− l2
e
0
⎤ ⎥
6EI
⎥ ⎥
l2 ⎥
2EI ⎥
l
⎥ ⎥
0⎥
⎥
6EI − l2 4EI
有限单元法及程序设计
绪论 1.力学分析方法:解析法,数值法 有限元法——实际结构形状和所受载荷比较复杂,大多用解析法很困难,因而数值法得到不断发 展,随着电子计算机的进步,而发展起来的一种新兴的数值分析方法. 2 基本步骤： （1）结构离散化：将结构从集合上用线或面划分为有限个单元。 （2）单元分析： 导出单元的节点位移和结点力之间的关系（单元刚度矩阵）。 （3）整体分析： 将各单元组成的结构整体进行分析，导出征个结构点位移与结点力之间的关 系。 3 程序设计的步骤： （1） 提出问题，拟定解决方案 （2） 构造数学模型 （3） 画出程序流程图 （4） 编写程序 （5） 编译调试程序 （6） 试算验证程序
+
k3 ii
k3 ji
0 ⎤⎡∆1 ⎤ ⎥⎢ ⎥
0
⎥ ⎥
⎢ ⎢
∆
2
⎥ ⎥
k3 ii
⎥ ⎥
⎢ ⎢
∆
3
⎥ ⎥
k
3 jj
⎥⎦
⎢⎣
∆
4
⎥⎦
（1-37）
或 其中
P0 = K0∆0
（1-38）
⎡ K11
K12
K13
K14 ⎤
⎡
k1 ii
k1 ij
⎢
⎥⎢
0 0⎤ ⎥
⎢K 21
K0 = ⎢
⎢ ⎢
K
31
K 22 K 32
即 T e −1 = T eT
第四节 单元未知量编码
为了便于编程计算，需要按一定规律对结点的位移分量编号。结构的节点位移有自由 结点位移和支座结点位移（亦称支座结点位移）之分。自由结点位移是未知量。建立结构 整体结构方程求解未知节点位移的方式有两种：“前处理法”和“后处理法”。
用后处理法分析改结构时，设所有点位移都是未知量，则结点位移列阵为（参看图 1-9））
[P= M 1 M 2
k13⎤
k
⎥ 23⎥
Biblioteka Baidu
k 33⎥⎦
M 3]
整体刚度矩阵 节点载荷列阵
（1-7）
∆ = [θ 1 θ 2 ]θ 3 T
3.有限元位移法分析连续梁需要考虑的问题
(1) 刚度集成法:
①将(1-3)K 扩阶,扩大的元素为 0,得到单元贡献矩阵
123
123
⎡kii kij 0⎤ 1
⎡0 0 0 ⎤ 1
在建立整体刚度方程式（1-38）时，假定所有点位移都是未知量，相当整体结构无知
座，因而在外力作用下，除了弹性变形外，还有可能发生刚体位移，此时，各结点位移不
能唯一确定。这说明式（1-39）为奇异矩阵，不能求逆，故利用（1-38）不能求结点位移。
实际在图 2-9a）所示的刚架中，结点 1、4 为固定点，因此结点位移是已知的，支承条
在此单元中，单元杆端力列阵和杆端位移列阵分别为
[ ]T
F = Xi Yi Mi X j Y j M j
[ ]T
δ = ui vi θi u j vj θ j
单元杆端力列阵 杆端位移列阵
为了导出一般单元杆端力与杆端位移之间的关系，我们分别考虑以下两种情况。
首先分析两个杆端轴力 X i、X j 与轴向位移 ui、u j 的关系。根据胡克定律，有
图 1-8 在两个坐标系中，力偶分量不变，即
Mi = Mi
b)
同理，对于单元 ○e j 端的杆端力可得
X j = X j cosα + Y j sinα ⎫
Yj
=
−X
j
sin α
+Yj
cosα
⎪⎪ ⎬
c)
Mj =Mj
⎪ ⎪⎭
将 a)、b)、c)合起来，并用矩阵形式表示，可得
⎡
Xi
e
⎤
⎡ cosα
sinα 0
构，分析时所划分的各单元的局部坐标系显然不同。因此在研究结构平衡条件和变形连续条件
时，必须选定一个统一的坐标系 xOy,称为整体坐标系。同时，还必须把在局部坐标系中建立的
单元刚度矩阵转换为整体坐标系下的单元刚度矩阵。
图 1-8a)、图 1-8b)分别表示单元○e 在局部坐标系 xOy 和整体坐标系 xOy 种的杆端力分
第一篇 杆件结构的有限单元法及程序设计
第一章 平面杆件单元的有限单元法 第一节 有限单元法的基本概念
1． 基本思路：先分后合（先单元分析，再整体分析）
2． 基本概念：整体号：节点端点号按自然数 1，2，3，……（在整体坐标系 xOy 下）
局部号：每一个单元始末用 i，j 标记（在单元的局部坐标 xyz 系下，方向与整
⎢ ⎢
Y
j
⎥ ⎥
⎢⎣M
j
⎥ ⎦
⎢⎣ 0
00 0
0
1⎥⎦
⎢⎣M
j
⎥ ⎦
（1-24）
此式即为两种坐标系中单元杆端力的变换式，亦可简写为
e
F
= T eF e
（1-25）
[ ] e
式中： F = X i Yi
Mi
Xj Yj
Mj T
局部坐标系中的单元杆端力列阵
[ ] F e = Xi Yi Mi Xj Yj Mj T 整体坐标系中的单元杆端力列阵
杆端位移分量等于 1 时,所引起的第 r 各杆端力分量值.
2 是对称矩阵,其元素 krs = ksr(r ≠ s) .
3 是奇异矩阵,它的元素行列式等于零,即 k = 0 .
4 具有分快性质.
3. 轴力单元:只考虑轴向杆端位移和杆端力的单元
第三节 单元刚度矩阵的坐标变换
上述单元刚度方程和单元刚度矩阵实在局部坐标系 xOy 中建立起来的，对于一般杆件结
Mi
=
6EI l2
vi
+
4EI l θi
−
6EI l2
vj
+
2EI l θj
⎫ ⎪ ⎪
M
j
=
6EI l2
vi
+
2EI l
θi
−
6EI l2
vj
+
4EI θ l
j
⎪ ⎪
12EI 6EI 12EI 6EI Yi = l 3 vi + l 2 θ i − l 3 v j + l 2 θ j
⎬ ⎪ ⎪
(b)
⎡ cosα sinα 0 0
0 0⎤
⎢
⎥
⎢− sinα cosα 0 0
0 0⎥
⎢
⎥
⎢0 Te = ⎢
01 0
0 0⎥ ⎥
⎢0 ⎢
0 0 cosα sinα 0⎥ ⎥
⎢0 ⎢
0 0 − sinα cosα 0⎥ ⎥
⎢⎣ 0
00 0
0 1⎥⎦
T 为正交矩阵,其逆矩阵等于其转置矩阵
单元坐标变换矩阵(1-26)
Yj
=
12EI − l3
vi
−
6EI l2 θi
12EI + l3
vj
−
6EI l2
θ
j
⎪ ⎪ ⎭
将(a)、(b)两式合在一起，并写成矩阵形式如下
⎡ EA
e
⎡ Xi ⎤ ⎢⎥
⎢ ⎢
Yi
⎥ ⎥
⎢ ⎢
M
i
⎥ ⎥
⎢ ⎢
X
j
⎥ ⎥
⎢ ⎢
Yj
⎥ ⎥
⎢⎣M
j
⎥ ⎦
=
⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢− ⎢ ⎢ ⎢ ⎢ ⎢
（1-40）
可以分为两组方程，一组是
⎡ P2 ⎢
⎤ ⎥
=
⎡K ⎢
22
⎢⎣P3 ⎥⎦
⎢ ⎣
K 23 K 33
K 24
⎥ ⎥
=
⎢k ⎢
1 ji
K
34
⎥ ⎥
⎢0 ⎢
k1 jj
+
k2 ii
k2 ji
k2 ij
k2 jj
+
k3 ii
0⎥
⎥
k3 ii
⎥ ⎥
（1-39）
⎢⎣K 41 K 42 K 43 K 44 ⎥⎦ ⎢⎣ 0
0
k3 ji
k
3 jj
⎥⎦
为结构的整体刚度矩阵，或称为结构的原始刚度矩阵。
0
6EI l2 4EI l
0
6EI −
l2 2EI
l
EA −
l 0
0 EA l 0
0
0
12EI − l3
6EI −
l2
0
12EI l3 6EI
− l2
0
⎤ ⎥
6EI
⎥ ⎥
l2 ⎥
2EI ⎥
l
⎥ ⎥
0⎥
⎥
6EI −
l2 4EI
⎥ ⎥ ⎥ ⎥
l⎦
（1-19）
2. 单元刚度矩阵的性质
1 每个元素代表单位杆端位移引起的杆端力,任一元素 krs (r、s 取 1 至 6)的物理意义是第 s 个
元的固端力矩之和,顺时针为正.
2 去掉附加约束(相当在各结点施加外力荷载 P e ,其大小与约束力矩相同,方向相反) 3 将两部分杆端弯矩叠加起来.
第二节 局部坐标系中的单元刚度矩阵 1. 一般单元
设单元○e 的弹性模量、截面惯性矩、截面积分别为 E、I、A,杆长为 l。单元的 i、j 端各有
三个杆端力、 X、Y和M (即轴力、剪力和弯矩)和与其相应的三个杆端位移 u、v、θ ，如图 1-7 所示。图中 xOy 为单元局部坐标系，取 i 点位于坐标原点， x 轴与杆轴重合，规定由 i 到 j 为 x 轴的正方向，由 x 轴顺时针旋转 90◦为 y 轴正方向。力和位移的正方向如图 1-7 所示。
体坐标系一致）。
3．Fe= keδe 其中：
ke
⎡kii
=
⎢ ⎣
kji
kij ⎤
kjj
⎥ ⎦
单元刚度矩阵，各元素为刚度系数
δe=
⎡θi⎤
⎢⎣θj
⎥ ⎦
单元杆段位移列阵
Fe=
⎡ ⎢ ⎣
Mi Mj
⎤ ⎥ ⎦
单元杆端力列阵
K ∆ =P
⎡k11 k12 K= ⎢⎢k 21 k 22
⎢⎣k 31 k 32
位移列阵
位移的正方向与结点力的正方向一致。
在求出各单元刚度方程之后，根据结点平衡条件和位移连续条件，可建立整个结构的
位移法方程
⎡ P1 ⎤
⎡
k
1 ii
⎢⎥⎢
⎢ P2 ⎢
⎥ ⎥
=
⎢ ⎢
k
1 ji
⎢ ⎢
P3
⎥ ⎥
⎢0 ⎢
⎢⎣P4 ⎥⎦ ⎢⎣ 0
k1 ij
k1 jj
+
k2 ii
k2 ji
0
0
k2 ij
k2 jj