第9章(动力学fem)

1、公式号、图号第九章 机械结构的动力有限元本章介绍介绍机械结构动力学的有限元建模方法，然后利用有限元求解机械结构系统的固有动特性，对典型工况下的动态响应进行分析，并介绍了几个典型应用实例。8.1 动力学方程的建立在动力学问题中，位移、速度、应变、应力和载荷都是时变的。弹性结构的动力学问题的基本方程如下平衡方程 (4.70)几何方程 (4.71)物理方程 (4.72)边界条件 , (在Su边界上)、, (4.73)以及初始条件 。式中，是密度，是阻尼系数，分别是对t的二阶导数和一阶导数，即加速度和速度。分别表示惯性力和阻尼力，作为体积力的一部分出现在平衡方程中。只对空间域进行离散，单元内位移的插值

2、函数分别为 (4.74)其中，u。平衡方程式(4.71)以及力边界条件的等效积分形式的Galerkin提法如下 (4.75)对上式第一项进行分部积分，并引入物理方程，则由上式可以得到 (4.76)将位移空间离散后的表达式(4.74) 代入上式，注意这里的相当于u, v, w，并根据变分原理，最终得到结构系统的动力学方程 (4.77)其中，和分别是结点加速度向量和速度向量，和分别是结构系统的质量矩阵、阻尼矩阵、刚度矩阵和结点载荷向量，分别由各自的单元矩阵和向量集成 ， , (4.78)其中，, 分别是单元的质量矩阵、刚度矩阵、阻尼矩阵和单元载荷向量。4.3.2 质量矩阵和阻尼矩阵式(4.79)所

3、表达的单元质量矩阵 （4.79）称为协调质量矩阵或一致质量矩阵，这是因为公式导出时采用了和导出刚度矩阵一致的位移插值函数，其质量分布是按照实际分布情况考虑的。式中的形函数矩阵N随单元类型而异。例如，对于平面三角形常应变单元，其位移插值函数为, 其中，I是单位矩阵，系数见前文，A是三角形单元面积。利用上式可以算得平面三角形常应变单元的协调质量矩阵具体表达式为 如果是等参元，设其形函数矩阵为，单元协调质量矩阵为 (4.80)为计算方便，还经常采用集中质量矩阵，即假设单元的质量集中在结点上，即把每个单元的分布质量按静力学中的平行分解原理，平均分配到每个结点上，形成一个阶数等于单元自由度数的对角线质量

4、矩阵，而非主角线元素均为0。例如，对于平面三角形常应变单元，其集中质量矩阵为 (4.81)对于式(4.78)所表达的单元阻尼矩阵 （4.82）称为协调阻尼矩阵，它是假设阻尼力正比于质点运动速度的结果。通常将介质阻尼简化为这种情况，这时单元阻尼矩阵比例于单元质量矩阵。除此之外，还有比例于应变速度的阻尼，例如由于材料内摩擦引起的结构阻尼通常可简化为这种情况，这时的阻尼力可表示成，可以得到单位阻尼矩阵 （4.83）此单元阻尼矩阵比例于单元刚度矩阵。由于系统的固有振型对于和具有正交性，因此固有振型对于比例于和的阻尼矩阵也具有正交性。所以这种阻尼矩阵称为比例阻尼矩阵或振型阻尼矩阵。通常允许将实际结构的阻

5、尼矩阵简化为和的线性组合。即 （4.84）其中是不依赖于频率的常数。这种振型阻尼称为Rayleigh阻尼。8.2 利用有限元分析机械结构的固有动特性4.4.1 机械结构的固有频率和固有振型不考虑阻尼影响的机械结构的自由振动方程是 (4.85)它的解可以假设为以下形式 (4.86)其中,是阶向量，是向量振动的频率，是时间变量，是由初始条件确定的时间常数。将式(4.86)代入式(4.85)，就得到一广义特征值问题 (4.87)求解以上方程可以确定和，结果得到个特征解，其中特征值，代表系统的个固有频率，并有.对于结构的每个固有频率，由式(4.87)可以确定出一组各节点的振幅值，它们互相之间保持固定的

6、比值，但绝对值可任意变化，所构成的向量称为特征向量，在工程上通常称为结构的固有振型。设特征向量，代表结构的个固有振型，它们的幅度可按以下要求规定 (4.88)这样规定的振型又称为正则振型，即所谓的固有振型。固有振型具有如下重要性质。将特征解代回方程式(4.86)，得到 (4.89)上式前一式两端前乘以，后一式两端前乘以，并由和的对称性推知 (4.90)所以可以得到 (4.91)由上式可见，当时，必有 (4.92)上式表明固有振型对于矩阵是正交的。和(4.88)式在一起，可将固有振型对于的正则正交性质表示为 (4.93)进而可得 (4.94)定义固有振型矩阵，则 (4.95)特征解的性质还可表示

