“有限元法基础及应用”补充讲义(一)-ok-m.doc

“有限元法基础及应用”补充讲义一、引子弹簧单元与弹簧系统目标：掌握离散结构直接刚度法分析的原理和形式。了解有限元位移法列式的形式和基本概念。1、弹簧单元分析弹簧的物理特性：弹簧单元描述：图 1-1图 1-22个节点：节点位移：节点力：单元自由度：2 已知弹簧力位移关系：弹簧刚度 弹簧伸长量弹簧力，拉伸为正（1-1）考虑弹簧变形平衡时的条件和弹簧物理特性，得到下列方程：（1-2）写成矩阵形式：（1-3）写成矩阵符号形式：单元节点位移列阵弹簧单元的刚度矩阵式（1-2）、（1-3）为弹簧单元的刚度方程，反映了单元特性：节点力与节点位移之间的关系。式中：单元节点力列阵（注意：单元节点力是节点对单元的作用力。） 弹簧单元的刚度矩阵为：弹簧单元刚度方程讨论：1）有何特点？对称、奇异、主对角元素恒正。2）中元素代表什么含义？刚度系数大小等于弹簧刚度；每列元素代表一端固定、另一端产生单位位移时加在弹簧单元上的节点力。3）上面单元刚度方程可以求解吗？为什么？ 不可以。刚度方程仅仅表征一个典型单元的弹性特性，单元水平上无法确定单元节点位移。只有把系统中所有单元特性集成后，在系统水平上才可能求出所有未知位移和反力。单元水平上，若已知单元的节点位移，可由刚度方程求出所有单元节点力分量。若节点力已知，单元节点位移不能确定，单元可作刚体运动（小位移）。这也是单元刚度矩阵奇异性的物理解释。2、弹簧系统整体分析原理 以右图的一个弹簧系统为例，研究如何由单元特性集成系统特性并建立对系统进行求解的控制方程。图 1-3 由前面得到的弹簧单元的刚度方程公式（1-2），分别写出2个弹簧单元的特性方程如下：（1-4）单元1（1-5）单元2（注：右端节点力分量的下标1，2为单元节点的局部编号，上标是单元号）下面按两个方法完成系统特性的装配和控制方程的建立。并在特定条件下求解。1）由节点平衡方程导出： 系统处于平衡时，考虑各节点（节点1，2，3）的平衡条件：（1-6） 节点受到的外载荷与节点受到与其连接的所有单元对其作用力（单元节点力的反作用力）之和等于零。因此有下列（节点）平衡方程（组）：（1-7）把单元特性（1-4），（1-5）代入（1-6）得到： 写成矩阵形式：（1-8）（1-9）矩阵符号形式：式（1-8），（1-9）就是系统平衡方程，该方程建立了离散系统的外载荷与节点位移之间的关系，是求解节点位移的控制方程。 弹簧系统结构总刚度矩阵 系统节点位移列阵 系统节点载荷列阵讨论： （1） K 有那些特点和性质？ （2）上述方程能求解吗？2)由单元刚度方程叠加导出将单元1，2的刚度方程（1-4），（1-5）进行增广（扩大到系统规模）：（1-10）（1-11） 上述两个矩阵方程叠加，得：（1-12）上式中代入节点力平衡关系（1-6），就得到与（1-8）相同的节点平衡方程。 上述两种方法都必须考虑：单元特性集成；离散结构的节点上外载荷（系统外力）与节点力（系统内力）的平衡。因此方程（1-8）的本质是节点上力平衡关系，左边是由节点位移表示的（总）节点力，右边是节点所受外载荷。3）给定载荷和约束条件下的求解 （1-13）设边界条件为：则节点平衡方程（1-8）变化为：（1-14）该方程组展开后分为2个部分： 第2，3个方程变化为：（1-15） （1-16） 第1个方程变化为： 先后解方程（1-15）、（1-16）得到：（1-17）（1-18）至此解出了系统的未知位移和未知反力，并可以进一步求弹簧力。3、例题图1-4所示一个3个弹簧的系统。图 1-4 求： （a） 系统总刚度矩阵 （b） 节点2，3的位移 （c） 节点1、4的反力 （d） 弹簧2中的力解：（a）:写出各单元刚度矩阵：应用叠加法直接得到系统总刚度矩阵：或：该矩阵具有如下特点：对称、奇异、稀疏、非零元素沿主对角线呈带状分布。（b）:参考前面做法（1-8）式）及已求出的总刚度矩阵，写出系统节点平衡方程：（1-19）考虑到位移边界条件：则平衡方程组（1-19）第2，3方程化为：求解上式得：（c）： 由（1-19）的方程1，4得：（d）：（拉力）弹簧2内力为：4、练习题对图示弹簧系统，试用叠加法求其总刚度矩阵。并根据节点平衡方程的含义，尝试由各单元刚度矩阵的元素直接写出总刚度矩阵。图 1-5二、杆单元目标：通过杆单元特性方程的建立，初步掌握有限元法单元分析的过程和原理。了解杆系结构分析的原理。1、等截面杆单元及其刚度矩阵研究2节点等截面杆单元。L 杆长A 截面积E 弹性模量杆单元位移杆单元应变杆单元应力单元的力学量和基本关系如下：（2-1）应变位移关系：（2-2）应力应变关系：单元节点位移：单元节点力：下面研究杆单元的单元特性。1）直接法导出杆单元特性（2-3）杆单元伸长量：应用材料力学基本知识对单元进行力学分析。（2-7）（2-5）（2-6）（2-4）（2-8）杆的轴向刚度：杆内力：杆应力：杆应变：由于轴向变形模式下，杆单元的行为与弹簧单元相同，因此可比照弹簧单元的刚度方程（1-2），考虑到（2-7），直接写出杆单元的刚度方程：（2-9） 写成符号形式：（2-10）杆单元刚度矩阵为：2）公式法导出杆单元特性 步骤如下：（1）在单元上定义近似位移函数（位移模式）把一个单元上的位移分布假设为简单多项式函数。有限元法中用插值法通过节点位移分量作为待定参数来构造单元位移函数。对图2-1的杆单元，方便起见引入局部坐标 由于该杆单元只有2个未知位移分量，因此单元上假设的简单位移函数采用一次多项式。故对单元的节点位移进行线性插值。 容易定义出节点的插值函数如下：（2-11）因此单元上近似位移函数的插值形式为：（2-12）该位移函数也称为单元的位移模式，这里是线性位移模式。式（2-11）中的插值函数又称为形状函数，简称形函数。（2-13）式（2-12）写成矩阵形式为：上式中称为单元的形函数矩阵。式（2-13）是有限元法中最重要的关系式之一，通过该式把单元上的近似位移分布函数用节点位移来表示，为进行单元层次上的分析打下了基础。（2）单元应变和单元应力（2-14）由杆一维变形的应变位移方程（2-1）和单元的位移函数（2-13）求出单元的应变分布和节点位移的关系：（2-15）式中：称为单元的位移应变转换矩阵，简称应变矩阵。（2-16）由一维杆的应力应变关系（2-2），得单元应力和单元节点位移的关系：（3）用弹性体的虚位移原理导出杆单元刚度方程变形体的虚位移：假想在弹性体上发生的，满足位移许可条件（内部连续，边界协调）的微小、任意位移场。可以理解为某个位移场的微小扰动（变分）。虚位移的特征：1）假想的，与真实位移无关；2）几何上是许可的（可能的）：连续、协调；3）微小、任意大小。虚位移原理：弹性体受力平衡时，若发生虚位移，则外力虚功等于弹性体的虚应变能（应力在虚应变上做的虚功）。下面把虚位移原理应用在所研究的杆单元上。单元虚位移：节点虚位移：定义杆单元的虚位移：节点虚位移单元虚位移单元虚应变（2-17）单元虚应变能：那么，节点力虚功为：单元虚应变：对杆单元应用虚功原理，那么上述节点力（外力）虚功等于虚应变能，因此有下列关系：（2-18）考虑到的任意性，从上式可以得到：（2-19）上式就是杆单元的刚度方程，杆单元的刚度矩阵为：这就是单元刚度矩阵的通式，其导出原理和计算方法可推广到其他类型的单元。具体计算如下：（2-20）显然，与前面直接法得到的单元刚度矩阵（2-8）式相同。3）杆单元讨论a. 只有拉伸、压缩变形的杆单元在局部坐标系下是一维问题，2节点单元只有2个节点位移分量单元有2个自由度，单元刚度方程、刚度矩阵为2阶。b. 单元刚度矩阵元素的物理意义：（2-21）设单元刚度方程为：（2-22）令：（2-23）上式代入（2-21）得到：上式表明，单元刚度矩阵第一列元素就是当单元节点位移满足式（2-22）时的单元节点力分量。如果能设法求出此时的节点力，就得到第一列的刚度系数。因此，一般地，单元刚度矩阵的第i（这里i=1,2）列元素等于当维持单元的第i个自由度位移为，其它自由度位移为0时，施加在单元上的所有节点力分量。可以用此方法直接导出杆单元的刚度矩阵元素，试练习。c. 显然，杆单元的刚度矩阵对称、奇异、主对角元恒正。4）例题例1求图示杆中的应力。解：结构划分为个杆单元，单元之间在节点铰接。根据杆单元刚度矩阵的公式分别写出两个单元的刚度矩阵为：图 2-2参照前面弹簧系统的分析方法，装配杆系统的有限元方程（平衡方程）如下：（2-24）（2-25）考虑图2-2中的约束和载荷情况后，方程（2-24）变化为：则上式第2个方程为：（2-26）（2-27） 求解该方程后得到系统的位移解：计算单元应力：单元1:单元2:提示：1） 本例中单元应力的计算采用了材料力学中的方法,与采用有限元单元应力公式的结果相同，请验证。2） 对锥形杆，单元截面积可以用平均值，从而转化为类似本题的问题求解。3） 求应力之前需要先求出节点位移，因此本方法称为有限元位移法。4） 如果杆上受连续分布的轴向载荷或单元内部受轴向集中载荷，分析时可以按照虚功相等的原则先把单元上的载荷等效移置到节点上。例2 如图2-3所示直杆结构有限元模型。已知：图 2-3求：杆两端的支反力。解：单元、节点的定义如图2-3。先检查杆右端（节点3）与墙壁是否接触。先计算右端的自由伸长：根据上式可以判断：结构受力平衡时，右端间隙将闭合，节点3与刚性墙壁接触。参考前面的讨论，直接写出系统平衡方程为：（2-28）引入下列载荷和位移边界条件：则有限元平衡方程（2-28）成为：（2-29）分离出第二个方程：（2-30）即：（2-31）解得：全部位移解为：（2-32）根据上式位移解，从系统平衡方程（2-28）的第1，3个方程分别求出支反力如下：解毕。2、 2-D和3-D空间中的杆单元（平面和空间桁架单元）1）2-D空间中的杆单元2-D空间中建立杆单元的基本思路是根据前面在杆的一维局部坐标系下建立的单元特性方程通过坐标变换，转换为2-D总体坐标系下的方程，同时得到坐标变换后的单元刚度矩阵。而系统整体分析的原理和方法与一维情况相同。图 2-4（1）变换图2-4为一个杆单元及其局部坐标系与2-D总体坐标的关系。节点的位移分量和节点力分量在2-D局部坐标系x-y下描述，杆节点i具有2个自由度：位移分量为，；节点力分量为，其中只有x方向的位移分量和节点力分量用来描述单元特性。节点上的位移和节点力向量在2-D局部坐标系与2-D总体坐标系下的坐标变换如下：（称为方向余弦。）上述变换的矩阵形式：（2-33）（2-34）符号形式： （2-35）其中（2-36）称为向量的坐标变换矩阵，是单元特性坐标变换的基本元素。显然是正交矩阵，即：因此，由（2-33）可得单元节点位移列阵的坐标变换式如下： （2-37）（2-38）或：（2-39）其中（2-40）比照（2-38）得到单元节点力的坐标变换式：（2）2-D空间刚度矩阵下面可以导出单元刚度方程和单元刚度矩阵的坐标变换式。已经知道杆的一维局部坐标系下的刚度方程为：（2-41）（2-42）把该方程扩充到2-D局部坐标系x-y下的4阶形式：（2-43）符号形式：（2-44）引入单元列阵变换式（2-38），（2-40）得：考虑到变换矩阵的正交性，得到：（2-46）（2-45）或：（2-47）其中：式（2-46）就是2-D总体坐标系下杆单元的刚度方程，是二维空间杆单元刚度矩阵，其计算式如下：（2-48）其中方向余弦由杆节点整体坐标求得。（3）单元应力 由单元应力计算公式（2-16）和位移向量变换得：即：（2-49）（4）思考与讨论 如何从二维空间总体坐标系下杆单元刚度方程（2-46），根据刚度矩阵元素的物理意义，直接导出总体坐标系下杆单元刚度矩阵计算式（2-48）？2）例题如图2-5，平面桁架由2根相同的杆组成，已知E，A，L。求：1）节点2位移2）每根杆应力解：图 2-5求出每个单元在总体坐标下的刚度矩阵：单元1，（i,j）=(1,2)：（2-51）（2-50）单元2，（i,j）=(2,3): 上面的刚度矩阵也可以按（2-48）直接计算。将上述总体坐标系下的单元刚度矩阵扩充到整体规模，装配得到系统有限元方程：（2-52）引入边界约束和载荷：（2-53）则系统有限元平衡方程凝聚为：（2-54）解得：（2-55）下面按公式（2-49