有限元第五章 有限元动力学基本原理

1、第五章第五章 有限元动力学分析基本原有限元动力学分析基本原理理 一、单元质量矩阵的计算一、单元质量矩阵的计算 1.1.单元一致质量矩阵单元一致质量矩阵 2.2.单元集中质量矩阵单元集中质量矩阵 3.3.常用单元的一致质量矩阵常用单元的一致质量矩阵二、单元阻尼矩阵二、单元阻尼矩阵 1.1.速度阻尼矩阵速度阻尼矩阵 2.2.应变阻尼矩阵应变阻尼矩阵三、机械结构的固有频率和振型三、机械结构的固有频率和振型 1.1.无阻尼自由振动方程无阻尼自由振动方程 2.2.矩阵迭代法矩阵迭代法 3.3.其他方法其他方法 四、机械结构的动力响应计算四、机械结构的动力响应计算 1.1.振型叠加法振型叠加法 2.2.直

2、接积分法直接积分法 在前面的介绍中，我们均假设作用在弹性体（或结在前面的介绍中，我们均假设作用在弹性体（或结构）上的载荷与时间无关，与此相应的，位移、应力构）上的载荷与时间无关，与此相应的，位移、应力及应变等也都和时间无关，即前面介绍的全部内容皆及应变等也都和时间无关，即前面介绍的全部内容皆称结构静力学有限元方法。但工程实际中还存在着另称结构静力学有限元方法。但工程实际中还存在着另外一类载荷与时间有关的动载荷作用于结构或弹性体，外一类载荷与时间有关的动载荷作用于结构或弹性体，此时，相应的位移、应力、应变等都与时间有关，而此时，相应的位移、应力、应变等都与时间有关，而且必须考虑惯性力和加速度等因

3、素，这类分析或问题，且必须考虑惯性力和加速度等因素，这类分析或问题，成为动力学分析。成为动力学分析。 对于质点对于质点弹簧系统的振动，大家比较熟悉，例如弹簧系统的振动，大家比较熟悉，例如一个自由度为一个自由度为n n的质点的质点弹簧振系，其动平衡方程为弹簧振系，其动平衡方程为 PKCM&第五章第五章 有限元动力学分析基本原有限元动力学分析基本原理理 上式中每一项的含义不同上式中每一项的含义不同 为弹性力K 为阻尼力C?&M对于单元体而言，可以得到类似的上述方程对于单元体而言，可以得到类似的上述方程 eeeeeeepkcm&第五章第五章 有限元动力学分析基本原有限元动力学分析基本原理理 单元质量

4、矩阵根据其形成过程分为一致质量阵和单元质量矩阵根据其形成过程分为一致质量阵和集中质量阵，各有自身的优点和缺点。集中质量阵，各有自身的优点和缺点。1.1.一致质量矩阵一致质量矩阵一、单元质量矩阵的计算一、单元质量矩阵的计算 在离散后的结构中，取出一个单元，根据达朗贝在离散后的结构中，取出一个单元，根据达朗贝尔原理，单位体积上作用的惯性力为：尔原理，单位体积上作用的惯性力为： 惯性力是分布力，按分布力向节点等效的原则和惯性力是分布力，按分布力向节点等效的原则和实施过程，有：实施过程，有： eetNNttq222222一、单元质量矩阵的计算一、单元质量矩阵的计算 1.1.一致质量矩阵一致质量矩阵于是

5、，令于是，令 eVTVeTTVeqdVNNdVtNNdVqNR 22 VTedVNNm一、单元质量矩阵的计算一、单元质量矩阵的计算 1.1.一致质量矩阵一致质量矩阵 的计算式是通式，并因为计算质量矩阵和刚度矩的计算式是通式，并因为计算质量矩阵和刚度矩阵使用的形状函数一致，因此被称为一致质量阵。阵使用的形状函数一致，因此被称为一致质量阵。em2.2.集中质量矩阵集中质量矩阵 在工程实际中，为了求解方便，有人把单元质量在工程实际中，为了求解方便，有人把单元质量平均分到单元的各个节点上，如平面三角形单元的平均分到单元的各个节点上，如平面三角形单元的质量可分配为：质量可分配为：VkjidVmmm3一、

