（原创，转载请注明出处）

1 概述

本系列文章研究成熟的有限元理论基础及在商用有限元软件的实现方式，通过

（1）   基础理论

（2）   商软操作

（3）   自编程序

三者结合的方式将复杂繁琐的结构有限元理论通过简单直观的方式展现出来，同时深层次的学习有限元理论和商业软件的内部实现原理。

有限元的理论发展了几十年已经相当成熟，商用有限元软件同样也是采用这些成熟的有限元理论，只是在实际应用过程中，商用CAE软件在传统的理论基础上会做相应的修正以解决工程中遇到的不同问题，且各家软件的修正方法都不一样，每个主流商用软件手册中都会注明各个单元的理论采用了哪种理论公式，但都只是提一下用什么方法修正，很多没有具体的实现公式。商用软件对外就是一个黑盒子，除了开发人员，使用人员只能在黑盒子外猜测内部实现方式。

                                            

一方面我们查阅各个主流商用软件的理论手册并通过进行大量的资料查阅猜测内部修正方法，另一方面我们自己编程实现结构有限元求解器，通过自研求解器和商软的结果比较来验证我们的猜测，如同管中窥豹一般来研究的修正方法，从而猜测商用有限元软件的内部计算方法。我们关注CAE中的结构有限元，所以主要选择了商用结构有限元软件中文档相对较完备的Abaqus来研究内部实现方式，同时对某些问题也会涉及其它的Nastran/Ansys等商软。为了理解方便有很多问题在数学上其实并不严谨，同时由于水平有限可能有许多的理论错误，欢迎交流讨论，也期待有更多的合作机会。

通用结构有限元软件iSolver介绍视频：

        

                                          http://www.jishulink.com/college/video/c12884

==第32篇：谐响应分析算法==

在系列文章30：谐响应分析原理中，我们介绍谐响应分析的原理，但原理怎么转换为计算机的算法并在程序中实现的过程没有提及，在Abaqus中，谐响应分析求解可以分为模态叠加法、直接法和子空间法，其中前两者最常用，在这章中我们将介绍谐响应分析的这两种求解算法，并对上一章的例子在iSolver中采用这两种算法后和Abaqus进行结果比对分析，证明iSolver的结果和Abaqus无任何误差。

1.1 谐响应分析算法方程

谐响应分析是对瞬态动力学的简化，对一般的瞬态动力学，由于有质量和阻尼响应，方程可以写为：

                                             

M、C、K分别表示质量阵、阻尼阵和刚度阵。

对任意一个外载荷R，譬如随时间的变换如下：

将曲线做FFT变换，外载荷可以表示为下述表达式：

可以得到F随频率变化的输入载荷曲线：

对于每个频率点载荷

，最终系统稳定后的响应也是频率的正弦函数，就是

所以，可以将瞬态动力学表示式中的位移、加速度等都表示为随时间正弦变化的函数。带入方程后，就只剩K和C随时间的表达式包括时间t变量，下面就是具体采用哪种方法得到这个K和C的问题。

1.2 直接法

谐响应分析都是假定系统变形量比较小，且应力应变变化曲线在弹性范围内，同时也不会在变形中产生碰撞，所以是个在某个初始位置的线性的小扰动的分析步，这也是为何谐响应分析在Abaqus设置中位于Linear Pertubation的分析组内的原因。如果这个初始位置是通过非线性分析过程得到的，那么从整体分析来说，也可以说是包括了非线性效应。

线性假设下，K可以认为常量，方程中只剩阻尼C的表达式，一种简单的表达式就是C认为是K和M的线性组合。

alpha和beta分别是质量和刚度正比的阻尼系数。Abaqus中有两种方法设置：

（1）   针对每种材料中分别设置damping，下方例子将采用这种方式。

（2）   针对所有材料设置全局的damping，必须在inp中采用关键词才能设置。Abaqus的CAE中直接法也应该在step设置界面中有类似模态叠加法设置一个全局damping，但不知为何没有实现。

去掉正弦项，那么整个方程可以完全转换为频率的方程。

上式可以直接求解出u的值，u表示振幅。这种方法称为直接法，顾名思义，就是在每个频率点通过对模型的原始方程直接积分计算系统的稳态谐波响应。当然，上式严格来说应该是复数形式，为了简单起见，上述公式约去了复数的表达式，同时C也不是只有K和M正比的阻尼，还有其它阻尼，具体请参考有限元理论相关书籍，我们在此不再累述。

