用有限元素法建立结构振动的数学模型

1、第三章 用有限元素法建立结构振动的数学模型31 引言【工程要求】：对于简单的连续结构，如单件的杆、板、梁，可以建立结构振动的偏微分方程，但对于杆、板、梁组成的复杂结构，仍然采用建立偏微分方程的方法则十分困难。如果用假设模态法（李兹方法），对实际工程结构假设出品质良好的整个结构的假设模态也十分困难。要对结构振动进行数值分析，必须建立振动的数学模型振动方程。工程结构振动分析中，要采用将结构离散为有限自由度系统的方法有限元素法，来建立结构的数学模型。【发展简况】有限元素法，是在上一世纪五十年代中期，经过M.T.Turner及J.H.Argyris等人的开拓性工作以及后来许多研究者的大量工作，发展起来

2、的一种结构分析的有效方法，上一世纪六十年代初，由J.S.Archer及J.H.Argyris等人引入到结构动力学分析中来。有限元素法发展到今天，已经非常成熟，而且与先进的计算机技术结合，已经形成了一个以有限元分析方法为基础的计算机辅助工程（CAE）的技术领域以及更进一步的虚拟产品设计（VPD）这样的先进概念。世界上著名的CAE分析软件商主要有MSC.software和Ansys等公司的产品。【有限元动力学分析的任务】在结构振动分析领域，有限元素法处理的问题主要是两类：结构固有振动特性计算和结构振动响应计算（包括频率响应分析与响应时间历程分析）。两类问题中，用有限元法建立振动数学模型是最基础的工

3、作。【有限元素法（分析结构振动问题）的特点】：原则上，有限元素法由于其对复杂边界的适应性，它可以处理任何复杂的结构。求解结果的精度可以根据需要不断改善，建模过程规范统一，计算形式适合于计算机求解。【存在的问题】：随着精度要求的不断提高，所要求的计算机容量和计算时间急剧增加，从而引出了大型特征值问题的快速求解方法、将大型结构振动问题转化为若干小型结构振动问题集合的子结构求解方法，以及结构振动问题的并行求解方法等问题的研究。【工程结构振动分析方法】从结构振动分析的发展历史看，经典的方法有：1集中质量法将质量分别集中在若干节点处，形成集聚质量阵。结构的刚度仍然连续分布，采用材料力学中求柔度的方法，求

4、出柔度系数，得到柔度矩阵，即用柔度法来形成刚度矩阵。集中质量法存在问题：对大型复杂结构，用材料力学的方法，进行柔度矩阵的求解显然是不现实的。2假设模态法以李兹法为基础，选择一组假设模态组成的模态矩阵，对结构进行离散，如第二章所述的方法。假设模态法存在的问题：1对几何形状复杂的结构，假设模态难以选择。2对整个系统用假设模态法得到的运动方程是高度耦合的，求解困难。3对不同的结构，要根据实际情况选取不同的假设模态，求解过程不规范统一。引入有限元素法的思想既解决了上述方法的缺点，又保留了它们的优点。【有限元法分析振动问题的基本原理】用有限元法分析结构动态问题的基本思想，与结构静态分析的思想是一样的。它

5、采用的方法仍然是：将结构分解为有限数目的单元，各元素间由节点相连，各单元内结构的变形用位移形函数（相当于元素级的假设模态）来表示，以节点位移作为控制变量（元素的广义坐标）。元素间的位移连续条件通过引入的形函数来满足，动态平衡条件通过最后导出的有限元方程来体现。由于节点数目是有限的，最后得到的方程是一个多自由度、离散的、线性的矩阵微分方程。32 运动方程的建立仍然采用熟悉的拉格朗日方程法建立其数学模型（运动方程）。对任一单元内部任一点的位移与节点位移的关系：（31）N称为假设的已知位移形函数（可以看成是单元的假设模态，一般仍采用静态变形函数）显然： （32）单元的动能：（33）称为单元质量矩阵，

