有限元单元法及程序设计

有限单元法及程序设计

绪论

1.力学分析方法:解析法,数值法

有限元法——实际结构形状和所受载荷比较复杂,大多用解析法很困难,因而数值法得到不断发展,随着电子计算机的进步,而发展起来的一种新兴的数值分析方法.

2基本步骤：

（1）结构离散化：将结构从集合上用线或面划分为有限个单元。

（2）单元分析：导出单元的节点位移和结点力之间的关系（单元刚度矩阵）。

（3）整体分析：将各单元组成的结构整体进行分析，导出征个结构点位移与结点力之间的关系。

3程序设计的步骤：

（1）提出问题，拟定解决方案（2）构造数学模型（3）画出程序流程图（4）编写程序

（5）编译调试程序（6）试算验证程序

4.根据国家标准（GB-1526-89）规定的程序流程图标准化符号及规定

：

a) b) c) d) e) f) 图表示程序流程图的起点和终点；图表示数据信息的输入和输出；

图表示数据进行系列运算之前要完成的数据预置；图表示判断条件；

图表示各种处理功能，如数学运算方式等；图表示流程的路径和指向。

第一篇

杆件结构的有限单元法及程序设计

第一章平面杆件单元的有限单元法第一节有限单元法的基本概念

1．基本思路：先分后合（先单元分析，再整体分析）2．基本概念：整体号：节点端点号按自然数1，2，3，……（在整体坐标系xOy 下）

局部号：每一个单元始末用i，j标记（在单元的局部坐标系下，方向与整

体坐标系一致）。3．Fe =k e δe 其中：

k e =⎢

⎡k ii

⎣k ji

e

k ij ⎤k jj ⎥⎦

单元刚度矩阵，各元素为刚度系数单元杆段位移列阵单元杆端力列阵

⎡θi ⎤

δ=⎢⎥

⎣θj ⎦

F e =⎢

⎡M i ⎤⎥⎣M j ⎦

K ∆=P

（1-7）

整体刚度矩阵

⎡k 11k 12k 13⎤⎢⎥K=k 21k 22k 23⎢⎥⎢⎣k 31k 32k 33⎥⎦

位移列阵P=[M 1M 2

∆=[θ1θ2θ3]T

M 3]

节点载荷列阵

3.有限元位移法分析连续梁需要考虑的问题(1)刚度集成法:

①将(1-3)K扩阶,扩大的元素为0,得到单元贡献矩阵

123⎡k ii k ij ⎢单元①:K ①=k ji k jj ⎢⎢⎣00⎡k ii

①②⎢1

K=K +K =k ji

⎢⎢⎣0

0⎤1

0⎥⎥20⎥⎦3

0⎤1

k ji 2⎥⎥2k jj 2⎥⎦3

123

⎡000⎤1⎢⎥单元②:K ②=0k ii k ij 2⎢⎥⎢⎣0k ji k jj ⎥⎦3

②将单元贡献矩阵想叠加,形成整体刚度矩阵

123

k ij 1

k jj 1+k ii 2

k ji 2

(2)两端支承条件的引入

先不考虑约束条件,得到整体刚度矩阵后,将其主对角线元素k ii 改为1,第i 行,第j 列其余元素改为0,对应的载荷元素也改为0. (3)非结点荷载的处理

利用等效结点荷载进行分析:

各结点(包括两端结点)加约束,阻止结点转动,其约束力矩分别为交于该结点的各相关单元的固端力矩之和,顺时针为正.

去掉附加约束(相当在各结点施加外力荷载P e ,其大小与约束力矩相同,方向相反)将两部分杆端弯矩叠加起来.

第二节局部坐标系中的单元刚度矩阵

1. 一般单元

e 的弹性模量、截面惯性矩、截面积分别为E、I、A,杆长为l 。单元的i 、j 端各有设单元○三个杆端力、(即轴力、剪力和弯矩) 和与其相应的三个杆端位移，如图

1-7所示。图中为单元局部坐标系，取i 点位于坐标原点，

