杆梁结构的有限元分析原理.ppt

第3章 杆系结构的有限元法杆、梁单元概述讨论杆梁单元和由它们组成的平面和空间杆梁结构系统.o从构造上来说其长度远大于其截面尺寸的一维构件o承受轴力或扭矩的杆件成为杆o杆梁问题都有精确解（且是唯一的）o承受横向力和弯矩的杆件称为梁o平面桁架 平面刚架 连续梁 空间刚架 空间桁架等o承受轴力或扭矩的杆件称为杆o将承受横向力和弯矩的杆件称为梁o变截面杆和弯曲杆件本章主要内容3.1有限元分析的完整过程E1=E2=2E7PaA1=A2=2cm2l1=l2=10cmP3为10N作用下二杆结构的变形。E1、A1E2、A2说明：说明：u1、u2、u2分别表示分别表示节点节点1、2、3的水平位移的水平位移1）用标准化的分段小单元来逼近原结构2）寻找能够满足位移边界条件的许可位移场3）基于位移场的最小势能原理来求解 基本变量为：节点位移内部各点位移应变应力(1)(3)(2)问题的解题思路完整的求解过程1）结构离散化 该构件由两根杆件做成，因此可以自然离散成2个杆单元。假定以这类单元位移的特征为两个端点位移，就这两个离散单元给出节点编号和单元编号。单元1：i=1，j=2 单元2：i=2，j=3 单元位移模式：u(x)=a0+a1x 单元节点条件：u(0)=u1，u(l)=u2 将式（b）代入式（a），从而得2）单元分析（a）（b）回代得 写成矩阵形式为其中Ni，Nj是形函数。形函数矩阵形函数矩阵说明：说明：u表示位移列阵表示位移列阵 ue表示单元位移表示单元位移根据几何方程可得单元应变的表达单元应变写成矩阵形式为简记为几何函数矩阵或者是应变转换矩阵几何函数矩阵或者是应变转换矩阵根据物理方程可得单元应力的表达单元应力写成矩阵形式为简记为单元应力矩阵或者是应力转换矩阵单元应力矩阵或者是应力转换矩阵节点位移列阵单元e势能的表达 说明说明 积分域，积分域，P1、P2、分别表示作用单元分别表示作用单元e上的节点在上的节点在u1、u2的力的力写成矩阵形式为单元单元e刚度矩阵刚度矩阵单元单元e节点力列阵节点力列阵在得到各个单元的势能表达式后，需要进行离散单元的装配，以求出整个系统的总势能，对于该系统，总势能包括两个单元部分3）离散单元的装配处理边界条件是获取可能位移场，将左端的约束条件，即u1=0代入上式可以得到简化的势能表达式4）边界条件的处理由于上式是基于许可位移场的表达的系统势能，这是由全部节点位移分段所插值出的位移场为全场许位移场，且基本未知量为节点位移，根据最小势能原理（即针对未知位移求一阶变分）有5）建立刚度方程将结构参数和外载荷代入上式有求解得（单位m）6）求解节点位移7）计算单元应变8）计算单元应力对于单元势能的表达，对其取极值有具体地对于单元1，有其中R1是节点1的支反力，P2是单元1的节点2所受的力，即单元2对该节点的作用力，将前面求得的节点位移代入上式可得支反力大小。9）计算支反力以上是一个简单结构有限元方法求解得完整过程，对于复杂结构，其求解过程完全相同，由于每一个步骤都具备标准化和规范性的特征，所以可以在计算机上编程而自动实现。讨论1：对于一个单元的势能取极值，所得到的方程为节点的位移和节点力之间的关系，也称为单元的平衡关系，由此可以求出每一个单元所受的节点力。讨论2：由前面的步骤，我们也可以直接将各个单元的刚度矩阵按照节点编号的对应位置来进行装配，即在未处理边界条件之前，先形成整体刚度矩阵。其物理意义是，表示在未处理边界条件前的基于节点描述的总体平衡关系。在对该方程进行位移边界条件的处理后就可以求解，这样与先处理边界条件再求系统势能的最小值所获得的方程完全相同。小结 有限元分析的基本步骤及表达式1、物体几何区域的离散化2、单元的研究（所有力学信息都用节点位移）来表达3、装配集成4、边界条件的处理并求解节点位移5、支反力的求取以及其它力学量（应力、应变及位移三大物理量）的计算有限元分析的基本步骤及表达式一拉压杆单元图2.1 拉压杆单元示意图设杆单元长度为 ，横截面面积为 ，单元材料的弹性模量为 ，在局部坐标系中杆端荷载分别为 和 ，杆端位移分别为 和 ，单元上的轴向分布荷载为 。3.2 局部坐标下的杆单元分析用结点位移表示单元上任意截面的位移。对拉压杆单元，可以取其位移为一次多项式，即 由位移的边界条件：可得系数 、为：这样，截面任意一点的位移 为：用矩阵表示为：其中 （3-1）（3-2）单元位移模式。根据材料力学中应变的定义，有 这里 为应变矩阵。由虎克定律，其应力为：（3-4）（3-3）进行应力、应变分析其中利用虚位移原理求单元刚度矩阵，设杆端i、j分别产生虚位移 、，则由此引起的杆轴截面任意位置的虚位移为：对应的虚应变为：根据虚位移原理虚功方程（力乘以虚位移得虚功、外力虚功等于变形虚功），有：将上式整理得：（3-5）（3-6）求单元刚度矩阵式中 ：为局部坐标系下单元结点荷载矩阵。记 则可以得到拉压杆单元的单元刚度方程为：这里 为局部坐标系下的单元刚度矩阵，为局部坐标系下等效结点荷载矩阵，但值得指出的是：分布荷载 中可以包含集中荷载。根据定义，可以进一步求得单元刚度矩阵为：（3-10）（3-7）（3-8）（3-9）等效结点荷载二扭转杆单元图2 扭转杆单元示意图设扭转杆单元的长度为 ，截面惯性矩为 ，剪切模量为 ，杆端扭矩分别为 、，杆端扭转角分别为 、，单元上的分布荷载集度为 ，则任意截面的扭转角为 位移函数求得如（一）式中 为局部坐标系下扭转杆单元的结点位移矩阵（3-11）由材料力学可知，截面扭矩为：式中：我们利用极小势能原理来进行单元分析，杆单元的势能用泛函表示为：内力势能外力势能其中 为局部坐标系下扭转杆单元的结点集中荷载矩阵（3-12）由极小势能原理，取上述泛函的变分 ，可得：或者写为：设：可得扭转杆单元的单元刚度方程为：可以看到，其形式与拉压杆单元的单元刚度方程完全一致。同样，由上式可以进一步求得其局部坐标系下得单元刚度矩阵为：（3-13a）（3-13b)（3-14、3-15）（3-16）（3-17）等效结点荷载三只计弯曲的杆单元设杆单元的长度为 ，截面惯性矩为 ，弹性模量为 ，杆端集中剪力为 、，杆端集中弯矩分别为 、，杆端横向位移为 、，杆端扭转角分别为 、，在单元上分布有荷载集度为 的竖向分布荷载和集度为 的分布力偶。根据梁的平截面假定可知平面纯弯梁单元的轴向应变为：这里利用平截面假设（这里利用平截面假设（变形后横截面仍保持平面，与纵线正交）如图：材料力学基础知识取挠曲线方程为 的三次多项式，即单元上任意一点的挠度为：根据单元的位移边界条件：时：时：可以得到式中的待定系数结点位移矩阵和结点荷载矩阵分别为（3-18）将系数a、b、c、d代入式，并将挠曲线方程用矩阵形式表示为：式中 为形函数矩阵，其中：上式为平面弯曲单元的形函数。（3-19）（3-20）根据式（2-19）确定的单元位移场，可得单元上某一点得曲率为：截面的弯矩为：这里：为平面弯曲杆单元的应变矩阵。根据虚位移原理。有：（3-21）则平面弯曲杆单元的单元刚度方程为：其中的单元刚度矩阵可由式（2-23）求得为：记：（3-23）（3-22）等效结点荷载（3-24）（3-25）四平面一般杆单元（考虑拉伸、弯曲、不考虑扭转）杆单元的长度为 ，截面面积为 ，截面惯性矩为 ，弹性模量为 ，单元的、端各有三个力为 、和 、，其对应的位移为 、和 、，建立如图所示的局部坐标系，各物理量的正向如图中所标。设单元上没有荷载作用，首先考虑轴向力的作用，由于杆端轴力 、只引起杆端轴向位移 、，根据拉压杆单元的单元刚度方程，有：则结点位移矩阵和结点荷载矩阵分别为：（3-26）（3-27）其次，杆端弯矩 、和杆端剪力 、只与杆端的转角位移 、和杆端的横向位移 、有关系，根据只计弯曲杆单元的单元刚度方程（注意，由于不考虑单元上的荷载作用，故方程式中的等效结点荷载 等于零）可得：结构力学相关知识这样，上述表达式合并在一起，写成矩阵形式如下：可以将上式简写为：（3-28）（3-29）其中单元刚度矩阵（3-30）五 空间受力杆件单元（考虑扭转、拉伸、弯曲）对空间杆件单元，除了杆端力和结点位移数目较平面单元多外，其分析方法与平面杆单元类似（包含拉伸、扭转、两个方向弯曲）设局部坐标系的 轴为单元的形心主轴，横截面的两个主轴分别为 轴和 轴（如图所示）。设杆横截面面积为 ，杆单元长度为 ，在 平面内抗弯刚度为 ，在 平面内的抗弯刚度为 ，杆件的抗扭刚度为 。空间刚架有6个位移分量和6个结点力分量，设局部坐标系下它们分别为o纯轴向拉压o纯扭转oxoy面内弯曲oxoz面内弯曲（3-31）其中的单元刚度矩阵可写为将式（2-31）写成矩阵的形式有（3-32）（3-33）六单元刚度矩阵的性质 单元刚度矩阵 为对称矩阵，其元素 单元刚度矩阵 中的每个元素代表单位杆端位移引起的杆端力。其中的任意元素 的物理意义是第 个杆端位移分量等于1（其余位移分量等于0）时，所引起的第 个杆端力的分量值。一般单元的单元刚度矩阵 是奇异矩阵，它的元素组成的行列式等于零，即 。根据奇异矩阵的性质，没有逆矩阵。也就是说，如果给定杆端位移 ，根据（2-29）或（2-31）式可以求出杆端力 的惟一解，但反过来，如果已知杆端力 ，则不能根据 来确定杆端位移 的惟一解。因为即使在杆端力已知的情况下，由于单元两端无任何约束，因此除出杆端自身变形外，还可以发生任意的刚体位移。举例来说，如果物体处于静止状态，我们可以说其处于平衡状态，但反过来，如果物体处于平衡状态，则我们不能说其一定处于静止。单元刚度矩阵 具有分块的性质，即可以用子矩阵表示 。用虚线把 分为四个子矩阵，把 和 各分为两个子矩阵，因此，又可以写为：这里：或 或 或 或 用子矩阵形式表示单元刚度矩阵和单元刚度方程，可以使其表达的物理意义更加明显。在单元刚度矩阵 中，其任意子矩阵 表示杆端力 和杆端位移 之间的关系。（3-34）3.3杆系结构的整体分析（整体坐标系）一平面问题坐标转换矩阵图 平面问题杆端力转换示意图一般情况下，在进行单元分析时是在局部坐标下完成的。对于某一单元而，如果局部坐标系与整体坐标系不一致，则有单元分析的物理量必须通过坐标转换到整体坐标系中，然后再进行整体坐标系下的分析 这里 表示由 轴到 轴的角，角度转动的正负由右手定则确定，本书中以顺时针方向转动为正。在两个坐标系中，力偶分量保持不变，即有：同理，对于 端的杆端力，有：根据力的投影定理，将整体坐标下的杆端力分别投影到局部坐标下，有如下关系（3-35）（3-36）（3-36）将这些式子用矩阵形式可表示为：上式可以简写成：这即为两种坐标系下单元杆端力的坐标变换式。其坐标转换矩阵为：（3-38）（3-39）（3-40）从坐标转换矩阵 的表达式可以看出，为正交矩阵，其逆矩阵等于其转置矩阵，即有：并且有：式中 为单位矩阵。同样的推导，可以得到两种坐标系下的杆端位移之间的转换关系为：这里 和 分别为局部坐标系和整体坐标系下的杆端位移矩阵，为前面介绍的转换矩阵。整体坐标下的杆端力为局部坐标下的杆端力为（3-42）（3-41）（3-43）（3-44）（3-45）因此可得：上式两边同乘以 ，可以得到：设 可得：上式即为整体坐标系下的单元刚度方程。（3-47）（3-46）刚度矩阵转换）二 空间问题的坐标变换（空间问题）考虑结点i在局部坐标下 的杆端力 与在整体坐标系 的杆端力 的关系 设 轴与x、y、z轴的方向余弦分别为：，则将杆端力 、向 轴投影，可以求得杆端力 ，即：同理可以求得：用矩阵形式可以表示为：上式即为结点i的杆端力在局部坐标系和整体坐标系下的转换关系，其中的矩阵称为关系矩阵。与上面的推导类似，同样可以推出以 、表示 、，以及对于结点j的相对应的转换关系，其中转换关系矩阵都是 。综上所述，整体坐标系下单元杆端力矩阵与局部坐标系下单元杆端力矩阵具有如下的关系表达式：（3-48）其中的 为：称为空间坐标系的单元转换矩阵，它是一个正交矩阵，即：对于杆端位移，同样可推导出在两种坐标系中的转换关系：这样，可得空间杆件单元在整体坐标系中单元刚度方程为：其中表示空间单元在整体坐标系中的单元刚度矩阵。（3-49）（3-50）（3-51）（3-52）（3-53）三杆系结构的整体分析对杆系结构进行单元分析，仅仅是有限元分析中的第一步。我们的目的是要对整个结构进行分析，研究结构的整体性能。因此，在对结构的各单元分析完成后，必须将单元分析的结果进行整合，对结构进行整体分析。整体分析的过程实际上是如何将单元分析的结果进行有效组合，建立整体刚度方程并求解结点位移的过程。根据对结点位移的编码方式，可以采用“先处理法”和“后处理法”来建立整体刚度方程。1后处理法所谓后处理法，就是由单元刚度矩阵形成整体刚度矩阵，建立刚度方程后再引入支承条件，进而求解结点位移的方法。运用这种方法时，假设所有结点位移均为未知量，按照顺序统一进行编码，如图所示的平面杆件单元。结点位移矩阵为：结点荷载矩阵为：求出各单元刚度方程后，根据平衡条件和位移连续条件，可以建立整个结构的位移法方程：简 写成：这里 为结构的整体刚度矩阵，有：注意，在建立方程的过程中，我们假设所有结点都有位移。因此整个结构在外力作用下，除了发生弹性变形外，还可能发生刚体平动位移，这样各结点位移不能唯一确定。这说明整体刚度矩阵为一奇异矩阵，不能求逆矩阵，即根据整体刚度方程可得到无穷多个解。实际上，在图所示刚架中，结点1和结点4均为固定端，其三个位移分量均为0，即有：这样，将上述支承条件引入到方程中，对整体刚度方程进行修改，可得：对上述方程进行化简，可以得到两组方程：和 这样，利用第1式可以求得结点位移 和 ，再根据第2式可以求得支座反力 和 。2先（前）处理法所谓先处理法，就是先引入支承条件，根据支承条件仅对未知的自由结点位移分量编号，得到的位移矩阵中不包含已知的约束位移分量，即可以直接得到方程求解自由结点位移分量。如图所示的平面刚架结构ABCD，由于在A处和D处均为固定端，其位移为0，故位移编码均为0，在C处为铰接，故BC杆在C端的角位移与DC杆在C端的角位移不相同，因此在C处编两个结点3和4，但结点3和4的横向位移和竖向位移相同，故采用相同的编号，各结点位移编码如图所示。图先处理法位移编码示意图3.4等效结点荷载与边界条件的处理非结点等效荷载和边界条件的处理是有限元分析中必须考虑的的两个重要方面。由于只考虑结点荷载，因此必须将作用于单元上的非结点荷载转换到节点上。有限元的结点荷载来自两部分（1）作用于结点处的集中力、集中力偶，前面多次提到的 直接叠加到结点上即可（整体坐标下）（2）非结点荷载首先须在局部坐标下等效到结点荷载 ，然后再转到整体坐标系下的 1非结点荷载的处理在前面的分析中，我们已经介绍了求等效结点荷载的方法，如3-7式、3-15式、3-22式分别可用来求不同情况下的等效结点荷载。此外，可以这样来考虑：第一步，在局部坐标系下求单元的固端力 。对于某个单元，我们假定单元的两端均固定，然后根据静力平衡求得固定端的反力。第二步，根据单元固端力求单元的等效结点荷载 。根据局部坐标系与整体坐标系单元杆端力的变换式，固端内力在两种坐标系下的变换形式可以写成：因此，整体坐标系下的等效结点荷载矩阵可以由下式计算：（3-7）（3-15）（3-22）（拉压杆）纯弯杆扭转杆平面刚架单元固端力（遵行结构力学）2边界条件的处理（1）铰结点）铰结点在杆系结构中，除了刚性结点外，通常会遇到一些杆件通过铰结点与其它杆件联结，如下图所示杆件系统，有4根杆件汇交于D点，其中BD杆在D端通过铰支座与其它杆件铰接，其余3根杆为刚性接触。对于这样的铰结点，具有如下的性质：铰结点上各杆具有相同的线位移，但截面的转角位移不相同；结点上具有铰接杆端不承受弯矩作用。如下图所示结构中，BD杆在D端的杆端弯矩为0，只有CD、ED、GD杆在结点D上与外弯矩保持平衡。对于这样的结点，我们在对其进行单元划分时，通常考虑在D处设置2个结点。按照先处理法，对图示结构进行位移编码，如图2b所示。图铰结点的处理示意图 ANSYS通过结点耦合实现2弹性支承点弹性支承点在实际工程中，有时会遇到弹性支承的情况（如图），这时一般将弹性支座看作是在结构约束点沿约束方向的一个弹簧，弹簧的刚度系数为 ，在数值上等于使弹簧支座沿约束方向产生单位位移时所需施加的力。ANSYS引入弹簧单元即可具体做法可以归结为：先解除弹性支承点约束，在i处给一个结点号，形成总刚度矩阵，然后在总刚度矩阵中将第i行的主元素加上弹性支承的刚度系数，此时第行变为：以上的分析也适用与角位移为弹性约束的情况。若有多个弹性支座，可同时引入，即只需将相应的主对角线元素加上相应的弹性刚度系数即可。例1 3.5杆系结构有限元计算实例1 结构的离散化与编号 2 各个单元的矩阵描述 结构包括有斜杆，所以必须在总体坐标下对节点位移进行表达，所推导的单元刚度矩阵也要进行变换3 建立整体刚度方程 1.将所得到的各个单元刚度矩阵按节点编号进行组装，可以形成整体刚度矩阵;2.同时将所有节点载荷也进行组装。4 边界条件的处理及刚度方程求解 5 各单元应力的计算 6 支反力的计算 将节点位移的结果代入整体刚度方程中对于单元2：取i=1，j=2，则 ，故 对于单元1：取i=3，j=1，则c=1，s=0，故 对于单元3：取i=2，j=3，则c=0，s=1，故 例 2整体编号，对号入座得总刚消除边界条件后有 求解左式即可得到相关位移、应力、应变等等例3其中E1=3e4MPa,E2=3e4MPa解(1)单元划分，建立局部坐标系和整体坐标系，并对数据进行整理，对单元和结点编号(2)求局部坐标系下的单元刚度矩阵，由于单元1和2的尺寸完全一样，因此其单元刚度矩阵一样，为(3)求整体坐标系下的单元刚度矩阵，由于单元1和整体坐标为-90角，单元2的坐标和整体坐标完全一样，单元的转换矩阵为整体坐标系下的单元刚度矩阵整体坐标系下的总刚度矩阵（4）求总结点荷载 （a)等效结点荷载：首先求各单元在局部坐标下的固端内力，对单元1（局部y坐标向左），其q=-12kN/m，l=5m，对于单元2有F=8kN，得到局部坐标下的结点力向量利用转换坐标矩阵，得到整体坐标下的单元等效结点荷载为（b)直接作用在结点上的荷载（C)整体坐标下的结点荷载向量（5）建立整体平衡方程（6）求解（7）计算各单元的杆端内力，首先从求出的结点位移中取出各单元在整体坐标下的杆端位移，有然后计算杆端力（杆端力一定要恢复到局部坐标下完成）（c）剪力图 （b）弯矩图（d）轴力图（a）位移曲线例 4(1)单元划分，建立局部坐标系和整体坐标系，其中局部坐标下的x方向如箭头所示，y方向向下，z轴依据右手法则（2)形成局部坐标下的单刚将以上数据代入单刚矩阵且有（3）单元位移分量：依据先处理法：即先引入支承条件，根据支承条件仅对未知的自由结点位移分量进行编号，得到位移分量中不包含已知约束的位移分量（4）整体刚度矩阵：依据位移分量的对应情况，可以得到整体刚度矩阵，叠加同位置的单刚（对应零位移分量去掉）依据位移分量的对应情况，可以得到坐标下的结点力（对应零位移分量去掉）