（原創，歡迎轉載，轉載請說明出處）

1 概述

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

（1） 基礎理論

（2） 商軟操作

（3） 自編程序

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

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

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

百度網盤鏈接: https://pan.baidu.com/s/10d6jHdZ01SBY2JxiS6bffw 提取碼: 6fdf

2 聲學分析

Abaqus雖然是結構計算軟件，但如今隨著系統的復雜，結構與其它物理場的耦合也越來越頻繁，因此，Abaqus也將自己的計算能力往其它物理場擴展，其中，熱學、聲學、流體和電磁學等都是Abaqus多物理場的擴展的典型應用。熱學主要計算熱膨脹時結構應力的變化和熱傳導等現象，本身和結構關系密切，所以Nastran、Abaqus等結構軟件都會包括熱分析功能；流體由于達索已經收購了業界兩款最好的基于最近大熱的LBM粒子算法的商業軟件PowerFlow和Xflow，可能處于產品線的發展考慮，把基于有限元方法的Abaqus/CFD模塊直接停止開發了，但保留了Abaqus中的CEL模塊，專門處理流固耦合問題；而聲學模塊主要處理聲腔內的模態分析和聲傳播問題；電磁模塊至今不涉及，在此不討論了。這幾個多物理場的擴展每一項都需要單獨的其它物理場的專業知識和結構專業的融合，也許大家都和我們一樣只熟悉結構的算法，對陌生的其它專業總覺得非常難，但其實靜下心來系統學習一下，其它專業的內容并不是想象的那么困難，結構有限元學好了，其它專業都是一通百通的。譬如下面的聲學有限元理論公式我們就可以完全按照結構有限元的流程推導。我們在實際研發CAE軟件過程中發現難的不是這些多物理場公式等的編寫，難的是如何在現有的結構CAE軟件基礎上在前后處理和求解器層面分別加入其它分析能力，同時又保持原有框架的可擴展性和易維護性，當然，拋開多物理場，就算僅研發結構CAE軟件如何在維持現有框架結構基礎上在軟件層面添加新的功能也是實際中我們需要花大力氣研究的一個關鍵點，這也是我們理解的在研發上理論上的有限元和工程上的有限元分析的最大區別。

現階段我們只討論聲學物理場，由于聲學涉及的問題較多，我們僅限于和實際工作中用到的和結構強相關的部分內容，水平有限，同時對聲學理論等理解的也不夠透徹，可能有很多問題，大家有疑問也歡迎多討論。我們分幾篇文章來說明，本篇是聲學分析（3）-濕模態。

==第44篇：濕模態分析==

2.1.1 濕模態概念

一個結構體在真空中的振動和在周圍流體介質中的振動有時差異巨大，必須分開考慮，此時有干濕模態之分，以下干濕模態概念選之：結構振動問題中的干模態和濕模態 - 知乎 (zhihu.com)

模態分析是結構動力學問題中最為普遍的一個問題，大家對這個概念也非常熟悉，它是結構的固有特性，在不考慮阻尼的影響下，只與結構的剛度和質量有直接關系。

本文所涉及的問題是有關干、濕模態的計算，屬于模態問題中經常遇到的問題。通常我們所計算的模態其實是干模態，主要是由于結構放置在空氣中而空氣這一流體對結構的影響幾乎可以忽略不計，所以基本上可以默認為就是干模態。

但是嚴格意義上講，受流體附帶阻尼及剛度的影響，這類的模態仍然是濕模態。所謂濕模態是考慮流體對結構的作用，也就是在通用的振動方程中加入了有關流體的附加質量及剛度矩陣（Kx、Mx），這塊相互作用對結構的模態存在一定的影響，尤其是涉及諸如油、水等液體的作用。舉例來講，司馬光砸缸的那口缸，在加水和不加水的情況下，砸缸的聲音肯定是不一樣的，一個明顯清脆寬廣，一個就顯得沉悶。再比如我們車載的油箱結構，在加滿汽油和未加滿汽油的情況下，兩者的模態肯定是不一樣的。另外還有諸如潛艇、船舶、油罐車等結構。所以，對于受到液體作用顯著的結構，我們在計算的時候需要研究其濕模態。

2.1.2 濕模態分析方法

當結構體內部或者外部和流體接觸時，結構受到流體的壓力作用而發生變形，而變形又會反過來對流體邊界有影響，形成一個流體壓力->結構變形->流體壓力變化的典型流固耦合系統。流固耦合的體系相當龐大，我們只討論流體載荷對一個振動結構的耦合作用。

有兩種可能的情況:

（1）一種可能是結構變形導致流體壓縮或膨脹，影響流體對結構的壓力，典型例子是船體或者浮動平臺越來越長后，會影響船體周圍的波浪形式。

