1平面杆件结构的有限单元法

第一篇

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

第1章

平面杆件结构的有限单元法1.1有限单元位移法的基本概念

所谓杆件是指从构造上来说其长度远大于其截面尺寸的一维构件。在结构力学上我们通常将承受轴力或扭矩的杆件称为杆，而将承受横向力和弯矩的杆件称为梁。在有限单元法中这两种情况的单元分别称为杆单元和梁单元。为方便起见，本书都称之为杆单元。杆系结构是最简单的一类结构，也是我们在工程上最常见的一类结构。一、有限单元法的基本思路

先分后合，既先将结构划分成各个单元，进行单元分析，然后再将各单元集合成结构整体，进行整体分析。二、举例说明

以简单的结构连续梁为例，介绍有限元位移法的基本概念和求解方法

图示简单的结构连续梁

1、划分单元取每一跨为一独立单元，单元编号为①、②，单元的线刚度为i1,i2

，结点编号为：1，2，3，（称为整体号）结点角位移为：

1

，

2

，3结点力偶：M1，M2，M3。

建立梁的整体坐标系xoy

根据右手坐标系规定，结点角位移和力偶以顺时针为正。

每个单元也建立坐标系

，称为局部坐标系其方向与整体坐标系一致。每一单元的始、末端分别记为i,j端，称为局部号。

每个单元也建立坐标系

，称为局部坐标系其方向与整体坐标系一致。每一单元的始、末端分别记为i,j端，称为局部号。

离散化后的杆件单元和结点，其内力只画出杆端弯矩，各杆端弯矩以顺时针为正。

2、建立数学模型

将结点1，2，3处的结点角位移为基本未知量。得出单元①，②的杆端弯矩和杆端转角关系式：

单元①

单元②根据结点1、2、3处位移连续条件，有代入式（a)、(b)：为本例的位移法方程根据结点1、2、3处平衡条件有式（d)、(e)代入(f)得：

3、采用矩阵形式

以表示单元序号，则上面的（a）、（b）式可统一写为：

或简写为：

式中：单元刚度矩阵单元杆端位移列阵单元杆端力列阵

或简写为：用矩阵形式表示得：

式中：整体刚度矩阵结点位移列阵结点荷载列三、几点说明

（1）、刚度集成法的应用

建立整体刚度矩阵，常用的方法：刚度集成法（直接刚度法），它是直接利用单元刚度矩阵的“叠加”来形成整体刚度矩阵。

①、将式（1－3）单元刚度矩阵扩阶单元①121233单元贡献矩阵

单元②121233

单元贡献矩阵相叠加，形成整体刚度矩阵121233121233对于具有n个结点和（n-1）个单元的多跨连续梁，其整体刚度矩阵如下：

（2）、两端支承条件的引入

以上在推导连续梁整体刚度矩阵时，没有涉及连续梁两端有固定端支座的情况。图示右端为固定端的两跨连续梁，先不考虑右端的约束条件，得出整体刚度矩阵与式（1－6）相同。

考虑右端转角为零的支承条件，求解基本未知量的基本方程为：

为了便于编写程序，希望引入支承条件后，矩阵的阶数和排列次序不变，，而又达到修正整体刚度方程的目的，将（k）式修改为如下形式：等效于基本方程式（1）和支承条件

连续梁两端支承条件的引入方法：将整体矩阵的主对角线元素改为1，第i行、i列的其余元素改为0，对应的荷载元素也改为0，其中i=1或n。

i＝1对应左端为固定支座；i＝n对应右端为固定支座。左右端同时为固定端支座时，应同时进行修改。

（3）、非结点荷载的处理首先在各结点加约束，阻止结点转动，各杆独立承担所受的荷载，杆端产生固端弯矩，记为：

各结点的约束力矩分别为交于该结点的各相关单元的固端力矩之和，以顺时针为正。

然后去掉这些附加的约束，这相当于在各结点施加一外力荷载，其大小与约束力矩相同，但方向相反。原非结点荷载的等效结点荷载

叠加既得非结点荷载作用下得各杆杆端弯矩。

非结点荷载作用下的杆端弯矩为两部分：1、在结点加阻止转动的约束条件下的固端弯矩；2、在等效结点荷载作用下的杆端弯矩。

例1.1应用有限元位移法求图示连续梁的内力

解（1）结构离散化，单元及结点编号如图所示（2）求固端力矩及等效结点荷载三个单元固端力矩为：

