第9章 动力学有限元

1问题描述如图所示有一工字形截面的外伸梁，外伸端长度为a=1m，跨度l=2m，外伸端受到W=10KN/m的均布载荷的作用。工字形截面的截面面积为A=45cm2，弹性模量E=200GPa，抗弯惯性矩Iz=5000cm4，求此外伸梁跨中的最大挠度。2

问题描述有一材料为钢的轴类零件，其结构如图所示，两端受50MPa的面载荷作用。已知钢的弹性模量是200GPa，泊松比为0.3，试分析该零件内部的应力分布情况。3

问题描述现有一个薄壁圆筒，如图所示。圆筒长度L为0.5m，壁厚t为5mm，内径R为0.2m，薄壁圆筒在其长度的中心处受一对沿着直径方向的压力F的作用，力的大小为1000N，求薄壁圆筒在受力点处的径向位移，圆柱的两端在边界处自由。已知薄壁圆筒的弹性模量为200GPa，泊松比为0.3。4梁单元板单元的应用长宽均为1m的厚度为0.05m的钢板，在两边和中间位置均焊接有加强筋，建立其有限元分析模型。第9章动态分析有限元法9.1引言9.2动力学有限元基本方程9.3质量矩阵和阻尼矩阵9.4结构的固有频率和固有振型9.5结构动力响应9.6动力响应算例9.1引言动力学问题中最经常遇到的是结构动力学问题，它有两类研究对象。一类是在运动状态下工作的机械或结构，例如，高速旋转的电机，往复运动的内燃机，以及高速运行的飞行器，如何保证它们运行的平稳性及结构的安全性是极为重要的研究课题。另一类是承受动力载荷作用的工程结，例如建于地面的高层建筑和厂房，正确分析和设计这类结构，在理论和实际上都是具有重要意义的。动力学研究的另一重要领域是波在介质中的传播问题。有限元方程（刚度方程）：静力学问题：[K]{δ}={F}静力问题：1)静止;2)匀速动力问题：外载随时间变化大动态分析的必要性：当产品受到随时间变化的动载荷时，需要进行动态分析，以了解产品动态特性。动载荷（又称动力分析）固有特性分析响应分析固有频率振型位移响应速度响应加速度响应动应变动应力固有特性：是一组模态参数构成，它由结构本身（质量与刚度分布）决定，而与外部载荷无关，但决定了结构对动载荷的响应；响应分析：是计算结构对给定动载荷的各种响应特性。以三维实体动力分析为例，用有限元法求解的基本步骤如下：（1）连续区域的离散化（2）构造插值函数由于只对空间域进行离散，所以单元内位移u,v,w的插值分别表示为:（9.1）其中（3）形成系统的求解方程（9.2）其中分别是系统的结点加速度向量和结点速度向量，M,C,K和Q(t)分别是系统的质量、阻尼、刚度和结点载荷向量。（4）求解运动方程（9.3）如果忽略阻尼的影响，则运动方程简化为如果上式的右端项为零，则上式进一步简化为（9.4）这是系统的自有振动方程，又称为动力特性方程。（5）计算结构的应变和应力结构动力学问题的有限元法的实质就是将一个弹性连续体的振动问题，离散为一个以有限个节点位移为广义坐标的多自由度系统的振动问题。其基本原理和分析方法类同静力学的有限元法，按杆梁、薄板等不同结构进行分析。不同的是，应用振动理论建立动力学方程时，在单元分析中除需形成刚度矩阵外，还需形成质量矩阵，阻尼矩阵；在整体分析中，不仅求动力响应，还有求解特征值问题（结构振动的固有频率及相应的振动型（或模态））从以上步骤可以看出，和静力分析相比，在动力分析中，由于惯性力和阻尼力出现在平衡方程中，因此引入了质量矩阵和阻尼矩阵，最后得到求解方程不是代数方程组，而是常微分方程组。其它的计算步骤和静力分析是完全相同的。关于二阶常微分方程组的解法有两类：直接积分法和振型叠加法。直接积分法是直接对运动方程积分。而振型叠加法是首先求解一无阻尼的自由振动方程，然后用解得的特征向量，即固有振型对运动方程式进行变换。动力分析的计算工作量很大，因此提高效率，节省计算工作量的数值方案和方法是动力分析研究工作中的重要组成部分。目前两种普遍应用的减缩自由度的方法是Guyan减缩法和动力子结构法。从静力学有限元法可知，有限元的基本思想是将弹性体离散成有限个单元，建立整体刚度平衡方程：关于静力问题和动力问题的区别，据达朗贝尔原理，动力学问题只要在外力中计入惯性力后，便可按静力平衡处理。考虑到动力问题中的载荷和位移均为时间的函数，上式可记为：由于动力载荷可为作用于弹性体上的动载荷，也可为弹性体的惯性力，也可为与速度相关的阻尼力，即：据惯性力定义表示为：如阻尼力正比与速度，则动力学基本方程：9.2振动基本方程的建立1、单元刚度阵任取一个单元，单元节点位移为，节点速度和加速度为：，则单元节点内任一点的位移[N]为形函数，与时间t无关，为X、Y、Z的函数，它与静力分析中一样；由于[N]与时间无关，则单元应变矩阵，应力矩阵仍与静力分析完全相同：则刚度矩阵同样与静力情况相同：9.3单元质量、阻尼、刚阵计算2、惯性力与单元质量阵设单元节点加速度为，则单元内任一点的加速度：设单元的质量密度为，则单位体积中的惯性力为：负号表示惯性力与加速度相反。显然，整个单元上惯性力即为上式的积分。如何将这个作用于单元上的惯性力移置到单元节点上，通常有两种方法：1）虚功原理法——求得一致质量矩阵2）直接分配法——即按重心不变原则分配，求得集中质量矩。在动态分析中，单元的质量矩阵通常采用以下两种形式。1、一致质量矩阵按形成的单元质量矩阵称为一致质量矩阵，因为它采用了和刚度一致的形函数。这种质量矩阵取决于单元的类型和形函数的形式。2、集中质量矩阵集中质量矩阵将单元的分布质量按等效原则分配在各个节点上，等效原则就是要求不改变原单元的质量中心，这样形成的质量矩阵称为集中质量矩阵。集中质量矩阵是一个对角阵，集中质量矩阵：是一个对角阵，因而可简化动态计算，减小存储容量。利用这种矩阵计算出的结构固有频率偏低。不过有限元模型本身比实际结构偏刚，两者相互补偿，计算出的固有频率反而更接近真实值。一致质量矩阵：由于分布较合理，因此可以求得更精确的振型，另外，整个模型的质量分布还受网格划分形式的影响。这里[M]为单元的一致质量矩阵。显然，对于不同的单元，因形函数不同，则质量矩阵也是不同的。1）虚功原理法设单元中发生虚位移为则单元惯性力作的虚功为：单元节点上节点惯性力所作的功为：将和代入可得平面常应变三角形单元的一致质量阵为：单元质量矩阵一般而言，一致质量较准确地反映了单元内质量分布的实际情况，集中质量精度不如前者，但不存在耦合，使计算大大简化，是工程中常用的方法。2）直接分配法将单元内分布质量按重心不变原则分配至单元节点上，所产生的质量矩阵是没有耦合项的对角矩阵。如六自由度的平面三角形单元，单元总质量为W/g，则平均分配至三个节点上的质量所形成的质量阵为：3、单元阻尼阵单元阻尼力主要指结构阻尼力，它是由结构内部材料内摩擦引起的阻尼。设结构阻尼系数为，则单位体积产生的阻尼力（即阻尼力密度）为：利用虚功原理同理可得：一旦单元刚阵、质量矩阵、阻尼矩阵求得，则动力学方程中的整体刚阵、质量阵等可类似静力分析的刚度矩阵组装得到：计算结构的固有频率和振型是结构动力学分析的主要内容，也是分析结构动力响应和其它动力特性问题的基础。由于一般结构阻尼对结构的固有频率和振型影响极小，所以，求结构的固有频率和振型时，直接用无阻尼的自由振动方程求解。即因任意弹性体的自由振动都可分解为一系列的简谐振动的迭加：即结构上各节点位移为为节点位移振幅向量（即振型），与时间t无关的位移幅值；为与该振型对应的频率。9.3结构无阻尼自由振动