6、质量阵是对称矩阵。整个结构的动能为： （34）是全结构的节点位移列阵，为全结构质量阵，代表对整个结构各单元的组集（Assemble）。单元的应变向量:， （35）几何矩阵，弹性矩阵，应力向量单元的势能为：（36）全结构的势能： （37） （38）作用在单元上的分布力的虚功： （39）单元节点力（广义力）（310）全结构的外力虚功：（311） （312）【阻尼的处理】采用粘性阻尼假定：阻尼力与运动速度成正比，方向与速度相反。单元中分布阻尼的耗散函数（瑞利耗散函数）： （313）耗散力（即瑞利耗散力）与耗散函数的关系为： （314）全结构的耗散函数： （315） （316）将全结构的动能、势能、耗

7、散函数和广义力代入非保守系统的拉格朗日方程， （317）得到： （318）【几个相关问题】：1进行实际结构的振动分析时，在各个单元的矩阵组集之前，还要对单元矩阵进行由单元的局部坐标系向结构的总体坐标系转换。记局部坐标系下节点位移向量向总体坐标系下节点位移向量的转换阵为，则坐标转换关系为： （319）其中为对单元矩阵组集时“对号入座”的定位矩阵。（320）（321）（322） （323）经过上述变换后的单元矩阵可以直接叠加得到结构总体矩阵。2为了求解结构的固有振动特性，需要求解无阻尼情况下结构的自由振动方程：(324)将固有振动的简谐运动形式（325）代入得到的结构的特征方程：或 （326）数学

8、上构成所谓的广义特征值问题。3振动分析中采用的质量阵问题在结构振动分析中，常采用的质量阵形成方法有：集中质量模型和一致质量模型。在采用集中质量模型时，一般是按照杠杆原理将单元质量向单元各个节点上进行分配，在局部转动效应显著时，还要考虑单元的转动惯量。集中质量模型得到的质量阵为对角矩阵。将单元内惯性分布视为与静力形函数同样规律的分布，导出的质量矩阵称为一致质量阵。即上面（）推导出的质量阵。这样的质量阵为满阵。注意：集中质量阵与一致质量阵都不是振动结构在实际上精确的质量分布模型。理论上，结构的动位移是与频率相关的。动位移在不同振动频率和振型下是不同的，即不同的频率对应有不同的惯性。所以严格地讲，质

9、量阵也是与频率相关的。下面以轴向振动的杆元为例，说明这个问题。作为连续体的二力杆元的振动偏微分方程（波动方程）为：xL12 （327）记称为波速，分别为杆的弹性模量和密度。方程（327）的解为：（328）故此时形函数为: （329）代入单元质量阵公式得：（330）显然，这个质量阵为满阵，且各元素为与频率相关的量。理论上这样的质量阵真实反映实际的惯性分布，计算得到的固有特性更精确，但却使计算大大复杂化，而且只有对简单的情况（如二力杆元）下，才能从偏微分振动方程得出（动力）形函数，因此，工程实际中很少采用这种方法。而是采用与推导刚度阵时一致的（静力）形函数。至于采用集中质量阵还是一致质量阵，得到的

10、结果更好，没有固定的规律和结论，要视具体情况而定。从一般经验上讲，在单元划分较细时，用集中质量阵较好，反之宜采用一致质量阵。根据特征值隔离定理知道，采用一致质量阵分析得到的是固有频率精确解的上界，随着分元的细化，计算结果单调地向精确解逼近，而集中质量阵给出的结果就不具备这种特性，可能偏低也可能偏高。从工程分析经验看，由于建立有限元模型的离散过程，已经使结构比实际结构的刚度增大，因此采用集中质量阵可能有时反而会得到误差较小的结果，但这需要经验和技巧。但采用集中质量阵得到的振型一般误差较大，对振型要求较高时，还是宜采用一致质量阵。集中质量阵是对角阵，在计算时可以节省计算时间。我们希望能获得一种优于