等效结点荷载列阵为：（3）求单元刚度矩阵及相应的单元贡献矩阵

（4）求整体刚度矩阵（5）引入支承条件本例支承条件为：得基本方程：（6）求解基本方程（7）计算各杆端弯矩连续梁的弯矩图1.2局部坐标系中的单元刚度矩阵单元刚度方程是单元的杆端力与杆端位移之间的关系式，单元刚度矩阵是单元的杆端位移－杆端力变换矩阵。1.2.1、一般单元

杆单元的长度为，截面面积为，截面惯性矩为，弹性模量为，单元的、端各有三个杆端力为、、，其对应的位移为、、单元杆端力列阵为：单元杆端位移列阵为：为了导出一般单元杆端力与杆端位移的关系，考虑以下两种情况。

考虑杆端弯矩、和杆端剪力、与杆端转角位移、和杆端的横向位移、的关系，可得：将式（a)(b)合在一起，得：简写为：其中单元刚度矩阵为：②单元刚度矩阵为对称矩阵，其元素

①

单元刚度矩阵中的每个元素代表单位杆端位移引起的杆端力。其中的任意元素的物理意义是第个杆端位移分量等于1（其余位移分量等于0）时，所引起的第个杆端力的分量值。1.2.2、单元刚度矩阵的性质

③一般单元的单元刚度矩阵是奇异矩阵，它的元素组成的行列式等于零，即。根据奇异矩阵的性质，没有逆矩阵。也就是说，如果给定杆端位移，根据（1-18）式可以求出杆端力的惟一解，但反过来，如果已知杆端力，则不能根据来确定杆端位移的惟一解。因为即使在杆端力已知的情况下，由于单元两端无任何约束，因此除出杆端自身变形外，还可以发生任意的刚体位移。

④

单元刚度矩阵具有分块的性质。用虚线把分为四个子矩阵，把和各分为两个子矩阵，因此，单元刚度矩阵方程（1－17）写为：

式中：

其任意子矩阵表示杆端力和杆端位移之间的关系。1.2.3、轴力单元只考虑轴向杆端位移和轴向杆端力的单元，称为轴力单元。为了便于坐标变化，上式写为：式中：1.3单元刚度矩阵的坐标变换

局部坐标系中的杆端力分量整体坐标系中的杆端力分量两个坐标系中，力偶分量不变，即：对于单元j端的杆端力：将(a)(b)(c)式子用矩阵形式可表示为：上式可以简写成：两种坐标系中单元杆端力的变换式局部坐标系中的单元杆端力列阵：单元坐标转换矩阵为：整体坐标系中的单元杆端力列阵：

为正交矩阵，其逆矩阵等于其转置矩阵，即有：

同样也适用于两种坐标系下的杆端位移之间的变换：

将式（1－25）、（1－28）代入式（1－18），得：上式两边同乘以，可以得到：令则得：两种坐标系中单元刚度矩阵的变换公式写成展开形式为：式中：以上推导方法和步骤完全适用于轴力单元。轴力单元在整体坐标系下的杆端力列阵和杆端位移列阵为：轴力单元中不需要考虑杆端转角

和杆端力M，坐标变换转换矩阵为：在整体坐标系下轴力单元中的单元刚度矩阵展开形式为：式中：

单元刚度矩阵中的每个元素代表单位杆端位移引起的杆端力。其中的任意元素的物理意义是第个杆端位移分量等于1（其余位移分量等于0）时，所引起的第个杆端力的分量值。为对称矩阵，可以用分块子矩阵表示。在整体坐标系下的单元刚度矩阵与局部坐标系下的单元刚度矩阵有类似的性质：对于一般单元，具有奇异性。1.4单元未知量编码建立结构整体刚度方程求解位置结点位移的方式有：“先处理法”和“后处理法”。1.4.1后处理法由单元刚度矩阵形成整体刚度矩阵，建立刚度方程后再引入支承条件，进而求解结点的未知位移的方法称为“后处理法”。图示一平面刚架，设所有结点位移都是未知量。结点位移列阵为：求出各单元刚度方程后，根据平衡条件和位移连续条件，可以建立整个结构的位移法方程：式中：结构的整体刚度矩阵应该注意到，在建立整体刚度矩阵方程式（1－38）时，我们假设所有结点都可能位移，相当于整体结构并无支座。因此在外力作用下，除了发生弹性变形外，还可能发生刚体位移，这样各结点位移不能唯一确定。这说明整体刚度矩阵式（1－39）为一奇异矩阵，不能求逆矩阵，故利用整体刚度方程式（1－38）也不能求出结点位移。

