（原創，轉載請注明出處）

1 概述

本系列文章研究成熟的有限元理論基礎及在商用有限元軟件的實現方式，通過

（1）   基礎理論

（2）   商軟操作

（3）   自編程序

三者結合的方式將復雜繁瑣的結構有限元理論通過簡單直觀的方式展現出來，同時深層次的學習有限元理論和商業軟件的內部實現原理。

有限元的理論發展了幾十年已經相當成熟，商用有限元軟件同樣也是采用這些成熟的有限元理論，只是在實際應用過程中，商用CAE軟件在傳統的理論基礎上會做相應的修正以解決工程中遇到的不同問題，且各家軟件的修正方法都不一樣，每個主流商用軟件手冊中都會注明各個單元的理論采用了哪種理論公式，但都只是提一下用什么方法修正，很多沒有具體的實現公式。商用軟件對外就是一個黑盒子，除了開發人員，使用人員只能在黑盒子外猜測內部實現方式。

                                            

一方面我們查閱各個主流商用軟件的理論手冊并通過進行大量的資料查閱猜測內部修正方法，另一方面我們自己編程實現結構有限元求解器，通過自研求解器和商軟的結果比較來驗證我們的猜測，如同管中窺豹一般來研究的修正方法，從而猜測商用有限元軟件的內部計算方法。我們關注CAE中的結構有限元，所以主要選擇了商用結構有限元軟件中文檔相對較完備的Abaqus來研究內部實現方式，同時對某些問題也會涉及其它的Nastran/Ansys等商軟。為了理解方便有很多問題在數學上其實并不嚴謹，同時由于水平有限可能有許多的理論錯誤，歡迎交流討論，也期待有更多的合作機會。

通用結構有限元軟件iSolver介紹視頻：

        

                                          http://www.yqgqt.org.cn/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 視頻講解和操作驗證演示

模態疊加法的算例演示可看下方視頻

http://www.yqgqt.org.cn/college/video/c12884 20理論系列文章30-諧響應分析原理

 

==總結==

本文首先分別介紹了直接法和模態疊加法，并對兩種方法從適用范圍、計算效率和求解精度三個方面進行了對比分析，最后通過算例比較iSolver和Abaqus求解結果，對iSolver進行驗證。

以往的系列文章：

1.6.1 ========第一階段========

第一篇：S4殼單元剛度矩陣研究。

http://www.yqgqt.org.cn/content/post/338859

第二篇：S4殼單元質量矩陣研究。

http://www.yqgqt.org.cn/content/post/343905

第三篇：S4殼單元的剪切自鎖和沙漏控制。

http://www.yqgqt.org.cn/content/post/350865

第四篇：非線性問題的求解。

http://www.yqgqt.org.cn/content/post/360565

第五篇：單元正確性驗證。

http://www.yqgqt.org.cn/content/post/373743

第六篇：General梁單元的剛度矩陣。

http://www.yqgqt.org.cn/content/post/403932

第七篇：C3D8六面體單元的剛度矩陣。

http://www.yqgqt.org.cn/content/post/430177

第八篇：UMAT用戶子程序開發步驟。

http://www.yqgqt.org.cn/content/post/432848

第九篇：編寫線性UMAT Step By Step。

http://www.yqgqt.org.cn/content/post/440874

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

http://www.yqgqt.org.cn/content/post/531029

1.6.2 ========第二階段========

第十一篇：自主CAE開發實戰經驗第一階段總結。

http://www.yqgqt.org.cn/content/post/532475

第十二篇：幾何梁單元的剛度矩陣。

http://www.yqgqt.org.cn/content/post/534362

第十三篇：顯式和隱式的區別。

http://www.yqgqt.org.cn/content/post/537154

第十四篇：殼的應力方向。

http://www.yqgqt.org.cn/content/post/1189260

第十五篇：殼的剪切應力。

http://www.yqgqt.org.cn/content/post/1191641

第十六篇：Part、Instance與Assembly。

http://www.yqgqt.org.cn/content/post/1195061

第十七篇：幾何非線性的物理含義。

http://www.yqgqt.org.cn/content/post/1198459

第十八篇：幾何非線性的應變。

http://www.yqgqt.org.cn/content/post/1201375

第十九篇：Abaqus幾何非線性的設置和后臺。

http://www.yqgqt.org.cn/content/post/1203064

第二十篇：UEL用戶子程序開發步驟。

http://www.yqgqt.org.cn/content/post/1204261

1.6.3 ========第三階段========

第二十一篇：自主CAE開發實戰經驗第二階段總結。

http://www.yqgqt.org.cn/content/post/1204970

第二十二篇：幾何非線性的剛度矩陣求解。

http://www.yqgqt.org.cn/content/post/1254435

第二十三篇：編寫簡單面內拉伸問題UEL Step By  Step

http://www.yqgqt.org.cn/content/post/1256835

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

http://www.yqgqt.org.cn/content/post/1261165

第二十五篇：顯式分析的穩定時間增量。

http://www.yqgqt.org.cn/content/post/1263601

第二十六篇：編寫線性VUMAT Step By Step。

http://www.yqgqt.org.cn/content/post/1266640

第二十七篇：Abaqus內部計算和顯示的應變。

http://www.yqgqt.org.cn/content/post/1273788

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

http://www.yqgqt.org.cn/content/post/1282956

第二十九篇：幾何非線性的T.L.和U.L.轉換關系

http://www.yqgqt.org.cn/content/post/1286065

第三十篇：諧響應分析原理

http://www.yqgqt.org.cn/content/post/1290151

第三十一篇：自主CAE開發實戰經驗第三階段總結

http://www.yqgqt.org.cn/content/post/1294824