1.3 模态叠加法

谐响应分析的目标之一是捕捉导致大振幅的输入频率，而产生大的振幅一般都是由于和结构的共振引起的，但直接法由于没有做模态分析，那么对与输入的频率只能人为的划分，很多时候无法一次性的精确捕捉到共振频率，同时，直接法的公式也决定了所有的自由度都参与了计算，这导致结果计算相对较慢。为了克服这些缺点，有限元中还可以使用模态叠加法来求解谐响应分析。

模态叠加法，在求解模态稳态响应之前先提取无阻尼系统的特征模态，然后通过变换解耦系统为一组使用模态振型表示的单自由度方程。求解各单自由度方程得到系统在模态振型下的稳态响应后，通过变换最后求得系统的稳态响应。譬如下方的悬臂梁，各个颜色的曲线表示不同阶次的模态振型。

那么，模态叠加法的本质是什么呢？在数学上，空间中的任意物理量可以在不同的基上进行分解。

基的形式只要满足两个条件：

（1）   正交性：构成这个空间的基是相互正交的。

（2）   完备性：空间任意一个量可以空间中由这组基线性表示

举个例子，譬如我们二维空间的存在一组基如下：

显然，满足上面两个条件。任意一个u可以表示为u1和u2的线性表示

其中，fi是线性系数，如果U=（1，1），那么

；我们可以取另一组基函数，譬如

譬如，下面板一个角点受激励的情况：

通过对比和检查，很容易得到在小变形振动下整个系统对这个载荷的响应主要受下方四个频率代表的基函数影响，所以，谐响应分析时我们也只要取这四个频率的基函数就足够了。

在Abaqus或者iSolver程序中，模态叠加法的谐响应分析也有两步：

（1）进行模态分析Freq，求广义特征值和特征向量。

（2）按用户输入的频率范围求出所需频率点，然后对各个频率点，求各模态的系数，对模态位移做线性叠加。

1.4 算法对比

1.4.1 适用范围

模态叠加法求得的系统动力响应本质上是系统各阶模态振型的线性组合，因此根据模态分析求解特点，可以知道模态叠加法在一些情况下并不适用，而直接法具有更广泛的适用范围，主要体现如下：

（1）   材料具有频变特性，由于各阶模态无法体现材料的频变特性，因此在此情况下模态叠加法不适用，直接法支持具有频率变化特性的粘弹性材料，如下Abaqus设置。

（2）   系统具有非对称刚度矩阵，模态分析无法求解系统各阶模态，因此在此情况下模态叠加法不适用；

（3）   系统具有除模态阻尼以外的其它阻尼，根据模态叠加法的原理可知模态叠加法仅能考虑模态阻尼，因此在此情况下模态叠加法不适用。

（4）   在声振耦合中，需要谐响应分析首先计算出结构的稳态频率响应，直接法和模态叠加法都适用。

1.4.2 计算效率

由于直接法在扫频范围内的所有频率点处都需要计算系统整体的刚度阵、质量阵和阻尼阵，显然效率上要远低于模态叠加法。同时，下面的例子也说明了，直接法不一定一次就能精确定位到共振峰，存在多次计算的可能性。

1.4.3 求解精度

直接法是在每个频率点对系统进行复杂积分运算的，因此具备更高的求解精度。同时，需要注意的是用户在使用模态叠加法时应合理选择模态求解阶数以保证足够精确描述系统的动力学特性。

1.5 算例

1.5.1 模型介绍

模型依然采用30章：谐响应分析原理的悬臂梁模型，模型参数如下：

几何：悬臂梁长1000，I型截面参数如下，厚度都是5：

材料：杨氏模量210000，泊松比0.3，密度7.85e-9。

网格：取S4R单元类型，网格大小为12.5。

左端约束：

1.5.2 直接法

设置直接法分析步，频率范围如下：

和模态叠加法不同，我们这儿设置材料的结构阻尼：0.05

设置载荷：另一端加Y方向载荷方向1

查看载荷节点在Y方向的谐响应振幅曲线如下：（左侧是Abaqus结果，右侧为iSolver结果，结果无任何误差），可发现在0-250Hz的频率范围内，只有一个共振峰114.4Hz。