实际上，在图示刚架中，结点1和结点4均为固定端，因此结点位移是已知的。支承条件为：

这样，将上述支承条件引入到方程中，对整体刚度方程进行修改，可得：

对上述方程进行化简，分为两个方程组：

和

这样，利用第（1－41）式可以求得结点位移和，再根据第（1－42）2式可以计算未知的支座反力。对于一般杆件结构，都可以按上述步骤进行分析。无论结构具有多少个结点位移分量，经过调整其排列次序，总可以将他们分成两组：一组包括所有未知结点位移分量（或称自由结点位移分量），以表示；另一组为支座结点位移分量，以表示。全部结点力分量也分为两组，与相应者为已知的结点力列阵，以表示；与相应者为支座节点力列阵，以表示。则整体刚度方程可写成下列形式:展开上式得:当已知及时，上式可用以求自由结点位移分量；后式可用来计算支座反力。当无支座移动，既时，以上两式为：修正的整体刚度方程后处理法是将单元刚度矩阵集成整体刚度矩阵后，再引入支承条件予以修改，最后建立用以求解自由结点位移的修正的整体刚度方程。后处理法在对单元、结点的位移、力等未知量进行编号时比较方便。但当结构复杂，例如结构为复式框架、组合结构且支承较为复杂时，应用后处理法解题有诸多不便，此时常采用“先处理法”。

采用“先处理法”解题时，先考虑支承条件，根据支承条件仅对未知的自由结点位移分量编号，得到的结点位移矩阵中不包含已知的约束结点位移分量，这样得到的整体刚度方程就是式（1－44）所表示的“修正的整体刚度方程”。用它可以直接求解自由结点的位移分量。1.4.2先处理法具有组合结点的刚架划分为3个单元，其编号为①、②、③，各杆之上的箭头表示局部坐标系的轴正方向。考虑各单元结点的位移分量编号。采用“先处理法”需做如下规定：①仅对独立的位移分量按自然顺序编号，称为位移号。若某些位移分量由于联结条件或直杆轴向刚性条件（忽略轴向变形）的限制彼此相等，则将他们编为同一位移号。②在支座处，由于刚性约束而使某些位移分量为0，此位移分量编号为0。杆段单元编号单元结点编号单元位移分量编号始端末端AB①120,0,01,2,3BC②231,2,34,5,6DC③540,0,04,5,7在计算机程序中，单元两端的结点号可采用二维数组JE（i，e）表示，称为“单元两端结点号数组”JE（i，1）＝单元e始端的结点号JE（i，2）＝单元e末端的结点号本例中：JE（1，1）＝1JE（2，1）＝2JE（1，2）＝2JE（2，2）＝3JE（1，3）＝5JE（2，3）＝4任意结点位移分量的位移号可采用二维数组JN（i，j）表示，称为“结点位移数组”JN（1，j）＝结点j沿x方向的位移号JN（2，j）＝结点j沿y方向的位移号JN（3，j）＝结点j角位移的位移号本例中，对第3结点而言：JN（1，3）＝4JN（2，3）＝5JN（3，3）＝6将单元的始端及末端的位移编码排成一行（始端在前)，此数码称为“单元定位数组”。本例中例1.2试对图示结构的结点和单元进行编号，并用单元两端的结点号数组表表示EF杆单元两端结点号；用结点位移号数表示F结点的位移分量编号；写出CD杆单元定位数组。分别考虑以下两种情况：（1）考虑各杆的轴向变形；（2）忽略各杆的轴向变形。解：建立整体坐标系，将结构划分为6个单元，用箭头表示各单元的局部坐标系轴的正方向，根据结构特征编为9个结点号。其中联结①、③、④单元的铰接处编有3个结点号，原因是三个单元在铰接点处具有相同的线位移，但角位移不相同。单元②、③、⑤相交点编有两个结点号，因为单元②、③在此处为刚性联结，具有相同的角位移和线位移。（1）考虑各杆的轴向变形时，各结点位移分量编号如图所示。EF杆为第⑥单元，单元两端结点号数组表示为：JE（1，6）＝8JE（2，6）＝9F点的结点号为9，结点位移号数组表示为：JN（1，9）＝14JN（2，9）＝15JN（3，9）＝16CD杆为第③单元，单元定位数组表示为：（1）忽略各杆的轴向变形时，各结点位移分量编号如图所示。EF杆为第⑥单元，单元两端结点号数组不变，仍为：JE（1，6）＝8JE（2，6）＝9F点的结点号为9，结点位移号数组表示为：JN（1，9）＝8JN（2，9）＝0JN（3，9）＝10CD杆为第③单元，单元定位数组表示为：1.5平面结构的整体刚度矩阵在进行了单元分析得出单元刚度矩阵之后，需要进行整体分析，以“先处理法”为例，将离散单元重新组合成原结构，使其满足结构结点的位移连续条件和力的平衡条件，从而得到修正的结构整体刚度方程，既前面的式（1－44）。修正的结构整体刚度矩阵自由结点位移分量列阵自由结点荷载分量列阵由单元刚度矩阵集成整体刚度矩阵通常采用“直接刚度法”。其计算过程可以分为两步：首先求出各单元的贡献矩阵，然后将它们叠加起来，得出整体刚度矩阵。然而在实际电算中，这种做法是不便采用的，因为在计算中需先将所有单元的贡献矩阵储存起来，而各单元贡献矩阵的阶数与整体刚度矩阵的阶数相同，因这样就要占用大量的计算机存储容量。故在实际中并不是采用贡献矩阵法，而是利用各单元的定位数组，采用“边定位，边累加”的方法。