1、固有频率和振型计算将节点位移代入动力方程，化简得广义特征值问题：由于结构自由振动时，各个节点的振幅不可能全为零，则称为结构的特征方程，即求结构的固有频率和振型归结为特征值问题。设计结构的自由度为n，则特征方程为的n次代数方程，其n个根称为特征值，记为它们的平方根称为系统的固有频率，即将这些固有频率从小到大依次排列为最低的频率称为基频，它是所有频率中最重要的一个。这个过程称之为正规化利用正规化，可得2、特征向量对应每个固有频率，可有方程由此求得一组节点振幅不全为0的向量称为特征向量，也称为振型或模态向量。由于上述方程为齐次方程，显然解不唯一，也就是说：振型的形状是唯一的，但其振幅不是唯一的；或一个特征值可对应有多个特征向量，但一个特征向量只对应一个特征值。实际中，常选特征向量使则对应所有的特征值问题:3、特征向量的性质正交性：任意两个特征值对应的特征向量关于质量矩阵或刚度矩阵正交。即设则有若将所有的特征值对应的特征向量组装成特征向量矩阵，即考虑到正规化:可进一步记为：可简记为矩阵形式：1、幂迭代法特点：用于计算最大（主）特征值十分有效。这里[D]称为动力矩阵，也即一个变换矩阵，它可将任一特征向量变换为一常数与其自身的乘积.9.4特征值问题的解法结构固有频率和振型的计算归结为求的特征值和特征向量。由于有限元法将结构离散为n个自由度，n一般相当大，故n次特征方程的直接求解十分困难，常求其近似解，常用的求解方法有幂迭代法、逆迭代法、子空间迭代法等。由于任两个特征值对应的特征向量是正交的，则n个特征向量可组成特征向量空间中的一个特征向量基，其特征向量空间中的任一特征向量可表示为基向量的线性组合。即存在任一向量：设这个向量被[D]变换后形成一新的特征向量为：类推，可得：由于所有的特征值排列为：即存在考虑到问题为齐次方程，特征向量前的系数可以略去，则上式在p趋近无穷时，其第一项就趋近实际计算，只需迭代有限次即可得精确解。幂法迭代格式1、选初始特征向量，如单位向量2、构造新特征向量，并归一化3、计算特征值近似值4、计算相邻两次迭代的特征值误差，检查是否收敛若需计算二阶、三阶等特征值，则需构造新的动力矩阵2、逆迭代法逆迭代法也称为反幂法，类似于幂法，特征值问题改写为：其具体迭代格式为：1）选初始向量如单位向量2）计算中间向量3）求解线性方程组4）归一化5）计算特征值近似值6）计算相邻两次迭代的特征值误差，检查是否收敛9.5动力响应的计算对于受迫振动，基本方程为求解此方程通常有两种数值方法：振型迭加法和逐次积分法1、振型迭加法振型迭加法的基本思想是利用结构固有振型的正交性，把结构的复杂振动分解为一组相互独立的单自由度振动（即解耦），从而求得结构的位移响应。设结构无阻尼自由振动的各阶固有频率和相应的固有振型为：则结构任意时刻的受迫振动产生的位移可认为是n个固有振型为基的线性组合，即为组合系数，是时间t的函数，也称为振形坐标广义质量阵广义阻尼阵广义刚度阵广义激振力上式可记为这里代入动力学方程：左乘据正交性可知，这些广义矩阵均为对角矩阵，即表示方程各个变量之间是没有耦合项的，从而动力方程转化为n个相互独立的单自由度振动的动力方程，即：分别求解这n个方程可求得从而求得动力方程的位移解:进而可求得速度、加速度。采用瑞利阻尼，即[C]=α[M]+β[K]2、逐次积分法基本思想：将时间t离散为n个区间，并假设在一个时间区间内，结构的加速度响应为线性变化，由此，对加速度积分，可得速度和位移，一旦所有区间计算完毕，则求出结构的动力响应。假设在至t的很小时间间隔内，加速度线性变化：对积分，并引入初始条件待定积分常数将代入t时刻的动力方程并整理后即可逐步求解各时刻的加速度，然后求出各时刻的速度和位移。3.直接积分法一、中心差分法