最后一个频率的最大响应位移分布如下，可见iSolver和Abaqus也是没有误差：

1.5.3 模态叠加法

我们把系列文章30：谐响应分析原理中模态叠加法的结果贴过来：

1.5.4 结果比较

可发现直接法只在求解频率域内有20个频率点，由于点太少，无法精确捕捉到共振峰，所以，只能再次计算，此时将频率设为200个点，由于没有模态频率的信息，也无需用bios=3来划分了。

把直接法和模态叠加法的结果放在一起，此时得到位移响应如下：

可以发现由于我们设置的阻尼是完全一致的structural阻尼，用两种方法计算得到的结果差异非常小。

1.6 视频讲解和操作验证演示

模态叠加法的算例演示可看下方视频

https://www.jishulink.com/college/video/c12884 20理论系列文章30-谐响应分析原理

 

==总结==

本文首先分别介绍了直接法和模态叠加法，并对两种方法从适用范围、计算效率和求解精度三个方面进行了对比分析，最后通过算例比较iSolver和Abaqus求解结果，对iSolver进行验证。

以往的系列文章：

1.6.1 ========第一阶段========

第一篇：S4壳单元刚度矩阵研究。

http://www.jishulink.com/content/post/338859

第二篇：S4壳单元质量矩阵研究。

http://www.jishulink.com/content/post/343905

第三篇：S4壳单元的剪切自锁和沙漏控制。

http://www.jishulink.com/content/post/350865

第四篇：非线性问题的求解。

http://www.jishulink.com/content/post/360565

第五篇：单元正确性验证。

https://www.jishulink.com/content/post/373743

第六篇：General梁单元的刚度矩阵。

https://www.jishulink.com/content/post/403932

第七篇：C3D8六面体单元的刚度矩阵。

https://www.jishulink.com/content/post/430177

第八篇：UMAT用户子程序开发步骤。

https://www.jishulink.com/content/post/432848

第九篇：编写线性UMAT Step By Step。

http://www.jishulink.com/content/post/440874

第十篇：耦合约束（Coupling constraints）的研究。

https://www.jishulink.com/content/post/531029

1.6.2 ========第二阶段========

第十一篇：自主CAE开发实战经验第一阶段总结。

http://www.jishulink.com/content/post/532475

第十二篇：几何梁单元的刚度矩阵。

http://www.jishulink.com/content/post/534362

第十三篇：显式和隐式的区别。

http://www.jishulink.com/content/post/537154

第十四篇：壳的应力方向。

https://www.jishulink.com/content/post/1189260

第十五篇：壳的剪切应力。

https://www.jishulink.com/content/post/1191641

第十六篇：Part、Instance与Assembly。

https://www.jishulink.com/content/post/1195061

第十七篇：几何非线性的物理含义。

https://www.jishulink.com/content/post/1198459

第十八篇：几何非线性的应变。

https://www.jishulink.com/content/post/1201375

第十九篇：Abaqus几何非线性的设置和后台。

http://www.jishulink.com/content/post/1203064

第二十篇：UEL用户子程序开发步骤。

https://www.jishulink.com/content/post/1204261

1.6.3 ========第三阶段========

第二十一篇：自主CAE开发实战经验第二阶段总结。

https://www.jishulink.com/content/post/1204970

第二十二篇：几何非线性的刚度矩阵求解。

http://www.jishulink.com/content/post/1254435

第二十三篇：编写简单面内拉伸问题UEL Step By  Step

http://www.jishulink.com/content/post/1256835

第二十四篇：显式求解Step By Step。

https://www.jishulink.com/content/post/1261165

第二十五篇：显式分析的稳定时间增量。

http://www.jishulink.com/content/post/1263601

第二十六篇：编写线性VUMAT Step By Step。

https://www.jishulink.com/content/post/1266640

第二十七篇：Abaqus内部计算和显示的应变。

https://www.jishulink.com/content/post/1273788

第二十八篇：几何非线性的T.L.和U.L.描述方法

https://www.jishulink.com/content/post/1282956

第二十九篇：几何非线性的T.L.和U.L.转换关系

https://www.jishulink.com/content/post/1286065

第三十篇：谐响应分析原理

https://www.jishulink.com/content/post/1290151

第三十一篇：自主CAE开发实战经验第三阶段总结

https://www.jishulink.com/content/post/1294824