杆件结构有限单元法

2.1结构离散2.2单元的刚度矩阵2.3坐标变换2.4结构刚度方程2.5支座约束处理2.6刚度方程求解及内力计算例1.桁架结构计算示例例2.刚架结构计算示例第2章杆件结构分析的有限单元法2.1结构离散离散化要点：杆件的转折点、汇交点、自由端、集中载荷作用点、支承点以及沿杆长截面突变处等均可设置成结点。结构中两个结点间的每一个等截面直杆可以设置为一个单元。变截面杆件可分段处理成多个单元，仍按等截面杆单元计算，截面刚度取自各段中点处截面。对曲杆结构，可细分，用两点之间的直线代替曲线。在有限元法计算中，载荷作用到结点上。当结构有非结点载荷作用时，应该按照静力等效的原则将其变换为等效结点载荷。建立单元坐标系和结构整体坐标系。杆系结构离散化示意图

2.1结构离散

2.1结构离散符号规定：线位移及相应力与坐标轴方向一致时为正；转角位移和力矩，按右手法则定出的矢量方向若与坐标轴正向相一致时为正。基本未知量：链杆单元，平面每个节点两个线位移，空间每个节点3个线位移梁单元，平面每个节点两个线位移和一个转角位移，空间每个节点3个线位移和三个转角位移节点力与节点位移要一一对应平面刚架的梁单元e：平面桁架的链杆单元e：2.1结构离散6单元刚度方程反映单元的节点力和节点位移之间的关系;单元刚度矩阵可根据刚度系数的物理意义，由力和变形之间的关系确定。单元刚度矩阵也可以由位移函数和虚功原理或最小势能原理计算。（具有一般性）2.2单元刚度矩阵返回7单元刚度矩阵（方法一）

对于一般平面刚架梁单元，单元刚度方程：单元的刚度方程：描述单元的结点位移δ(e)与结点力

F(e)之间的关系.结点位移结点力对于两结点平面梁单元：单元的刚度矩阵K(e)为6*6的方阵对于平面桁架的链杆单元，K(e)为4*4的方阵。8平面刚架两结点梁单元的刚度方程根据单元刚度系数的物理意义,由梁单元受力和变形及等截面直杆的刚度方程可以给出。 注意：所有力与位移的正负号均决定于坐标系方向。平面刚架两节点梁单元：F(e)K(e)第i列元素的物理意义:

代表ui=1单独作用于基本结构引起的六个杆端力大小。第krs元素的物理意义:

代表us=1单独作用于基本结构引起的第r个杆端力Fr大小。δ(e)ui=1

vi=1

θi=1

uj=1

vj=1

θj=1

ui=1平面梁单元的单元刚度矩阵单独作用下返回ui=1

vi=1

θi=1

uj=1

vj=1

θj=1

平面梁单元的单元刚度矩阵vi=1单独作用下返回ui=1

vi=1

θi=1

uj=1

vj=1

θj=1

平面梁单元的单元刚度矩阵单独作用下返回ui=1

vi=1

θi=1

uj=1

vj=1

θj=1

uj=1平面梁单元的单元刚度矩阵单独作用下返回ui=1

vi=1

θi=1

uj=1

vj=1

θj=1

vj=1平面梁单元的单元刚度矩阵单独作用下返回ui=1

vi=1

θi=1

uj=1

vj=1

θj=1

平面梁单元的单元刚度矩阵单独作用下15平面一般梁单元的单元刚度方程为:杆端力

向量F

(e)单元杆端位移

向量δ(e)单元刚度矩阵

K(e)16梁单元的单元刚度矩阵为:单元刚度矩阵常用子块形式表示:其中每个都是3×3的方阵，子块

Kij(e)

表示杆端j作用一单位位移时,杆i端引起的杆端力。返回17平面桁架链杆单元

平面桁架单元只有轴向变形,杆端力只有轴力；

矩阵表示：返回1）单元位移函数选用准则单元位移函数的项数，至少应等于单元的自由度数，它的阶数至少包含常数项和一次项。位移函数中包含单元的刚体位移和常应变状态。单元的位移函数应保证在单元内连续，以及相邻单元之间的位移协调性。单元刚度矩阵（方法二）

根据结构分析类型选取单元类型和位移场，根据单元类型采用对应的单元位移模式，位移函数能反映真实结构的位移分布规律，保证计算精度及解的收敛性。选用的单元位移函数应当满足下列要求:

形函数在其对应节点的值为1，在其他节点的值为0。2）轴向拉压杆（链杆）单元的位移的函数由单元结点位移，代入位移函数中确定待定系数项：轴向位移形函数：轴向位移形函数：单元刚度矩阵（方法二）

梁单元挠度函数：3）梁单元平面弯曲的位移函数梁单元平面弯曲包括四个结点位移分量：由单元结点位移，代入位移函数中确定待定系数项：单元刚度矩阵（方法二）

称为形函数矩阵。

