1、,第9章 动力问题有限元法,张 洪 伟,2,动力学问题,第2节 质量矩阵和阻尼矩阵,第1节 引言,第3节 直接积分法,第4节 振型叠加法,第5节 解的稳定性,第6节 大型特征值问题的解法,第7节 减缩系统自由度的方法,第8节 小结,3,第1节 有限元动力学方程的建立,动力学问题中最经常遇到的是结构动力学问题，它有两类研究对象。一类是在运动状态下工作的机械或结构，例如，高速旋转的电机，往复运动的内燃机，以及高速运行的飞行器，如何保证它们运行的平稳性及结构的安全性是极为重要的研究课题。另一类是承受动力载荷作用的工程结构，例如建于地面的高层建筑和厂房，核电站的安全壳和热交换器，这些结构的破裂、倾覆和
2、坍塌等破坏事故的发生，将给人民的生命财产造成巨大损失。正确分析和设计这类结构，在理论和实际上都是具有重要意义的。 动力学研究的另一重要领域是波在介质中的传播问题。,4,三维弹性动力学的基本方程是：,平衡方程,几何方程,物理方程,边界条件,初始条件,（在V域内）,（在V域内）,（在V域内）,（在Su域内） （在S域内）,（1.1）,（1.2）,（1.3）,（1.4） （1.5）,（1.6）,在动载荷作用下，对于任一瞬时，设单元节点发生虚位移 ，则单元 内也产生相应的虚位移 和虚应变 。单元内产生的虚应变能为:,单元除受动载荷外，还有加速度和速度引起的惯性力 和阻尼力 ，其中为材料密度，v是线性阻
3、尼系数。外力所做的虚功为：,式中，Pv、Ps、Pc分别为作用于单元上的动态体力、动态面力和动态集中力；V为单元面积；A为单元面积。,动力学方程建立：,且形函数仅为坐标x、y、z的函数，与时间无关，因此有,根据虚位移原理，有,代入经整理，可得单元运动方程为,由于,式中,分别称为单元的刚度矩阵、质量矩阵和阻尼矩阵，它们就是决定单元动态性能的特性矩阵。,称为单元节点动载荷列阵，它是作用在单元上的体力、面力和集中力向单元节点移置的结果。在动态分析和静力分析中，单元的刚度矩阵是相同的，外部载荷的移置原理也一样。,8,动力学有限元分析基本步骤如下：,（1）连续区域的离散化 （2）构造插值函数 由于只对空间
4、域进行离散，所以单元内位移u,v,w的插值分别表示为:,（1.7）,其中,9,（3）形成系统的求解方程,（1.8）,其中,分别是系统的结点加速度向量和结点速度向量，,M,C,K和Q(t)分别是系统的质量、阻尼、刚度和结点载荷向量。,10,（4）求解运动方程,（1.9）,如果忽略阻尼的影响，则运动方程简化为,如果上式的右端项为零，则上式进一步简化为,（1.10）,这是系统的自有振动方程，又称为动力特性方程。,（5）计算结构的应变和应力,11,从以上步骤可以看出，和静力分析相比，在动力分析中，由于惯 性力和阻尼力出现在平衡方程中，因此引入了质量矩阵和阻尼矩 阵，最后得到求解方程不是代数方程组，而是
5、常微分方程组。其 它的计算步骤和静力分析是完全相同的。 关于二阶常微分方程组的解法有两类：直接积分法和振型叠加法。,直接积分法是直接对运动方程积分。而振型叠加法是首先求解一 无阻尼的自由振动方程，然后用解得的特征向量，即固有振型对 运动方程式进行变换。,动力分析的计算工作量很大，因此提高效率，节省计算工作量的 数值方案和方法是动力分析研究工作中的重要组成部分。目前两 种普遍应用的减缩自由度的方法是减缩法和动力子结构法。,12,第2节 质量矩阵和阻尼矩阵,一、协调质量矩阵和集中质量矩阵,单元质量矩阵,称为协调质量矩阵。,集中质量矩阵假定单元的质量集中在结点上，这样得到的质量矩 阵是对角线矩阵。以
6、下分实体单元和结构单元进行讨论。,1.实体单元 介绍两种常用方法 （1）第一种方法,其中，ne是单元的结点数。该式的力学意义是：Mle每一行的主元 素等于Me中该行所有元素之和，而非主元素为零。,（2.1）,13,第2节 质量矩阵和阻尼矩阵,（1）第二种方法,该式的力学意义是：Mle每一行的主元素等于Me中该行主元素乘 以缩放因子a，而非主元素为零。,（2.2）,14,第2节 质量矩阵和阻尼矩阵,例1 计算平面应力（应变）单元的协调质量Me矩阵和集中质量矩阵Mle 。单元采用3结点三角形单元。,（1）协调质量矩阵 位移插值函数是,（2.3）,其中I是22单位矩阵。,（2.4）,Me,Ce,Ke
7、和Qe分别是单元的质量、阻尼、刚度和载荷矩阵。,15,第2节 质量矩阵和阻尼矩阵,算得单元的协调质量矩阵,（2.5）,其中，WtA是单元的质量，t是单元的厚度。,16,第2节 质量矩阵和阻尼矩阵,（2）集中质量矩阵 按第一种方法计算，得到集中质量矩阵为,（2.6）,此式的力学意义是：在单元的每个结点上集中1/3的质量。,17,第2节 质量矩阵和阻尼矩阵,按第二种方法计算，得到集中质量矩阵与第一种方法结果一样。,注：对于8结点矩形单元，两种方法得到的集中质量矩阵不同。,在实际分析中，更多的是推荐用第二种方法来计算集中质量矩阵。,2.结构单元,2结点经典梁单元、协调质量矩阵和集中质量矩阵如下所示：
8、,（1）协调质量矩阵 位移插值函数是,（2.7）,其中,18,第2节 质量矩阵和阻尼矩阵,计算得单元的协调质量矩阵为,（2.8）,其中，l是单元长度，WlA是单元的质量，A是截面面积。,（2）集中质量矩阵,（2.9）,此式的力学意义是在每个结点上集中1/2的单元质量。,19,第2节 质量矩阵和阻尼矩阵,需要指出，虽然质量矩阵M在理论上是正定的，但通常需要在计 算中对,进行精确积分才能保证此性质。如果计,算中采用低阶的积分，则M可能是奇异的，这将使后续的动力分 析发生困难，因此在选择Me的积分阶次时应予注意。,20,第2节 质量矩阵和阻尼矩阵,二、振型阻尼矩阵,它是假定阻尼力正比于质点运动速度的
9、结果，通常均将介质阻尼 简化为这种情况。这时单元矩阵比例于单元质量矩阵。,在以后的讨论中，将知道系统的固有振型对于M和K是具有正交 性的，因此固有振型对于M和K的阻尼矩阵C也是具有正交性的。 所以这种阻尼矩阵称为比例阻尼或振型阻尼。,在实际分析中要精确地决定阻尼矩阵是相当困难的，通常允许将 实际结构的阻尼矩阵简化为M和K的线性组合。这种振型阻尼称 为Rayleigh阻尼。,求解方法,求解运动方程,直接积分法,模态叠加法,隐式积分,显式积分,完整矩阵法,缩减矩阵法,完整矩阵法,缩减矩阵法,逐步积分法按是否需要联立求解耦联方程组，可分为两大类：隐式方法：逐步积分计算公式是偶联的方程组，需联立求解，
10、计算工作量大，通常增加的工作量与自由度的平方成正比，例如Newmark法、Wilson 法。 显式方法：逐步积分计算公式是解偶的方程组，无需联立求解，计算工作量小，增加的工作量与自由度成线性关系，如中心差分方法。重点介绍两种常用的时域逐步积分法中心差分法和Newmark法。,23,第3节 直接积分法,一、中心差分法,在中心差分法中，加速度和速度可以用位移表示，即,（3.2）,（3.1）,中心差分法的递推公式,（3.3）,上式是求解各个离散时间点解的递推公式，这种数值积分方法又 称为逐步积分法。,24,第3节 直接积分法,需要指出，此算法有一个起步问题，为此利用(3.1)，(3.2)得到。,将利
11、用中心差分法逐步求解运动方程的算法步骤归结如下：,1.初始计算,形成刚度矩阵K、质量矩阵M和阻尼矩阵C。 给定 选择时间步长t， t tcr，并计算积分常数,计算 形成有效质量矩阵 三角分解,25,第3节 直接积分法,2.对于每一时间步长（t0， t ，2 t ）,计算时间t的有效载荷,求解时间t t的位移,如果需要，计算时间t的加速度和速度,26,第3节 直接积分法,关于中心差分法还需要着重指出一下几点：,中心差分法是显式算法。 中心差分法是条件稳定算法。 显式算法用于求解由梁、板、壳等结构单元组成的系统的动态响应时如果对角化后的质量矩阵M中已略去了与转动自由度相关的项，则M的实际阶数仅是对
12、于位移自由度的阶数。 中心差分法比较适合于由冲击、爆炸类型载荷引起的波传播问题的求解。,对于结构动力学问题，一般说，采用中心差分法就不太适合。,中心差分法的精度和数值稳定性以上给出的中心差分逐步积分公式，是收敛的；具有2阶精度，即误差O(t2) ；是有条件稳定，稳定条件tTn/；具有较高的计算效率。,28,第3节 直接积分法,二、Newmark方法,在tt t的时间区域内， Newmark积分法采用下列的假设,（3.4）,（3.5）,其中和是按积分精度和稳定性要求决定的参数。另一方面， 和取不同数值则代表了不同的数值积分方案。,Newmark方法中的时间t t的位移解答a t t是通过满足时间
13、t t的运动方程的。,29,第3节 直接积分法,计算a t t的两步递推公式,（3.6）,将利用Newmark法逐步求解运动方程的算法步骤归结如下：,1.初始计算,形成刚度矩阵K、质量矩阵M和阻尼矩阵C。 给定,30,第3节 直接积分法,选择时间步长t 及参数和，并计算积分常数。这里要求： 0.50， 0.25(0.5+)2,形成有效刚度矩阵 三角分解,31,第3节 直接积分法,2.对于每一时间步长（t0， t ，2 t ）,计算时间t t的有效载荷,求解时间t t的位移,如果需要，计算时间t的加速度和速度,32,第3节 直接积分法,关于Newmark法还需要着重指出一下几点：,Newmark
14、法是隐式算法。 关于Newmark法的稳定性。证明，当0.50， 0.25(0.5+)2时，算法是无条件稳定的。 Newmark法适合于时程较长的的系统瞬态响应分析。 Newmark法的其它表达形式。,和是按积分精度和稳定性要求决定的参数。,=1/2和 =1/4，平均常加速度法。 =1/2和 =1/6，线性加速度法。 =1/2和 =0，中心差分法。,34,第4节 振型叠加法,振型叠加法在积分运动方程以前，利用系统自由振动的固有振型 将方程转化为n个相互不耦合的方程，对这种方程可以解析或数 值地进行积分。当采用数值方法时，对于每个方程可以采取各自 不同的时间步长，即对于低阶振型可采用较大的时间步
15、长。,这两者结合起来相当于直接积分法时很大的优点，因此当实际分 析的时间历程较长，同时只需要少数较低阶振型的结果时，采用 振型叠加法将时十分有利的。,35,第4节 振型叠加法,一、求解系统的固有频率和固有振型,此计算步骤是求解不考虑阻尼影响的系统自由振动方程，即,它的解可以假设为以下形式,（4.1）,其中，是n阶向量，是向量的振动频率，t是时间变量，t0 是由初始条件确定的时间常数。,36,第4节 振型叠加法,解方程确定和。特征向量1， 2， n代表系统的n个 固有振型。它们的幅度可按以下要求规定,这样规定的固有振型又称为正则振型，今后所用的固有振型，只 指这种正则振型。固有振型对于矩阵M是正
16、交的。,在有限元分析中，特别是动力分析中，方程的阶数很高而求解 的特征解又相对较少的特征值问题，称为大型特征值问题。,（4.2）,n自由度的振动系统，具有n个固有频率和与之对应的n阶主振型。且这些主振型之间存在着关于质量矩阵和刚度矩阵的正交性。,对应于,相减,关于正则振型,表明，对应于不同固有频率的主振型之间，即关于质量矩阵相互正交，又关于刚度矩阵相互正交，这就是主振型的正交性。还可以证明，零固有频率对应的主振型也必定与系统的其它主振型关于质量矩阵和刚度矩阵正交。,Ki称为第i阶主刚度或第i阶模态刚度；Mi称为第i阶主质量或第i阶模态质量。,可见，由于主振型的正交性，不同阶的主振动之间不存在动
17、能的转换，或者说不存在惯性耦合。同样可以证明第i阶固有振动的广义弹性力在第j阶固有振动的微小位移上的元功之和也等于零，因此不同阶固有振动之间也不存在势能的转换，或者说不存在弹性耦合。对于每一个主振动来说，它的动能和势能之和是个常数。在运动过程中，每个主振动内部的动能和势能可以互相转化，但各阶主振动之间不会发生能量的传递。因此，从能量的观点看，各阶主振动是互相独立的，这就是主振动正交性的物理意义。,以各阶主振型矢量为列，按顺序排列成一个nn阶方阵，称此方阵为主振型矩阵或模态矩阵，即,根据主振型的正交性，可以导出主振型矩阵的两个性质,主质量矩阵,主刚度矩阵,使Mr由对角阵变换为单位阵,将主振型矩阵
18、的各列除以其对应主质量的平方根，即,这样得到的振型称为正则振型。,正则振型的正交关系是,以各阶正则振型为列，依次排列成一个nn阶方阵，称此方阵为正则振型矩阵，即,由正交性可导出正则矩阵两个性质,在一般情况下，具有有限个自由度振动系统的质量矩阵和刚度矩阵都不是对角阵。因此，系统的运动微分方程中既有动力耦合又有静力耦合。对于n自由度无阻尼振动系统，有可能选择这样一组特殊坐标，使方程中不出现耦合项亦即质量矩阵和刚度矩阵都是对角阵，这样每个方程可以视为单自由度问题，称这组坐标为主坐标或模态坐标。由前面的讨论可知，主振型矩阵U与正则振型矩阵 ，均可使系统的质量矩阵和刚度矩阵转换成为对角阵。因此，可利用主
19、振型矩阵或正则振型矩阵进行坐标变换，以寻求主坐标或正则坐标。,44,第4节 振型叠加法,二、系统动力响应分析,1.位移基向量的变换,引入变换,（4.3）,此变化的意义是a(t)看成是i(i=1,2,n)的线性组合，i可以看 成是广义的位移基向量，xi是广义的位移值。从数学上看，是将 位移向量a(t)从以有限元系统的结点位移为基向量的n维空间转换 到以i为基向量的n维空间。,通常在实际分析中，需要求解的自由度方程数远小于系统的自 由度数n,45,第4节 振型叠加法,2.求解单自由度系统振动方程,单自由度系统振动方程的求解，通常采用杜哈美积分，又称为叠 加积分。这个方法的基本思想是将任意激振力ri
20、(t)分解为一系列 微冲量的连续作用，分别求出系统对每个微冲量的响应，然后根 据线性系统的叠加原理，将它们叠加起来。得到系统对任意激振 的响应。,杜哈美积分的结果是,其中,ai，bi是由起始条件决定的常数。,（4.4）,46,第4节 振型叠加法,3.振型叠加得到系统的响应,在得到每个振型的响应后，将它们叠加起来就是系统响应。,对振型叠加法的一些性质和特点：,振型叠加法中，将系统的位移转换到以固有振型为基向量的空 间这对系统的性质并无影响，而是以求解广义特征值为代价，得到n个单自由度系统的运动方程。 振型叠加法中对于n个单自由度系统运动方程的积分，比联立方程组的直接积分节省计算时间。 对于非线性
21、系统通常必须采用直接积分法。,47,第4节 振型叠加法,例3 以三自由度系统为例，现在用振型叠加法求解。 此时应求解的广义特征值问题是,（1）,按照一般的线性代数方法可以得到（1）式的解答为,（2）,48,第4节 振型叠加法,利用（2）式，可以转换为以1，2和3为基向量的 3个互不耦合的运动方程，即：,（3）,原系统的初始条件是,经转换后为,（4）,49,第4节 振型叠加法,利用无阻尼情形的杜哈美积分可以得到（3）式的精确解为：,（5）,最后利用振型叠加得到系统的位移为,（6）,2018/10/5,50,第4节 振型叠加法,根据（6）式计算得到每一时间步长的位移值如下：,对于tT3/100.3
22、63时，算得位移值：,t5T318.14时，算得位移值：,51,此结果是系统响应的精确解，可以用来检验中心差分法和 Newmark方法的结果。对于t0.363的情况，三者的比较见右图,由图可见，由于t较小，两种 直接积分法的结果都相当好。 而对于t18.14的情况，由于 t已相当大，虽然此时 Newmark方法的解仍然保持稳 定，但误差较大。,52,第5节 解的稳定性,解的稳定性定义是：如果在任何时间步长t条件下，对于任何 初始条件的解不是无限制的增长，则称此积分方法是稳定的； 如果t必须小于某个临界值，上述性质才能保持，则称此积分 方法是有条件稳定的。,讨论解的稳定性实质上是讨论误差引起的响
23、应。要讨论的方程是：,（5.1）,53,第5节 解的稳定性,一、中心差分法,利用中心差分法对（5.1）进行积分，可以写出,（5.2）,假定解的形式为,（5.3）,代入解出方程的根,（5.4）,54,第5节 解的稳定性,的根关系到解的性质，分析解稳定性的条件。,(1)真正解在小阻尼情况下应具有振荡特性，因此必须是复数 这就要求,即,（5.5）,因为,所以上式可以得到,（5.6）,55,第5节 解的稳定性,(2)真正解不应无限的增长，这就要求,（5.7）,表示无阻尼的自由振动。,为了保持解的稳定性，中心差分法的时间步长必须服从一下条件,（5.8）,其中tcr是临界时间步长，Tn是系统的最小固有周期
24、。,56,第5节 解的稳定性,二、Newmark方法,方程的根为,（5.10）,（5.9）,57,第5节 解的稳定性,分析解稳定性的条件,(1)真正解在小阻尼情况下应具有振荡特性，因此必须是复数 这就要求,即,（5.11）,当pi很大时，即t不受限制时，仍要求上式成立，必须是,（5.12）,58,第5节 解的稳定性,(2)真正解不应无限的增长，这就要求,即,（5.13）,当pi很大时，仍要求上式成立，必须是,（5.14）,（5.15）,59,第5节 解的稳定性,综合以上分析可以得到Newmark方法无条件稳定的条件是,（5.16）,如果不满足上述条件，要得到稳定解，时间步长t必须满足,（5.1
25、7）,60,第6节 大型特征值问题的解法,反迭代法 算法简单比较适合于只要求得到系统的很少数目特征值的情况,子空间迭代法 求解大型矩阵特征值问题的最常用且有效的方法，它适合于求解 部分特征值，被广泛应用于结果动力学的有限元分析中。,里兹向量直接叠加法 直接产生一组里兹向量，对运动方程进行缩减，求解缩减了的运 动方程，进而得到原系统方程的特征解。,Lanczos方法 直接产生一组Lanczos向量，对运动方程进行缩减，求解缩减了 的运动方程，进而得到原系统方程的特征解。,61,第7节 减缩系统自由度的方法,Guyan缩减法 又称为主从自由度法，通常不宜分析高阶的频率和振型。,动力子结构法 又称为
26、模态综合法，它能够大幅度地缩减动力分析的规模。大 型复杂系统分析如果采用动力子结构方法，计算效率将成量级 的提高。现今大型动力系统分析中广泛采用的就是该方法。,旋转周期分析方法 在理论上，分析中未引进自由度缩减方法所带来的近似性，因此 可以得到和整体结构分析时相同的精度。但它有局限性，只能用 于具有旋转周期的结构，不如子结构法应用范围广泛。,62,第8节 小 结,在结构动力学有限元求解方程的解法中，关于二阶常微分方程组的直接积分法，分别以中心差分法和Newmark 法为代表讨论了显式算法和隐式算法的各自算法步骤、特点、稳定性条件及其适合使用的情况。 振型叠加法也是动力分析中一种成熟而被广泛应用的方法。它的核心内容是动力特性方程的求解。将振型叠加法推广于有限元分析，本章讨论的反迭代法和子空间迭代法是现行最常用的方法。 关于系统自由度的减缩方法，现今大型动力系统分析中广泛采用的是动力子结构法，它能够大幅度地缩减动力分析的规模。,