6、单元质量矩阵的计算一、单元质量矩阵的计算 2.2.集中质量矩阵集中质量矩阵单元质量矩阵为：单元质量矩阵为：kkjjiiemmmmmmdiagm3.3.常用单元的一致质量矩阵常用单元的一致质量矩阵一次杆单元一次杆单元 21126222121212121ldxAdxAdxNNAmllTle一、单元质量矩阵的计算一、单元质量矩阵的计算 3.3.常用单元的一致质量矩阵常用单元的一致质量矩阵二次杆单元二次杆单元 1688841814304224222122212121222121AldxAdxNNAmTlTle一、单元质量矩阵的计算一、单元质量矩阵的计算 3.3.常用单元的一致质量矩阵常用单元的一致质量

7、矩阵三次梁单元三次梁单元2222422313221561354313422135422156420llllllllllllAlme一、单元质量矩阵的计算一、单元质量矩阵的计算 3.3.常用单元的一致质量矩阵常用单元的一致质量矩阵三角形三角形平面问题单元平面问题单元20210201021010201010212称对tme一、单元质量矩阵的计算一、单元质量矩阵的计算 3.3.常用单元的一致质量矩阵常用单元的一致质量矩阵矩形平面问题单元矩形平面问题单元4042040204102040102042010204020102049称对abtme二、单元阻尼矩阵的计算二、单元阻尼矩阵的计算 阻尼矩阵非常复杂

8、，主要是阻尼本身的复杂性引阻尼矩阵非常复杂，主要是阻尼本身的复杂性引起的，一般均为假设，如阻尼力正比于单元的运动起的，一般均为假设，如阻尼力正比于单元的运动速度，此时得到的阻尼矩阵正比于单元质量矩阵；速度，此时得到的阻尼矩阵正比于单元质量矩阵；也可以假设阻尼力正比于单元的应变速度，此时得也可以假设阻尼力正比于单元的应变速度，此时得到的阻尼矩阵则正比于单元刚度矩阵，还有一些其到的阻尼矩阵则正比于单元刚度矩阵，还有一些其他类型的假设，如上述两者的组合，分别有：他类型的假设，如上述两者的组合，分别有： eeeeeeekmckcmc二、单元阻尼矩阵的计算二、单元阻尼矩阵的计算 对于组合阻尼，如已知结构