梁单元挠度函数：梁单元位移模式：3）梁单元平面弯曲的位移函数单元刚度矩阵（方法二）

4）平面刚架梁单元的应力应变

平面刚架梁单元内任一点的轴向线应变由两部分组成。则[B]—平面刚架梁单元的应变转换矩阵。轴向应变：弯曲应变：平面刚架梁单元应变：

单元刚度矩阵（方法二）

5）平面刚架梁单元的刚度矩阵

梁单元的i，j结点发生虚位移为

单元内相应的虚应变应为：

由虚功原理有

由于结点虚位移的任意性，可得到单元的刚度方程：单元的刚度矩阵：单元刚度矩阵（方法二）

横截面积

：横截面对主惯性轴z的惯性矩：

通过积分同样可以得到平面刚架梁单元的刚度矩阵：

平面刚架梁单元刚度矩阵25在整体分析时，要在结构整体坐标下进行杆端力（结点力）的叠加;通过坐标变换使所有单元的杆端力和杆端位移都变换到结构整体坐标下。在进行单元分析时，使用的是单元坐标系，各单元杆端力和杆端位移的排列顺序和符号要参照单元坐标系。力和位移均为矢量，方向不同不能代数相加。2.3杆件单元的坐标变换返回设两个坐标系之间的夹角为α，以整体到单元坐标系的转向与x到y轴的转向一致为正。26杆端力变换：杆端位移变换：平面链杆单元:平面杆件单元的坐标变换为正交矩阵:返回单元坐标变换矩阵：平面梁单元:整体坐标描述的单元刚度矩阵：

2.4结构刚度方程主要内容：231返回结构的刚度方程结构刚度矩阵的装配结构的荷载向量结点位移向量Δ：所有结点的位移按顺序排成一列；刚架每个结点都有3个位移，桁架每个结点都有2个位移。结点力向量P:作用在结点上的力按结点位移的顺序排成一列;单元内部的荷载等效到节点上。返回结构刚度方程：反映结点力向量P与结点位移向量Δ之间的关系，即：KΔ=PK为结构的总刚度矩阵，且为对称方阵。由个单元刚度矩阵装配叠加而成。结点力向量P与结点位移Δ一一对应，也按结点位移的顺序排列；以上所有的量都要用整体坐标表示。结构的刚度方程29具体做法：把整体坐标下的单元刚度矩阵根据结点编码把各子块送到总刚度矩阵K对应的位置中去。装配过程：“子块搬家，对号入座”。

整体坐标下单元刚度矩阵：如果i,j对应的结点编号为g,h，则单元刚度矩阵的各子块在总刚中的位置分别为：结构刚度矩阵的装配123ghn123ghn总刚度矩阵的集成子块列子块行每个单元的刚度矩阵都经过如上扩展和对号入座后，总刚度矩阵的各个子块经过简单的叠加即可得到最终的总刚度矩阵。“子块搬家，对号入座”如：图示平面刚架的总刚度矩阵的集成各单元结点编码如图。用矩阵记为：如果求出各单元整体坐标下的单元刚度矩阵：123456123456根据结点编号把各单元的子块搬入总刚K中的对应位置。同一位置各子块的对应元素相加；空位补0。如果修改各单元结点编号。则用矩阵：123456123456对角线下方划红色的元素需要存储。单元的子块搬入总刚度矩阵中的位置，完全取决于结构结点编码。对同一结构，如果改变了结点的编码，则总刚度矩阵完全不同。33总刚度矩阵的特点总刚度矩阵是一个对称矩阵；处于主对角线对称位置的两个元素是相等的，即kij=kji

。总刚度矩阵是一个稀疏的矩阵；大片的区域都是零元素，它的非零元素只分布在主对角线两侧的带状区域内。

最大半带宽d=（c+1)*m;其中：c——各单元两端结点编号差的最大值；

m——每个结点的自由度数；不相关结点对应的刚度子块均为0。总刚度矩阵是一个奇异矩阵；当没有引进支座约束条件的情况下，总刚度矩阵不存在逆矩阵。返回判断图示结构总刚矩阵的最大半带宽。d=(3+1)*3=1234d=(6+1)*3=21总刚矩阵中元素的排列与结点的顺序直接相关。d=(7+1)*2=16返回最大相关结点差3最大相关结点差6最大相关结点差735结构荷载向量不考虑约束反力，只由外荷载引起的结点力排成的向量则称为荷载向量。如图1，结构的结点力向量:结构的荷载向量:荷载向量P的构成：直接结点荷载Pd等效结点荷载PE：单元内的非结点荷载（如分布荷载，温荷载，惯性力等）等效移置到结点上得到的。

1234等效移置的方法:首先求出基本结构在非结点荷载作用下引起的固端力。最后将各固端力反向作用到结构的结点上去，即为该结点的等效结点荷载。1234等效结点荷载向量:荷载向量:2312荷载等效图基本结构:两端固定梁返回例4：计算图示结构的荷载向量。

返回荷载等效图

荷载等效图