轴与杆轴重合，规定由i 到j 为轴的正方向，由轴顺时针旋转90◦为轴正方向。力和位移的正方向如图1-7所示。

在此单元中，单元杆端力列阵和杆端位移列阵分别为

=F =X i

i

Y i M i X

j

Y

j

i i j j j

]

M

j

]

T

单元杆端力列阵杆端位移列阵

T

为了导出一般单元杆端力与杆端位移之间的关系，我们分别考虑以下两种情况。首先分析两个杆端轴力i 、j 与轴向位移i 、j 的关系。根据胡克定律，有

EA EA EA ⎫(i −j ) =i −j

⎪⎪l l l （a）⎬

EA EA EA ⎪

j =−(i −j ) =−i +j ⎪

l l l ⎭

其次考虑杆端弯矩i 、j 与杆端剪力i j 与杆端转角i 、j 和横向位移i 、j 的关系。

i =

根据结构力学位移法的转角位移方程，并按照本节规定的符号和正负号，可得

6EI 4EI 6EI 2EI ⎫

+−+i i j j ⎪l l l 2l 2

⎪

6EI 2EI 6EI 4EI j =2i +i −2j +j ⎪

⎪l l l l

⎬

12EI 6EI 12EI 6EI

i =3i +2i −3j +2j ⎪

⎪l l l l

12EI 6EI 12EI 6EI ⎪j =−3i −2i +3j −2j ⎪

l l l l ⎭

i =

(b)

将(a)、(b)两式合在一起，并写成矩阵形式如下

⎡i ⎤

⎢⎥⎢i ⎥⎢⎥⎢i ⎥⎢⎥⎢⎥⎢j ⎥⎢⎥⎢j ⎥⎢⎥⎣j ⎦

e

⎡EA

0⎢l

⎢12EI ⎢0

l 3⎢

6EI ⎢0⎢2

=⎢

⎢−0⎢l

12EI ⎢0−⎢l 3⎢6EI ⎢0

l 2⎣

06EI

l 24EI 0−6EI l 22EI l

e

EA −l 00

0−12EI l 36EI −2

012EI l 36EI −2

l

l 00

e

e

⎤0⎥6EI ⎥

⎥2

l ⎥2EI ⎥⎥⎥0⎥

⎥6EI ⎥−2⎥l 4EI ⎥

⎥l ⎦

e

⎡i ⎤⎢⎥⎢i ⎥⎢⎥⎢i ⎥⎢⎥⎢j ⎥⎢⎥⎢⎥⎢j ⎥⎢⎥⎣j ⎦

e

(1-17)

上式可简写成

其中单元刚度矩阵为

F =k δ

012EI l 36EI 20−12EI l 36EI l 2

−06EI l 24EI 06EI l 22EI l

−

(1-18)

⎡EA ⎢l ⎢⎢0⎢⎢⎢0e

k =⎢

⎢−⎢l ⎢⎢0⎢⎢0⎣

EA l 00

0−12EI l 36EI −2

012EI l 36EI −2

l

l 00

⎤⎥6EI ⎥

⎥2

l ⎥2EI ⎥⎥⎥0⎥

⎥6EI ⎥−2⎥l 4EI ⎥

⎥l ⎦0

（1-19）

2. 单元刚度矩阵的性质

每个元素代表单位杆端位移引起的杆端力,任一元素k rs (r、s取1至6)的物理意义是第s 个杆端位移分量等于1时,所引起的第r 各杆端力分量值.是对称矩阵,其元素k rs =k sr (r ≠s ) .

是奇异矩阵,它的元素行列式等于零,即具有分快性质.

k =0.

3. 轴力单元:只考虑轴向杆端位移和杆端力的单元

第三节单元刚度矩阵的坐标变换

上述单元刚度方程和单元刚度矩阵实在局部坐标系中建立起来的，对于一般杆件结构，分析时所划分的各单元的局部坐标系显然不同。因此在研究结构平衡条件和变形连续条件时，必须选定一个统一的坐标系xOy,称为整体坐标系。同时，还必须把在局部坐标系中建立的单元刚度矩阵转换为整体坐标系下的单元刚度矩阵。

