有限单元法基础第2版 课件 第11章 动力学问题

第11章动力学问题结构动力学问题运动构件运行的平稳性及其结构的安全性。冲击动力学问题承受动力载荷的工程结构。爆炸和高速冲击问题，波传播问题。特征值问题结构固有频率和模态。10.1引言悬臂梁动力响应简支梁模态分析弹性体表面波传播1弹性动力学基本方程1）控制方程2）边界条件3）初始条件11.2有限元动力方程2利用Galerkin法推导有限元方程1）插值函数构造2）单元平衡方程平衡方程和等效积分形式的Galerkin提法的任意性单元积分形式相邻单元作用力相邻单元作用于当前单元的等效结点载荷质量矩阵阻尼矩阵刚度矩阵等效结点载荷（时间函数）（3）系统有限元方程忽略阻尼的影响：无外载荷：自由振动方程，动力特性方程3利用哈密尔顿原理推导有限元方程哈密尔顿原理第一项进行分部积分在t1和t2时刻的结点位移值ae给定，则在该两个时刻的虚位移应为零的任意性可得由与利用Galerkin方法得到的完全相同协调质量矩阵平面问题11.3质量矩阵和阻尼矩阵1一致质量矩阵平面问题四边形等参单元三角形3结点单元三角形单元对称阵单元重量轴对称问题四边形单元三角形单元三角形3结点单元对称阵三维问题六面体等参单元四面体单元四面体4结点单元一维轴力单元2结点Lagrange单元平面桁架单元2结点平面桁架单元平面框架单元2结点Euler梁单元梁单元局部坐标系下的单元质量矩阵2集中质量矩阵将质量集中施加在单元的结点上。对角矩阵，可以提高计算效率。将一致质量转换为集中质量矩阵的方法a）实体单元方法一：行加法各行的主元素（矩阵的对角元素）等于各行所有元素之和，非主元素均为零。所有方法的核心是保证质量不变方法二：对角元素放大法各行主元素等于该行主元素乘以缩放因子a。a根据质量守恒原则确定。方法三：积分法在积分一致质量矩阵时，将结点取为积分点，采用标准形状函数时积分得到的质量阵就是对角阵。Ⅰ：行加法；Ⅱ：放大系数法；Ⅲ：结点积分法几种二维单元质量矩阵对角化三角形3结点单元b）结构单元通常忽略对应于转动自由度的元素。对于与位移自由度相关的元素则采用和实体单元相同的处理方法。实际计算结果表明，采用两种质量矩阵的计算结果相差不大。但采用集中质量矩阵可以简化计算，特别是对于显示积分方案，在对阻尼矩阵对角化处理后可以大大提高时间积分的效率。不考虑转动自由度质量时Euler梁单元的一致质量矩阵协调阻尼矩阵通常假设阻尼正比于质点的运动速度，或正比于应变率，由此可得：比例阻尼，或振型阻尼实际中经常采用Rayleigh阻尼3阻尼矩阵11.3系统特征值问题及其求解1无阻尼自由振动实特征值标准特征值问题非零解的条件特征方程展开得n次多项式n个特征对A为正定矩阵，特征值为正实数A为半正定矩阵零特征值个数等于刚体模态个数自由振动响应无阻尼自由振动方程代入振动方程广义线性特征值问题非零解的条件广义特征值问题广义特征值问题的特征方程对应于固有振型（模态）特征向量若规定正则振型若K和M为n×n的矩阵，可求得n个特征值固有频率位移约束足够为静定结构，无刚体位移，K正定，特征值为正实数没有足够的位移约束，K半正定，会出现固有频率为零的情况特征解的性质将特征解代入特征方程M、K对称①②由①和②可得固有振型对于矩阵M是正交的由①可得定义固有振型矩阵固有频率矩阵无论是求解系统动力特性，还是动力响应，实际需求解的特征解的个数通常远小于系统的自由度。2质量和刚度矩阵奇异的特征值问题结构上没有足够的约束条件，K奇异，结构自由振动响应不具备唯一性，结构有刚体模态，对应的固有频率为零。任意常数非奇异特征值有一个平移量移轴法梁和板壳单元忽略转动自由度，集中质量矩阵对角元mii=0特征对为如果M有r个零元素，系统具有r个无穷大特征值，相应的特征向量为r个单位向量，各单位向量中的1出现在与矩阵M中零元素相对应的位置。3特征值和特征向量的算法广义特征问题参考专门介绍特征值求解方法的著作向量迭代法变换法多项式迭代法Sturm序列法应用最广的方法矩阵反迭代法、子空间迭代法、里兹向量直接迭代法、Lanczos向量直接迭代法。子空间迭代法和Lanczos向量迭代法是公认的高效算法，已用于许多著名的商用有限元软件系统中。4算例：简支梁特征值计算简支梁的长度为l=40m，矩形截面考虑两种情况：宽度a=1m，高度b=2m和a=1m，高度b=1m。材料的杨氏模量E=3×104Pa，泊松比为ν=0，密度ρ=0.1kg/m3。简支Eular梁的固有频率理论解为简支梁固有频率有限元计算结果及理论解求解方法f1(Hz)f2(Hz)f3(Hz)a=1mb=2ma=1mb=1ma=1mb=2ma=1mb=1ma=1mb=2ma=1mb=1m8结点单元0.30930.15511.22410.61872.70731.38592结点梁单元0.30940.15511.22510.61892.71221.3864Eular梁理论0.31050.15531.24180.62092.79411.3971三阶简支梁模型及用8结点单元计算得到的前三阶模态二阶一阶简支梁模型11.4动力学方程求解方法当实际分析的时间历程较长，同时只需要少数较低阶振型的结果时，采用振型叠加法十分有利。振型叠加法计算分两步：1）求解系统的固有频率和振型（模态）2）求解系统的动力响应1模态叠加法求解系统的动力响应1）位移基向量的变换作变换：代入方程如果C为振型阻尼，则由Φ的正交性可得M和K的正交性则变换后的方程成为n个相互不耦合的二阶常微分方程：与单自由度系统的振动方程一致。2）求解单自由度系统振动方程（a）Duhamel（杜哈美）积分忽略阻尼（b）数值积分方法求解3）振型叠加得到系统的响应（a）通常高阶频率成份对系统的实际响应影响很小，只需求较低阶频率。（d）有限元求得的高阶频率精度较差。（e）非线性系统必须采用直接积分法，因为K＝K(t)，系统的特征解也随时间变化，无法利用振型叠加法。（b）模态叠加法的效率取决于在响应分析中所包含的模态数p，而p的取