如图所示的结构有3个单元，5个结点，7个独立的位移分量，这样其整体刚度矩阵应为阶矩阵，即：

各元素在整体刚度矩阵中的位置为：单元单元刚度矩阵：单元②的刚度矩阵为：

中各元素在整体刚度矩阵中的位置为：

单元③的刚度矩阵为：

中各元素在整体刚度矩阵中的位置为：

这样，按照以上所讲的定位方法，将、、中的相关元素累加到整体刚度矩阵对应的元素上，可以得到整体刚度矩阵为：

在实际电算程中，采用“边定位，边累加”的方法，过程如下：

例1.3求图示刚架的整体刚度矩阵，设各杆截面尺寸相同。

解（1）整理数据并进行编号（2）求局部坐标系中单元刚度矩阵（3）求整体坐标系中单元刚度矩阵123000123004（4）形成整体刚度矩阵1.6非结点荷载的处理为分析平面结构而建立的整体刚度方程，反映了结构的结点荷载与结点位移之间的关系。作用在结构上的荷载除了直接作用在结点上的荷载之外，还有作用在杆件上的分布荷载、集中荷载等。这些非结点荷载应转换成等效结点荷载。将与相叠加，可得综合结点荷载。综合结点荷载亦称“总结点荷载”其下标c通常可略去不写，即：

等效结点荷载的计算步骤如下：第一步，在局部坐标系下求单元的固端力。简图剪力弯矩QABQBAMABMBA表1－2

两端固定梁的固端力

第二步，求单元的等效结点荷载。

单元杆端力的变换式：固端力在两种坐标系下的变换式：有：等效结点荷载列阵可以由下式写出：展开，得到：

第三步，求整体结构的等效结点荷载。

求得单元等效结点荷载之后，利用单元结点位移分量编号，将中的各分量叠加到结构等效荷载列阵中去。因为中的各元素也是按结点位移分量编号排列的，中的6个元素也与结点位移分量编号一一对应，所以也按对号入座方法，将其逐一累加到相应的位置上去。

例1.4求图示刚架的等效结点荷载和综合结点荷载。解：单元局部坐标系及结点位移分量编号如图所示。（1）：求各单元在局部坐标系中的固端内力单元①：由表1－2可得：单元②：因此有：（2）求整体坐标系下各单元的等效结点荷载：单元①，＝900单元②，＝00利用式（1－49）或（1－50）可得：（3）求刚架等效结点荷载列阵，结果为：（4）求直接作用在结点上的荷载，（5）求总结点荷载列阵，应用有限元位移法对杆系结构进行分析，如果按照“先处理法”的思想，其步骤可以归结为如下：（1）整理原始数据，对单元和结点进行编号，确定每个单元的局部坐标系整个结构的整体坐标系。（2）计算局部坐标系中的单元刚度矩阵。（3）按结点位移分量编号，“对号入座，同号叠加”，形成结构整体刚度矩阵。

1.7平面结构分析算例（5）求总的结点荷载列阵。在这一步里，首先必须将非结点荷载转换成等效结点荷载，再与对应的结点荷载叠加，形成总的结点荷载矩阵。