11、集中质量阵的对角化质量阵。例如，对于梁的弯曲振动，可以按下式来计算对角化的非一致质量阵的单元： （331）而非对角元全部置零，是梁的四个形函数。由于这种质量矩阵在一定程度上反映了单元形函数的特性，因此可以给出精度较好的计算结果。33 典型结构单元的有限元建模结构有限元模型的自由度数节点数×节点位移数一、 纵向振动杆元杆元的任意截面处位移由两节点位移插值得到。 （332）应满足边界条件： （333）杆元受节点力时的静态方程为： （334）其解为： （335）代入边界条件得到形函数： （336）从而杆元的质量阵和刚度阵为： （337） （338）二、 弯曲振动梁元梁元的任意截面处位移由两

12、节点位移插值得到。（339）为梁元的节点位移。均匀梁元受常值节点力作用时的挠曲线偏微分方程为： （340）从而 （341）应满足边界条件： （342）代入边界条件得到： （343）代入位移表达式，整理后，得到形函数为： （344）形函数阵： （345）梁元的质量阵和刚度阵为： （346） （447）上面给出的是平面梁元的特性矩阵，对于空间梁元特性矩阵形成过程完全相同。可以参考任何一本有限元素法方面的专著或教材。这里不再赘述。【单元的坐标变换阵】：单元局部坐标与总体坐标不一致时，需要进行坐标变换，记局部坐标与总体坐标之间夹角为，则局部坐标下节点位移列阵与总体坐标下位移列阵间的变换关系为： (34

13、8)则坐标变换阵为： (349)三面内振动的平板元以三节点三角元为例，取三角元三个顶点的位移为单元广义坐标， （350）单元内任一点的位移为：（351）三个形函数应满足边界条件： （352）设满足此条件的形函数为： （353）代入边界条件可得到： （354） （355）为三角形单元的面积。由上式按下图轮换下标求得：（356a）几何矩阵：（356） （357）对平面应力问题： (358a) （358） （359）质量阵和刚度矩阵为： （360）记板厚为，且由于都是常数矩阵，则 （361）的分块形式为 （362）节点刚度阵（363）单元质量阵（364）四弯曲振动的平面板元仍采用三节点三角元，引入局

14、部坐标系，原点在单元的1节点，轴沿单元的一条边，轴垂直于单元平面，如图所示。单元的三个节点共有九个节点位移条件，故设单元内任意一点的离面垂直位移表达式为： （365）或 （366）其中 （367）（368）单元的节点位移边界条件为： （369） （365）代入（369）式得到： （370） （371）显然，各节点的离面位移垂直位移（挠度）、绕轴的转角、绕轴的转角，就是单元特性分析时所用的广义坐标。节点的位移列阵为： （372）从而有：， （373）因此，弯曲板元的形函数为： （374）按薄板弯曲问题的基本假定（直法线假定），板内各点线位移为： （375）应变列阵： （376） （377） （3

15、78）其中 （379）为微分算子列阵。记（380）则 （381）从而得到弯曲板元的几何矩阵：（382）单元刚度阵 （383）（384）为弹性矩阵。单元质量矩阵为： （385）（具体的计算结果略）积分后单元的刚度矩阵和质量矩阵表达式可参见有限元素法的专著。其它形式的单元的特性矩阵，均可按上述过程推导出。当然，在实际结构振动有限元分析中，如果自己编程，不需要自行推导，可以查阅有关参考书籍。应该承认，在各种关于结构有限元静力分析书籍中，对有限元的介绍非常详尽，也给出了各种高精度的单元形函数和刚度矩阵，但一般有限元素法的书籍中没有给出单元质量矩阵。如果采用商用有限元分析软件，则在选择单元类型时，要根据

