Subroutine Elem_Stiff(···)
说明
Stiff=0.0 ! 单元刚度清零
Select Case (Type)
Case (1)
平面杆系结构单元
Case (2)
空间杆系结构单元
Case Default
出错信息
End Select
End Subroutine Elem_Stiff
3.4 杆系结构单元分析子程序
3.4.1 单元刚度总体设计 3.4.2 说明部分设计
Integer,Intent (in),,· · ·
入口整型参数
Real(8),Intent(in),,· · ·
入口实型参数
Real(8),Intent(out),,· · ·
出口实型参数
Real(8),,Work1,· · ·
Integer,,i,j,k,· · ·
实型和整型工作变量
3.4 杆系结构单元分析子程序
3.4.3 平面杆系结构设计
Select Case (Plane)
Case (1)
平面桁架元素赋值
Case (2)
平面梁柱元素赋值
Case (3)
· · · · · ·
Case Default
出错信息
End Select
3.4.4 空间杆系结构设计
Select Case (Space)
Case (1)
空间桁架元素赋值
Case (2)
空间梁柱元素赋值
Case (3)
交叉梁元素赋值
Case Default
出错信息
End Select
3.4 杆系结构单元分析子程序
3.4.5 有关单元等效结点荷载设计和进一步的考虑
1) 单元等效结点荷载设计同仿单元刚度。
2) 从各类单元刚度元素的计算，可看到要用到长度、
单元弹性特性、单元截面特性等数据。因此，要
确定存放它们的数据结构。要将它们作为出口。
3) 为计算单元等效结点荷载元素，首先要建立各种
荷载情况等效荷载表达式，它们可由积分或载常
数表得到。然后要解决荷载信息的存放结构，也
要将它们作为出口量。
4) 单元刚度矩阵、等效结点荷载矩阵都应先清零。
4.1 杆系结构整体分析
首先就全刚结点平面刚架进行讨论，然后推广。
4.1.1 总的思路
在单元特性搞清后，将单元拼装回去。在结点处
位移自动协调基础上，如果全部结点平衡，则求得
的结点位移将是实际结构的解。因此，整体分析就
是设法建立结点平衡方程。
4.1.2 坐标转换
?
x
y
y x
xF
yF
xF
yF
yd
xd
yd
xd
组成结构的杆件可以各个
方向，单元分析对局部坐标，
因此，必须将物理量转为统
一坐标 —— 整体坐标。
1) 力的转换关系
?
?
?
?
?
?
?
?
?
?
?
? ?
?
?
?
?
?
?
?
y
x
y
x
c o ss i n
s i nc o s
F
F
F
F
??
??
4.1 杆系结构整体分析
2) 位移转换关系
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
y
x
y
x
c o ss i n-
s i nc o s
d
d
d
d
??
??
3) 转换矩阵
? ?
?
?
?
?
?
?
?
?
?
?
??
100
0
0
CS
SC
? ? ?
? ? ? ?
? ? ? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
0
0
s i n
c o s
T
S
C
转换矩阵是正交矩阵。
?
x
y
y x
xF
yF
xF
yF
yd
xd
yd
xd
? ? ? ?T1 TT ??
4.1 杆系结构整体分析
4) 杆端力转换
? ? ? ? ? ? ? ? ? ? )(T eeee PFTPF ???
5) 杆端位移转换
? ? ? ?? ?ee dTd ?
6) 刚度方程的转换
?
x
y
y x
xF
yF
xF
yF
yd
xd
yd
xd
? ? ? ? ? ? ? ? ? ?
? ?? ? ? ? ? ?? ? ? ? ? ? eeee
eeee
dTkTdkT
PFTPF
T
)(
??
????
? ? ? ? ? ? ? ?eeee dkPF ???如果记 称为整体单元刚度矩阵? ? ? ? ? ? ? ?TkTk ee T?
则
? ? ? ? ? ? ? ?eeee dkPF ??
这就是整体坐标下的单元刚度方程。
本节以后的讨论认为
都是对整体坐标的
4.1 杆系结构整体分析
4.1.3 结点平衡方程的建立
1) 一简单例子（如图）
图中有两套编号,红 的
是单元杆端编号,黑 的是
结构整体编号。
1-1) 结点示意
1
2 3
41
2
1 2 2
1
①
②
③
图中蓝色的表示结点荷载（已知）,红 色的表示
杆端力（未知的）,,分别 1,2单元杆端力
子矩阵。对 1,4结点“荷载”含有未知反力。
? ?2F ? ?1F
2
? ?2dP
? ?12F
? ?21F
1-2) 结点平衡
21 orj ?
由示意图可见，结构
结点的平衡方程为 ? ? ? ???
各杆交 i
ji FP
4.1 杆系结构整体分析
从例图可见，其全部结点
平衡方程为
1
2 3
41
2
1 2 2
1
①
②
③
若记
? ? ? ?
? ? ? ? ? ?
? ? ? ? ? ?
? ? ? ?
314d
32223d
21122d
111d
FP
FFP
FFP
FP
?
??
??
?
? ? ? ? ? ? ? ? ? ?? ?TT4dT3dT2dT1dd PPPPP ?
? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? TT32T31T22T21T12T11 FFFFFFF ?
2
? ?2dP
? ?12F
? ?21F
? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ?F
I
II
II
I
P
?
?
?
?
?
?
?
?
?
?
?
?
?
00000
0000
0000
00000
d
4.1 杆系结构整体分析
式中 [I],[0]分别为单位和零矩阵。
? ? ? ?? ?FAP ?d
1
2
3
若引入矩阵记号,则结点平衡方程可改写作
这一结论虽然是由一个例子得到的，但是显然
对一切结构都是成立的。问题在于不同结构, [A]
矩阵是不同的。
? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
00000
0000
0000
00000
I
II
II
I
A
4.1 杆系结构整体分析
4.1.4 杆端位移用结点位移来表示
1
2 3
41
2
1 2 2
1
①
②
③仍以简单例子来说明若记
? ? ? ? ? ?? ?TT4T1 ??? ????
? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? TT32T31T22T21T12T11 dddddd??
由结点、杆端位移的协调条件，可得 [?],[? ]
的对应关系为
? ? ? ? ? ??? TA?
式中 [A]T是前面力关系 [A]的转臵，因此 [A]T称
为 位移转换矩阵 。
4.1 杆系结构整体分析
4.1.5 整体刚度方程 —— 结点平衡
1
2 3
41
2
1 2 2
1
①
②
③若记
? ? ? ? ? ?? ?31 kkd i a gK ????
? ? ? ?? ? ? ? ? ?? ? ? ?eqd PKF;FAP ??? ??
? ? ? ? ? ? ? ?? ?TT3T2T1eq PPPP ?
引入位移转换关系，则
? ? ? ? ? ?? ? ? ?)( eqd PKAP ??? ?
? ? ? ?? ?eqeq PAP ?记
? ? ? ? ? ?? ?? ? ? ??Teqd AKAPP ??
这就是整体刚度方程，它的物理实质是结点平
衡。 [K]称作结构刚度矩阵（或整体刚度矩阵），
[R]称作综合等效结点荷载矩阵，它由两部分组成。
? ? ? ? ? ? ? ? ? ?? ? ? ?KAKA;RPP ??? Teqd记
? ? ? ?? ??KR ?则
单元个数
? ?? ?? ? ? ? ? ? ? ??
?
?
3
1
TT
i
iii AkAAKA
321
4.1 杆系结构整体分析
4.1.6 整体刚度矩阵的建立
1
2 3
41
2
1 2 2
1
①
②
③若将 [A]按单元分成图示
三个子矩阵
? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
00000
0000
0000
00000
I
II
II
I
A
则
由此可见，整体刚度矩阵可由各单元整体刚度
矩阵 装配累加 得到。为说明如何装配，先将单元
刚度矩阵进行分割
? ? ? ? ? ?? ? ? ?
i
i kk
kk
k ?
?
?
?
?
?
?
2221
1211
r
s
sr
整体结点码
则由矩阵乘法可证明,[A]i[k]i[A]iT的结果是，将
刚度矩阵子矩阵按整体结点码 r, s 送整体刚度矩阵
相应位臵。这一装配规则称为“对号入座”。
4.1 杆系结构整体分析
1) 任意结构情况
上面结论是通过具体例子（全刚结点平面刚架）
得到的，由虚位移原理或势能原理进行整体分析
（见讲义），可得任意结构其结论同此例。
2) 结点位移编号
如果按结点顺序，对结点
非零位移进行依次编号，这一
序号称作结点位移码。为便于
计算机处理并减少结构刚度矩
阵的阶次，将零位移的号码变
为零。
对图示三铰刚架，当仅
用一种单元（梁柱自由是
单元）时结点位移编号如
图所示。
3) 单元定位向量
按单元局部结点码顺序，
将结点位移码排成的向量，称作单元的 定位向量 。
(0,0,1) (0,0,2)
(3,4,5) (6,7,8) (6,7,9)
( 1 0,1 1,1 2 )
1 2
3 4 5
6
① ②
③ ④
(0,0,1) (0,0,2)
(3,4,5) (6,7,8) (6,7,9)
( 1 0,1 1,1 2 )
1 2
3 4 5
6
① ②
③ ④
4.1 杆系结构整体分析
对图示刚架各单元的定位
向量为
① (0,0,1,3,4,5)
② (0,0,2,10,11,12)
③ (3,4,5,6,7,8)
④ (6,7,9,10,11,12)
① 2 3)
②,6,7,8)
③ 1 2 3 4 5)
④ 4 5 6 7,8)
4) 按单元定位向量集装刚度矩阵和综合荷载
前面说明的是分块子矩阵集装,下面说明如何按定
为向量来集装,
如果如图所是采用各种不同的单元 （一端有铰
）,则定位向量为
(0,0,1) (0,0,2)
(3,4,5) (6,7,8) (6,7,9)
( 1 0,1 1,1 2 )
1 2
3 4 5
6
① ②
③ ④
(0,0) (0,0)
(1,2,3) (4,5) (6,7,8)
1 2
3 4 5
① ②
③ ④
如何获得
带铰的单元刚
度矩阵和等
效荷载矩阵
定位向量
(0,0,1) (0,0,2)
(3,4,5) (6,7,8) (6,7,9)
( 1 0,1 1,1 2 )
1 2
3 4 5
6
① ②
③ ④
(0,0) (0,0)
(1,2,3) (4,5) (6,7,8)
1 2
3 4 5
① ②
③ ④
4.1 杆系结构整体分析
4-1) 刚度集装
以 3 单元为例来说明
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
5554535251
4544434241
3534333231
2524232221
1514131211
3
k
87654
定位向量
单元局部位移码
54321
5
4
3
2
1
8
7
6
5
4
根据单元
局部位移码
和定位向量
的对应关系
用定位向量
位移码送元
素。
根据单元局部位移码和定位向量的对应关系用定
位向量位移码送元素，定位向量元素为零时不送。
(0,0,1) (0,0,2)
(3,4,5) (6,7,8) (6,7,9)
( 1 0,1 1,1 2 )
1 2
3 4 5
6
① ②
③ ④
(0,0) (0,0)
(1,2,3) (4,5) (6,7,8)
1 2
3 4 5
① ②
③ ④
4.1 杆系结构整体分析
4-2) 荷载集装
以 4 单元为例来说明
定位向量
8
7
6
0
0
? ?
45
4
3
2
1
4
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
P
P
P
P
P
P
5
4
3
2
1
局
部
位
移
码 此结论同样适用于刚度集装
4.1 杆系结构整体分析
4.1.7 整体分析总结
1) 对局部坐标和整体坐标不一致的单元，要对刚
度、荷载进行坐标转换。
2) 需对“结构”进行结点、位移的局部和整体编
号。
4) 集装所得整体刚度矩阵是对称、带状稀疏矩阵，
当支撑条件能限制刚体位移时，矩阵非奇异。
3) 根据单元局部位移码和定位向量的对应关系用
定位向量位移码送元素，定位向量元素为零时
不送。据此可集装、累加得到整体刚度矩阵。
5) 综合荷载由两部分组成，因此首先要将直接作
用结点的荷载按结点位移码送入，如果还有单
元等效荷载，再按定位向量集装、累加。
4.1 杆系结构整体分析
8) 如果有某位移码方向弹性支撑，需进行将弹簧
刚度送入位移码对应的对角线元素位臵累加。
9) 如果有某位移码方向已知支撑位移，需进行将
“边界条件处理”。具体做法以后介绍。
7) 整体刚度方程实质是全部结点的平衡条件。
6) 刚度矩阵带状稀疏，其带宽取决于结点、位移
编码。
最大半带宽 =定位向量中最大元素差 +1。
4.1 杆系结构整体分析
2.5.8 边界条件的处理
10) 当用虚位移或势能原理作整体分析时（势能为
例），应变能为单元应变能之和，外力势能为
单元外力势能之和 +结点外力势能。全部杆端力
的外力势能彼此抵消。
1) 乘大数法
设第 i 个位移为已知值 a, N =108 或更大的数。
乘大数法是 将刚度矩阵 Kii改为 N?Kii，将 Ri改为
N?a。
请考虑
为什麽这样做能使
边界条件得到满足？
2) 臵换法（划零臵 1）
设第 i 个位移为已知值 a 。
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
n
i
n
i
nnnjninn
inijiiii
nji
nji
R
R
R
R
KKKKK
KKKKK
KKKKK
KKKKK
2
1
2
1
21
21
2222221
1111211
?
?
?
?
4.1 杆系结构整体分析
上述臵换工作量大一些，显然可看出边界条件得
到精确满足。
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
aKR
a
aKR
aKR
KKKK
KKKK
KKKK
nin
i
i
n
i
nnnjnn
nj
nj
22
11
2
1
21
222221
111211
0
00100
0
0
?
?
?
?
4.1 杆系结构整体分析
★ 3) 关于斜边界的处理
如图示意的斜支座情况，有多种处理方案。
3-1) 通过单元的坐标转换来处理
?
x
y
?
xy
图示有斜支座单元,r 结点处
以倾角 ?- ? 来进行坐标转换，也
即在 r 结点处整体坐标为图示 xy 。
r
3-2) 通过增加一个单元来处理
图示有斜支座单元,r 结点处沿 y 方向增加一个
刚结的单元，此单元有“无穷大”的抗拉刚度、但
没有抗弯刚度。单元长度可任意。
3-3) 对整体刚度矩阵进行处理（参见匡文起教材）
最大
半带宽
2.5 杆系结构整体分析
2.5.10 总刚度矩阵的存储与对应解法
因为总刚度矩阵对称、带状稀疏，利用这一特
点可减少存储、提高解算效率。
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nnnin
iniii
ni
KKK
KKK
KKK
1
1
1111
零元素
零元素非零元素
最大
半带宽
主对角线元素
等
半
带
存
储
零元素
1
3
6
10
14
53
60
66
72
非
零
元
素
变
带
宽
一
维
存
储
带
宽
是
变
的
到 P.55
2.5 杆系结构整体分析
目前一般都用变带宽存储，下面结合程序说明
存储和解法。首先介绍一些 F 90的语法。
定义导出类型
导出类型 —— 结点
type,,typ_Joint
real,,x,y ！ 坐标
integer,,GDOF(3)！ 整体位移码
end type typ_Joint
1) 有关 F90 语法
?导出类型 新特性
用结点导出类型作为成员导出单元类型,
type,,typ_Element
integer,,JointNo(2) ！ 结点编号
type(typ_Joint),,Node(2) ！ 结点
integer,,GlbDOF(6) ! 定位向量
real,,EA,EI
end type typ_Element
type (typ_Element),,Elem(5) ！ 定义 5个单元类型
…
! 对单元 i的端点 j的 x,y,GDOF(1:3)的赋值
Elem(i)%Node(j) = Joint(Elem(i)%JointNo(j))
…
由导出类型定义新类型
由导出类型定义变量
real,,A(5),B(5,10),C(5)
B=0.0 !对 B清零
A=1.0 !对 A赋 1,A(i)=1.0,i=1,5
C=A+2 !数组与标量运算, A(1:5)+(/ 2,2,2,2,2 /)
A=C+A !数组与数组运算 （ 同形 ）
C=sqrt(A) !数组的函数运算, C(i)=sqrt(A(i),i=1,5
数组内部函数,
dot_product(vector_a,vector_b) !点积
如,dot_product((/1,2,3/),(/2,3,4/)) 的值为 20 （ 待续 ）
有关 F90 语法
数组运算与赋值,
matmul(matrix_a,matrix_b) 矩阵相乘
如,locEDisp = matmul(T,glbEDisp)
transpose(matrix) 矩阵转置
如,glbEDisp = matmul(transpose(T),locEDisp)
size(array,dim) 求数组第 dim维的长度
dim为可选变元, size(a,dim=2)
若 array为一维时, 可不用 dim。
sum(array,dim,mask) 数组元素求和
dim,mask为可选变元； mask=条件表达式
sum(a(1:10)) 对 a的 1到 10元素求和
sum(a(1:10),mask=a>0) 对 a(1:10)中大于 0的元素求和
? ? ? ?? ?ee dTd ?
（续）
有关 F90 语法
?where结构 新特性
例
where (C>0)
C = 0
A = B*D
end where
where (C>0)
A = B
end where
定义
where (数组关系表达式 )
数组赋值语句
…
elsewhere
数组赋值语句
...
end where
规则, 1)同形数组 ; 2)不许嵌套 ; 3)最多二分叉
有关 F90 语法
?cycle和 exit语句 新特性
用在 do循环中
cycle —— 作下一个循环步
exit —— 跳出循环，执行 end do后一条语句
等效例
do
...
if(.not.cond1)then
...
if(cond2) goto 5
...
end if
end do
5,..
用法
do
...
if (cond1) cycle
...
if (cond2) exit
...
end do
...
有关 F90 语法
?数组构造函数 spread
语法
spread(数组名,dim,ncopies)
将数组沿 dim维方向复制 ncopies形成新数组
dim,ncopies — 整型, 位置变元, 关键字变元
（ 若按位置引用, 可略关键字 ）
例, （仅限一维数组）
1） spread(one,dim=1,ncopies=3)
spread(one,1,3)
spread(one,ncopies=3,dim=1)
[1,1,1]
或
[1,1,1]T
2） ELocVec(1:6)=(/ 1,0,3,4,5,0 /)
spread(ELocVec,dim=1,ncopies=3)
?
?
?
?
?
?
?
?
?
?
054301
054301
054301
3）
spread(A(2,2:),dim=1,ncopies=2)
?
?
?
?
?
?
?
?
?
?
?
332313
322212
312111
A
...
...
...
??
?
??
?
3222
3222
..
..
如果 dim=2呢?
有关 F90 语法
?指针 pointer
pointer是变量的属性, 可以指向相同类型的变量 ;
被指向的目标必须具有 target属性或 pointer属性
可以将指针变量理解为别名, 称号
real,target,,a,b,EDisp(6) ! 可被指针所指
real,pointer,,p1,p2 ! 称号：班长, 课代表
！ p1,p2是指针, 可以指向实型数据
real,pointer,,G(:) ! 先进集体
！ G是指针, 可以指向一维实型数组
? 指针是一种, 称号,, 上述声明语句建立了, 称
号,, 但并未, 授予, 某个变量这个称号, 因此
是指向, 空,, 并未占用内存
a = 3.0
p1 => a ! p1指向 a
! 称号 p1授予 a,a的数据有两个名：固定名 a和流动
名 p1; 既可用 p1也可用 a（ p1 — 班长, a — 张三 ）
a = 4.0 ! a的值变为 4.0,p1也变为 4.0
p1 => b ! 班长换人了
G => EDisp ! 先进集体有了得主
！ EDisp(:)的长度就是 G(:)的长度, 用 G和用 EDisp
同样效果
又如：
real,target,,a,b
real,pointer:,p,q
a = 3.14
b = 2.0
p => a ! p = a = 3.14
q => b ! q = b = 2.0
q = p
! ( q指向的数据 b ) = ( p指向的数据 a )
！ 即,所有 = 3.14
? 指针可以指向有名的数据区, 也可以指向无名的
数据区
real,pointer:,b(:,:)
integer,,n
read(*,*) n
allocate (b(n,n)) !指针指向了一个刚开辟的数组
! 以下可以当作常规数组用
b(1,1) = 1.1
b(1,2) = 1.2
...
deallocate (b)
有关 F90 语法
?用指针建立动态数组
指针与 allocatable数组的区别
? 具备 allocatable数组的所有功能
? 还可以用在导出类型中, 例如整体刚度矩阵的变带宽存
储：
type,,typ_Kcol !整体刚度矩阵 K的列
real(8),pointer,,row(:) !该列的行元素
end type
...
type (typ_Kcol),allocatable,,Kcol(:)
…
allocate (Kcol(NGlbDOF)) ! 分配了 NGlbDOF列
...
allocate (Kcol(5)%row(3:5)) ! 第 5列只用 3至 5行
(1) NElem,NJoint,NGlbDOF,NJLoad,NELoad
单元数 结点数 总自由度数 结点荷载数 单元荷载数
[Joint - 结点 ] … NJoint行
(2) Joint%X,Joint%Y,GDOF(1:3)
X 坐标 Y坐标 结点位移码
[Elem - 单元 ] … NElem行
(3) JointNo1,JointNo2,EA,EI
起点号 终点号 刚度
[JLoad - 结点荷载 ] … NJLoad行
(4) JointNo,LodDOF,LodVal
作用点号 局部位移码 荷载值
[ELoad - 单元荷载 ] … NELoad行
(5) ElemNo,Indx,Pos,LodVal
单元号 类型号 位置 荷载值
Indx 类型 pos
1 均布 长度
2 集中 位置
3,..
2) 某程序输入数据说明
3,5,7,1,1
0,0,0,0,0
0,4,1,2,3
4,4,4,5,6
4,4,4,5,7
4,0,0,0,0
1,2,1.0E9,16
2,3,1.0E9,24
5,4,1.0E9,12
2,1,10.0E3 ! 结点 2，自由度 1，值为 10E3
2,1,4,-4.0E3 ！单元 2，均布，长 4米，值为 -4E3
2-1) 数据文件例子：
(2)
(1) (3)
2 4
1
3
5
i = 6
i = 4 i = 3
10 kN
4 kN/m
4m
4m
EA= 109 N
(1)
(2)
坐标
位移码
(3)
(4)
(5)
结点号 EA,EI
read(5,*) NElem,NJoint,NGlbDOF,NJLoad,NELoad
allocate (Joint(NJoint))
allocate (Elem(NElem))
allocate (JLoad(NJLoad))
allocate (ELoad(NELoad))
...
read(5,*) (Joint(i),i=1,NJoint)
read(5,*) (Elem(ie)%JointNo,Elem(ie)%EA,&
ELem(ie)%EI,ie=1,NElem)
if (NJLoad>0) read(5,*) (JLoad(i),i=1,NJLoad)
if (NELoad>0) read(5,*) (ELoad(i),i=1,NELoad)
2-2) 程序读入计算所需数据：
2-3) 求始行码和分配带宽子程序
!==================================subroutine SetMatBand (Kcol,Elem) ! 接口简单
!==================================
type(typ_Kcol),intent(in out),,Kcol(:) ！总刚列
type(typ_Element),intent(in)),,Elem(:) ！单元
integer(ikind),,minDOF,ELocVec(6)
integer(ikind),,Row1(size(Kcol,dim=1))
！ Row1为自动数组，子程序结束后自动释放。
！这样做可使接口简单，减少了数组的控制变量。
integer(ikind),,ie,j,NGlbDOF,NElem
NGlbDOF = size(Kcol,dim=1) ！使接口简单
NElem = size(Elem,dim=1)
Row1 = NGlbDOF ! 先设所有始行码同终行码
! 确定各列始行码向量
do ie=1,Nelem ! 对单元循环
! 确定定位向量
ELocVec(:) = Elem(ie)%GlbDOF(:)
! 寻找定位向量中大于零的最小值
minDOF = minval(ELocVec,mask=ELocVec>0)
! 屏蔽定位向量中小于等于零的
where (ELocVec > 0)
! 寻找 Row1(ELocVec)和 minDOF中的最小值
Row1(ELocVec) = min(Row1(ELocVec),&
minDOF)
end where
end do
! 为各列的带宽分配空间
do j=1,NGlbDOF
! 对总自由度数循环
allocate (Kcol(j)%row(Row1(j):j)
! 给 Kcol(j)%row 从 Row1(j) 到 j 个空间
Kcol(j)%row = zero
! 总刚元素清零
end do
return
end subroutine SetMatBand
621,.,,,,j,i ?
IJ)ie(ij Kk ? ??
??
?
?
?
??
??
0)(
0)(
)ie(
j
)ie(
i
J
I
?
?
3) 整体刚度矩阵的集成
do ie=1,NElem
… ! 计算单刚 EK(6,6),ELocVec(6)
do j=1,6 ! 对单元逐列集成
JGDOF = ELocVec(j) ! 取出位移码
if (JGDOF = 0) cycle ! 作下一循环步
where (ELocVec > 0, And, ELocVec <= JGDOF)
! 位移码非零同时小于第 j个局部码对应的位移码
Kcol(JGDOF)%row(ELocVec) = &
Kcol(JGDOF)%row(ELocVec) + EK(:,j) ! 集成
end where
end do
end do
局部码 位移码
4) 变带宽矩阵的分解求解
4-1) LDLT分解法求解 [A]{x} = {b}
若 ? ?
nn?A
对称正定，则可分解为
? ? ? ?? ?? ?TLDLA ? =
单位上三角阵
对角阵
原方程变为
求解步骤：
1,分解：
2,解 y,
3,解 x,
LDLT分解法 Gauss消去
? ?? ?? ? ? ? ? ?? ? ? ?byLxLDL ??T
? ? ? ?? ?? ?TLDLA ?
? ?? ? ? ?byL ?
? ?? ? ? ? ? ?yxLD ?T
前消去
处理右端项
回代
(不同的 b只做一次分解 )
? ? ? ?? ?? ? ? ?? ?ULLDLA ?? T
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nnnn
ijii
nj
nj
aaa
aaa
aaaa
aaaa
??
???
?
????
??
??
21
21
222221
111211
存放：
主对角 —
上三角 —
iid
jiij ll ?T
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1
1
1
1
21
3231
21
0
??
?
nn
ij
ll
l
ll
l
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nn
jiii
n
n
d
ld
ldldd
ldldldd
?
?
?
0
222322222
1113111211111
jiiiij ldu ?
ija
ijl
? ?
? ?
??
n
k
n
k
jkkkikkjikij ldlula
1 1
?
?
?
?
?
?
?
?
?
?
ikl
ik
ik
l
ik
ik 1
0
j,ik ?? 时不必求和
? ?
? ?
??
)j,i(
k
i
k
kkjkikij dlla
m i n
1 1
)(?
(上三角, i < j )
?
1111 da ? 1111 ad ?
211112 lda ?
121121 adl ?
222211122 dlda ?? 221112222 ldad ??
不动
存在 处
12a
存在 处
22a
(第 j 列系数 )
iiji
i
k
kkjkikij dldlla ?? ?
?
?
1
1
上三角, i < j
4-2) 分解一般情况,逐列分解
ii
i
k
kkjkikijji ddllal ??
?
?
???
? ?? ??
?
1
1
121 ?? j,,,i ?
对角线, i = j > 1
jj
j
k
kkjkjj ddla ?? ?
?
?
1
1
2
n.,,,,,j 32?
?
?
?
??
1
1
2
j
k
kkjkjjjj dlad
(第 j 列系数 )
第 i 列中第 1个非零元素的行码为：
4-3) 变列宽存贮的分解修正：
ii
i
ik
kkjkikijji ddllKl ?
?
?
?
???
?
?? ?
?
?
1
1
121 ?? j,,,ii ?
11 ?? ik
?
?
?
??
1
2
1
j
ik
kkjkjjjj dlKd
第 j 列中第 1个非零元素的行码为：
ii1
ji1
则
)m a x ( 111 ji i,ii ?
n.,,,,,j 32?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
77
66
575655
465444
36534333
2752423222
413111
0
0000
K
K
KKd
Kld
Klld
Kllld
lld
分解顺序
4-4) F90 实现
! 三角分解
diag(1:ncol) = (/ (Kcol(j)%row(j),j=1,ncol) /)
do j = 2,ncol
row1 = lbound(Kcol(j)%row,1) ! i_1 j
do i = row1,j-1
row_1 = max(row1,lbound(Kcol(i)%row,1)) ! i_1
s = sum(diag(row_1:i-1) *Kcol(i)%row(row_1:i-1) * &
Kcol(j)%row(row_1:i-1)) ! 求和部分
Kcol(j)%row(i) = (Kcol(j)%row(i)-s)/diag(i)
end do ！第 j列系数分解完毕
s = sum(diag(row1:j-1)*Kcol(j)%row(row1:j-1)**2)
diag(j) = diag(j) - s ! 第 j列的主对角元素
end do
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nn
nn
ij
b
b
b
y
y
y
ll
l
ll
l
??
??
?
2
1
2
1
21
3231
21
1
1
1
1
04-5) 一般情况解 y 公式
22121 byyl ??
iik
i
k
ik byyl ???
?
?
1
1
k
i
k
ikii ylby ?
?
?
??
1
1
n,,,i ?21?
12122 ylby ??
11 by ?
可不动
可存在 处2b
4-6) 变列宽存储解 y的修正
k
i
ik
ikii ylby ?
?
?
??
1
1
上三角的第 i 列
从第 行元素开始1i
! 回代步骤 1,自上而下
do i = 2,ncol
row1 = lbound(Kcol(i)%row,dim = 1)
GP(i) = GP(i) - sum( Kcol(i)%row(row1:i-1) &
* GP(row1:i-1) )
end do
求出后,不再用，可将 存在 处iy
iyib ib
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nnnn
n
n
dy
dy
dy
x
x
x
ll
lll
??
??
?
222
111
2
1
232
13121
1
1
1
1
0
4-7) 一般解 x 的 公式
nnnn dyx ?
??? 1nx
自下而上：
? ? ? ? ? ? ? ?yDxL 1T ??
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nnd
d
d
D
1
1
1
22
11
1
?
iii
n
ik
kkii dyxlx ?? ?
?? 1
?
??
??
n
ik
kkiiiii xldyx
1
?
这样做的缺欠：
nni)n(i)n(i xlyy ?? ? 1
自下而上逐行计算，行中遇到 0元素需跳过，不方便！
4-8) 改为自右向左逐列计算！
记
iiini dyy ?? )1(
nx
解出
)1( ?? n
nn yx
)1(
11
?
?? ?
n
nn yx
对 y 向量预处理
再解出
?
第 n 列乘上 后移到
右边去，修正 y 向量nx
1?nx
??
?
?
?
???
????
??
?
21
1
)(
11
1
)1()(
,,n,njyx
j,,iixlyy
j
jj
jji
j
i
j
i
?
?
! 回代步骤 2,自右向左
GP(:) = GP(:)/diag(:)
do j = ncol,2,-1
row1 = lbound(Kcol(j)%row,dim = 1)
GP(row1:j-1) = GP(row1:j-1) - GP(j) * &
Kcol(j)%row(row1:j-1)
end do
! 现在 GP就是解
4-9) 小结：
? 无须一维存贮，无须行列码转换定位
? 公式与程序直接对应翻译，直截了当
? 求和采用内部函数进行数组运算
? 动态内存，用多大、开多大
? 数据封装性好，接口简单,
subroutine SetMatBand (Kcol,Elem)
subroutine VarBandSolv (Disp,Kcol,GP)
? 通用性强
2.6 杆系结构内力计算
2.6.1 杆端内力计算公式
解方程的结果可以得到结点位移，有了位移就
可以进一步求各单元的内力。
解算步骤：
? ? ? ?? ? ? ?PdkF ??
? ?? ?? ? ? ?PdTk ??
从 [?]根据定位向量取出单元杆端位移。
由单元倾角确定是算 还是 。? ?? ?? ?dTk ? ?? ?dk
减去等效结点荷载或加上固端力矩阵。
2.6 杆系结构内力计算
2.6.2 杆中任意截面内力计算公式
注意事项：
由图示隔离体图，列杆段的平衡方程即可得到
任意截面的内力计算公式。 请自行写出。
按上述隔离体图所求得的内力是精确的。
M
N
Q
1F
2F3F
求得单元结点位移后，代入单元位移场、求应
变、求内力，这样所得的内力一般不满足平衡条
件。只是近似解。
说明
Stiff=0.0 ! 单元刚度清零
Select Case (Type)
Case (1)
平面杆系结构单元
Case (2)
空间杆系结构单元
Case Default
出错信息
End Select
End Subroutine Elem_Stiff
3.4 杆系结构单元分析子程序
3.4.1 单元刚度总体设计 3.4.2 说明部分设计
Integer,Intent (in),,· · ·
入口整型参数
Real(8),Intent(in),,· · ·
入口实型参数
Real(8),Intent(out),,· · ·
出口实型参数
Real(8),,Work1,· · ·
Integer,,i,j,k,· · ·
实型和整型工作变量
3.4 杆系结构单元分析子程序
3.4.3 平面杆系结构设计
Select Case (Plane)
Case (1)
平面桁架元素赋值
Case (2)
平面梁柱元素赋值
Case (3)
· · · · · ·
Case Default
出错信息
End Select
3.4.4 空间杆系结构设计
Select Case (Space)
Case (1)
空间桁架元素赋值
Case (2)
空间梁柱元素赋值
Case (3)
交叉梁元素赋值
Case Default
出错信息
End Select
3.4 杆系结构单元分析子程序
3.4.5 有关单元等效结点荷载设计和进一步的考虑
1) 单元等效结点荷载设计同仿单元刚度。
2) 从各类单元刚度元素的计算，可看到要用到长度、
单元弹性特性、单元截面特性等数据。因此，要
确定存放它们的数据结构。要将它们作为出口。
3) 为计算单元等效结点荷载元素，首先要建立各种
荷载情况等效荷载表达式，它们可由积分或载常
数表得到。然后要解决荷载信息的存放结构，也
要将它们作为出口量。
4) 单元刚度矩阵、等效结点荷载矩阵都应先清零。
4.1 杆系结构整体分析
首先就全刚结点平面刚架进行讨论，然后推广。
4.1.1 总的思路
在单元特性搞清后，将单元拼装回去。在结点处
位移自动协调基础上，如果全部结点平衡，则求得
的结点位移将是实际结构的解。因此，整体分析就
是设法建立结点平衡方程。
4.1.2 坐标转换
?
x
y
y x
xF
yF
xF
yF
yd
xd
yd
xd
组成结构的杆件可以各个
方向，单元分析对局部坐标，
因此，必须将物理量转为统
一坐标 —— 整体坐标。
1) 力的转换关系
?
?
?
?
?
?
?
?
?
?
?
? ?
?
?
?
?
?
?
?
y
x
y
x
c o ss i n
s i nc o s
F
F
F
F
??
??
4.1 杆系结构整体分析
2) 位移转换关系
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
y
x
y
x
c o ss i n-
s i nc o s
d
d
d
d
??
??
3) 转换矩阵
? ?
?
?
?
?
?
?
?
?
?
?
??
100
0
0
CS
SC
? ? ?
? ? ? ?
? ? ? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
0
0
s i n
c o s
T
S
C
转换矩阵是正交矩阵。
?
x
y
y x
xF
yF
xF
yF
yd
xd
yd
xd
? ? ? ?T1 TT ??
4.1 杆系结构整体分析
4) 杆端力转换
? ? ? ? ? ? ? ? ? ? )(T eeee PFTPF ???
5) 杆端位移转换
? ? ? ?? ?ee dTd ?
6) 刚度方程的转换
?
x
y
y x
xF
yF
xF
yF
yd
xd
yd
xd
? ? ? ? ? ? ? ? ? ?
? ?? ? ? ? ? ?? ? ? ? ? ? eeee
eeee
dTkTdkT
PFTPF
T
)(
??
????
? ? ? ? ? ? ? ?eeee dkPF ???如果记 称为整体单元刚度矩阵? ? ? ? ? ? ? ?TkTk ee T?
则
? ? ? ? ? ? ? ?eeee dkPF ??
这就是整体坐标下的单元刚度方程。
本节以后的讨论认为
都是对整体坐标的
4.1 杆系结构整体分析
4.1.3 结点平衡方程的建立
1) 一简单例子（如图）
图中有两套编号,红 的
是单元杆端编号,黑 的是
结构整体编号。
1-1) 结点示意
1
2 3
41
2
1 2 2
1
①
②
③
图中蓝色的表示结点荷载（已知）,红 色的表示
杆端力（未知的）,,分别 1,2单元杆端力
子矩阵。对 1,4结点“荷载”含有未知反力。
? ?2F ? ?1F
2
? ?2dP
? ?12F
? ?21F
1-2) 结点平衡
21 orj ?
由示意图可见，结构
结点的平衡方程为 ? ? ? ???
各杆交 i
ji FP
4.1 杆系结构整体分析
从例图可见，其全部结点
平衡方程为
1
2 3
41
2
1 2 2
1
①
②
③
若记
? ? ? ?
? ? ? ? ? ?
? ? ? ? ? ?
? ? ? ?
314d
32223d
21122d
111d
FP
FFP
FFP
FP
?
??
??
?
? ? ? ? ? ? ? ? ? ?? ?TT4dT3dT2dT1dd PPPPP ?
? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? TT32T31T22T21T12T11 FFFFFFF ?
2
? ?2dP
? ?12F
? ?21F
? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ?F
I
II
II
I
P
?
?
?
?
?
?
?
?
?
?
?
?
?
00000
0000
0000
00000
d
4.1 杆系结构整体分析
式中 [I],[0]分别为单位和零矩阵。
? ? ? ?? ?FAP ?d
1
2
3
若引入矩阵记号,则结点平衡方程可改写作
这一结论虽然是由一个例子得到的，但是显然
对一切结构都是成立的。问题在于不同结构, [A]
矩阵是不同的。
? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
00000
0000
0000
00000
I
II
II
I
A
4.1 杆系结构整体分析
4.1.4 杆端位移用结点位移来表示
1
2 3
41
2
1 2 2
1
①
②
③仍以简单例子来说明若记
? ? ? ? ? ?? ?TT4T1 ??? ????
? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? TT32T31T22T21T12T11 dddddd??
由结点、杆端位移的协调条件，可得 [?],[? ]
的对应关系为
? ? ? ? ? ??? TA?
式中 [A]T是前面力关系 [A]的转臵，因此 [A]T称
为 位移转换矩阵 。
4.1 杆系结构整体分析
4.1.5 整体刚度方程 —— 结点平衡
1
2 3
41
2
1 2 2
1
①
②
③若记
? ? ? ? ? ?? ?31 kkd i a gK ????
? ? ? ?? ? ? ? ? ?? ? ? ?eqd PKF;FAP ??? ??
? ? ? ? ? ? ? ?? ?TT3T2T1eq PPPP ?
引入位移转换关系，则
? ? ? ? ? ?? ? ? ?)( eqd PKAP ??? ?
? ? ? ?? ?eqeq PAP ?记
? ? ? ? ? ?? ?? ? ? ??Teqd AKAPP ??
这就是整体刚度方程，它的物理实质是结点平
衡。 [K]称作结构刚度矩阵（或整体刚度矩阵），
[R]称作综合等效结点荷载矩阵，它由两部分组成。
? ? ? ? ? ? ? ? ? ?? ? ? ?KAKA;RPP ??? Teqd记
? ? ? ?? ??KR ?则
单元个数
? ?? ?? ? ? ? ? ? ? ??
?
?
3
1
TT
i
iii AkAAKA
321
4.1 杆系结构整体分析
4.1.6 整体刚度矩阵的建立
1
2 3
41
2
1 2 2
1
①
②
③若将 [A]按单元分成图示
三个子矩阵
? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
? ? ? ? ? ? ? ? ? ? ? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
00000
0000
0000
00000
I
II
II
I
A
则
由此可见，整体刚度矩阵可由各单元整体刚度
矩阵 装配累加 得到。为说明如何装配，先将单元
刚度矩阵进行分割
? ? ? ? ? ?? ? ? ?
i
i kk
kk
k ?
?
?
?
?
?
?
2221
1211
r
s
sr
整体结点码
则由矩阵乘法可证明,[A]i[k]i[A]iT的结果是，将
刚度矩阵子矩阵按整体结点码 r, s 送整体刚度矩阵
相应位臵。这一装配规则称为“对号入座”。
4.1 杆系结构整体分析
1) 任意结构情况
上面结论是通过具体例子（全刚结点平面刚架）
得到的，由虚位移原理或势能原理进行整体分析
（见讲义），可得任意结构其结论同此例。
2) 结点位移编号
如果按结点顺序，对结点
非零位移进行依次编号，这一
序号称作结点位移码。为便于
计算机处理并减少结构刚度矩
阵的阶次，将零位移的号码变
为零。
对图示三铰刚架，当仅
用一种单元（梁柱自由是
单元）时结点位移编号如
图所示。
3) 单元定位向量
按单元局部结点码顺序，
将结点位移码排成的向量，称作单元的 定位向量 。
(0,0,1) (0,0,2)
(3,4,5) (6,7,8) (6,7,9)
( 1 0,1 1,1 2 )
1 2
3 4 5
6
① ②
③ ④
(0,0,1) (0,0,2)
(3,4,5) (6,7,8) (6,7,9)
( 1 0,1 1,1 2 )
1 2
3 4 5
6
① ②
③ ④
4.1 杆系结构整体分析
对图示刚架各单元的定位
向量为
① (0,0,1,3,4,5)
② (0,0,2,10,11,12)
③ (3,4,5,6,7,8)
④ (6,7,9,10,11,12)
① 2 3)
②,6,7,8)
③ 1 2 3 4 5)
④ 4 5 6 7,8)
4) 按单元定位向量集装刚度矩阵和综合荷载
前面说明的是分块子矩阵集装,下面说明如何按定
为向量来集装,
如果如图所是采用各种不同的单元 （一端有铰
）,则定位向量为
(0,0,1) (0,0,2)
(3,4,5) (6,7,8) (6,7,9)
( 1 0,1 1,1 2 )
1 2
3 4 5
6
① ②
③ ④
(0,0) (0,0)
(1,2,3) (4,5) (6,7,8)
1 2
3 4 5
① ②
③ ④
如何获得
带铰的单元刚
度矩阵和等
效荷载矩阵
定位向量
(0,0,1) (0,0,2)
(3,4,5) (6,7,8) (6,7,9)
( 1 0,1 1,1 2 )
1 2
3 4 5
6
① ②
③ ④
(0,0) (0,0)
(1,2,3) (4,5) (6,7,8)
1 2
3 4 5
① ②
③ ④
4.1 杆系结构整体分析
4-1) 刚度集装
以 3 单元为例来说明
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
5554535251
4544434241
3534333231
2524232221
1514131211
3
k
87654
定位向量
单元局部位移码
54321
5
4
3
2
1
8
7
6
5
4
根据单元
局部位移码
和定位向量
的对应关系
用定位向量
位移码送元
素。
根据单元局部位移码和定位向量的对应关系用定
位向量位移码送元素，定位向量元素为零时不送。
(0,0,1) (0,0,2)
(3,4,5) (6,7,8) (6,7,9)
( 1 0,1 1,1 2 )
1 2
3 4 5
6
① ②
③ ④
(0,0) (0,0)
(1,2,3) (4,5) (6,7,8)
1 2
3 4 5
① ②
③ ④
4.1 杆系结构整体分析
4-2) 荷载集装
以 4 单元为例来说明
定位向量
8
7
6
0
0
? ?
45
4
3
2
1
4
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
P
P
P
P
P
P
5
4
3
2
1
局
部
位
移
码 此结论同样适用于刚度集装
4.1 杆系结构整体分析
4.1.7 整体分析总结
1) 对局部坐标和整体坐标不一致的单元，要对刚
度、荷载进行坐标转换。
2) 需对“结构”进行结点、位移的局部和整体编
号。
4) 集装所得整体刚度矩阵是对称、带状稀疏矩阵，
当支撑条件能限制刚体位移时，矩阵非奇异。
3) 根据单元局部位移码和定位向量的对应关系用
定位向量位移码送元素，定位向量元素为零时
不送。据此可集装、累加得到整体刚度矩阵。
5) 综合荷载由两部分组成，因此首先要将直接作
用结点的荷载按结点位移码送入，如果还有单
元等效荷载，再按定位向量集装、累加。
4.1 杆系结构整体分析
8) 如果有某位移码方向弹性支撑，需进行将弹簧
刚度送入位移码对应的对角线元素位臵累加。
9) 如果有某位移码方向已知支撑位移，需进行将
“边界条件处理”。具体做法以后介绍。
7) 整体刚度方程实质是全部结点的平衡条件。
6) 刚度矩阵带状稀疏，其带宽取决于结点、位移
编码。
最大半带宽 =定位向量中最大元素差 +1。
4.1 杆系结构整体分析
2.5.8 边界条件的处理
10) 当用虚位移或势能原理作整体分析时（势能为
例），应变能为单元应变能之和，外力势能为
单元外力势能之和 +结点外力势能。全部杆端力
的外力势能彼此抵消。
1) 乘大数法
设第 i 个位移为已知值 a, N =108 或更大的数。
乘大数法是 将刚度矩阵 Kii改为 N?Kii，将 Ri改为
N?a。
请考虑
为什麽这样做能使
边界条件得到满足？
2) 臵换法（划零臵 1）
设第 i 个位移为已知值 a 。
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
n
i
n
i
nnnjninn
inijiiii
nji
nji
R
R
R
R
KKKKK
KKKKK
KKKKK
KKKKK
2
1
2
1
21
21
2222221
1111211
?
?
?
?
4.1 杆系结构整体分析
上述臵换工作量大一些，显然可看出边界条件得
到精确满足。
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
aKR
a
aKR
aKR
KKKK
KKKK
KKKK
nin
i
i
n
i
nnnjnn
nj
nj
22
11
2
1
21
222221
111211
0
00100
0
0
?
?
?
?
4.1 杆系结构整体分析
★ 3) 关于斜边界的处理
如图示意的斜支座情况，有多种处理方案。
3-1) 通过单元的坐标转换来处理
?
x
y
?
xy
图示有斜支座单元,r 结点处
以倾角 ?- ? 来进行坐标转换，也
即在 r 结点处整体坐标为图示 xy 。
r
3-2) 通过增加一个单元来处理
图示有斜支座单元,r 结点处沿 y 方向增加一个
刚结的单元，此单元有“无穷大”的抗拉刚度、但
没有抗弯刚度。单元长度可任意。
3-3) 对整体刚度矩阵进行处理（参见匡文起教材）
最大
半带宽
2.5 杆系结构整体分析
2.5.10 总刚度矩阵的存储与对应解法
因为总刚度矩阵对称、带状稀疏，利用这一特
点可减少存储、提高解算效率。
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nnnin
iniii
ni
KKK
KKK
KKK
1
1
1111
零元素
零元素非零元素
最大
半带宽
主对角线元素
等
半
带
存
储
零元素
1
3
6
10
14
53
60
66
72
非
零
元
素
变
带
宽
一
维
存
储
带
宽
是
变
的
到 P.55
2.5 杆系结构整体分析
目前一般都用变带宽存储，下面结合程序说明
存储和解法。首先介绍一些 F 90的语法。
定义导出类型
导出类型 —— 结点
type,,typ_Joint
real,,x,y ！ 坐标
integer,,GDOF(3)！ 整体位移码
end type typ_Joint
1) 有关 F90 语法
?导出类型 新特性
用结点导出类型作为成员导出单元类型,
type,,typ_Element
integer,,JointNo(2) ！ 结点编号
type(typ_Joint),,Node(2) ！ 结点
integer,,GlbDOF(6) ! 定位向量
real,,EA,EI
end type typ_Element
type (typ_Element),,Elem(5) ！ 定义 5个单元类型
…
! 对单元 i的端点 j的 x,y,GDOF(1:3)的赋值
Elem(i)%Node(j) = Joint(Elem(i)%JointNo(j))
…
由导出类型定义新类型
由导出类型定义变量
real,,A(5),B(5,10),C(5)
B=0.0 !对 B清零
A=1.0 !对 A赋 1,A(i)=1.0,i=1,5
C=A+2 !数组与标量运算, A(1:5)+(/ 2,2,2,2,2 /)
A=C+A !数组与数组运算 （ 同形 ）
C=sqrt(A) !数组的函数运算, C(i)=sqrt(A(i),i=1,5
数组内部函数,
dot_product(vector_a,vector_b) !点积
如,dot_product((/1,2,3/),(/2,3,4/)) 的值为 20 （ 待续 ）
有关 F90 语法
数组运算与赋值,
matmul(matrix_a,matrix_b) 矩阵相乘
如,locEDisp = matmul(T,glbEDisp)
transpose(matrix) 矩阵转置
如,glbEDisp = matmul(transpose(T),locEDisp)
size(array,dim) 求数组第 dim维的长度
dim为可选变元, size(a,dim=2)
若 array为一维时, 可不用 dim。
sum(array,dim,mask) 数组元素求和
dim,mask为可选变元； mask=条件表达式
sum(a(1:10)) 对 a的 1到 10元素求和
sum(a(1:10),mask=a>0) 对 a(1:10)中大于 0的元素求和
? ? ? ?? ?ee dTd ?
（续）
有关 F90 语法
?where结构 新特性
例
where (C>0)
C = 0
A = B*D
end where
where (C>0)
A = B
end where
定义
where (数组关系表达式 )
数组赋值语句
…
elsewhere
数组赋值语句
...
end where
规则, 1)同形数组 ; 2)不许嵌套 ; 3)最多二分叉
有关 F90 语法
?cycle和 exit语句 新特性
用在 do循环中
cycle —— 作下一个循环步
exit —— 跳出循环，执行 end do后一条语句
等效例
do
...
if(.not.cond1)then
...
if(cond2) goto 5
...
end if
end do
5,..
用法
do
...
if (cond1) cycle
...
if (cond2) exit
...
end do
...
有关 F90 语法
?数组构造函数 spread
语法
spread(数组名,dim,ncopies)
将数组沿 dim维方向复制 ncopies形成新数组
dim,ncopies — 整型, 位置变元, 关键字变元
（ 若按位置引用, 可略关键字 ）
例, （仅限一维数组）
1） spread(one,dim=1,ncopies=3)
spread(one,1,3)
spread(one,ncopies=3,dim=1)
[1,1,1]
或
[1,1,1]T
2） ELocVec(1:6)=(/ 1,0,3,4,5,0 /)
spread(ELocVec,dim=1,ncopies=3)
?
?
?
?
?
?
?
?
?
?
054301
054301
054301
3）
spread(A(2,2:),dim=1,ncopies=2)
?
?
?
?
?
?
?
?
?
?
?
332313
322212
312111
A
...
...
...
??
?
??
?
3222
3222
..
..
如果 dim=2呢?
有关 F90 语法
?指针 pointer
pointer是变量的属性, 可以指向相同类型的变量 ;
被指向的目标必须具有 target属性或 pointer属性
可以将指针变量理解为别名, 称号
real,target,,a,b,EDisp(6) ! 可被指针所指
real,pointer,,p1,p2 ! 称号：班长, 课代表
！ p1,p2是指针, 可以指向实型数据
real,pointer,,G(:) ! 先进集体
！ G是指针, 可以指向一维实型数组
? 指针是一种, 称号,, 上述声明语句建立了, 称
号,, 但并未, 授予, 某个变量这个称号, 因此
是指向, 空,, 并未占用内存
a = 3.0
p1 => a ! p1指向 a
! 称号 p1授予 a,a的数据有两个名：固定名 a和流动
名 p1; 既可用 p1也可用 a（ p1 — 班长, a — 张三 ）
a = 4.0 ! a的值变为 4.0,p1也变为 4.0
p1 => b ! 班长换人了
G => EDisp ! 先进集体有了得主
！ EDisp(:)的长度就是 G(:)的长度, 用 G和用 EDisp
同样效果
又如：
real,target,,a,b
real,pointer:,p,q
a = 3.14
b = 2.0
p => a ! p = a = 3.14
q => b ! q = b = 2.0
q = p
! ( q指向的数据 b ) = ( p指向的数据 a )
！ 即,所有 = 3.14
? 指针可以指向有名的数据区, 也可以指向无名的
数据区
real,pointer:,b(:,:)
integer,,n
read(*,*) n
allocate (b(n,n)) !指针指向了一个刚开辟的数组
! 以下可以当作常规数组用
b(1,1) = 1.1
b(1,2) = 1.2
...
deallocate (b)
有关 F90 语法
?用指针建立动态数组
指针与 allocatable数组的区别
? 具备 allocatable数组的所有功能
? 还可以用在导出类型中, 例如整体刚度矩阵的变带宽存
储：
type,,typ_Kcol !整体刚度矩阵 K的列
real(8),pointer,,row(:) !该列的行元素
end type
...
type (typ_Kcol),allocatable,,Kcol(:)
…
allocate (Kcol(NGlbDOF)) ! 分配了 NGlbDOF列
...
allocate (Kcol(5)%row(3:5)) ! 第 5列只用 3至 5行
(1) NElem,NJoint,NGlbDOF,NJLoad,NELoad
单元数 结点数 总自由度数 结点荷载数 单元荷载数
[Joint - 结点 ] … NJoint行
(2) Joint%X,Joint%Y,GDOF(1:3)
X 坐标 Y坐标 结点位移码
[Elem - 单元 ] … NElem行
(3) JointNo1,JointNo2,EA,EI
起点号 终点号 刚度
[JLoad - 结点荷载 ] … NJLoad行
(4) JointNo,LodDOF,LodVal
作用点号 局部位移码 荷载值
[ELoad - 单元荷载 ] … NELoad行
(5) ElemNo,Indx,Pos,LodVal
单元号 类型号 位置 荷载值
Indx 类型 pos
1 均布 长度
2 集中 位置
3,..
2) 某程序输入数据说明
3,5,7,1,1
0,0,0,0,0
0,4,1,2,3
4,4,4,5,6
4,4,4,5,7
4,0,0,0,0
1,2,1.0E9,16
2,3,1.0E9,24
5,4,1.0E9,12
2,1,10.0E3 ! 结点 2，自由度 1，值为 10E3
2,1,4,-4.0E3 ！单元 2，均布，长 4米，值为 -4E3
2-1) 数据文件例子：
(2)
(1) (3)
2 4
1
3
5
i = 6
i = 4 i = 3
10 kN
4 kN/m
4m
4m
EA= 109 N
(1)
(2)
坐标
位移码
(3)
(4)
(5)
结点号 EA,EI
read(5,*) NElem,NJoint,NGlbDOF,NJLoad,NELoad
allocate (Joint(NJoint))
allocate (Elem(NElem))
allocate (JLoad(NJLoad))
allocate (ELoad(NELoad))
...
read(5,*) (Joint(i),i=1,NJoint)
read(5,*) (Elem(ie)%JointNo,Elem(ie)%EA,&
ELem(ie)%EI,ie=1,NElem)
if (NJLoad>0) read(5,*) (JLoad(i),i=1,NJLoad)
if (NELoad>0) read(5,*) (ELoad(i),i=1,NELoad)
2-2) 程序读入计算所需数据：
2-3) 求始行码和分配带宽子程序
!==================================subroutine SetMatBand (Kcol,Elem) ! 接口简单
!==================================
type(typ_Kcol),intent(in out),,Kcol(:) ！总刚列
type(typ_Element),intent(in)),,Elem(:) ！单元
integer(ikind),,minDOF,ELocVec(6)
integer(ikind),,Row1(size(Kcol,dim=1))
！ Row1为自动数组，子程序结束后自动释放。
！这样做可使接口简单，减少了数组的控制变量。
integer(ikind),,ie,j,NGlbDOF,NElem
NGlbDOF = size(Kcol,dim=1) ！使接口简单
NElem = size(Elem,dim=1)
Row1 = NGlbDOF ! 先设所有始行码同终行码
! 确定各列始行码向量
do ie=1,Nelem ! 对单元循环
! 确定定位向量
ELocVec(:) = Elem(ie)%GlbDOF(:)
! 寻找定位向量中大于零的最小值
minDOF = minval(ELocVec,mask=ELocVec>0)
! 屏蔽定位向量中小于等于零的
where (ELocVec > 0)
! 寻找 Row1(ELocVec)和 minDOF中的最小值
Row1(ELocVec) = min(Row1(ELocVec),&
minDOF)
end where
end do
! 为各列的带宽分配空间
do j=1,NGlbDOF
! 对总自由度数循环
allocate (Kcol(j)%row(Row1(j):j)
! 给 Kcol(j)%row 从 Row1(j) 到 j 个空间
Kcol(j)%row = zero
! 总刚元素清零
end do
return
end subroutine SetMatBand
621,.,,,,j,i ?
IJ)ie(ij Kk ? ??
??
?
?
?
??
??
0)(
0)(
)ie(
j
)ie(
i
J
I
?
?
3) 整体刚度矩阵的集成
do ie=1,NElem
… ! 计算单刚 EK(6,6),ELocVec(6)
do j=1,6 ! 对单元逐列集成
JGDOF = ELocVec(j) ! 取出位移码
if (JGDOF = 0) cycle ! 作下一循环步
where (ELocVec > 0, And, ELocVec <= JGDOF)
! 位移码非零同时小于第 j个局部码对应的位移码
Kcol(JGDOF)%row(ELocVec) = &
Kcol(JGDOF)%row(ELocVec) + EK(:,j) ! 集成
end where
end do
end do
局部码 位移码
4) 变带宽矩阵的分解求解
4-1) LDLT分解法求解 [A]{x} = {b}
若 ? ?
nn?A
对称正定，则可分解为
? ? ? ?? ?? ?TLDLA ? =
单位上三角阵
对角阵
原方程变为
求解步骤：
1,分解：
2,解 y,
3,解 x,
LDLT分解法 Gauss消去
? ?? ?? ? ? ? ? ?? ? ? ?byLxLDL ??T
? ? ? ?? ?? ?TLDLA ?
? ?? ? ? ?byL ?
? ?? ? ? ? ? ?yxLD ?T
前消去
处理右端项
回代
(不同的 b只做一次分解 )
? ? ? ?? ?? ? ? ?? ?ULLDLA ?? T
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nnnn
ijii
nj
nj
aaa
aaa
aaaa
aaaa
??
???
?
????
??
??
21
21
222221
111211
存放：
主对角 —
上三角 —
iid
jiij ll ?T
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
1
1
1
1
21
3231
21
0
??
?
nn
ij
ll
l
ll
l
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nn
jiii
n
n
d
ld
ldldd
ldldldd
?
?
?
0
222322222
1113111211111
jiiiij ldu ?
ija
ijl
? ?
? ?
??
n
k
n
k
jkkkikkjikij ldlula
1 1
?
?
?
?
?
?
?
?
?
?
ikl
ik
ik
l
ik
ik 1
0
j,ik ?? 时不必求和
? ?
? ?
??
)j,i(
k
i
k
kkjkikij dlla
m i n
1 1
)(?
(上三角, i < j )
?
1111 da ? 1111 ad ?
211112 lda ?
121121 adl ?
222211122 dlda ?? 221112222 ldad ??
不动
存在 处
12a
存在 处
22a
(第 j 列系数 )
iiji
i
k
kkjkikij dldlla ?? ?
?
?
1
1
上三角, i < j
4-2) 分解一般情况,逐列分解
ii
i
k
kkjkikijji ddllal ??
?
?
???
? ?? ??
?
1
1
121 ?? j,,,i ?
对角线, i = j > 1
jj
j
k
kkjkjj ddla ?? ?
?
?
1
1
2
n.,,,,,j 32?
?
?
?
??
1
1
2
j
k
kkjkjjjj dlad
(第 j 列系数 )
第 i 列中第 1个非零元素的行码为：
4-3) 变列宽存贮的分解修正：
ii
i
ik
kkjkikijji ddllKl ?
?
?
?
???
?
?? ?
?
?
1
1
121 ?? j,,,ii ?
11 ?? ik
?
?
?
??
1
2
1
j
ik
kkjkjjjj dlKd
第 j 列中第 1个非零元素的行码为：
ii1
ji1
则
)m a x ( 111 ji i,ii ?
n.,,,,,j 32?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
77
66
575655
465444
36534333
2752423222
413111
0
0000
K
K
KKd
Kld
Klld
Kllld
lld
分解顺序
4-4) F90 实现
! 三角分解
diag(1:ncol) = (/ (Kcol(j)%row(j),j=1,ncol) /)
do j = 2,ncol
row1 = lbound(Kcol(j)%row,1) ! i_1 j
do i = row1,j-1
row_1 = max(row1,lbound(Kcol(i)%row,1)) ! i_1
s = sum(diag(row_1:i-1) *Kcol(i)%row(row_1:i-1) * &
Kcol(j)%row(row_1:i-1)) ! 求和部分
Kcol(j)%row(i) = (Kcol(j)%row(i)-s)/diag(i)
end do ！第 j列系数分解完毕
s = sum(diag(row1:j-1)*Kcol(j)%row(row1:j-1)**2)
diag(j) = diag(j) - s ! 第 j列的主对角元素
end do
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nn
nn
ij
b
b
b
y
y
y
ll
l
ll
l
??
??
?
2
1
2
1
21
3231
21
1
1
1
1
04-5) 一般情况解 y 公式
22121 byyl ??
iik
i
k
ik byyl ???
?
?
1
1
k
i
k
ikii ylby ?
?
?
??
1
1
n,,,i ?21?
12122 ylby ??
11 by ?
可不动
可存在 处2b
4-6) 变列宽存储解 y的修正
k
i
ik
ikii ylby ?
?
?
??
1
1
上三角的第 i 列
从第 行元素开始1i
! 回代步骤 1,自上而下
do i = 2,ncol
row1 = lbound(Kcol(i)%row,dim = 1)
GP(i) = GP(i) - sum( Kcol(i)%row(row1:i-1) &
* GP(row1:i-1) )
end do
求出后,不再用，可将 存在 处iy
iyib ib
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
??
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nnnn
n
n
dy
dy
dy
x
x
x
ll
lll
??
??
?
222
111
2
1
232
13121
1
1
1
1
0
4-7) 一般解 x 的 公式
nnnn dyx ?
??? 1nx
自下而上：
? ? ? ? ? ? ? ?yDxL 1T ??
? ?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
nnd
d
d
D
1
1
1
22
11
1
?
iii
n
ik
kkii dyxlx ?? ?
?? 1
?
??
??
n
ik
kkiiiii xldyx
1
?
这样做的缺欠：
nni)n(i)n(i xlyy ?? ? 1
自下而上逐行计算，行中遇到 0元素需跳过，不方便！
4-8) 改为自右向左逐列计算！
记
iiini dyy ?? )1(
nx
解出
)1( ?? n
nn yx
)1(
11
?
?? ?
n
nn yx
对 y 向量预处理
再解出
?
第 n 列乘上 后移到
右边去，修正 y 向量nx
1?nx
??
?
?
?
???
????
??
?
21
1
)(
11
1
)1()(
,,n,njyx
j,,iixlyy
j
jj
jji
j
i
j
i
?
?
! 回代步骤 2,自右向左
GP(:) = GP(:)/diag(:)
do j = ncol,2,-1
row1 = lbound(Kcol(j)%row,dim = 1)
GP(row1:j-1) = GP(row1:j-1) - GP(j) * &
Kcol(j)%row(row1:j-1)
end do
! 现在 GP就是解
4-9) 小结：
? 无须一维存贮，无须行列码转换定位
? 公式与程序直接对应翻译，直截了当
? 求和采用内部函数进行数组运算
? 动态内存，用多大、开多大
? 数据封装性好，接口简单,
subroutine SetMatBand (Kcol,Elem)
subroutine VarBandSolv (Disp,Kcol,GP)
? 通用性强
2.6 杆系结构内力计算
2.6.1 杆端内力计算公式
解方程的结果可以得到结点位移，有了位移就
可以进一步求各单元的内力。
解算步骤：
? ? ? ?? ? ? ?PdkF ??
? ?? ?? ? ? ?PdTk ??
从 [?]根据定位向量取出单元杆端位移。
由单元倾角确定是算 还是 。? ?? ?? ?dTk ? ?? ?dk
减去等效结点荷载或加上固端力矩阵。
2.6 杆系结构内力计算
2.6.2 杆中任意截面内力计算公式
注意事项：
由图示隔离体图，列杆段的平衡方程即可得到
任意截面的内力计算公式。 请自行写出。
按上述隔离体图所求得的内力是精确的。
M
N
Q
1F
2F3F
求得单元结点位移后，代入单元位移场、求应
变、求内力，这样所得的内力一般不满足平衡条
件。只是近似解。