e 在局部坐标系和整体坐标系xOy 种的杆端力分图1-8a)、图1-8b)分别表示单元○

量。

为了导出整体坐标系中杆端力X i 、Y i 、M i 和局部坐标系中i 、i i 之间的关系，将X i 、Y i 分别向轴上投影，可得

i =X i

cos α+Y i sin α⎫

⎬

i =−X i sin α+Y i cos α⎭

式中，α表示由x 轴到轴之间的夹角，以顺时针为正。

a)

图

在两个坐标系中，力偶分量不变，即e j 端的杆端力可得同理，对于单元○

1-8

b)

i =M i

j =X j cos α+Y j sin α⎫

⎪⎪

j =−X j sin α+Y j cos α⎬

⎪

j =M j ⎪⎭

将a)、b)、c)合起来，并用矩阵形式表示，可得

c)

⎡X i ⎤⎡cos α

⎢⎥⎢⎢Y i ⎥⎢−sin α⎢⎥⎢⎢M i ⎥⎢0⎢⎥=⎢⎢X j ⎥⎢0⎢⎥⎢⎢Y ⎥⎢0j ⎢⎥⎢⎢M ⎥⎢0

⎣⎣j ⎦

e

sin αcos α0000

0010

000cos α

000sin αcos α0

0−sin α0

0⎤⎡X i ⎤

⎥⎥⎢0⎥⎢Y i ⎥

⎥⎥⎢0⎥⎢M i ⎥

⎥⎥⎢0⎥⎢X j ⎥

⎥⎥⎢0⎥⎢Y j ⎥

⎥⎥⎢⎢⎥1⎥⎦⎣M j ⎦

e

（1-24）

此式即为两种坐标系中单元杆端力的变换式，亦可简写为

F =T e F e

式中：F =i

e

e

（1-25）

i i j j j Y j M j ]00010

00cos α

T

]

T

局部坐标系中的单元杆端力列阵整体坐标系中的单元杆端力列阵