382.5支座约束处理在形成结构刚度方程时没有考虑支座等约束条件;总刚度矩阵是一个奇异矩阵；刚度方程P=K没有唯一解，方程中包含任意大小的刚体位移。必须引进约束条件，消除刚体位移，才能得到唯一解。支座条件引进的目的就是使：约束结点，位移已知0,约束反力未知。

自由结点，荷载已知,位移是未知量。结构刚度方程：K11(1)K12(1)0000K21(1)K22(1)+K22(2)0K24(4)0000K33(3)K34(2)0000K43(2)K44(2)+K44(3)+K44(4)0K46(4)0K52(4)00K55(5)K56(5)000K64(4)K65(5)K66(4)+K66(5)123456789101112131415161718P1P2P3P4P5P6P7P8P9P10P11P12P13P14P15P16P17P18对总刚度方程P=KΔ初等变换--行列交换K11(1)00K12(1)000K33(3)00K34(2)000K55(5)K52(4)0K56(5)K21(1)00K22(1)+K22(2)K24(4)00K43(2)00K44(2)+K44(3)+K44(4)K46(4)00K65(5)0K64(4)K66(4)+K66(5)123789131415456101112161718P1P2P3P7P8P9P13P14P15P4P5P6P10P11P12P16P17P18Px

未知支座反力P1已知节点

荷载Δ1支座已知0位移Δx求解的未知量Kx1KxxK11K1x把初等变换后的总刚度方程KΔ=

P可写成将方程式展开得:已知节点力（节点荷载）未知节点力（支座反力）

未知节点位移求解的未知量已知节点位移（支座0位移)求解未知节点位移Δx计算未知约束反力Px对于刚性支座：支座反力对位移的计算没有影响，但位移决定支座反力。有限元法是面向计算机的方法，边界处理不要改变刚度方程的阶数，消除刚体位移且不影响自由结点的位移。常用的处理方法有：1.主元素置1法：要修正刚度矩阵和荷载向量；作法：0位移对应的刚度矩阵主元素置为1，相应的副元素置0；荷载对应项置为0。2.主元素乘大数法：仅修正刚度矩阵；作法：仅把0位移对应的刚度矩阵主元素乘以一大数（如1030）。返回支座约束条件的引进主元素置1法若第i个自由度（位移）为0，应保证刚度方程解得Δi=0。

将总刚度矩阵K中第i行的主元素（第i行的主对角线元素）改为1，即令K（i，i）=1。将第i行、i列的所有副元素都改为零。即令K（i，j）=0，K（j，i）=0(i≠j)。将荷载向量中的第i元素置为零，即令Pi=0。

经过这三步改动后，便可实现Δi=0。

第3步第2步第1步第i行方程支座约束条件的引进123456789101112131415161718主元素置1法引进支座条件：

P1P2P3P4P5P6P7P8P9P10P11P12P13P14P15P16P17P18000P4P5P6000P10P11P12000P16P17P18根据0位移修改刚度矩阵；根据0位移修改荷载向量。(1)(2)(3)(4)(5)123456789101112131415161710P11P12P13P14P15P16P17P18如改变单元结点的编号；如图引入边界条件后的总刚方程为：主元素置1法引进支座条件：

返回形成了荷载向量P，集成了总刚度矩阵K并且引进支座条件后，便可由刚度方程K△=

P求解结点位移△。这就转化为求解大型线性代数方程组问题。线性代数方程式组的解法：直接法和迭代法。直接解法：如高斯消去法，及其派生的LU、LDLT三角分解法。迭代法：塞德尔迭代法。2.6刚度方程求解及内力计算一、刚度方程求解返回47首先从整体的结点位移向量△中取出该单元的结点位移。设单元e左右两端的结点编号分别为m和n，则该单元的整体坐标表示的结点位移向量为：然后将进行坐标变换，换成用单元坐标表示。最后代入单元刚度方程,便可求各单元的杆端力如果单元上还作用非结点荷载，则需要叠加由非结点荷载引起的固端内力，得到真正的杆端力：二、单元杆端内力计算返回例1.用有限单元法计算图示桁架内力。EA=常量。

1.结构离散如图所示2.单元坐标描述的各单元刚度矩阵结点坐标：1（0,1），2（1.732,1），3（0,0）；单元结点：1（1,2）,2（3,2）3.整体坐标描述的各单元刚度矩阵a=00a=300返回4.按结点编码

装配整体刚度矩阵单刚坐标变换：1231

2

300123212325.结构荷载向量6.根据边界位移条件修正结构的刚度方程：5.结构荷载向量6.根据边界位移条件修正结构的刚度方程：主元素置1法5.结构荷载向量6.根据边界位移条件修正结构的刚度方程：主元素乘大数法7.求解结构的刚度方程组得位移：8.各单元的杆端位移向量：整体坐标下：转换到单元坐标下：9.单元坐标下各单元的杆端力向量：单元内力：返回55例2.应用有限元法计算各杆端