1、v第一节第一节 有限单元法的基本概念有限单元法的基本概念v第二节第二节 局部坐标系中的单元刚度矩阵局部坐标系中的单元刚度矩阵 v第三节第三节 单元刚度矩阵的坐标变换单元刚度矩阵的坐标变换 v第四节第四节 单元未知量编码单元未知量编码 v第五节第五节 平面结构的整体刚度矩阵平面结构的整体刚度矩阵 v第六节第六节 非结点荷载处理非结点荷载处理 v第七节第七节 平面结构分析算例平面结构分析算例第四讲 杆件有限元3.有限元位移法分析连续梁需要考虑的问题(1)刚度集成法:将(1-3)K 扩阶,扩大的元素为 0,得到单元贡献矩阵 单元:K=32132100000jjjiijiikkkk 单元:K=3213
2、2100000jjjiijiikkkk 将单元贡献矩阵想叠加,形成整体刚度矩阵 K=K+K=321321002222111jjjijiiijjjiijiikkkkkkkk (1)两端支承条件的引入 先不考虑约束条件,得到整体刚度矩阵后,将其主对角线元素 kii改为 1,第 i 行,第j 列其余元素改为 0,对应的载荷元素也改为 0.(2)非结点荷载的处理 利用等效结点荷载进行分析:各结点(包括两端结点)加约束,阻止结点转动,其约束力矩分别为交于该结点的各相关单元的固端力矩之和,顺时针为正.去掉附加约束(相当在各结点施加外力荷载 Pe,其大小与约束力矩相同,方向相反)将两部分杆端弯矩叠加起来.第
3、二节第二节 局部坐标系中的单元刚度矩阵局部坐标系中的单元刚度矩阵 1.一般单元 设单元 e 的弹性模量、截面惯性矩、截面积分别为E、I、A,杆长为 l。单元的 i、j 端各有三个杆端力、MYX和、(即轴力、剪力和弯矩)和与其相应的三个杆端位移、vu，如图 1-7 所示。图中yOx为单元局部坐标系，取 i 点位于坐标原点，x轴与杆轴重合，规定由 i 到 j 为x轴的正方向，由x轴顺时针旋转 90?为y轴正方向。力和位移的正方向如图 1-7 所示。在此单元中，单元杆端力列阵和杆端位移列阵分别为 TjjjiiiMYXMYXF 单元杆端力列阵 Tjjjiiivuvu 杆端位移列阵 为了导出一般单元杆端
4、力与杆端位移之间的关系，我们分别考虑以下两种情况。首先分析两个杆端轴力jiXX、与轴向位移jiuu、的关系。根据胡克定律，有 jijijjijiiulEAulEAuulEAXulEAulEAuulEAX)()(（a）其次考虑杆端弯矩jiMM、与杆端剪力jiYY、与杆端转角ji、和横向位移jivv、的关系。根据结构力学位移法的转角位移方程，并按照本节规定的符号和正负号，可得 jjiijjjiiijjiijjjiiilEIvlEIlEIvlEIYlEIvlEIlEIvlEIYlEIvlEIlEIvlEIMlEIvlEIlEIvlEIM23232323222261261261261246262646
5、 (b)将(a)、(b)两式合在一起，并写成矩阵形式如下 ejjjiiieevuvulEIlEIlEIlEIlEIlEIlEIlEIlEAlEAlEIlEIlEIlEIlEIlEIlEIlEIlEAlEA460260612061200000260460612061200000222323222323jjjiiiMYXMYX (1-17)上式可简写成 eeekF (1-18)其中单元刚度矩阵为 lEIlEIlEIlEIlEIlEIlEIlEIlEAlEAlEIlEIlEIlEIlEIlEIlEIlEIlEAlEAke4602606120612000002604606120612000002223
6、23222323 （1-19）1.单元刚度矩阵的性质 每个元素代表单位杆端位移引起的杆端力,任一元素rsk(r、s 取 1 至 6)的物理意义是第 s个杆端位移分量等于 1 时,所引起的第 r 各杆端力分量值.是对称矩阵,其元素)(srkksrrs.是奇异矩阵,它的元素行列式等于零,即0k.具有分快性质.第三节第三节 单元刚度矩阵的坐标变换单元刚度矩阵的坐标变换 上述单元刚度方程和单元刚度矩阵实在局部坐标系yOx中建立起来的，对于一般杆件结构，分析时所划分的各单元的局部坐标系显然不同。因此在研究结构平衡条件和变形连续条件时，必须选定一个统一的坐标系 xOy,称为整体坐标系。同时，还必须把在局部
7、坐标系中建立的单元刚度矩阵转换为整体坐标系下的单元刚度矩阵。图 1-8a)、图 1-8b)分别表示单元 e 在局部坐标系yOx和整体坐标系 xOy 种的杆端力分量。为了导出整体坐标系中杆端力 Xi、Yi、Mi和局部坐标系中iiiZYX、之间的关系，将 Xi、Yi分别向yx、轴上投影，可得 cossinsincosiiiiiiYXYYXX a)式中，表示由 x 轴到x轴之间的夹角，以顺时针为正。图 1-8 在两个坐标系中，力偶分量不变，即 iiMM b)同理，对于单元e j 端的杆端力可得 jjjjjjjjMMYXYYXXcossinsincos c)v单元坐标变换矩阵(1-26)1000000
8、cossin0000sincos0000001000000cossin0000sincoseTeTeTTT1于其转置矩阵为正交矩阵，逆矩阵等 第四节第四节 单元未知量编码单元未知量编码 为了便于编程计算，需要按一定规律对结点的位移分量编号。结构的节点位移有自由结点位移和支座结点位移（亦称支座结点位移）之分。自由结点位移是未知量。建立结构整体结构方程求解未知节点位移的方式有两种：“前处理法”和“后处理法”。前处理法前处理法:仅对未知的自由结点位移分量编号,得到的结点位移列阵中不包含已知的约束结点位移分量.后处理法后处理法:由单元刚度矩阵形成整体刚度矩阵,建立刚度方程后在引入支承条件,进而求解结点
9、的位置位移的方法.图(1-9a)为一平面刚架。其中 1、2、3、4 为结点号；、为单元杆号。XOy 为选定的整体坐标系。该结构承受任意荷载如图 1-9a)所示。用后处理法分析改结构时，设所有点位移都是未知量，则结点位移列阵为（参看图 1-9）TTvuvuvuvu4443332221114321 Pix、Piy、Pi分别代表作用在结点 i（i=1，2，3，4）上的水平力、竖向力和力偶。规定，结点力 Pix、Piy的正方向与整体坐标系 x、y 的正方向相同，Pi以顺时针指向为正；结点位移的正方向与结点力的正方向一致。在求出各单元刚度方程之后，根据结点平衡条件和位移连续条件，可建立整个结构的位移法方
10、程 43213333222211114321000000jjjiiiiijjjiijiijjjiijiikkkkkkkkkkkkPPPP （1-37）或 000 KP （1-38）其中 333322221111444342413433323124232221141312110000000jjjiiiiijjjiijiijjjiijiikkkkkkkkkkkkKKKKKKKKKKKKKKKKK （1-39）为结构的整体刚度矩阵，或称为结构的原始刚度矩阵。在建立整体刚度方程式（1-38）时，假定所有点位移都是未知量，相当整体结构无知座，因而在外力作用下，除了弹性变形外，还有可能发生刚体位移，此时，
11、各结点位移不能唯一确定。这说明式（1-39）为奇异矩阵，不能求逆，故利用（1-38）不能求结点位移。实际在图 2-9a）所示的刚架中，结点 1、4 为固定点，因此结点位移是已知的，支承条件全为零。将该支承条件引入到整体刚度方程，得 0032444342413433323124232221141312114321KKKKKKKKKKKKKKKKPPPP （1-40）可以分为两组方程，一组是 323332232232KKKKPP （1-41）它可以求结未知结点位移2、3。令一组是 3243222100KKPP （1-42）称为反力方程。利用式（1-41）求出结点位移2、3并代入上式后，便可计算未知
12、的支座反力。对于一般杆件结构，都可以按上述步骤进行分析。无论结构有多少个结点位移分量，经过调整其排列顺序，总可以将它分为两组：一组包括所有的未知结点的位移分量，以F表示；另一组为支座结点位移分量，以R表示。相应的，将全部结点分为两组，与F相应者为已知的结点力列阵，以 PF 表示；与R相应者为支座结点力列阵，以 PR表示。于是有 RFRFPPP00,与以上分析方法相配合，将整体刚度矩阵 K0 中的各元素重新排列，则 K00=P0 可写成 FRRRRFFRFFPPKKKKRF （1-43）展开上式得 RRRRFRFFRFRFFFPKKPKK 式（1-44）、式（1-45）式（1-44）为“修正的整
13、体刚度方程”，它与式（1-38）的区别在于引进了支承条件。后处理法:由单元刚度矩阵形成整体刚度矩阵,建立刚度方程后在引入支承条件,进而求解结点的位置位移的方法.前处理法:仅对未知的自由结点位移分量编号,得到的结点位移列阵中不包含已知的约束结点位移分量.束结点位移分量.图 1-10 所示的具有组合结点的刚架划分 为三个单元，其编号为、，各 杆之间的箭头表示局部坐标系的正方向，刚架结构编号为 15。下面考虑各单元 结点位移分量编号。采用“先处理法”作如下规定：仅对独立的位移分量按自然数编号，称为位移号。若某些位移分量由于联结 条件或直接杆轴向刚性条件的限制彼此 相等，则编号相同。相等，则编号相同。
14、(2)在支座处，由于刚性约束而使某些位移分量为零时，此位移分量编号为零。因此图 1-10 编号如下 单元结点编号 杆数 单元编号 始端 末端 单元位移分量编号 AB 1 2 000 123 BC 2 3 123 456 DC 5 4 000 457 在计算程序中，单元两端结点号可采用二维数组 JE(i,e)表示，称为“单元两端结点号数组”。JE(1,e)=单元始端的结点号 JE(2,e)=单元末端的结点号 在本例中 JE(1，1)=1，JE(2，1)=2 JE(1，2)=2，JE(2，2)=3 JE(1，3)=5，JE(2，3)=4 任意结点位移分量的位移号可用二维数组 JN(i,j)表示，称
15、为“结点位移号数组”JN（1，j）=结点 j 沿 x 方向的位移号 JN（2，j）=结点 j 沿 y方向的位移号 JN（2，j）=结点 j 角位移的位移号 在本例中，对第三结点而言，JN（1，3）=4 JN（2，3）=5 JN（3，3）=6 将单元 e 始端及末端得位移号排成一行（始端在前），此数码为“单元定位数组”，利用它可方便的形成杆端位移与相应结点位移间的协调条件。它的展开式为 edemmmm),.,(21 式中，d 个元素 m1md分别是单元 e 的两端位移分量所对应的位移号数值。在本例中，)754000()654321()321000(321mmm 第五节 平面结构的整体刚度矩阵 在
16、进行了单元分析得出单元刚度矩阵之后，需要进行整体分析。以“先处理法”为例，将离散单元重新组合成原结构，使其满足结构结点的位移连续条件和力的平衡条件，从而得到修正的结构刚度方程，即前面给出的式（1-44）KFFF=PF 式中：KFF称为修正的结构整体刚度矩阵；F、PF分别为自由结点位移与自由结点荷载列阵。当已知计算对象为自由结点位移分量而不至引起误解时，式（1-44）也常称为整体刚度方程，KFF简称为整体刚度矩阵，F、PF分别简称为结点位移列阵与荷载列阵。为了书写方便，下标常常略去。有单元刚度矩阵集成整体刚度矩阵，通常采用“直接刚度法”，把计算步骤分为两步，首先求出各单元贡献矩阵，然后将它们叠加
17、起来，得出整体刚度矩阵。然而在实际电算中，不便采用，原因是在计算中需要先将所有单元的贡献矩阵 Ke都保存起来，Ke的阶数与整体刚度矩阵 K 相同，这就占用了大量存储容量，因此实际运算中，采用“边定位，便累加”的方法。其原理没有变，而且结果相同。图 1-10 中单元的单元矩阵为 32100032100016665646361615655545352514645444342413635343332312625242322211615141312111kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 式中单元刚度矩阵的上面和右侧标记了单元结点位移分量编号。因为整体刚度矩阵个
18、元素是按位移分量编号排列的，按先处理法，单元刚度矩阵中对应于分量编号为零的元素不进入整体刚度矩阵，非零编号指明了其余各元素在整体刚度矩正中的行、列号。所以 311642115411144KkKkKk 321652215512145KkKkKk 331662315613146KkKkKk 单元的刚度矩阵为 65432165432126665646361615655545352514645444342413635343332312625242322211615141312112kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 各单元在 K 中的位置为 6,.,2,16,.
19、,2,12jiKkijij 单元的刚度矩阵为 75400075400036665646361615655545352514645444342413635343332312625242322211615141312113kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 各元素在 K 中的位置为 743645435444344KkKkKk 753655535545345KkKkKk 773765635646346KkKkKk 按以上定位方法，将三个单元的刚度矩阵有关元素一道整体刚度矩阵对应位置，得到 366365364266265264263262261356256211
20、14421114425325225134624634524534424424324224123623523423316623216523116422622522422315622215522115421621521421314621214521114400000000kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkK 在实际电算时，采用 将 K 置零,这是 K=077;将 k 中的相关元素,按照“对号入座”，累加到 K；将 k 中的相关元素，继续按照“对号入座”，累加到 K；将 k 中的相关元素，继续按照“对号入座”，累加到 K，整体
21、刚度矩阵最后完成。例1-3 求图 1-12 所示刚架的整体刚度矩阵 K。设各杆截面尺寸相同。A=0.5m2 I=1/24 m2 E=3104MPa 解（1）整理数据并进行编号。4342441012121030610100410300lEAlEIlEAlEA（2）求局部坐标系中单元刚度矩阵ek。由于单元、的尺寸完全相同，故有21kk，可直接利用公式（1-19）求得 1003005030030120301200030000300503001003003012030120003000030021kk（3）求整体坐标系中单元刚度矩阵ek。按式（1-32）和式（1-33），求得各单元在整体坐标系中的单元
22、刚度矩阵，并将单元结点位移分量编号标记于上面与左侧。411000032110003050030030000300030012300125003010003003000030003001230012000321k 4211040032110030050300301203012000300003005030010030030120301200030000300400321 kk（4）成整体刚度矩阵 K。采用本节中介绍的方法建立结构整体刚度矩阵如下 41010050300502003030303031200300312K 式中，子向量fjfiFF和分别为单元 e 在端点 i、j 的固端内力。几种非结点荷载作用下的单元固端力列于表 1-2 表表 1 1-2 2 两端固定梁的固端力两端固定梁的固端力 剪 力 弯 矩 简 图 QAB QBA MAB MBA )22(23233alallqa )2(233allqa )386(122222alallqa )34(1223allqa )2(32allPb )2(32bllPa 22lPab 22lbPa 36labm 36labm 2)3(llabm 2)3(llbam