F e =[X i Y i M i X j

⎡cos αsin α⎢

⎢−sin αcos α⎢⎢00e

⎢T =

⎢00⎢⎢00⎢⎢00⎣

即T

000sin αcos α0

0−sin α0=T e

T

0⎤

⎥0⎥⎥0⎥⎥0⎥⎥0⎥⎥1⎥⎦

单元坐标变换矩阵(1-26)

T 为正交矩阵,其逆矩阵等于其转置矩阵

e −1

第四节单元未知量编码

为了便于编程计算，需要按一定规律对结点的位移分量编号。结构的节点位移有自由

结点位移和支座结点位移（亦称支座结点位移）之分。自由结点位移是未知量。建立结构整体结构方程求解未知节点位移的方式有两种：“前处理法”和“后处理法”

。

用后处理法分析改结构时，设所有点位移都是未知量，则结点位移列阵为（参看图1-9））

∆=[∆1=[u 1v 1θ1

∆2∆3

∆4]

T

u 2v 2θ2u 3v 3θ3

u 4v 4θ4]

T

P ix 、P iy 、P i θ分别代表作用在结点i（i=1，2，3，4）上的水平力、竖向力和力偶。规定，结点力P ix 、P iy 的正方向与整体坐标系x、y的正方向相同，P i θ以顺时针指向为正；结点

位移的正方向与结点力的正方向一致。

在求出各单元刚度方程之后，根据结点平衡条件和位移连续条件，可建立整个结构的位移法方程

⎡P 1⎤⎡k ii ⎢⎥⎢⎢P 2⎥⎢k 1ji

⎢⎢⎥=⎢P 3⎥⎢0⎢⎥⎢

⎢⎢⎣P 4⎥⎦⎣0

或其中

1

1

k ij

k 1jj +k ii 2

k

2ji

2jj

k ij 2k +k k 3ji

3

ii

0⎤⎡∆1⎤⎥⎢⎥0⎥⎢∆2⎥⎥⎢⎥3⎥⎢∆3⎥k ii

⎥⎢⎥3⎥k jj ⎦⎢⎣∆4⎥⎦

（1-38）

（1-37）

P 0=K 0∆0K 12K 22K 32K 42

K 13K 23K 33K 43

1

K 14⎤⎡k ii

⎥⎢K 24⎥⎢k 1ji

⎥=⎢K 34⎥⎢0⎥⎢⎢K 44⎥⎦⎣0

⎡K 11

⎢⎢K 21

K 0=⎢

⎢K 31⎢⎢⎣K 41

1

k ij

k 1jj +k ii 2

k

2ji

2jj

k ij 2k +k k 3ji

3

ii

0⎤⎥0⎥⎥3⎥k ii ⎥⎥k 3jj ⎦

（1-39）

为结构的整体刚度矩阵，或称为结构的原始刚度矩阵。

在建立整体刚度方程式（1-38）时，假定所有点位移都是未知量，相当整体结构无知座，因而在外力作用下，除了弹性变形外，还有可能发生刚体位移，此时，各结点位移不能唯一确定。这说明式（1-39）为奇异矩阵，不能求逆，故利用（1-38）不能求结点位移。

实际在图2-9a）所示的刚架中，结点1、4为固定点，因此结点位移是已知的，支承条件全为零。将该支承条件引入到整体刚度方程，得

⎡P 1⎤⎡K 11⎢⎥⎢⎢P 2⎥⎢K 21⎢⎥=⎢⎢P 3⎥⎢K 31⎢⎥⎢⎢⎣P 4⎥⎦⎢⎣K 41⎡P 2⎤⎡K 22⎢⎥=⎢

K 32⎢⎣P 3⎥⎦⎢⎣⎡P 1⎤⎡K 22

⎢⎥=⎢

0⎢⎣P 2⎥⎦⎢⎣

K 12K 22K 32K 42

K 13K 23K 33K 43

K 14⎤⎡0⎤

⎥⎢⎥K 24⎥⎢∆2⎥⎥⎢⎥K 34⎥⎢∆3⎥⎥⎢⎥K 44⎥⎣0⎥⎦⎦⎢

（1-40）

可以分为两组方程，一组是

K 23⎤⎡∆2⎤⎥⎢⎥K 33⎥⎢∆3⎥

⎦⎣⎦0

⎤⎡∆2⎤⎥⎢⎥⎥⎢∆3⎥⎦⎣⎦

（1-41）

它可以求结未知结点位移Δ2、Δ3。令一组是

（1-42）

K 43

称为反力方程。利用式（1-41）求出结点位移Δ2、Δ3并代入上式后，便可计算未知的支座反力。

对于一般杆件结构，都可以按上述步骤进行分析。无论结构有多少个结点位移分量，

经过调整其排列顺序，总可以将它分为两组：一组包括所有的未知结点的位移分量，以ΔF 表示；另一组为支座结点位移分量，以ΔR 表示。相应的，将全部结点分为两组，与ΔF 相应者为已知的结点力列阵，以PF 表示；与ΔR 相应者为支座结点力列阵，以P R 表示。于是有

⎡∆F ⎤∆0=⎢⎥,

⎢⎣∆R ⎥⎦⎡K FF

⎢⎢⎣K RF

展开上式得

⎡P F ⎤

P 0=⎢⎥

⎢⎣P R ⎥⎦

（1-43）

与以上分析方法相配合，将整体刚度矩阵K0中的各元素重新排列，则K0Δ0=P0可写成

K FR ⎤⎡∆F ⎤⎡P R ⎤

⎥⎢⎥=⎢⎥K RR ⎥⎦⎢⎣∆R ⎥⎦⎢⎣P F ⎥⎦

K FF ∆F +K FR ∆R =P F

K RF ∆F +K RR ∆R =P R

式（1-44）、式（1-45）

式（1-44）为“修正的整体刚度方程”，它与式（1-38）的区别在于引进了支承条件。

后处理法:由单元刚度矩阵形成整体刚度矩阵,建立刚度方程后在引入支承条件,进而求解结点的位置位移的方法.前处理法:仅对未知的自由结点位移分量编号,得到的结点位移列阵中不包含已知的约束结点位移分量.

图1-10所示的具有组合结点的刚架划分为三个单元，其编号为①、②、③，各杆之间的箭头表示局部坐标系的正方向，刚架结构编号为1～5。下面考虑各单元结点位移分量编号。采用“先处理法”作如下规定：

仅对独立的位移分量按自然数编号，称为位移号。若某些位移分量由于联结条件或直接杆轴向刚性条件的限制彼此相等，则编号相同。

在支座处，由于刚性约束而使某些位移分量为零时，此位移分量编号为零。因此图1-10编号如下

单元结点编号

杆数单元编号单元位移分量编号

始端末端

AB 12000123①BC 23123456②DC 54000457③

在计算程序中，单元两端结点号可采用二维数组JE(i,e)表示，称为“单元两端结点号数组”。

JE(1,e)=单元始端的结点号JE(2,e)=单元末端的结点号在本例中

JE(1，1) =1，JE(2，1) =2JE(1，2) =2，JE(2，2) =3JE(1，3) =5，JE(2，3) =4

任意结点位移分量的位移号可用二维数组JN(i,j)表示，称为“结点位移号数组”

JN （1，j ）=结点j 沿x 方向的位移号JN （2，j ）=结点j 沿y 方向的位移号JN （2，j ）=结点j 角位移的位移号在本例中，对第三结点而言，JN （1，3）=4JN （2，3）=5JN （3，3）=6

e 始端及末端得位移号排成一行（始端在前）将单元○，此数码为“单元定位数组”，利用

它可方便的形成杆端位移与相应结点位移间的协调条件。它的展开式为

m e =(m 1, m 2,..., m d ) e

e 的两端位移分量所对应的位移号数值。在本例中，式中，d 个元素m 1～m d 分别是单元○

m 1=(000123) m 2=(123456) m 3=(000457)

第五节

平面结构的整体刚度矩阵

在进行了单元分析得出单元刚度矩阵之后，需要进行整体分析。以“先处理法”为例，

将离散单元重新组合成原结构，使其满足结构结点的位移连续条件和力的平衡条件，从而得到修正的结构刚度方程，即前面给出的式（1-44）

K FF △F =PF

式中：K FF 称为修正的结构整体刚度矩阵；△F 、P F 分别为自由结点位移与自由结点荷载列阵。当已知计算对象为自由结点位移分量而不至引起误解时，式（1-44）也常称为整体刚度方程，K FF 简称为整体刚度矩阵，△F 、P F 分别简称为结点位移列阵与荷载列阵。为了书写方便，下表常常略去。

有单元刚度矩阵集成整体刚度矩阵，通常采用“直接刚度法”，把计算步骤分为两步，首先求出各单元贡献矩阵，然后将它们叠加起来，得出整体刚度矩阵。然而在实际电算中，不便采用，原因是在计算中需要先将所有单元的贡献矩阵K e 都保存起来，K e 的阶数与整体刚度矩阵K 相同，这就占用了大量存储容量，因此实际运算中，采用“边定位，便累加”

的方法。其原理没有变，而且结果相同。

图1-10中单元①的单元矩阵为

k =

1

0⎡k 11⎢⎢k 21⎢⎢k 31⎢⎢k 41⎢⎢k 51⎢⎢k ⎣61

0k 1201k 13k 14

2

k 15k 25k 35k 45k 55k 65

k 22k 32k 42k 52k 61

k 23k 33k 43k 53k 63

k 24k 34k 44k 54k 64

31k 16⎤⎥k 26⎥⎥k 36⎥⎥k 46⎥⎥k 56⎥⎥k 66⎥⎦

000123

式中单元刚度矩阵的上面和右侧标记了单元结点位移分量编号。因为整体刚度矩阵个元素

是按位移分量编号排列的，按先处理法，单元刚度矩阵中对应于分量编号为零的元素不进入整体刚度矩阵，非零编号指明了其余各元素在整体刚度矩正中的行、列号。所以

1k 44→K 111k 54→K 211k 64→K 31

1

k 45→K 121k 55→K 221k 65→K 32

1

k 46→K 131k 56→K 231k 66→K 33

单元②的刚度矩阵为

k =

2

1

⎡k 11⎢⎢k 21⎢⎢k 31⎢⎢k 41⎢⎢k 51⎢⎢k ⎣61

2k 1234k 13k 14

5

k 15k 25k 35k 45k 55k 65

2ij

k 22k 32k 42k 52k 61

k 23k 33k 43k 53k 63

k 24k 34k 44k 54k 64

62

k 16⎤⎥k 26⎥⎥k 36⎥⎥k 46⎥⎥k 56⎥⎥k 66⎥⎦

123456

各单元在K 中的位置为

⎛i =1, 2,..., 6⎞

⎟k →K ij ⎜

⎜j =1, 2,..., 6⎟⎝⎠

单元③的刚度矩阵为

k =

3

0⎡k 11⎢⎢k 21⎢⎢k 31⎢⎢k 41⎢⎢k 51⎢⎢k ⎣61

0k 1204k 13k 14

57

k 15k 25k 35k 45k 55k 65

k 22k 32k 42k 52k 61

k 23k 33k 43k 53k 63

k 24k 34k 44k 54k 64

k 16⎤⎥k 26⎥⎥k 36⎥⎥k 46⎥⎥k 56⎥⎥k 66⎥⎦

3

000457

各元素在K 中的位置为

3

k 44→K 443k 54→K 543k 64→K 74

3

k 45→K 453k 55→K 553k 65→K 75

3

k 46→K 463k 56→K 563k 76→K 77

按以上定位方法，将三个单元的刚度矩阵有关元素一道整体刚度矩阵对应位置，得到

12

⎡k 44+k 11⎢12⎢k 54+k 21⎢12⎢k 64+k 31⎢

2

K =⎢k 41

⎢⎢k 2

51

⎢⎢k 2

61

⎢⎢⎣0

12k 45+k 1212k 55+k 2212k 65+k 32

2

k 422k 522k 62

12k 46+k 1312k 56+k 2312k 66+k 33

2

k 432k 532k 63

2k 142k 242k 3423k 44+k 4412k 44+k 11

2

k 643k 64

2k 152k 252k 3523k 45+k 4512k 44+k 11

2

k 653k 65

2k 162k 262k 362k 462k 562k 66

000

0⎤⎥0⎥⎥0⎥⎥3⎥k 46⎥3⎥k 56⎥0⎥⎥3⎥k 66⎦

在世机电算时，采用

将K 置零,这是K=07×7;

将k ①中的相关元素,按照“对号入座”，累加到K；

将k ②中的相关元素，继续按照“对号入座”，累加到K；

将k ③中的相关元素，继续按照“对号入座”，累加到K，整体刚度矩阵最后完成。

例1-3求图1-12所示刚架的整体刚度矩阵K 。设各杆截面尺寸相同。

A=0.5m2I=1/24m 2

E=3×104MPa

解

（1）整理数据并进行编号。

4EA

=100×104l 12EA

=12×1043l

e

（2）求局部坐标系中单元刚度矩阵。由于单元

12

①、②的尺寸完全相同，故有=，可直接利用公式

（1-19）求得

EA

=300×104l 6EI

=30×1042l

00−30000⎤⎡300

⎢⎥⎢012300−1230⎥⎢⎥⎢0301000−3050⎥

12

⎥==⎢

⎢−3000030000⎥⎢⎥⎢0−12−30012−30⎥⎢⎥⎢030500−30100⎥⎣⎦

（3）求整体坐标系中单元刚度矩阵k 。按式（1-32）和式（1-33），求得各单元在整体坐

标系中的单元刚度矩阵，并将单元结点位移分量编号标记于上面与左侧。

e

1230001⎡120−30−120−30⎤⎢⎥2⎢030000−3000⎥⎢⎥3⎢−[1**********]⎥1

⎥×104k =⎢

0⎢−1203012030⎥⎢⎥0⎢0−300003000⎥⎢⎥0⎢050300100⎥⎣−30⎦

1230041⎡30000−30000⎤⎢⎥2⎢012300−1230⎥⎢⎥3⎢0301000−3050⎥

12

⎥×104==⎢

0⎢−3000030000⎥⎢⎥0⎢0−12−30012−30⎥⎢⎥4⎢30500−30100⎥⎣0⎦

（3）形成整体刚度矩阵K。采用本节中介绍的方法建立结构整体刚度矩阵如下

0−300⎤⎡312

⎢⎥⎢03123030⎥

⎥×104K =

⎢

⎢−303020050⎥⎢⎥⎢3050100⎥⎣0⎦

第六节非结点荷载处理

为分析平面结构而建立的整体刚度方程，反映了结构的结点荷载与结点位移之间的关系。

作用在结构上的荷载除了直接作用在结点上的荷载P d 之外，还有作用在杆件上的分布荷载等。这些非结点荷载应转换成等效结点荷载P e 。将P d 和P e 叠加，课的综合结点荷载（总结点荷载）P c ,其下标c 通常可略去不写，即

P=Pd +Pe (1-46)

直接作用在结点上的荷载，可按其作用方位直接加入P 之中，而等效结点荷载的计算步骤如下。

e 的固端力f e 。第一步：在局部坐标系下，求单元○

f

e

⎡fi e ⎤

⎥=fi =⎢e ⎢fj ⎥⎣⎦

[

fi fi fj fj

fj ]

eT

（1-47）

e 在端点i、j的固端内力。几种非结点荷载作用下的式中，子向量fi 和fj 分别为单元○单元固端力列于表1-2

表1-2

简

图

两端固定梁的固端力

剪力Q AB Q BA

弯M AB

矩

M BA

qa 3(2l −32l

2la 2+a 3) −

−

qa

(2l −a ) 32l

3

qa 2−(6l 2−212l 8la +3a 2)

qa 3

(4l −3a ) 12l 2

Pb 2

−3(l +2a ) l Pa 2

−3(l +2b ) l

Pab 2−2

l Pa 2b l 2

m

6ab l 3

−m

6ab l 3

m

b (3a −l ) l 2

m

a (3b −l ) l 2

e 的等效结点荷载P e e 。第二步：求单元○

仿照局部坐标系与整体坐标系中单元杆端力的变换式（1-25）

e =T e F e

固端内力在两种坐标系下的变换公式，可写成

e =T e F f

e

e e

F f =T eT f

e

（1-48）

e

因此，等效结点荷载列阵P e 可以由下式求出

P e =−T eT f

将该式展开得

（1-49）

⎡P e 1⎤⎡−fi cos α+fi sin α⎤

⎥⎢⎥⎢

⎢P e 2⎥⎢−fi sin α−fi sos α⎥

⎥⎢⎥⎢

⎥⎢P e 3⎥⎢−fi

⎥（1-50）⎢⎥=⎢

⎢P e 4⎥⎢−fj cos α+fj sin α⎥

⎥⎢⎥⎢

⎢P ⎥⎢−fj sin α−fj sos α⎥

⎥⎢e 5⎥⎢

⎥⎢P ⎥⎢−fj ⎣e 6⎦⎣⎦

当α=0时，P e =−f

e

e

第三步：求整体结构的等效结点荷载P e

求得单元等效结点荷载P e 之后，利用单元结点位移分量编号，就可以将P e 中的各分量叠加到结构等效荷载列阵P e 中去。因为P e 中的各元素是按结点位移分量编号排列的，P e 中的6个元素也与结点位移分量编号一一对应，所以也按对号入座方法，将其之一累加到P e 中相应的位置上去。当直接作用在结点上的荷载等于零时，即P d =0时，由式（1-46）可知P =P e

在“后处理法”中，P、Pd 和P e 应分别与自由解点位移相对应，其表达式与式（1-46）相同。

e

e

e