在中心差分法中，加速度和速度可以用位移表示，即（3.2）（3.1）中心差分法的递推公式（3.3）上式是求解各个离散时间点解的递推公式，这种数值积分方法又称为逐步积分法。直接积分法需要指出，此算法有一个起步问题，为此利用(3.1)，(3.2)得到。将利用中心差分法逐步求解运动方程的算法步骤归结如下：1.初始计算形成刚度矩阵K、质量矩阵M和阻尼矩阵C。给定选择时间步长△t，△t<△tcr，并计算积分常数计算形成有效质量矩阵三角分解直接积分法2.对于每一时间步长（t＝0，△t

，2△t…）计算时间t的有效载荷求解时间t＋△t的位移如果需要，计算时间t的加速度和速度直接积分法关于中心差分法还需要着重指出一下几点：中心差分法是显式算法。中心差分法是条件稳定算法。显式算法用于求解由梁、板、壳等结构单元组成的系统的动态响应时如果对角化后的质量矩阵M中已略去了与转动自由度相关的项，则M的实际阶数仅是对于位移自由度的阶数。中心差分法比较适合于由冲击、爆炸类型载荷引起的波传播问题的求解。对于结构动力学问题，一般说，采用中心差分法就不太适合。直接积分法二、Newmark方法

在t～t＋△t的时间区域内，Newmark积分法采用下列的假设（3.4）（3.5）其中α和δ是按积分精度和稳定性要求决定的参数。另一方面，α和δ取不同数值则代表了不同的数值积分方案。Newmark方法中的时间t＋△t的位移解答at＋△t是通过满足时间t＋△t的运动方程的。直接积分法计算at＋△t的两步递推公式（3.6）将利用Newmark法逐步求解运动方程的算法步骤归结如下：1.初始计算形成刚度矩阵K、质量矩阵M和阻尼矩阵C。给定直接积分法选择时间步长△t及参数α和δ，并计算积分常数。这里要求：δ≥0.50，α≥0.25(0.5+δ)2形成有效刚度矩阵三角分解第3节直接积分法2.对于每一时间步长（t＝0，△t