9、的阻尼比及结构的固对于组合阻尼，如已知结构的阻尼比及结构的固有频率，其计算方法有：有频率，其计算方法有：2222)(2)(2ijiijjjiijijji如果如果jijiji22ji则则1.1.结构无阻尼自由振动的运动方程结构无阻尼自由振动的运动方程三、机械结构固有频率与振型三、机械结构固有频率与振型 机械结构的振动固有频率和振型问题，在有限元方机械结构的振动固有频率和振型问题，在有限元方法求解释，对应的数学问题既是矩阵的特征值和特法求解释，对应的数学问题既是矩阵的特征值和特征向量问题。关于矩阵的特征值及特征向量问题，征向量问题。关于矩阵的特征值及特征向量问题，是矩阵理论中比较热门的研究领域，下

10、面我们仅简是矩阵理论中比较热门的研究领域，下面我们仅简单地罗列以下常见方法的名称，具体的方法求解步单地罗列以下常见方法的名称，具体的方法求解步骤，可以参考有关书籍，有大量的软件保重均包含骤，可以参考有关书籍，有大量的软件保重均包含求解特征值和特征向量的软件程序。求解特征值和特征向量的软件程序。 结构在无外力作用时，得到的是自由振动，此时结构在无外力作用时，得到的是自由振动，此时阻尼影响不大，结构的自由振动可简化为：阻尼影响不大，结构的自由振动可简化为： 0KM 1.1.结构无阻尼自由振动的运动方程结构无阻尼自由振动的运动方程三、机械结构固有频率与振型三、机械结构固有频率与振型设结构作简谐运动设

11、结构作简谐运动 tsin0代入无阻尼振动方程，可得代入无阻尼振动方程，可得 002MK上式解存在的条件为上式解存在的条件为 02MK这是计算方法中最典型的特征值问题。这是计算方法中最典型的特征值问题。2.2.矩阵迭代法矩阵迭代法 这种方法用于求解基频或最高阶频是很有效的，并这种方法用于求解基频或最高阶频是很有效的，并且能得到相应的特征向量。且能得到相应的特征向量。 020MK将无阻尼自由振动方程改写将无阻尼自由振动方程改写为为三、机械结构固有频率与振型三、机械结构固有频率与振型 0201KM即有即有 020S迭代步骤迭代步骤 )0(令令 )1 (2)1 ()0(S代入代入2.2.矩阵迭代法矩阵

12、迭代法 ）（和12)1 (求得求得 )2(2)2(1）（S再代入再代入 )1(2)1(iiiS）（以此类推以此类推 )1(kk）（收敛条件收敛条件2.2.矩阵迭代法矩阵迭代法三、机械结构固有频率与振型三、机械结构固有频率与振型 330352023300020001KM例题：已知一振动系统的质量矩阵、刚度矩阵用迭例题：已知一振动系统的质量矩阵、刚度矩阵用迭代法计算其最高阶固有频率和振型。代法计算其最高阶固有频率和振型。3/10002/100011M解：解：2.2.矩阵迭代法矩阵迭代法三、机械结构固有频率与振型三、机械结构固有频率与振型 1105 . 15 . 210231SMS 在开始迭代时，需

13、选取初始迭代向量，可以按经验在开始迭代时，需选取初始迭代向量，可以按经验估计，也可以用静力学特性的位移值，选得合适可估计，也可以用静力学特性的位移值，选得合适可以减少迭代时间。先假设：以减少迭代时间。先假设： T111)0(于是有于是有 0011111105 . 15 . 21023)0(S2.2.矩阵迭代法矩阵迭代法三、机械结构固有频率与振型三、机械结构固有频率与振型 T0011)1(2)1(推得推得 031130011105 . 15 . 21023)1(S继续迭代继续迭代 T03/113)2(2)2(推得推得 1 . 05 . 016 . 303/111105 . 15 . 21023)