7、成 (4.96)式中，和分别称为固有振型矩阵和固有频率矩阵。由此，原特征值问题还可以表示成 (4.97)4.4.2 固有频率和固有振型的求解方法结构振动的固有频率和固有振型求解是振动模态分析的关键问题。求解固有频率和振型的方法有振型截断法、矩阵逆迭代法、李兹法、广义雅可比法等。对于一个连续体的结构，其固有频率有无限多阶，在有限元中，结构被离散成细小的单元，对于大型复杂的结构，单元的数目可能数以万计，由这些单元形成的振动方程组规模是很庞大的，由此而得到的特征方程，其矩阵的阶次通常很高。有限元中经常用来求解结构的低阶模态。另外，同样规模的特征值问题, 其计算量比静力问题的计算量要高出几倍，因此，如

8、何降低特征值问题的计算规模、减少计算量也是一个重要的课题。如下介绍振型截断法即所谓的Guyan缩聚的基本原理。对于一个机械结构，设其静力问题的总体刚度方程为 （4.98）其中, r为总自由度。将其写成分块矩阵的形式，有 （4.99）其中, 。上面的分块是根据节点位置的重要程度来划分的，一般情况下，将对应于结构关键位置的节点位移划分为，叫做“主自由度”（master DOF），而剩下的节点位移划分为，叫做“从自由度”（slave DOF）。假定进行分块时考虑到的特征（比如结构内部无载荷，或对应于内部自由度的载荷为零），同时还可以定义矩阵对角线的相对刚度系数来确定“主自由度”，即计算 （4.100

9、）其中，分别为系统的刚度矩阵和质量矩阵主对角线上的元素。当相对刚度系数超过某一临界值时，所对应的节点自由度选为“主自由度”，即 （4.101）其中，为临界相对刚度系数值，可以根据需要缩聚的自由度数来确定。方程（4.99）可以写为 (4.102)由上式中的第二式，有 (4.103)将q写成 (4.104)其中转换矩阵T为 (4.105)将该转换关系（4.105）代入到求动力学问题的虚功方程中，可以得到动力学系统的方程 (4.106)即有 (4.107)其中 (4.108)考虑无阻尼自由振动，则 (4.109)该方程的解为 (4.110)其中为系统的固有频率，将式（4.110）代入到式（4.109

10、）中，有 (4.111)方程（4.111）就是缩聚后的自由振动方程。求出“主自由度” 后，由（4.104）求出所有节点的位移。以上的缩聚是基于静力问题中对应于“从自由度”上无外载的情形下推导的，即得到缩聚关系式（4.103）。对于动力学方程，同样可以写出完全的分块矩阵方程，即 (4.112)对于（4.112）中的第二组方程，有 （4.113）其中为惯性力。与式(4.102)中的第二个方程进行比较，式（4.112）中由于存在惯性力，因为它不为零，由式（4.113）得到的为 (4.114)在获得“主自由度”的振型后，由（4.113）式来求取关于“从自由度”的振型和总的振型为。由于此时的转换关系（4

11、.114）式与原转换关系（4.103）式有差别，计算出的自然频率和振型会有一定误差，高阶振型的误差更大。4.4.3 机械结构固有频率和振型的计算举例对机械结构进行有限元动特性分析时，首先在结构的几何模型的基础上，建立有限元模型。建立模型的原则是在保证准确反映结构模态特性的前提下进行必要的简化，在获得必要的精度的同时，尽可能降低计算量与计算复杂性。例如一矩形钢板，长500mm，宽100mm，高12mm，一端固定在某支架上，可以视为一悬臂梁结构。进行有限元模态分析，我们获得了其前四阶弯曲模态的固有频率，38.5, 241.7, 678.9, 1335.6 Hz，这四阶弯曲模态对应的振型如图4-5所