，2△t…）计算时间t＋△t的有效载荷求解时间t＋△t的位移如果需要，计算时间t的加速度和速度直接积分法关于Newmark法还需要着重指出一下几点：Newmark法是隐式算法。关于Newmark法的稳定性。以后将证明，当δ≥0.50，α≥0.25(0.5+δ)2时，算法是无条件稳定的。Newmark法适合于时程较长的的系统瞬态响应分析。Newmark法的其它表达形式。Newmark法的另一种以为未知量的两步递推公式（3.7）直接积分法Newmark法的以为未知量的三步递推公式（3.7）其中Newmark法的两步递推公式和三步递推公式中，令α＝0，δ＝1/2，就可以得到中心差分法的两步和三步递推公式。这样一来，这两种时间积分公式就采用了统一的表达形式，便于程序编制，特别时便于应用在隐式－显式混合时间积分方案种。直接积分法例2

考虑一个三自由度系统。它的运动方程是（1）初始条件：当t＝0时，已知此系统的固有频率是：相应的振动周期是：T1＝1089，T2＝4.444，T3＝3.628。直接积分法（1）用中心差分法求解系统响应时间步长分别取△t＝T3/10＝0.363和△t＝5T3＝18.14进行计算。对于t＝0，可以计算得到然后按中心差分法所列步骤进行计算。△t＝T3/10＝0.363时c0＝7.589，c1＝1.377，c2＝15.178，c3＝6.588e－2直接积分法对于每一时间步长，先计算有效载荷（2）在从下列方程计算t＋△t时间的位移at＋△t（3）第3节直接积分法由上式得到的每一时间步长的位移结果如下：该结果将在后续内容中与精确解进行比较。直接积分法△t＝5T3＝18.14时，按相同的步骤计算，所得结果如下:在计算下去，位移将继续无限增大，这是不步稳定的典型表现。直接积分法（2）用Newmark法求解系统的响应时间步长分别取△t＝T3/10＝0.363和△t＝5T3＝18.14进行计算。对于t＝0，可以计算得到然后按Newmark法所列步骤进行计算。给定α＝0.25及δ＝0.5。△t＝T3/10＝0.363时c0＝30.356，c1＝5.510，c2＝11.019，c3＝1.0c4＝1.0，c5＝0.0，c6＝0.1815，c7＝0.1815第3节直接积分法对于每一时间步长计算有效载荷然后求解时间t＋△t的位移at＋△t并计算时间t＋△t的加速度和速度直接积分法得到的每一时间步长的位移结果如下：该结果将在后续内容中与精确解进行比较。△t＝5T3＝18.14时，按相同的步骤计算，所得结果如下:振型叠加法振型叠加法在积分运动方程以前，利用系统自由振动的固有振型将方程转化为n个相互不耦合的方程，对这种方程可以解析或数值地进行积分。当采用数值方法时，对于每个方程可以采取各自不同的时间步长，即对于低阶振型可采用较大的时间步长。这两者结合起来相当于直接积分法时很大的优点，因此当实际分析的时间历程较长，同时只需要少数较低阶振型的结果时，采用振型叠加法将时十分有利的。振型叠加法一、求解系统的固有频率和固有振型此计算步骤是求解不考虑阻尼影响的系统自由振动方程，即它的解可以假设为以下形式（4.1）其中，φ是n阶向量，ω是向量φ的振动频率，t是时间变量，t0是由初始条件确定的时间常数。振型叠加法解方程确定φ和ω。特征向量φ1，φ2，…φn代表系统的n个固有振型。它们的幅度可按以下要求规定这样规定的固有振型又称为正则振型，今后所用的固有振型，只指这种正则振型。固有振型对于矩阵M是正交的。在有限元分析中，特别是动力分析中，方程的阶数很高而求解的特征解又相对较少的特征值问题，称为大型特征值问题。（4.2）振型叠加法二、求解系统动力响应1.位移基向量的变换引入变换（4.3）此变化的意义是a(t)看成是φi(i=1,2,…n)的线性组合，φi可以看成是广义的位移基向量，xi是广义的位移值。从数学上看，是将位移向量a(t)从以有限元系统的结点位移为基向量的n维空间转换到以φi为基向量的n维空间。通常在实际分析中，需要求解的但自由度方程数远小于系统的自由度数n振型叠加法2.求解单自由度系统振动方程单自由度系统振动方程的求解，通常采用杜哈美积分，又称为叠加积分。这个方法的基本思想是将任意激振力ri(t)分解为一系列微冲量的连续作用，分别求出系统对每个微冲量的响应，然后根据线性系统的叠加原理，将它们叠加起来。得到系统对任意激振的响应。杜哈美积分的结果是其中ai，bi是由起始条件决定的常数。（4.4）振型叠加法3.振型叠加得到系统的响应在得到每个振型的响应后，将它们叠加起来就是系统响应。对振型叠加法的一些性质和特点：振型叠加法中，将系统的位移转换到以固有振型为基向量的空间这对系统的性质并无影响，而是以求解广义特征值为代价，得到n个单自由度系统的运动方程。振型叠加法中对于n个单自由度系统运动方程的积分，比联立方程组的直接积分节省计算时间。对于非线性系统通常必须采用直接积分法。振型叠加法例3

仍以例2中三自由度系统为例，现在用振型叠加法求解。此时应求解的广义特征值问题是（1）按照一般的线性代数方法可以得到（1）式的解答为（2）振型叠加法利用（2）式，可以将原文体转换为以φ1，φ2和φ3为基向量的3个互不耦合的运动方程，即：（3）原系统的初始条件是经转换后为（4）振型叠加法利