14、2(S继续迭代继续迭代2.2.矩阵迭代法矩阵迭代法三、机械结构固有频率与振型三、机械结构固有频率与振型 如此继续迭代，经过如此继续迭代，经过1010次迭代，可得次迭代，可得 20467. 069300. 01386. 42047. 0693. 011105 . 15 . 21023)10(S 2047. 06930. 0120467. 069300. 01)10()11(推得推得 T205. 0693. 01386. 4)11(2)11(2于是于是2.2.矩阵迭代法矩阵迭代法三、机械结构固有频率与振型三、机械结构固有频率与振型 得到的固有频率是最高阶频率，因为振型的变化是：得到的固有频率是最高

15、阶频率，因为振型的变化是：205. 0693. 01符号变化两次，振系是符号变化两次，振系是3 3自由度，因此，得到的是第自由度，因此，得到的是第3 3阶频率和振型。阶频率和振型。 在工程实际中，人们一般关心的主要是结构的低阶在工程实际中，人们一般关心的主要是结构的低阶频率。因此，在进行迭代过程中作适当的变换，使矩频率。因此，在进行迭代过程中作适当的变换，使矩阵不按阵不按 为特征值进行迭代，而是按为特征值进行迭代，而是按 为特征为特征值进行迭代，从而得到值进行迭代，从而得到 的最大值，也是的最大值，也是 的的最小值。最小值。22/12/12 KMMK01 MK1两边同左乘两边同左乘 ，得到，得

16、到2.2.矩阵迭代法矩阵迭代法三、机械结构固有频率与振型三、机械结构固有频率与振型 在计算过程中，引入参数在计算过程中，引入参数21 将其代入无阻尼自由振动运动方程，则有将其代入无阻尼自由振动运动方程，则有 1K T MKT1令三、机械结构固有频率与振型三、机械结构固有频率与振型2.2.矩阵迭代法矩阵迭代法依次类推依次类推 采用前述的迭代步骤，用采用前述的迭代步骤，用 代替代替 ，即可得到，即可得到 值值 T S ) 1 () 1 ()0(T ) 1(1)(iiiT）（直到直到 ）（1)(-kk停止迭代停止迭代）（ 121i得到得到 ) 1( i此时为低阶特性此时为低阶特性三、机械结构固有频率

17、与振型三、机械结构固有频率与振型2.2.矩阵迭代法矩阵迭代法 330352023300020001KM例题：已知一振动系统的质量矩阵、刚度矩阵用迭例题：已知一振动系统的质量矩阵、刚度矩阵用迭代法计算其最高阶固有频率和振型。代法计算其最高阶固有频率和振型。 6/115 . 115 . 15 . 111111K解：解：三、机械结构固有频率与振型三、机械结构固有频率与振型2.2.矩阵迭代法矩阵迭代法 5 . 5315 . 4313211MKT于是于是 6 . 14 . 1161115 . 5315 . 431321)0(T仍选仍选 T111)0(三、机械结构固有频率与振型三、机械结构固有频率与振型2

18、.2.矩阵迭代法矩阵迭代法6287. 14433. 11773. 8629. 1443. 11768. 8628. 1442. 116 . 8继续迭代继续迭代从而得到从而得到 T629. 1443. 11114. 01)4(42）（三、机械结构固有频率与振型三、机械结构固有频率与振型3.3.用滤波法计算最低用滤波法计算最低n n阶特征对阶特征对 工程中关心的不仅是最低阶特征对，而是最低阶的工程中关心的不仅是最低阶特征对，而是最低阶的n n阶特征对，这是仅用迭代法不行，可用滤波法。阶特征对，这是仅用迭代法不行，可用滤波法。4.4.行列式搜索法行列式搜索法 这一方法利用特征值分离定理，通过对称矩阵

19、的三这一方法利用特征值分离定理，通过对称矩阵的三角分解计算矩阵的行列式值，用加速割线法求出靠角分解计算矩阵的行列式值，用加速割线法求出靠近下一个未知特征值的移动，然后用移位逆迭代求近下一个未知特征值的移动，然后用移位逆迭代求特征向量。特征向量。三、机械结构固有频率与振型三、机械结构固有频率与振型5.5.广义雅克比法广义雅克比法 广义雅克比法通过广义雅克比旋转矩阵把刚度广义雅克比法通过广义雅克比旋转矩阵把刚度矩阵和质量矩阵同时变换成对角矩阵，然后求得矩阵和质量矩阵同时变换成对角矩阵，然后求得特征值和特征向量，当矩阵阶数不高时，求解速特征值和特征向量，当矩阵阶数不高时，求解速度较快。度较快。6.6

20、.子空间迭代法法子空间迭代法法 子空间迭代法是瑞利子空间迭代法是瑞利- -李兹法和同时逆迭代法结李兹法和同时逆迭代法结合的产物，用于仅求解工程问题的低阶固有特征合的产物，用于仅求解工程问题的低阶固有特征对，求解速度非常快。对，求解速度非常快。三、机械结构固有频率与振型三、机械结构固有频率与振型7.7.兰索斯法兰索斯法 兰索斯法也是迭代法的一种，这是目前求解低兰索斯法也是迭代法的一种，这是目前求解低阶特征值和特征向量速度最快的一种，有兴趣的阶特征值和特征向量速度最快的一种，有兴趣的同学可以参阅同学可以参阅振动与冲击振动与冲击杂志杂志19871987年第年第3 3期期上吴立系老师的文章上吴立系老师

21、的文章“求解大型稀疏对称矩阵广求解大型稀疏对称矩阵广义特征值问题的义特征值问题的LanczosLanczos方法及通用程序方法及通用程序”。8.8.奇异刚度矩阵的处理奇异刚度矩阵的处理 采用移轴技术，在弹性位能中加入部分给定的采用移轴技术，在弹性位能中加入部分给定的动能，以便使刚度矩阵成为正定矩阵，关键在移动能，以便使刚度矩阵成为正定矩阵，关键在移轴系数的确定。轴系数的确定。四、机械结构动力响应的计算四、机械结构动力响应的计算 机械结构的动力响应计算是结构动力学的另一个主机械结构的动力响应计算是结构动力学的另一个主要问题，最常用的方法有振型叠加法和直接积分法。要问题，最常用的方法有振型叠加法和

22、直接积分法。1.1.振型叠加法振型叠加法 设结构的运动方程为设结构的运动方程为 )(tPKCM 并设已求得其无阻尼自由振动的频率和振型，记并设已求得其无阻尼自由振动的频率和振型，记 为第为第i i阶固有振型，则有振型的线性叠加来表示运动阶固有振型，则有振型的线性叠加来表示运动状态的结构位移为状态的结构位移为 i0 )()()(0220110tztztznn四、机械结构动力响应的计算四、机械结构动力响应的计算1.1.振型叠加法振型叠加法令令 n020100 )()()(21tztztzzn则则 z0 z 0z&0 )(000tPzKzCzM )(0000000tPzKzCzMTTTT IMT00

23、 2222100nTdiagK nnTdiagC222221100四、机械结构动力响应的计算四、机械结构动力响应的计算1.1.振型叠加法振型叠加法于是，振动方程解耦为于是，振动方程解耦为 )(2)(2)(202202222222101211111tpzzztpzzztpzzzTnnnnnnnTT&M&)()()(21tztztzn，依次用求解常微分方程的方法可以得到解：依次用求解常微分方程的方法可以得到解：再代回再代回 z0可以得到问题的解。可以得到问题的解。四、机械结构动力响应的计算四、机械结构动力响应的计算1.1.振型叠加法振型叠加法 通常情况下，高阶振型对动力响应影响较小，因通常情况下，

24、高阶振型对动力响应影响较小，因此，只取最低的此，只取最低的3 35 5阶（阶（1010）阶振型就可以得到）阶振型就可以得到满意精度的动力响应。满意精度的动力响应。2.2.直接积分法直接积分法对于动力学方程对于动力学方程 )(tPKCM 假设已知其初始条件假设已知其初始条件 000 如果将求解时间如果将求解时间T T等分成等分成n nT T个时间区间个时间区间 ，并能，并能通过前若干个时刻的解来确定下一时刻的解，就可通过前若干个时刻的解来确定下一时刻的解，就可获得问题的解，直接积分法可以解决这一问题。获得问题的解，直接积分法可以解决这一问题。t四、机械结构动力响应的计算四、机械结构动力响应的计算

25、2.2.直接积分法直接积分法 直接积分中对动力学方程是逐步地进行数值积分的，直接积分中对动力学方程是逐步地进行数值积分的，“直接直接”的意思是指：进行数值积分前没有进行把方的意思是指：进行数值积分前没有进行把方程变为另一种形式的变换。程变为另一种形式的变换。 直接积分有两个假设，一个是动力学方程的解只在直接积分有两个假设，一个是动力学方程的解只在相隔相隔 的一些离散时间区间上满足方程，而不要求的一些离散时间区间上满足方程，而不要求在任意时刻都满足方程；另一个是假设位移速度和加在任意时刻都满足方程；另一个是假设位移速度和加速度在每一个时间区间速度在每一个时间区间 内按一定的规律变化。内按一定的规

26、律变化。tt 直接积分法有中心差分法、直接积分法有中心差分法、HouboltHoubolt法、法、Wilson- Wilson- 法、法、NewmarkNewmark和龙格和龙格- -库塔法。库塔法。四、机械结构动力响应的计算四、机械结构动力响应的计算2.2.直接积分法直接积分法中心差分法中心差分法对于加速度对于加速度 ttttttt212 上式的误差为上式的误差为 高阶小量高阶小量 2t对于速度对于速度 tttttt21方程在方程在t t时刻为时刻为 ttttPKCM 将速度和加速度的差分格式代入方程，可以得到将速度和加速度的差分格式代入方程，可以得到四、机械结构动力响应的计算四、机械结构动

27、力响应的计算2.2.直接积分法直接积分法 ttttttCtMtMtKPCtMt)211()2()211(222中心差分法中心差分法显然，求解显然，求解 ，需要，需要2 2个初始条件个初始条件 和和 tt tt t我们有初始条件我们有初始条件 ，根据，根据 000 、 tttttt-00202121 四、机械结构动力响应的计算四、机械结构动力响应的计算2.2.直接积分法直接积分法中心差分法中心差分法可以解得可以解得 02002 ttt注意，有注意，有 02002 ttt于是，可以求得方程的解。归纳起来后，中心差分于是，可以求得方程的解。归纳起来后，中心差分法的计算步骤为法的计算步骤为2 2大步，

28、大步，9 9小步。小步。四、机械结构动力响应的计算四、机械结构动力响应的计算2.2.直接积分法直接积分法中心差分法中心差分法第一步：初始计算第一步：初始计算230212012211aaaatata1.1.形成形成 CKM、2.2.计算初始值计算初始值 000 、3.3.选取时间步长选取时间步长 , ,并计算积分常数并计算积分常数 crttt ,4.4.计算计算 0300- att5.5.形成有效质量阵形成有效质量阵: : CaMaM106.6.做三角分解做三角分解: : TLDLM 四、机械结构动力响应的计算四、机械结构动力响应的计算2.2.直接积分法直接积分法中心差分法中心差分法第二步：对每

29、一时间步长第二步：对每一时间步长 tttttCaMaMaKPP)()(1027.7.计算在时刻计算在时刻t t的有效载荷的有效载荷8.8.求解在时刻求解在时刻 的位移的位移 tt tttTPLDL9.9.如果需要，计算时刻如果需要，计算时刻 的速度和加速度的速度和加速度 tt )()2(10tttttttttttaa 四、机械结构动力响应的计算四、机械结构动力响应的计算2.2.直接积分法直接积分法中心差分法中心差分法说明：说明：1. 1. 的取值的取值t 是结构最高频率对应的周期，是最小周期。是结构最高频率对应的周期，是最小周期。crtncrTtt2.2.求解方程是显式差分，对于对角线质量阵和

30、无求解方程是显式差分，对于对角线质量阵和无阻尼振动求解很方便，可以直接求解，不需进行阻尼振动求解很方便，可以直接求解，不需进行三角分解。三角分解。四、机械结构动力响应的计算四、机械结构动力响应的计算2.2.直接积分法直接积分法中心差分法中心差分法例题：一个二自由度振动系统，其动力方程为：例题：一个二自由度振动系统，其动力方程为：100422610022121&已知该系统的自由振动周期为已知该系统的自由振动周期为 ，试用中心差分，试用中心差分法求解步长为法求解步长为 和和 时，方程的解。时，方程的解。8 . 22T102Tt 210Tt 解：取解：取1212个步长的系统响应，假设个步长的系统响应

31、，假设 0,000100004226100221 计算计算 0 四、机械结构动力响应的计算四、机械结构动力响应的计算2.2.直接积分法直接积分法中心差分法中心差分法考虑考虑 ，有，有 0392. 01,5 .25279. 128. 021,8 .12)28. 0(12302120aaaaaa 1000 28. 0t因而有因而有 392. 001000392. 00028. 000t四、机械结构动力响应的计算四、机械结构动力响应的计算2.2.直接积分法直接积分法中心差分法中心差分法 8 .12005 .25000079. 110028 .12M ttttP8 .12005 .255 .212245100对每个时刻步长求解方程对每个时刻步长求解方程 tttP8 .12005 .25时间时间t2t3t4t5t6t7t8t9t10t11t12t100.030.170.491.021.702.402.913.072.772.041.0220.391.452.834.145.025.264.904.173.372.782.542.60四、机械结构动力响应的计算四、机械结构动力响应