12、示。图4-5 有限元分析得到的某悬臂梁振型图如图4-6所示是某发动机压气机二级叶片利用有限元计算得到的固有频率和振型图。 (a) 5阶2747Hz (b) 7阶4019Hz (c) 10阶6112Hz图4-6 有限元计算得到的某压气机二级叶片的模态振型图8.3 机械结构的动态响应分析通常把结构动力分析分为结构的瞬态响应分析和结构的基础响应分析。动力响应问题就是求解动力学方程(4.77)，即在的作用下，求出作为时间函数的，和. 根据所用方法不同，前者又有振型叠加法和逐步积分法；根据结构的基础加速度性质不同，后者又有频率响应分析和谱分析。4.5.1 直接积分法直接积分是将时间的积分区间进行离散化，

13、计算每一段时刻的位移数值。通常的直接积分法是从两个方面解决问题，一是将在求解域内的任何时刻都应满足运动方程的要求，代之以仅在一定条件下近似地满足运动方程，即，将域内每点都满足的微分平衡方程转化为只在每个节点处满足的节点平衡方程。例如可以仅在相隔的离散时间点满足运动方程。二是以在单元内分片连续的已知变化规律的位移函数，代替全域内连续的未知函数从而，将通过微分平衡方程求全域内连续的未知函数问题转化为通过节点平衡力求节点未知位移的问题。在以下的讨论中，假定时间的位移，速度，加速度已知，并假设时间求解域被等分为个时间间隔。在讨论具体算法时，假定时刻的解已经求得，计算的目的在于求时刻的解，由此建立求解所

14、有离散时间点解的一般算法步骤。中心差分法是一种显式算法，就是由上一时刻的已知计算值来直接递推下步结果，在给定的时间离散步中，逐步求解各个时间离散点的值。其中，加速度和速度可以用位移表示 （4.115） （4.116）时间的位移解是，可由下面时间的运动方程应得到满足而建立，即 （4.117）为此将式(4.115)和(4.116)带入上式，得到 （4.118）如已经求得和，则从上式可以进一步解出。所以上式是求解各个离散时间点解的递推公式，这种数值积分方法又称逐步积分法。但是当时，为了计算，除了从初始条件已知的外，还需要知道，所以必须用一种专门的起步方法。为此利用式(4.115)、(4.116)可以

15、得到 （4.119）上式中可从给定的初始条件得到，而则可以利用时的运动方程(4.77)得到。中心差分法求解运动方程的算法具体步骤如下：a. 初始计算：(1) 形成刚度矩阵，质量矩阵M和阻尼矩阵C；(2)给定，和；(3)选择时间步长，并计算积分常数，；(4)计算；(5)形成有效质量矩阵 ；(6)三角分解 ；b. 对于每一时间步长(1)计算时间的有效载荷（2）求解时间的位移（3）如果需要，计算时间的加速度和速度除了上述显式积分算法之外，常用的 Newmark积分方法则是一种隐式算法，这面加以介绍。首先假设 (4.120) (4.121)其中,和是按积分精度和稳定性要求而设定的参数。当和时，(4.1

16、20)和(4.121)式对应于线性加速度法，此时它们可以从下面时间间隔内线性假设的加速度表达式的积分得到 （4.122）式中。Newmark方法是从常平均加速度法这种无条件稳定积分方案提出的，要求和。这时，内的加速度为 (4.123)不同于中心差分法，Newmark方法中的时间的位移解是通过满足时间的运动方程 (4.124)得到的。为此, 首先从式(4.121)解得 （4.125）将上式代入（4.120）式，然后再一并代入(4.124)式，则得到从,计算的公式（4.126）采用Newmark方法逐步求解运动方程的算法具体步骤如下：a.初始计算:（1）形成刚度矩阵，质量矩阵和阻尼矩阵；（2）给定

17、给定，和；(3) 选择时间步长，参数和，并计算积分常数:， (4)形成有效的刚度矩阵：；(5)三角分解：；b. 对每一个时间步长:(1)计算时间的有效载荷(2)求解时间的位移c. 计算时间的加速度和速度:从Newmark方法的循环求解方程式(4.126)可见，有效刚度矩阵中包含了，而一般情况下总是非对角矩阵，因此在求解时，的求逆是必须的(而在线性分析中只需计算一次)。这是由于导出(4.126)时，利用了时刻的运动方程。因此，这种算法称为隐式算法。当时，Newmark方法是无条件稳定的，即时间步长的大小不影响解的稳定性。此时的选择主要根据解的精度，具体说可根据对结构响应有主要贡献的若干基本振型的

