维普资讯 http://www.cqvip.com 技术讲座Ｔｅｃｈｎｉｃａｌ Ｌｅｓｓｏｎ 第七讲 动力学问题的有限元方法 Ｌｅｓｓｏｎ ７ Ｔｈｅ Ｆｉｎｉｔｅ Ｅｌｅｍｅｎｔ Ｍｅｔｈｏｄ ｆｏｒ Ｄｙｎａｍｉｃｓ 陈乐生 （福州大学机械工程及自动化学院） 动力学问题是一个带普遍性的问题，工程机 构或结构系统的设计和开发离不开系统的动力学 特性的分析，这是控制系统设计的基础和前提。 由于动力学模型的高阶和非线性以及耦合问题的 存在，解决动力学问题一般情况下都是采用数值 解法，有限元分析不失为一种有效的具有较高精 度的数值方法。这一节主要介绍采用有限元方法 求解动力学问题的基本原理和步骤，讨沦如何建 立弹性结构的离散化运动微分方程、包括质量矩 阵和阻尼矩阵的形成，讨论运动微分方程求解的 直接积分法和振型叠加法。动力学分析中除了有 限元本身的特点外，更多地存在数学问题，因此， 更加深入的研究请读者参考有关文献资料。 １运动微分方程的建立 １．１单元内任意一点的位移、速度和加速度 在前面几讲中，我们已经知道，对于第ｅ个单 元，单元内任意一点的位移｛“（ｔ）｝可以用单元的 节点位移｛６¨ （ｔ）｝表示为： ｛“（ｔ）｝＝［ ］ （ｔ）｝ （１．１） 式中，［ ］为插值函数或者称之为形函数，显 然，该函数与时问ｔ无关。 于是，单元内任意～一点的速度和加速度分别 为： ｛“ （ｔ）｝＝［Ｈ］｛ （ｔ）｝ （１．２） ｛“ （ｔ）｝＝［Ｈ］｛６ ‘ （ｔ）｝ （１．３） １．２单元虚功和虚应变能 根据弹性理论，第ｅ个单元虚应变能为： ６Ｕ ’＝』』』（ ） （Ｏ－）ｄＹ （１．４） 其中，单元应变矩阵（占）＝［Ｂ］（ ’，［Ｂ］为 几何矩阵，（ ’为单元节点位移向量； 单元应力矩阵（Ｏ－）＝［Ｄ］［Ｂ］，［Ｄ］为弹性矩 阵。 引入单元刚度矩阵，上式可以改写为： ４２ ６Ｕ‘ ＝｛６‘ （ｔ）｝ ［ ‘ ］｛６ （ｔ）｝ （１．５） 其中， （ｔ）｝为单元节点位移的变分，即 节点的虚位移； ［ ］＝』』』［Ｂ］７１［Ｄ］［Ｂ］ｄ 为单元刚度矩 阵。 单元外力的虚功为： ６ ‘ ＝｛６‘ （ｔ）｝ ｛Ｑ‘ （ｔ）｝ （１．６） 在单元载荷向量｛ｐ ’（ｔ）｝中，包括体积力、表 面力，惯性力，如果考虑阻尼作用，还要计人与速 度有关的阻尼力，正如前面几讲中提到的，这些外 力全部等效到节点上。有关这些力的计算这里不 再详细介绍。 １．３运动微分方程的建立 根据虚功原理，外力的虚功将全部转化为虚 应变能。将外力、惯性力以及阻尼力代人虚功方 程，可得运动微分方程： ∑（｛６‘ （ｔ）） 【 ‘ ］｛６‘ （ｔ））） ＝ ［｛６ （ ）｝ （｛Ｑ （ ）｝＿［ｍ］｛６触 （￡）｝一 ［ｃ］－｛１６引 （ｔ）｝］ （１．７） 将所有单元集成可得： ｛６Ｕ（ｔ）｝ Ｋ］｛Ｕ（ｔ）｝＝｛６Ｕ（ｔ）｝ （｛Ｒ（ｔ）｝＿ ［ ］｛６眦 ｝一［ｃ］｛ （ｔ）｝） （１．８） 式中，｛Ｕ（ｔ）｝为总体位移矢量；［Ｋ］为总体刚 度矩阵；［ ］为总体质量矩阵；［ｃ］为总体阻尼矩 阵；｛尺（ｔ）｝为总体载荷向量。 由虚位移的任意性，上式可进一步简化为： ［ ］｛ （ｔ）｝＋［ｃ］｛ （ｔ）｝＋［ ］｛Ｕ（ｔ）｝ ＝｛Ｒ（ ）｝ （１．９） 上式即为弹性结构离散化的运动微分方程， 是动力学问题有限元分析的基本方程。这是一个 二阶线性微分方程组，在求解之前必须引入支撑 条件、即根据位移边界条件进行修改。比如边界 零位移，弹性结构存在刚性位移，机构的基础存在 ＷＯＯＤＷＯＲＫＩＮＧ ＭＡＣＨＩＮＥＲＹ 维普资讯 http://www.cqvip.com
Ｔｅｃｈｎｉｃａｌ Ｌｅｓｓｏｎ技术讲座 运动等等情况的修改请参阅有关文献。 ２质量矩阵与阻尼矩阵 运动微分方程的建立，除了我们前面几讲中 着重讨论过的总体刚度矩阵外，还涉及到总体质 量矩阵和总体阻尼矩阵，下面简单介绍这两类矩 阵的形成。 ２．１单元质量矩阵的计算方法 在通用程序中，单元质量矩阵采用两种计算 方法：一致质量矩阵和集中质量矩阵。 ２．１．１一致质量矩阵 一致质量矩阵的形式为： ［ ］＝ｆ ｆ ｆＰ［Ｈ］ ［Ｈ］ｄＶ （２．１） 式中，Ｐ为密度；［日］为插值函数矩阵或称形 函数矩阵。该矩阵是对称和正定的。但是根据该 矩阵的定义可以看出，在组集整体质量矩阵时，非 零元素会占据很多存储空间并且耗费机时。ｌＩ ２．１．２集中质量矩阵 Ｌ 一一致质量矩阵是考虑单元内每个质点都存在 １ Ｏ Ｏ Ｏ Ｏ Ｏ 质量，如同单元外力的处理，集中质量是将单元的 Ｏ １ Ｏ Ｏ Ｏ Ｏ 质量近似地集中分配到各个节点，这样得到的单 元质量矩阵为一列 角阵： ｍｌ】０ ０ ０ ｍ２２ ０ ［ ］＝ Ｍ Ｍ Ｍ （２．２） ０ ０ ｍｎ 由于质量近似分配，计算结构振型时要比采 用一致质量会差些。 ２．１．３若干种单元的质量矩阵 （１）三维杆单元 对于密度为Ｐ、横截面积为Ａ、单元长度为 的 三维杆单元，通常采用集中质量矩阵，即将单元质 量 等分为二集中作用于单元两端，则有： 『 ０。。０。０ ＝ １ 三 三 （２．３） ０ ０ ０ ０ ｌ Ｌ０ ０ ０ ０ ０ （２）三节点三角形平面单元 木工机床２００７ Ｎｏ．４ 假设单元均质等厚，密度为Ｐ，厚度为ｔ，单元 面积为△，则集中质量矩阵为整个单元质量等分到 三个节点上： ０ ０ ０ ０ ０ ０ ０ ０ １ ０ ０ ０ ０ １ ０ ０ （２．４） ０ ０ １ ０ ０ ０ ０ １ 建立其一致质量矩阵时利用三角形单元的面 积坐标与形函数的关系： Ｌ ：Ｈ （ ，Ｙ） ｉ＝１，２，３ 代人（２．１）式，可得： ／２ ０ １／４ ０ １／４ ０］ ０ １／２ ０ １／４ ０ １／４ ｆ ／４ ０ １／２ ０ １／４ ［ｍ］＝ ０ １／４ ０ １／２ ０ ／４ ０ １／４ ０ １／２ ０ １／４ ０ １／４ ０ （２．５） （３）八节点四逆 等参数单元 假设四边形单 的面积为 ，厚度为ｔ，根据定 义，一致质量矩阵为： ［ｍ］＝ｐ ［Ｈ］ ［Ｈ］ｔｄｘｄｙ （２．６） 但是，Ｈ （ ，ｔ ７）是局部坐标系的函数，因此，将 上式改在母单元内积分： ［ｍ］＝ｐｌ ［Ｈ］ ［Ｈｉｔ（ｄｅｔ．，）晦咖（２．７） 式中，ｄｅｔＪ为雅克比行列式，雅克比矩阵反映 了整体坐标系与母单元局部坐标系的关系，形为： ｒ堕 ］ ．，：ｌ荔 ｊ若 若ｌ Ｊ｛ 若采用集中质量矩阵，则应将单元质量平均 分配到８个节点上，即： ［ｍ］＝÷Ｐ Ｌ ｍ ｊ Ｊ一１ Ｊｆ Ｉ ｉ ｔ（ｄｅｔ Ｊ）ｄ４ｄｎ １ ０ ＾０ Ｌ ０ １ ＾０ ０ Ｍ Ｍ ｏ Ｍ Ｍ （２．８） ０ ０ ＾ １ ０ ０ ０ ＾０ １ ４３ 维普资讯 http://www.cqvip.com 技术讲座Ｔｅｃｈｎｉｃａｌ Ｌｅｓｓｏｎ ２．２ 阻尼矩阵 ［ ］｛ 船（ｔ）｝＋［Ｋ］｛Ｕ（ｔ）｝＝０ ｛Ｕ（ｔ）｝：｛ ｝ｓｉｎ￣ｔ 有频率。 （３．１） （３．２） 工程结构的阻尼力主要来自两个方面，一是 由结构周边的粘性介质产生的粘性阻尼，二是材 此为常系数线性齐次微分方程组，通解为 这里，｛ ｝称为结构的振型， 称为结构的固 料内部分子问的磨擦产生的结构阻尼；前者与构 件的运动速度成正比，后者与弹性体内部应变的 速率成正比。 ２．２．１粘性阻尼力 粘性阻尼力定义为 ｛｛｝＝ Ｐ。｝＝ Ｐ ｑｏ ｛｝｛｝＝ｑｏ［Ｈ］Ｕ ｑｏ ｌ Ｊ｛Ｄ一｝ ｝ （２．９・ ） 【，‘ 如果将该阻尼力视为体积力，其等效的单元 节点阻尼力矢量为 ｛Ｆ‘ Ⅲ［Ｈ Ｐ ｝ｄＶ：［Ｃ ］｛６ ｝ （２．１０） 由此可得： ［Ｃ ］： Ⅲ［Ｈ］ ｌｐ［Ｈ］ｄＶ：ＯＺ［ｍ］ （２．１１） ２．２．２结构阻尼力 结构阻尼力定义为 ｛Ｐ ｝＝ ［Ｄ］导｛【，‘ ｝： ［Ｄ］［Ｂ］｛６酬 ｝ （２．１２） 其等效的单元节点阻尼力矢量为 ｛Ｆ‘ ｝＝Ⅲ［Ｂ］ ｛Ｐ ｝ｄＶ：［Ｃ ］｛６酬 ｝ （２．１３） 由此可得： ［Ｃ ］： ［Ｂ］ ’［Ｄ］［Ｂ］ｄＶ＝ ［ｋ ］（２．１４） 在一般情况下，认为整体阻尼矩阵可以写作 整体刚度矩阵和整体质量矩阵的线性组合，称为 瑞利阻尼： ［Ｃ］＝ＯＺ［ｍ］＋ ［Ｋ］ （２．１５） 以上各式中，比例系数 和 可以通过测定 结构自由振动的衰减率经过换算得到。实际上， 阻尼的作用在一般情况下对结构的动态响应影响 并不大，因此，上述处理方法会取得运算上的方 便。 ３弹性结构的无阻尼自由振动 弹性结构的动力学分析一般都要计算其固有 频率的振型，即在动力学方程（１．９）中，令载荷向 量为零，并忽略阻尼的影响，于是得到结构自由振 动的微分方程为： 将上述通解代人微分方程（３．１）式，可得： （［Ｋ］一Ａ［ ］）｛ ｝＝０ （３．３） 式中，记Ａ： ，Ａ和｛ ｝在数学描述中分别 称为广义特征值和广义特征向量。该齐次线性方 程组具有非零解的条件是系数行列式为零，即： ［Ｋ］一Ａ［Ｍ］：０ （３．４） ４基本方程的求解 在满足初始条件下，通过积分运算求解动力 学方程（１．９）式，以得到不同时刻弹性结构的动态 响应，包括位移、速度、加速度和力等。由于整体 质量矩阵、整体阻尼矩阵和整体刚度矩阵阶次都 较高，一般认为下面两种方法求解是较为有效。 ４．１直接积分法 直接积分法也叫逐步积分法，该方法的要点： 一是把时间历程等分为０，Ａｔ，２△￡，…， ，只要求相 隔△￡的各个离散时刻满足微分方程，而不要求所 有时刻都满足；二是在每个时问问隔△￡内，假定位 移、速度、加速度按照某种规律变化。有限元分析 中通常采用两种变化规律，即威尔逊（ｗｉｌｓｏｎ） 法 和纽马克（ｎｅｗｍａｒｋ）法。 ４．１．１威尔逊（ｗｉｌｓｏｎ） 法 威尔逊 法假定在ｔ到ｔ＋ｅａｔ内加速度呈线 性变化，并且证明了要达到无条件稳定，通常取 ：１．４。其数学描述为： （￡＋丁）｝＝ （￡）｝＋ （｛ 触（￡＋ ｅａｔ）｝一｛Ｕ础（ｔ）｝ （４．１） 式中，７－为时间增量，０≤７－≤ｅａｔ。 逐级积分上式可以得到ｔ＋△￡时刻单元内任 意一点的加速度、速度和位移。 ４．１．２纽马克（Ｒ，ｃｗｍａＦ￣）法 纽马克法同威尔逊 法一样，也是一种线性加 速度法的推广。该方法基于以下两个假设： ｛Ｕ （ｔ＋△￡）｝＝｛Ｕ （ｔ）｝＋［（１一 ｛ 触 （ｔ）｝＋６｛Ｕ （ｔ＋△￡）｝］△￡ （４．２） ＷＯＯＤＷＯＲＫｌＮＧ ＭＡＣＨｌＮＥＲＹ 维普资讯 http://www.cqvip.com
Ｔｅｃｈｎｉｃａｌ Ｌｅｓｓｏｎ技术讲座 ｛ｕ （ｔ＋△￡）｝：｛Ｕ（ｔ）｝＋｛ｕ （ｔ）｝△￡＋［（１ —｛（ｆ）｝ （ ［Ｍ］＋ ［Ｋ］）｛（ｆ）｝， ｆ２ 。（ｉ： ） ２ｏｒ）｛ｕ雌（￡）｝＋２ｏｒ｛Ｕ雌（￡＋△￡｝］ （４．３） 其中，系数６和 根据积分的精度和稳定性 的要求来确定。 已经证明，６ｉ＞０．５时， ｉ＞０．２５（０．５＋６） 纽马 ｌ ０ （ｉ￣ｊ）’ （ｔ）：｛（ｆ）｝ Ｒ（ｔ）｝，｛Ｒ（ｔ）｝为总载荷向量。 上式为一二阶常微分方程，可以采用积分方 克法是无条件稳定的。由上面两个假设以及运动 法求解。振型叠加法的特点就是在开始逐步积分 前从有限元的坐标系统变换成广义特征问题［ ］ ｛（ｆ）｝：Ａ［Ｍ］｛（ｆ）｝的特征矢量系。将系统位移转 微分方程可以求得ｔ＋△￡时刻位移向量、速度向量 和加速度向量。 有关以上两种方法的详细求解原理、过程以 及在有限元分析中的实施请参阅有关文献。 ４．２振型叠加法 直接积分法的基本要点是将整个动力学响应 的历程划分为若干时间步长，而为了求解的精度 要求，每个步长要求得很小，至少分为十步，每一 步都要对线性方程组进行求解，因此要耗费相当 多的机时，所以从一般意义上讲，计算短时间的动 力学响应，使用直接积分法比较有效，若长时间的 动力学响应采用该方法分析就不可取了。 振型叠加法是在逐步求解之前，现将运动微 分方程进行一些改造。考虑到弹性结构无阻尼自 由振动［ ］｛ ｝－Ａ［Ｍ］｛ ｝时的振型矢量具有对 质量矩阵和刚度矩阵的正交性，即有 Ｍ］ ：６ …） 【｛ Ｋ］｛ ｝：Ａ ６ 式中，６ｉｊ＝｛ 这样，如果采用振型矩阵｛（ｆ）｝作为变换矩阵， 该矩阵与时间ｔ无关，可将结构的整体位移写作结 构振型的线性组合，形如： ｛Ｕ（ｔ）｝＝ （ｔ）｛（ｆ）ｌ｝＋ （ｔ）｛（ｆ） ｝＋＾＋ （ｔ）｛（ｆ） ｝ （４．５） 由于振型可以在求解结构自由振动时得到， 因此，计算位移向量就转化为计算向量｛ （￡）｝。 将上式代人运动微分方程并经过整理，可得如下 关系式： ＆（ｔ）＋２ 。 ｌ＆（ｔ）＋ ２。 （ｔ）：Ｙ （ｔ） ｉ：１，２，＾， （４．６） 式中，０９ 为结构的固有频率； ，为第ｉ阶振型的阻尼比， 木工机床２００７ Ｎｏ．４ 换到以固有振型为基向量，这样处理对系统的性 质并无影响，而是以求解广义特征值为代价得到 非耦合的 个单自由度系统的运动方程以达到提 高计算效率的目的。但是，由于高阶广义特征值 的计算很困难也很费机时，而且，有限元分析的自 由度毕竟有限，对于高频振型的计算精度较差，因 此，通常只考虑部分低频及其相应振型、而略去其 余高频的影响。 ４．３特征值问题 由前面的分析我们已经知道，运动微分方程 的求解、包括固有频率、振型的确定以及振型叠加 法等，均涉及到特征值和特征向量的问题。但是 标准特征值问题得到的刚度矩阵和质量矩阵失去 了原有的带状分布和稀疏性，因此，一般通用程序 都是从广义特征方程直接求解。广特征值和 特征向量的求解通常可以采用广义雅克比法、行 列式搜索法和子空间迭代法等，这些方法的原理 和有限元分析中的实施请参考相关文献，这里不 再讨论。 参考文献 １．《有限元法概论》，龙驭球编，高等教育出版社， １９８０再版 ２．《弹性和塑性力学中的有限单元法》，谢贻权、何福 保编，机械工业出版社，１９８１第一版 ３．《固体力学有限元分析》，蒋维城编著，北京理工大 学出版社，１９９３再版 ４．《有限单元法基本原理和数值方法》，王勖成、邵敏 编著，清华大学出版社，１９９７年第二版 作者：陈乐生，福州大学机械工程及自动化学院教授 邮编：３５０００２ （收稿日期．２００７．１１．２９） ４５