值与结构和载荷的频率成分和空间分布有关。（c）对于爆炸或冲击载荷则需要考虑更多的振型，对于n自由度系统，可

能需要考虑高达2n/3阶振型。4）模态截断误差及修正利用前p阶模态叠加得到的系统响应为与由所有n阶模态叠加得到的响应之间的误差采用p阶模态模态分析中丢弃的部分载荷由振动理论可知，高频载荷在结构中产生的动态响应很小，不会超过静态响应，因此载荷产生的静态响应可以用来修正模态叠加法的结果5）Rayleigh阻尼系数确定方法由固有模态对质量矩阵和刚度矩阵的正交性若已知两个特定频率下的阻尼比通常假设应用于两个控制频率的阻尼比相同取多自由度系统的基频在对动力响应有显著贡献的高阶模态中选取2时间积分法的稳定性时间直接积分法将连续的时间离散成0,t1,t2,…tN时刻，利用0,t1,t2,…tn时刻的解，求tn+1时刻的解，又称为逐步积分法。递推形式：包含位移和速度等的向量rn：载荷A：递推矩阵或放大矩阵H：载荷系数初始条件X0对矩阵A进行特征分解（谱分解）P为由A的特征向量组成的矩阵，D为矩阵A的约当标准形，其对角元为矩阵A的特征值。当矩阵A有重特征值时，D中的非对角元可能为1。如果D对称固有频率矩阵振型矩阵若A为二阶非对称矩阵，且λ为其重特征值，则其约当标准形为其中α为0或者1若α=1Dn有界必须若α=0只要Dn有界若A没有重根，其约当标准形为只要Dn有界A的特征值模的最大值，谱半径。直接积分法的稳定性准则如果递推矩阵A没有重特征值，要求如果递推矩阵A存在重特征值，它们的模必须小于1讨论：积分方法存在数值阻尼人工阻尼越小时An+1趋于0的速度越快，积分方法的数值阻尼越大。谱半径除了影响算法的稳定性外，还反映了算法的数值阻尼。一般的结构动力学问题，响应主要受低阶振型控制，高阶振型贡献很小。受离散化的影响，有限元法对系统的高阶振型的近似精度很差。一个好的积分方法，应该在高频段具有一定的可控的数值阻尼，以有效地滤掉虚假的高频振型对系统振动响应的影响，同时在低频段的数值阻尼应尽可能地小，以保证结果的精度。谱半径在低频段应尽可能地接近于1，而在高频段应逐步光滑地减小，并趋于某个给定值。3中心差分法1）算法中心差分法t时刻速度近似时刻位移利用Tayler展开式时刻位移利用Tayler展开式忽略高阶小量时刻的运动方程：代入加速度和速度表达式由2）起步法：3）算法步骤初始计算形成刚度矩阵K、质量矩阵M和阻尼矩阵C给定确定选择时间步长，计算积分常数令n=0计算形成有效质量矩阵三角分解时间步循环计算时间t时刻的有效载荷求时刻的位移计算t时刻的加速度和速度4)解的稳定性利用中心差分法单自由度自由振动方程将上式和恒等式联立，写成递推公式（a）（b）无阻尼系统递推矩阵A的特征方程为稳定性条件：特征值的模要小于等于1（c=1）不稳定稳定（c）为复数稳定算法稳定条件：无阻尼单自由度系统中心差分法稳定性条件无阻尼多自由度系统中心差分法稳定性条件Tm是系统的最小固有振动周期是系统的最高阶固有频率阻尼系统中心差分法临界时间步长可以证明，系统最小固有振动周期总是大于或等于最小尺寸单元的最小固有振动周期。网格中最小尺寸的单元将决定中心差分时间步长的选择。单元尺寸越小，要求的时间步长也越小。确定最小固有振动周期的方法：划分网格后，找出尺寸最小的单元，形成单元的特征方程求解其最大特征值ωm，从而得到划分网格后，找出最小单元的最小边长L，可以近似估计其中C是声波传播速度4）讨论（a）采用集中质量矩阵且阻尼矩阵为对角阵时，方程求解无需求系数矩阵的逆，中心差分法是显式算法。忽略阻尼：求解方程效率高，在非线性分析中更有意义。（b）是条件稳定算法。稳定性条件：系统最小固有振动周期系统最高阶固有振动频率（d）适用于由冲击、爆炸类型载荷引起的波传播问题的求解。（g）结构动力响应中通常低频成分是主要的，容许采用较大的步长，通常采用无条件稳定的隐式算法。（c）如有对角元素等于或接近于零，如梁和板壳结构单元，不能采用中心差分法。