18、周期来确定，例如可选择(对应若干基本振型中的振动周期中的最小者)的若干分之一。一般说比结构系统的最小振动周期大得多，所以无条件稳定的隐式算法以求逆为代价，但与有条件稳定显式算法相比，可以采用大得多的时间步长，而且采用较大的还可滤掉高阶不精确特征解对系统响应的影响。4.5.2 振型叠加法结构在随时间变化的结点力作用下，由于系统中可能存在着各种阻尼(材料阻尼、滑移阻尼、介质粘性阻尼等等)，各结点产生有阻尼的强迫振动。因此，与系统初始条件有关的自由衰减振动，总要随时间增长而消失，只保留稳态的强迫振动。求得结构系统稳态的强迫振动解或稳态响应，并进一步算出动应力响应，是动力问题有限元法的核心内容之一。振

19、型叠加法在一定条件下比直接积分法的计算效率更高。在积分运动方程以前，利用系统自由振动的固有振型将几何坐标下的方程组转换为个正则坐标下的相互不耦合的方程，对这种方程可以解析或数值积分。当采用数值方法时，对于每个方程可以采取各自不同的时间步长，即对于低阶振型可采用较大的时间步长。这两者结合起来，相对于直接积分法是有很大的优点，因此当实际分析的时间历程较长，同时只需要少数较低阶振型的结果时，采用振型叠加法是十分有利的。如前所述，在求得系统的固有频率和固有振型、并将运动方程转换到正则振型坐标系后，进行如下位移基向量的变换 (4.127)其中，。此变换的意义是将看成的线性组合，可以看成是广义的位移基向量

20、，是广义的位移值。从数学上看，是将位移向量从以有限元系统的结点位移为基向量的维空间转换到以为基向量的维空间。将此变换代入运动方程，两端前乘以，并注意到的正交性，得到新基向量空间内的运动方程 (4.128)初始条件也相应的转换成 (4.129)阻尼矩阵如果是振型比例阻尼矩阵，也可从的正交性得到 (4.130)亦即 (4.131)其中，是第阶振型阻尼比。在此情况下，(4.128)式就成为个互相不耦合的二阶常微分方程 (4.132)上式每一个方程相当于一个单自由度系统的振动方程，可以方便地求解。式中，是载荷向量在振型上的投影。若是按一定的空间分布模式随时间变化，即 (4.133)则有 (4.134)

21、上式中引入符号表示空间坐标，表示在上的投影，是一常数。单自由度系统的振动方程(4.132)的求解，在一般情况下可采用上节讨论的直接积分方法。但在振动分析中常常采用杜哈梅积分，又称为叠加积分。这个方法的基本思想是将任意激振力分解为一系列微冲量的连续作用，分别求出系统对每个微冲量的响应，然后根据线性系统的叠加原理，将它们叠加起来，得到系统对任意激振的响应。杜哈梅积分的结果是 (4.135)其中, ，是由起始条件决定的常数。上式右端前一项代表引起的系统强迫振动项，后一项代表在一定起始条件下的系统自由振动项。当阻尼很小，即时，这时杜哈梅积分的结果是 (4.136)杜哈梅积分(4.135)或(4.136

22、)式，在一般情况下，也需要利用数值积分法计算，但是对于少数简单情况，则可得到解析的结果。在得到每个振型的响应以后，按(4.137)式将它们叠加起来，就得到系统的响应，亦即每个结点的位移值 (4.137)在叙述了振型叠加法的算法步骤以后，下面再对此方法的一些性质和特点作一定的分析。首先应看到将系统位移转换到以固有振型为基向量，对系统的性质并无影响，而是以求解广义特征值问题为代价，得到非耦合的个单自由度系统的运动方程，以达到提高计算效率的目的。对于个单自由度系统运动方程的积分，比对联立方程组的直接积分节省计算费用。另外，通常只要对非耦合运动方程中的一小部分进行积分。例如只要得到对应于前个特征解的响

23、应，就能很好地近似系统的实际响应。这是由于高阶的特征解通常对系统的实际影响较小，且有限元法得到的高阶特征解和实际相差也很大(因为有限元的自由度有限，对于低阶特征解近似性较好，而对于高阶则较差)，因此求解高阶特征解的意义不大。而低阶特征解对于结构设计则常常是必要的。但是采用振型叠加法需要增加求解广义特征值问题的计算费用，所以在实际分析中究竟采用哪种方法，应根据具体情况确定。应当指出，如果在振型叠加法中，对于个单自由度系统的运动方程都进行积分，且采用和直接积分法相同的积分方案和时间步长，则最后通过振型叠加得到的和直接积分法得到的结果在积分方案的误差和计算机舍入误差的范围内将是一致的。此外，对于非线性系统通常必须采用直接积分法。因此此时，这样