16、具体结构、对分析结果的精度要求以及网格划分的情况，来选择最合适的单元。这些问题在有限元著作中都有详细论述，这里也不再赘述。34结构振动分析有限元建模中的几个问题【阻尼矩阵的工程处理】在对结构进行振动响应分析时，必须计及阻尼的影响，然而由于实际结构的阻尼特性比较复杂，即使我们可以借助瑞利粘性阻尼假设，得到带粘性阻尼结构的振动方程（386）但是，结构的阻尼特性分布参数无法从理论上确定，因而就无法用理论方法来计算矩阵中的元素。在工程振动分析中，常常是采用瑞利比例阻尼假设，即假定阻尼阵是质量阵和刚度阵的线性组合： （387）比例常数通过对结构或类似结构，进行模态试验，测得其第阶和第阶模态频率和阻尼比来

17、确定。理论上可以从下面方程解出：（388）一般是测定多阶模态参数，通过最小二乘法求解出。【刚度矩阵奇异性的处理】由于采用有限单元法，通过对各单元刚度矩阵的装配，得到的全结构总体刚度阵一般是无约束结构的刚度矩阵，它是奇异的，即不能求逆。因此需要按照结构实际的约束边界条件进行处理，在对刚度矩阵进行处理时，有四种方法：1 第一种方法：当给定的边界约束条件足以使结构成为静定或超静定结构时，只需在进行全结构的总刚度阵装配时，抽掉刚度矩阵中对应于结构节点位移为零的行和列，即得到一个非奇异的降阶的刚度矩阵，对质量阵进行相应的处理。这样处理的缺点是不能进一步求出约束反力。2 第二种方法：用无限大或计算机允许的

18、最大数来代替刚度阵中对应于零位移节点的主对角元素，从而使刚度阵成为非奇异的可逆矩阵，并使对应的节点位移具有一个非常小、十分接近于零的解。这样处理的缺点是未能利用边界节点约束降低结构矩阵阶数而不能减小计算工作量。3 第三种方法：对刚度矩阵和质量矩阵按节点位移是否为零进行分块，如对结构的自由振动方程 （389）为约束反力列阵，下表，分别代表自由节点和约束节点。将矩阵方程分块展开得到： （390）解第一式的特征值问题得到及，再代入第二式得到约束反力。显然由这样的方法可以得到结构的全部信息，但是增加了编程计算的难度。4 当结构属于自由自由结构，或结构具有刚体运动模态时，我们不可能根据结构的边界约束条件

19、来消除刚度矩阵的奇异性，从而给特征值问题的求解带来困难。为了消除刚度矩阵的奇异性，可以利用特征值问题的移轴特性，将刚度矩阵转换成正定的矩阵。如下所示： （391）为一个任意的、适当的正常数，上方程中是一个正定矩阵，可以按通常的方法求解，得到的频率与原来结构的固有频率间具有关系： （392）而振型不变，就是由上方程解得的。值得注意的是，由于结构具有个刚体模态，对应的固有频率为零，所以，上方程的解也包括个刚体模态，对应着个重频（这里用了二阶小量符号，因为一般数值求解的结果不一定得到理论上的零），在求解时可以用“扫模法”将这个刚体模态扫去后再用一般的求解特征值方法（程序）求解。【自由度的静力缩聚法】为了提高计算精度，在结构的有限元建模时，对结构可以采用网格细分的方法得到较精细的有限元网格；也可以采用较粗的网格划分，而采用带有内节点的高阶元素。在采用这种高阶元素时，内节点参数与单元间的位移连续条件是无关的，与单元本身以外的节点参数也没有直接联系。因此可以用某种方法在组成全结构矩阵前，将其消去以减少结构的总自由度数。还有一种情况，如对梁弯曲振动问题进行分析时，往往只保留节点横向位移参数，而略去节点的转角参数，也需要在组成全结构的矩阵后，用节点横向位移参数来表示节点转角参数，从而达到特征值问题降阶的目的。此外，许