一般認為當無量綱

不是一個小量時，會出現這個情況。其中l是結構體的特征長度，是流體中載荷的波長，也就是結構體的特征長度和流體中載荷的波長可比擬。個人理解工程上可以采用兩種方式求解：

a. 流體用邊界元形式，但邊界考慮了流固耦合，屬于氣彈性（aero-elasticity）和水彈性（hydro-elasticity）的范疇，一直覺得照原理來說應該比下面的虛擬質量法來算濕模態更準，和做水彈性的專家聊過，實際上水彈性都還是基于干結構的模態基展開來計算結構動響應，所以輸入的是干結構的模態分析后處理文件，程序內部應該沒有計算濕模態，水彈性和氣彈性還是只計算流體載荷，不是算結構振動的。

b. 流體用聲學有限元形式，然后用Nastran的glue或者Abaqus的Tie和結構有限元綁定，此時流體的壓力會作為外載荷作用到結構上，耦合面結構的加速度也會作為聲學有限元的邊界，即前面聲學有限元章節中的Inward volume acceleration。這種方法從原理上來說也比下面的虛擬質量方法更加正確，但建模相對復雜，譬如如果是水域的話，見過外面的流體域大概取到結構體的6倍以上共振頻率解才能穩定，這樣導致網格量大大上升，所以在工程上反而不夠實用。

（2）另一種結構的變形不會導致流體載荷的變化，譬如海浪中一塊小木板的變形不會對波浪載荷有影響，即單向耦合。也就是當kl<<1，此時結構體相對于流體載荷中的波長可忽略。此時流體假定不可壓、無粘，可以用勢流理論的邊界元表示流體，然后流體的作用等效為虛擬質量加到結構質量上，而流體對結構的剛度陣的影響忽略，至于和上面的氣彈性或者水彈性的邊界元+結構有限元的方法的不同，有可能僅在耦合邊界的處理上，這邊只有單向的流體對結構的壓力傳遞，而上面的是雙向的，僅猜測，不一定對。這種方法雖然理論上不如上面的流固雙向耦合精確，但流體只要建一層邊界就行，所以工程上相對更實用。在結構商軟中，Nastran的MFluid和ACMODL分別模擬虛擬質量法和流固耦合模態分析，其中MFluid虛質量分析功能是由戴姆勒.奔馳公司提供資金支持開發的一個功能模塊，目的是計算油箱等流固耦合部件在高頻響應。而Ansys、Abaqus暫時沒有虛擬質量方法，都是直接采用基于聲學有限元的流固耦合來求濕模態。

本章只介紹基于虛擬質量的濕模態計算。

2.2 基于虛擬質量的濕模態計算理論

無粘狀態下，不考慮剪切應力，所以只有一個正壓力未知量，流體正壓力中包括兩部分：一部分是靜壓，不隨時間變化，靜壓對結構會產生一個初始應力和位移，結構強度足夠的話，這種小變形導致的振動影響非常小，另一部分是動壓，只有動壓才會對振動有大的影響，無粘狀態下這種在穩定的流體場基礎上疊加的動壓就是聲壓。既然是聲壓，那么就可以用前面統一的Helmholtz方程：

在結構振動過程中，流體域中我們假定沒有聲源，這樣依然能直接使用上一章的邊界元理論得到的近似方程。

上式是聲壓和速度的關系，需要轉換為力和加速度的關系。做以下變換

? Pi是流體對邊界元的壓力，朝結構內部法向，但速度邊界還是超流體法向，所以需要將Pi變號，使得力和速度方向一致。

? 壓強需要轉換到節點力，轉換方式和結構有限元的面載荷轉到集中力載荷一致，設為：

? 速度轉換為加速度，此時需要考慮方向

最終得到：

（2）因為P可以由直接獲得，所以P和Vn都是同類型單元更合理，但如果都采用線性單元，那么B的計算量相對要多很多，因此還是都是常單元效率更高些，此時：

其實和上面單元相比只是C的形函數從積分內移動到了積分外，沒有本質的區別，結果也非常相近，無論哪種都不足以解釋后面我們將談到的和Nastran的精度誤差。

最后來說說Green函數G。因為一般來說聲音是可壓縮的，G是復數，M正常來說也會是一個復數形式，這樣就會得到復數形式的模態，但前面說過，在這兒，因為假設不可壓，那么相當于聲速無窮大，也就是k->0，這也和k*l<<1的假設是一致的。

則

也就是G最終簡化為實數，使得M也是實數，相當于原始Helmholtz方程簡化為了Laplace方程：

2.3 基于iSolver流程的濕模態的工程實現

對標Nastran的MFluid功能，基于iSolver的濕模態工程實現上還做了以下幾個關鍵步驟：

