机械振动数值分析.ppt

参考书目 有限单元法，王勖成，清华大学出版社 计算动力学，王天舒，清华大学出版社 有限元分析的概念与应用，Robet D.Cook，西安交通大学出版社,求解域的离散 将连续体求解域离散为若干个子域，通过子域边界协调条件构造连续体； 利用近似函数求解 每个子域内用近似函数分片表示全域内待定的未知变量； 原问题的等效 利用变分或加权余量法，建立基本变量代数方程组或常微分方程组。,概 述,有限元法要点,复杂几何构型的适应性； 各种物理问题的可应用性； 建立于严格理论基础上的可靠性； 适合计算机实现的高效性。,有限元法特性,单元类型和形式； 有限元法的理论基础和离散格式； 有限元方程的求解方法； 有限元分析程序开发。,有限元法的发展和现状,标准化 任意复杂问题 模块化分解，单元建模 有限种类模块化单元 规范化 几何建模 力学建模 求解 后处理分析 通用化 形成标准模块化程序 应用规模化、普及性 求解问题规模庞大，易于为工程技术人员掌握,有限元法的优势,连续性假设 变形体内部处处连续 均匀性假设 变形体内部物质分配均匀 各向同性假设 物质在各方向上特性相同 线弹性假设 变形与外力作用的关系为线性 小变形假设 变形量远小于物体本身尺寸,基本假设,第一章 线弹性动力学变分原理,运动微分方程,应变方程,本构方程（物理方程）,边界条件,初始条件,变形体域内任意一点在任意时刻均满足运动微分方程。 变形体边界上任意一点在任意时刻均满足边界条件。,问题的精确解的特点,加权余量法的特点,变形体域内和边界上任意一点在任意时刻均近似满足运动微分方程。,残余力方程,加权余量法允许运动平衡方程和边界条件在若干局部存在残差，但要求残余力在域内和边界上的加权积分为零，即,当上式对任意权函数均满足，则称为式（7）为微分方程（1）和边界条件（4）的等效积分形式。,一般情况可选择近似解,将式（8）代入式（7），通过确定系数强迫残余力在域内和边界上在某种平均意义下为零。下面的讨论中假设近似函数完全满足边界条件，只考虑域内残差问题。,权函数可以选N个函数的线性组合，即,将式（9）代入式（7），得,配点法 取Dirac函数为权函数 子域法 权函数在N个子域内取1，在子域外取零，即,几类常用权函数,最小二乘法 调整近似函数中的参数，使余量均方和最小，即,几类常用权函数,迦辽金法 取试探函数为权函数 迦辽金法的特点,几类常用权函数,余量方程相当于虚功； 求解方程系数矩阵有对称性； 当存在泛函时与变分法有等效结果。,例1 用各种加权余量法计算图示弹性基础梁的挠度,弹性基础梁的基本微分方程和边界条件为,取试探解为无弹性基础时的精确解,则近似解可表示为,容易发现，式（16）严格满足边界条件,残差方程可写为,达朗伯拉格朗日原理,在式（7）中取权函数为真实位移的变分，可得与运动微分方程和边界条件的等效积分形式。,对方程中的第一项进行分部积分，得,式（19）为运动微分方程和边界条件等效积分的“弱”形式,阶连续性函数,n-1导数连续，且第n阶导数仅有有限个可积间断点的函数,哈密顿原理,对式（18）在任意时间间隔内积分,对给定时刻，方程中的第一项可转化为,将式（21）代入式（20），得普遍意义下的哈密顿原理,式（22）说明，对真实运动，系统动能变分和内、外力虚功在任意时间间隔内对时间的积分为零。,考虑粘滞力后,式（23）中考虑了粘滞力的虚功。,考虑到应力应变关系，外力虚功可改写为,将式（25）代入式（22），可得,式（26）说明，完整有势系统在任意时间间隔内满足几何关系和给定位移边界条件的所有可能运动中，真实运动使哈密顿作用量取驻值,例2 用哈密顿原理推导受均布动载荷的等截面悬臂梁的振动微分方程,梁内任意一点位移可表示为,任意一点的应变可表示为,系统总势能可表示为,系统总动能可表示为,对系统势能取变分,对系统动能取变分,代入哈密顿方程得,拉格朗日乘子法,约束泛函,罚函数法,引入附加泛函后，原泛函的有附加条件的驻值问题转化为修正泛函无附加约束条件的驻值问题。,拉格朗日乘子法修正泛函的变分,以离散结构为例，需要满足位移边界条件,修正泛函为,驻值条件,由此得,拉格朗日乘子法修正泛函的变分,拉格朗日乘子法的特点 方程组的阶数增加； 拉格朗日乘子有明确物理意义，相当于力； 导出的系数矩阵存在零对角元。,罚函数法修正泛函的变分,仍以离散结构为例，修正泛函为,驻值条件,由此得,罚函数法的特点 附加条件近似满足； 不增加方程阶数； 罚参数相当于刚度系数； 解的精度与罚参数有关。,罚函数法修正泛函的变分,例3 如图所示弹簧系统，在第一个自由度上施加外力，在第二个自由度上施加强制位移。用拉格朗日乘子法和罚函数法计算第一个自由度的位移,未施加强制位移时可写出系统平衡方程,将强制位移边界代入，可解得精确解,为满足强制约束条件施加在第二个自由度上的载荷,拉格朗日乘子法,可解得,为约束反力的负值,罚函数法,广义变分原理,将场函数应满足的附加条件引入泛函，将的变分原理转化为无附加约束条件的变分原理。,利用拉格朗日乘子法将哈密顿原理附加条件引入泛函，得H-W变分原理泛函为,其中,广义变分原理,由于所有变分项均独立，可推导出泛函取驻值条件,广义变分原理,拉格朗日乘子的物理意义为,由此得H-W变分原理泛函为,习题2 如图所示弹簧系统，在第一个和第三个自由度上施加外力，在第二个自由度上施加强制位移。用拉格朗日乘子法和罚函数法计算第一个和第三个自由度的位移，并计算第二个自由度上的作用力,几何关系 应力应变关系 ,杆单元的刚度矩阵和质量矩阵,第二章 有限元离散,泛函,杆单元的刚度矩阵和质量矩阵,取位移插值函数，用单元结点位移将杆内任意点位移表示为,将式（2）代入式（1），得,取变分，得,杆单元的刚度矩阵和质量矩阵,记,将式（8）代入式（7），结合哈密顿原理，并考虑到变分的任意性可解得,杆单元的刚度矩阵和质量矩阵,2结点杆单元的刚度矩阵和质量矩阵如式（10）所示,杆单元的刚度矩阵和质量矩阵,若采用集中质量矩阵，则为,考虑在平面中变形的杆单元，如图所示,平面杆单元的坐标转换,总体坐标系内的杆单元,则单元刚度矩阵、质量矩阵和载荷向量转换到总体坐标系中为,平面杆单元的坐标转换,对空间杆系，转换原理与平面杆系原理相同，转换矩阵可表示为,考虑如图所示杆系结构,平面杆单元例子,平面杆系结构,总体坐标系下的单元刚度矩阵和单元质量矩阵为,平面杆单元例子,平面杆单元例子,有限元离散的一般流程,几何体的离散逼近； 自由度选择与形函数的确定； 单元应变场的表示； 单元应力场的表示； 哈密顿原理或达朗伯-拉格朗日原理； 系统泛函变分； 单元内部积分； 单元坐标系到总体坐标系的转换； 总体刚度矩阵和质量矩阵的组装； 约束条件的处理或引入。,形函数的确定,形函数的阶数与单元形式和自由度数有关,以平面三结点三角形单元为例,2、结点位移表示待定系数,1、假设域内位移的插值函数,形函数的确定,3、解出系数，并合并关于结点位移的项次，得形函数,其中,注意下标轮换,形函数的确定,4、确定由结点位移表示的域内位移场，并表示为矩阵形式,形函数的确定,形函数的几个性质,结点上的插值函数满足 单元中任意点处形函数和应有 选取多项式时，常数项和一次项要完备 为表达统一起见，可采用参数变换,单元内部积分,单元刚度矩阵和单元质量矩阵的一般形式,单元等效结点载荷的一般形式,单元内部积分,定积分的数值积分（Guass积分）,核心思想 通过若干个点的函数值的加权组合计算,单点积分 两点积分,单元内部积分,Guass积分点位置及权系数,单元内部积分,二维和三维问题,约束条件的处理,直接代入法 特点 方程阶数降低,约束条件的处理,对角元素置1法 特点 适用于零位移情况,约束条件的处理,对角元素置大数法 特点 相当于罚函数法,等参变换,等参变换的基本思想就是局部系向总体系的映射,局部坐标系向物理坐标系的映射 类似于构造形函数，可构造物理坐标和局部坐标之 间的关系。仍然以三角形平面单元为例,等参变换,需要注意的地方,等参变换,等参变换,对面积分来说,三维问题？,等参变换,等参变换后的单元刚度矩阵、单元质量矩阵格式,等参变换,单元等效结点载荷的一般形式,第三章 大型系统特征值问题,系统运动方程,对特征值问题，可建立系统运动方程如下,自由振动解,将式（2）代入式（1），得,系统运动方程,由式（3）得特征方程,解得n个特征根及相应的特征向量,n次代数方程,另一种形式,特征向量正交性,取,结合式（6），得,特征向量正交性,进一步得正交性关系,一个重要的关系,刚度矩阵为半正定矩阵 质量矩阵为半正定矩阵,半正定矩阵的特征值问题,特征向量 特征值,2自由度系统的例子,Rayleigh商,对给定的振型向量，可定义Rayleigh商,Rayleigh商的上下界,考虑式（12），可将式（15）表示为,Rayleigh商,将式（16）展开，得,物理意义 假设模态相当于增加了约束，使得系统刚度提高,误差估计,设 和 是标准特征值问题的近似解，定义残差,考虑下式,得,误差估计,定义向量 和 矩阵的2范数,对式（21）两边取2范数，整理得,误差估计,由于,将式（24）代入式（23）得,式（23）说明，只要残差的2范数很小，给定的 就接近系统的某一阶特征值,逆迭代法,基本假定 刚度矩阵非奇异 迭代方程 考虑到任意向量都可由特征向量表示为,逆迭代法,综合式（26）和式（27），得,逆迭代法,从式（30）中容易看出，当 考虑Rayleigh商,逆迭代法,结合式（27）、（28）和（29）可将式（33）化为,逆迭代法,逆迭代法的一般步骤 给定初始迭代向量； 若已计算出前i-1阶模态，根据已得到的特征向量进行正交化处理；,逆迭代法,逆迭代法的一般步骤 解方程； 计算特征值,逆迭代法,几个要点 正交化处理； 归一化处理； 当刚度矩阵奇异时的处理； 系统包含近频时的处理,移 轴 法,带移轴的逆迭代法,考虑特征值问题,其特征值与原系统特征值之间满足如下关系,由式（26）给出的迭代关系可得,怎么定？,Rayleigh商迭代法,Rayleigh商迭代法 以Rayleigh商作为移轴量,计算流程,变换法,变换法的基本思想是通过迭代的方法逐次构造变换矩阵 ，使矩阵 和 趋近对角形，其一般迭代思路如下式所示,当 时，应有,变换法,并且有,两类应用广泛的变换法 雅可比变换法 Householder-QR法,雅可比变换法,适用范围 标准特征值问题，即 考虑标准特征值问题,取正交矩阵,使得,雅可比变换法,雅可比法中选择的变换矩阵为旋转矩阵，在一次变换中使待变换矩阵中的一个非对角元变为零。 例,雅可比变换法,将 变换为零，即 ，由此得,雅可比变换法,定义阀值 收敛条件,雅可比变换法,收敛性检验 考虑特征值问题 ，其中,由式（47）可计算得,雅可比变换法,经第一次循环，对非对角元（1，2）消零，得,由式（52）可计算得,雅可比变换法,类似的，分别对非对角元（1，3）和（2，3）消零，得,雅可比变换法,雅可比变换法的一般步骤,给定当前循环的阀值，一般取第m次迭代的阀值为10-2m； 对上三角块的非对角元按式（48）计算，若计算结果大于给定的阀值，则按式（47）进行旋转变换； 检查式（49）式是否满足，若满足则由迭代最终步骤计算出的矩阵给出特征值，计算特征向量，否则转向第1步，开始新一轮循环。,广义雅可比变换法,对雅可比变换矩阵进行修改，得转换矩阵,广义雅可比变换法,系数的选取应使刚度矩阵和质量矩阵在变换后其第i行j列的元素同时为零，即,若 ，可取 ，否则,广义雅可比变换法,其中,广义雅可比变换法,例1，用雅可比法求解图示系统的固有频率,广义雅可比变换法,广义雅可比变换法,精确解,Householder-QR变换法,Householder-QR变换法可大致分为三个步骤,Householder变换； 特征值的QR迭代计算； 特征向量计算。,Householder-QR变换法,Householder变换的目的在于将矩阵变换为三对角阵，加速QR迭代收敛速度。 QR迭代 考虑到实对称矩阵可由正交矩阵和上三角矩阵分解为,Householder-QR变换法,由此定义迭代格式,当 时， ，,Rayleigh-Ritz法,若假设振型与系统的第 i 阶主振型相差一阶小量，则Rayleigh商与系统的第 i 阶特征值相差二阶小量。 假设近似振型为,其中,里兹基向量,里兹坐标,Rayleigh-Ritz法,定义Rayleigh商,取上式的驻值，得,Rayleigh-Ritz法,由式（65）可以计算出新系统的特征值和特征向量,满足,则原系统的特征值和特征向量为,给出的近似特征值前一半精度较高，因此，若需要求解前 r 阶模态，要先定义出 2r 阶假设振型； 经变换后，系统规模大大降低； Rayleigh-Ritz法计算结果依赖于里兹基的选择； 里兹基可通过静力问题求得。,Rayleigh-Ritz法,Rayleigh-Ritz法的特点,Rayleigh-Ritz法,例2，用Rayleigh-Ritz法计算下述广义特征值问题的近似解。,假设载荷为,Rayleigh-Ritz法,由新系统特征值问题解得,精确解,Rayleigh-Ritz法,如果假设载荷为,由新系统特征值问题解得,Rayleigh-Ritz法,进一步得,子空间迭代法,子空间迭代法假设 r 个初始向量同时进行迭代，计算前 p 个特征值和特征向量。子空间迭代法可以看为是逆迭代法的推广，算法步骤类似，其具体迭代步骤可表示为,1、初始计算 （1） 形成刚度矩阵 和质量矩阵 ； （2） 按所给定的边界条件修正刚度矩阵 ； （3） 三角分解 ，得,子空间迭代法,2、给定初始迭代向量矩阵 3、对每次迭代 ，按下面步骤执行 （1） 赋值： （2） 求解方程组： （3） 赋值： （4） 赋值： （5） 赋值： （6） 解广义特征值问题：,子空间迭代法,（7） 检验特征值是否满足精度要求： 若满足精度要求，则转到第 4 步，若不满足，则 （8） 赋值： （9） 令： ，并返回步骤（2）； 4、赋值： 5、输出结果：,第四章 系统动力学响应,系统运动方程,对一般的振动响应，考虑初值问题可建立系统运动方程如下,求解方式,振型叠加法 直接积分法,振型叠加法,引入坐标变换,其中,广义位移,由式（2）和式（1）得,振型叠加法,若阻尼矩阵可解耦，则方程可变换为,利用Duhamel积分，得,直接积分法,直接积分法不对运动方程进行任何变换，直接对方程进行积分求解。其原理基于两个基本概念,将求解域上任意时刻均需满足运动方程的要求更改为在一定条件下近似的满足运动方程，如，在若干离散的时间点上满足方程； 在一定数目的时间间隔区域内，假设位移、速度和加速度的近似插值函数。,误差来源,截断误差 近似表达式中略去高阶小量； 舍入误差 超过计算机字长的四舍五入。,直接积分法解的稳定性,无条件稳定 对任意初始条件，无论时间步长怎么选取，积分结果都不会无界增大； 条件稳定 为保证积分结果的有界性，要求时间步长必须小于某个临界时间步长。,考虑直接积分法递推格式,在讨论稳定性时，可令 ，由式（7）得,给定的初始条件,递推矩阵,直接积分法解的稳定性,对递推矩阵进行谱分解，得,递推矩阵为二阶非对称矩阵，且有重特征根，其约当标准形为,式中 为递推矩阵的约当标准形，其对角元为递推矩阵的特征值。,当 时，只有 ，才能使 在 时有界,当 时，只要 ，就能使 在 时有界,直接积分法解的稳定性,若递推矩阵为二阶非对称矩阵，没有重特征根，其约当标准形为,易知，只要满足,时， 就有界,递推矩阵的谱半径,直接积分法解的稳定性,稳定性准则,若递推矩阵没有重特征根，要求其谱半径不大于1； 若递推矩阵存在重特征根，要求其特征值的模小于1。,数值阻尼（人工阻尼） 谱半径小于 1.0 时，随迭代增加，矩阵,中心差分法,将位移函数按Taylor级数展开，得,由式（13）可得,中心差分法,将式（14）代入任意时刻的运动方程,得,其中,中心差分法,起始步的计算,定义积分常数,中心差分法,计算步骤,1、初始计算 1） 形成刚度矩阵 和质量矩阵 和阻尼矩阵 ； 2） 给定初始条件 和 ，并计算 ； 3） 给定时间步长，计算积分常数； 4） 计算起步位移： ； 5） 计算有效质量矩阵： ； 6） 对有效质量矩阵进行三角分解： ；,中心差分法,2、对每一个时间步 1） 计算时刻 的有效载荷 ； 2） 计算 时刻的位移： ； 3） 如有需要，计算时刻 的加速度和速度,中心差分法,中心差分法的数值稳定性,对任意时刻的运动方程,考虑其中心差分格式，有,将式（21）写为递推矩阵格式有,中心差分法,其中,对无阻尼系统，由式（23）计算出递推矩阵特征值为,中心差分法,根据收敛条件，可知,从而得,对有阻尼系统，中心差分法的临界时间步长为,中心差分法,由于对任意一个解耦后的运动方程来说均应满足式（27），因此，满足收敛条件的时间步长由系统中的最高阶模态对应的固有频率决定。容易看出，当满足式（25）时，谱半径为 1 ，此时中心差分法没有数值阻尼。,Newmark法,在位移和速度的Taylor级数中只保留一阶导数项，则可得,若在式（28）中将导数的平均值代替 时刻的导数值，则可将式（28）改写为,Newmark法,将式（29）中的第二式代入第一式得,从式（29）中可以看出，加速度近似取平均加速度，因此该方法又称为平均加速度法。若在积分过程中对加速度取不同近似格式，可得到各种不同的积分格式。,Newmark引入如下速度、位移关系,未知量,Newmark法,由式（31）中的第二式可解出,再将式（32）代入式（31）中的第一式可解出,考虑 时刻的系统运动方程,Newmark法,将式（32）和式（33）代入式（34），得,简记为,Newmark法,令,并记,得,计算步骤,1、初始计算 1） 形成刚度矩阵 和质量矩阵 和阻尼矩阵 ； 2） 给定初始条件 和 ，并计算 ； 3） 给定时间步长，计算积分常数； 4） 形成有效刚度矩阵 ； 5） 对有效刚度矩阵进行三角分解 ；,Newmark法,2、对每一个时间步 1） 计算 时刻的有效载荷 ； 2） 计算 时刻的位移： ； 3） 计算 时刻的加速度和速度,Newmark法,Newmark法的数值稳定性,仍然讨论解耦后的无阻尼模态运动方程,将Newmark法格式的加速度近似表达式代入，得,利用式（31），并考虑到如下关系,Newmark法,得,从式（42）中可联立解得,由于时间的任意性，可将上式改写为,Newmark法,将上式代入式（40），得,写成递推矩阵格式，有,Newmark法,其中,将上式代入式（40），得,Newmark法,递推矩阵的特征值为,要求Newmark法稳定，即要求,Newmark法,由式（49）的第一式和第二式可分别得到,从式（51）可解出，Newmark法无条件稳定的条件为,Newmark法,若不满足式（52）所示条件，则算法是条件稳定的，要得到稳定的解，要求计算时间步长小于临界时间步长，即,