（e）涉及多体接触和大变形等复杂非线性问题，采用显式算法往往具有更好的收敛性。（f）可用于模拟拟静态问题，采用质量放大（Massscaling）增大临界步长。

4Newmark法假设平均加速度法中心差分法线性加速度法1）算法代入：2）算法步骤初始计算形成刚度矩阵K、质量矩阵M和阻尼矩阵C给定确定选择时间步长及参数，计算积分常数形成有效刚度矩阵三角分解时间步循环计算时间t时刻的有效载荷求时刻的位移计算时刻的加速度和速度3）解的稳定性将Newmark法得循环计算公式用于方程上两式写成递推公式A的特征方程（a）（b）不稳定综上，Newmark法无条件稳定条件若不满足上式条件，要得到稳定解（c）数值阻尼Newmark法的谱半径无数值阻尼有数值阻尼。为了使数值阻尼在高频段最大，以有效地滤除系统的虚假高频响应，应使β取最小值，从而使c取最小值。由稳定性条件最小极限谱半径Newmak法具有数值阻尼，且数值阻尼随的增大光滑地增大。4）讨论（a）方程的求解需要求系数矩阵的逆，是隐式算法。（b）时是无条件稳定算法。时间步长的选择取决于精度要求。（c）一般的结构动力学问题中，低频分量是结构响应的主要成分，因此只需要精确求解系统前p阶固有振型的贡献，这里p的取值和载荷的频率及其空间分布有关。一般地，可取时间步长大约Tp/20，Tp为系统第p阶固有周期。（d）适用于较长时间的瞬态分析。较大的时间步长可以滤掉高阶不精确特征解对系统响应的影响。5广义α方法利用如下平衡方程1）算法αf和αm为参数谱半径用于调节数值耗散的参数广义α法是无条件稳定直接积分法的谱半径在低频域内尽可能接近1，而随着谱半径的值应光滑减小，以达到数值阻尼逐渐增大的目的。的增大，广义α法对高频衰减最大为了尽可能地减小算法在低频段的数值阻尼，应满足如下条件Newmark法Hilber-Hughes-Taylor-α法HHT-α法或HHT法Wood-Bossak-Zienkiewic-α法WBZ-α法或Bossak法弹性体长95m，高40m。其密度为300kg/m3，杨氏模量为2.0×108Pa，泊松比为0.3。底边竖直方向约束，左右水平位移约束。上表面中点受向下三角形集中冲击力载荷作用，计算模拟表面波的传播过程。采用四边形4结点单元，减缩积分方案，网格尺寸设置为1m。利用ABAQUS/Explicit显示积分模块，采用时间自动增量步，模拟得到的该弹性体的变形和Mises应力响应如图中所示。有限元模型冲击载荷6算例：波传播及梁弯曲动力响应问题t=0.10st=0.061799st=