（1）自由面處理：Nastran的MFluid可以處理包括流體自由面（Free Surface）的情況，這也是非常常見的一個情況，一個典型問題是譬如單面法向向上的平板，自由面在板的上方，當自由面越來越接近平板時，理論上應該會得到真空中模態結果。在iSolver中，我們采用了鏡像濕面元的形式求解流體自由面的情況，同時在iSolver界面中可設置有無自由面及自由面高度。

（2）界面流程：Nastran的MFluid沒有界面，只能在bdf關鍵詞中添加，譬如，

MFLUID,2,0,10.0,1025.0,3,,N,N

ELIST,3,110001,THRU,116000

表示海水密度是1025，10m高的自由面，對110001到116000號單元中位于自由面下方的單元施加虛擬流體質量。MFLUID的關鍵詞很容易編輯，但ELIST單元號對不連續的濕表面單元手工很難列舉。在算法層面自主軟件很難超越商軟，只能在這些前后處理的個別小點上做些更加人性化的設置，在iSolver中可以直接在三維圖形上選取施加濕面元的單元，同時，不要求濕單元連續編號，后續我們還將把iSolver前處理產生的虛擬質量關鍵詞和濕單元信息自動輸出到bdf文件中，為那些想調用Nastran來求濕模態的用戶解決點Patran無法實現的功能。

（3）精度的考核：iSolver實現后，我們發現和Nastran的MFluid的部分結果存在部分差異，基頻普遍略大，而高階頻率差異較小，具體如下面示例展示，研究了很久，對上面我們貼的算法做了大量的調整，依然得不到和Nastran的MFluid差異非常小的濕模態結果，我們把和Nastran的對比結果貼在下方，還需做過這方面工作的高手指點。

2.4 矩形板濕模態模型驗證

2.4.1 算例介紹

單位采用SI制，模型為一個10X5大小的矩形板，厚度t=0.1,材料楊氏模量2.01e11,泊松比=0.31，密度7920。四周鉸支。

2.4.2 iSolver虛擬質量結果

iSolver中將矩形板劃分為0.5大小網格，單面加無限深水時，iSolver基頻結果：

? *網格大小影響不是很大，如果劃分成0.25，頻率基本不變

前10階如下

2.4.3 Nastran虛擬質量結果

在iSolver中將模型一鍵輸出bdf到Nastran，由于Nastran的MFluid關鍵詞中沒找到不設置自由面的方法，所以直接設了一個108的深度，對10X5米的板來說自由面影響基本可忽略，得到基頻6.57Hz：

2.4.4 Abaqus結構+聲學有限元濕模態結果

Abaqus用聲學有限單元模擬流體結果，水密度設置為1025，用Tie將板結構和聲學單元綁定，為排除水域尺寸的影響，參考相關文件，取12倍結構尺寸的水域，Abaqus中分析后三階頻率分別為5.9085，11.605，21.349

2.4.5 SAM、Nastran和Abaqus單側加水濕模態比較

將上述三個軟件得到的濕模態的結果匯總到一個表中，可發現iSolver和Nastran的MFluid相比基頻較小，但和Abaqus基于結構和聲學有限元的基頻更接近。2、3階頻率基本一致。我們一直沒找到單面加水的矩形板的模態理論值或者試驗值，雖然我們的基頻結果更接近Abaqus，但實際使用模態分析的時候，用戶更關心的是和Nastran的對標，只可惜一致沒找到Nastran的MFluid的后臺真實的理論修正方法，如果你恰好也做過基于虛擬質量的濕模態，可以嘗試一下這個算例結果，或者下載我們的軟件對比一下，看看是不是也是這種情況，有問題我們可以一起校核。

2.5 以往的系列文章

2.5.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

2.5.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

2.5.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

2.5.4 ========第四階段========

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

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

第三十二篇：諧響應分析算法

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

第三十三篇：線性瞬態動力學

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

第三十四篇：非線性瞬態分析

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

第三十五篇：接觸求解算法

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

 第三十六篇：DLOAD用戶子程序開發步驟

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

 第三十七篇：梁單元差異（1）-理論基礎

https://jishulink.com/content/post/1872208

 第三十八篇：梁單元差異（2）-梁截面方向

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

 第三十九篇：梁單元差異（3）-剪力和彎矩

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

 第四十篇：梁單元差異（4）-形心、剪心和偏置

http://www.yqgqt.org.cn/post/1888000

1.5.5 ========第五階段========

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

http://www.yqgqt.org.cn/post/1905963

 第四十二篇：聲學分析(1)-有限元

http://www.yqgqt.org.cn/post/1912491

 第四十三篇：聲學分析(2)-邊界元

http://www.yqgqt.org.cn/post/1923936