图1所示的平面桁架,其所受到的外荷载已在图中标出,已知各杆件的材料物理属性为EA=2×108 N, 试求各杆件内力.

1)建立单元刚度矩阵.根据图1中各节点的几何坐标,采用式(1)推导各单元刚度矩阵 kei. 以单元6为例,由单元6的整体坐标位置可得单元长度L=4 m,单元倾角θ=36.87°, 由式(1)得到单元6的单元刚度矩阵(为简化这里略去单位,下同)为

Ke6=[256 192 -256 -192

192 144 -192 -144

-256 -192 256 192

-192 -144 192 144]×102(14)

2)建立结构各单元的过渡矩阵.建立与整体刚度矩阵同阶的全零矩阵作为过渡矩阵(本例模型有6个结点,每结点2个自由度,共12阶),根据单元节点的连接关系建立结构关联向量(表2).以结点6为例,其关联向量(行向量)为[3 4 7 8],构造一个与关联向量同阶的全1行向量[1 1 1 1],求这两个行向量的直积[15] (也叫Kronnecker张量积),得到所有16个元素的行位置为

[3 3 3 3 4 4 4 4 7 7 7 7 8 8 8 8].

交换两个向量的顺序,再次求得两个行向量的直积,得到所有16个元素的列位置为

[3 4 7 8 3 4 7 8 3 4 7 8 3 4 7 8].

把单元6的刚度矩阵各元素按行排成一个行向量,根据得到的行列指标定位,将此行向量的各个元素置于过渡矩阵 kt6的相应位置,得到相应单元的过渡矩阵为

kt6=[0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0

0 0 256 192 0 0 -256 -192 0 0 0 0

0 0 192 144 0 0 -192 -144 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0

0 0 -256 -192 0 0 256 192 0 0 0 0

0 0 -192 -144 0 0 192 144 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0](15)

3)建立结构整体刚度矩阵.将所求得的各单元过渡矩阵加起来得到结构整体刚度矩阵为

KG=∑12i=1Kti=[756 -192 -256 192 -500 0 0 0 0 0 0 0

-192 144 192 -144 0 0 0 0 0 0 0 0

-256 192 1 012 0 0 0 -256 -192 -500 0 0 0

192 -144 0 955 0 -667 -192 -144 0 0 0 0

-500 0 0 0 1 256 -192 -500 0 -256 192 0 0

0 0 0 -667 -192 811 0 0 192 -144 0 0

0 0 -256 -192 -500 0 1 012 0 0 0 -256 192

0 0 -192 -144 0 0 0 955 0 667 192 -144

0 0 -500 0 -256 192 0 0 1 256 -192 -500 0

0 0 0 0 192 -144 0 667 -192 0 0 0

0 0 0 0 0 0 -256 192 -500 0 756 -192

0 0 0 0 0 0 192 -144 0 0 -192 144]×105(16)

值得指出的是,因为还没有考虑边界约束,此步骤得到的整体刚度矩阵是一个奇异矩阵.

4)施加边界条件,由于结点1在水平和竖直方向的位移被约束,即在结点1上水平和竖直方向的位移为0,结点4在竖直方向的位移被约束.即在结点4上竖直方向的位移为0,采用划0置1法[16]将这些边界条件应用到建立的整体刚度矩阵 KG中,得到修改后的整体刚度矩阵 K为

K=[1×10-5 0 0 0 0 0 0 0 0 0 0 0

0 1×10-5 0 0 0 0 0 0 0 0 0 0

0 0 1 012 0 0 0 -256 -192 -500 0 0 0

0 0 0 1×10-5 0 0 0 0 0 0 0 0

0 0 0 0 1 256 -192 -500 0 -256 192 0 0

0 0 0 0 -192 811 0 0 192 -144 0 0

0 0 -256 0 -500 0 1 012 0 0 0 -256 192

0 0 -192 0 0 0 0 955 0 667 192 -144

0 0 -500 0 -256 192 0 0 1 256 -192 -500 0

0 0 0 0 192 -144 0 667 -192 0 0 0

0 0 0 0 0 0 -256 192 -500 0 756 -192

0 0 0 0 0 0 192 -144 0 0 -192 144]×105(17)

5)施加荷载并建立结构刚度方程.由于荷载均施加在结点6的竖直向下方向上,因此由图1直接可得

f=[0 0 0 0 0 0 0 0 0 0 0 -10]T(18)

由此建立结构刚度方程为

KD=f(19)

6)计算单元内力.将式(17)和式(18)代入式(19)中求解结构刚度方程,得到各结点的位移,图2中的虚线即结点位移后的结构变形图.根据得到的结点位移可求得桁架杆单元的内力,见图2中数值.由于本研究主要讨论有限元建模过程中的整体矩阵集成方法,这里不再赘述.

图2 单元内力(单位:N)Fig.2 Diagram of element internal force(unit:N)