有限元理論基礎及Abaqus內部實現方式研究系列43：聲學分析(2)-邊界元

SnowWave02

2023年11月18日 13:00

瀏覽：3812 評論：2 收藏：6

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

1 概述

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

（1）基礎理論

（2）商軟操作

（3）自編程序

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

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

有限元理論基礎及Abaqus內部實現方式研究系列43：聲學分析(2)-邊界元的圖1

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

通用結構有限元CAE軟件iSolver操作視頻：

有限元理論基礎及Abaqus內部實現方式研究系列43：聲學分析(2)-邊界元的圖2

http://www.yqgqt.org.cn/college/video/c12884

==第43篇：聲學分析==

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

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

由于聲腔區域很多都是均勻的，這個區域很容易只要表達出邊界情況就能確定解，所以邊界元在聲學中經常用到，iSolver僅僅聚焦有限元，對于聲學有限元，我們團隊自己將結構有限元改造，加入了基于有限元的聲學分析模塊，但對于聲學邊界元，我們沒有做實際的開發，也不打算在邊界元領域深耕，只是希望和其它團隊合作在iSolver上二次開發，當然，很多人連自主CAE軟件直接操作都懶得試一下，更不要說愿意投入大把時間基于自主CAE軟件二次開發，但只要有人有意愿在iSolver平臺上二次開發，我們都會鼎力支持，現在的iSolver平臺已完全支持基于python和C++的二次開發。另一個團隊在iSolver前后處理平臺上二次開發完成聲學邊界元模塊，因為涉及界面的邊界元整體流程和底層的數據接口、聲振耦合傳遞等，所以我們對聲學邊界元略有了解，但比起真正做邊界元的還是差很多，僅記錄下我們理解的邊界元公式和開發難點。

1.1 聲學邊界元方程推導

聲學邊界元的標準推導可以直接看理論的書，我們只懂有限元，所以照著邊界元的理論書盡量按上一篇有限元的強方程-->弱方程-->有限元近似的流程重新簡要推導了一下。

1.1.1 強方程

聲學邊界元方程的推導最主要的是如何把一個體內域V的計算直接轉為只要知道邊界上的物理量就行。

還是從Helmholtz聲學波動方程出發。

這個方程有個特解，就是當一個源強為1的點源位于r0點時：

1.1.2 積分方程

一般，邊界元解決的工程問題都是結構表面振動或者流體表面脈動壓力導致聲音在外空間均勻介質傳播的問題，譬如汽車蓋在空氣中由于發動機振動導致的噪聲、潛艇在無限水域的輻射噪聲。此時聲源只在表面，輻射區域沒有聲源，對V內的任意一點。

其中 ${\mathrm{r}}_{\alpha }$ 為輻射體域內任意一個場點位置。

再假設這個輻射面是封閉的，這類問題也稱為直接邊界元問題。設這個輻射域由兩個面組成，一個是結構表面Sa，另一個應該是無窮大的球面SR2。

因為聲源只在表面，所以（2）式變為：

其中表示在r處有源強為1的聲源，有限元理論基礎及Abaqus內部實現方式研究系列43：聲學分析(2)-邊界元的圖20 還是為輻射體域內任意一個場點位置。用一個小球SR1包圍r點的聲源。在Sa（去掉SR1相交的部分，當SR1無窮小時可忽略）、SR1（去掉和Sa相交的部分）和SR2包圍的輻射域Va內。

如果r位于角上，則需要修正，否則結果和商軟將差異較大。

最終轉換為只有邊界上的積分：

粗看這個方程神奇的把k平方去掉了，但實際只不過轉換到了G函數中。

1.1.3 邊界元近似

上述方程個人覺得不能稱為有限元中的弱方程，因為不是像有限元一樣強方程乘以權重函數來形成積分形式的，但這個積分形式的離散可以和有限元類似，最終目點還是但滿足以下兩點：

嚴格的系統域內積分可以表達為另一種形式的函數，這些函數的未知量只有這些有限點上的物理量，這樣可以采用計算機求解。
當這些有限點的數目足夠密時，嚴格的系統域內積分就可以近似為另一種形式的函數了。

將系統邊界Sa劃分為有限的邊界單元，可以認為原有邊界域就和有限的邊界單元的面積累加起來近似。

在實際編程中，如果不是做純有限元，而是結構或者流體導致的外場問題，那么法向盡量取成從結構或者流體指向聲腔，因為一般結構或者流體體單元時，體上的面默認都是指向體的外面的。

我們還是以符號n為指向流體域的方向，得到：

結構或者流體導致的外場問題，一般都已知法向速度Vn，求聲壓P（r）。劃分為網格后，設網格點上的坐標為rn，只要將網格點rn對應的聲壓求出來就行。

上述還需解決單元內未知量P(ra)，和有限元類似，對一個單元內的聲壓P(ra)采用邊界元節點表示，等參形式下，任意點的聲壓P可以表示網格節點聲壓Pi的線性組合。

這樣，單元內的未知量就直接轉換為了節點的聲壓未知量組合了。

將r=rb，b=1…N（N為邊界上節點總數），遍歷，得到N個線性方程組，且包括N個P(rb)未知量。

上式直接求解可得Pi。

1.2 基于iSolver結構流程的聲學邊界元的程序實現的關鍵

邊界元的原理相對有限元來說并不復雜，但要做出一個能帶操作界面且精度效率可靠的工程軟件還是非常難的。合作團隊和平臺都投入了大量時間來共同實現基于iSolver的聲學邊界元的代碼開發，我們理解邊界元分析存在以下幾個關鍵技術：

1.2.1 大模型的存儲和計算問題

邊界元第一個問題是求解精度和內存消耗，其實拋開內存談效率和精度沒有意義，由上述最終的方程可知，待求的和所有單元的有關，導致得到的P(r)前系數陣是個滿陣，滿陣的存儲內存和計算都是個問題。在iSolver中實測1G內存大概只能存下1500個邊界單元，40G內存大概存下1w個邊界單元，基本上內存開銷和單元個數成平方比，對于計算效率，雖然邊界元涉及的矩陣均為非稀疏矩陣，但實際上我們發現現在成熟矩陣數學庫一般已經能兼容稀疏陣和非稀疏陣，如果工程使用的話直接調用成熟庫就行，當然，有能力也可自己編寫更快更適合自己實際情況的矩陣求解。現在商用聲學軟件都采用FMM快速多級子通過減少節點間的相互關系來減少內存和計算加速，缺點是也會損失些精度，而且算法復雜編程實現和維護相對也要花費更多時間。

1.2.2 聲學邊界元和結構有限元耦合分析問題

第二個問題是前后處理的流程：工程實際的噪聲源一般都是結構CAE或者流體CFD計算得到的邊界振速，對于聲學量會影響噪聲源的情況，需要雙向耦合，而對于影響較小的，可以直接將這些CAE軟件得到的邊界傳遞到邊界元的網格上，這些都涉及到兩套網格之間的匹配問題，網格不同需要通過形函數等插值保持力和力矩平衡。實際工程問題可以不單獨構建邊界元網格，而是直接在結構或者流體的CAE網格抽取耦合面的面網格作為邊界元網格，在Virtual.Lab等聲學邊界元商軟中很多工程就是這么做的，缺點是CAE網格一般比較密，導致抽取的邊界元網格相對較多，容易導致前面大模型導致的存儲和計算問題。就算是通過抽取表面的簡單方式來耦合，iSolver在實際前后處理編程中也為如何操作耦合考慮了很久，畢竟Abaqus等結構商軟沒有邊界元的現成流程可參考，同時Virtual.Lab也是直接導入CAE后處理開始的，最終iSolver前后處理的流程實現邊界元有3個關鍵步驟：

（1）直接從一個part模型上抽取表面形成邊界元

（2）類似穩態動力學的聲學邊界元的Step創建

（3）采用子模型的設置中來關聯邊界元和結構后處理數據。

1.2.3 融入到有限元分析流程中的邊界元求解過程

上面的是前后處理的關鍵流程，還有求解器的問題：對邊界元求解流程如何融入到有限元分析流程也是一個難點。上一章講過聲學有限元只要加入聲學單元，求出類似的剛度陣和非平衡力，就很容易嵌入到基于增量迭代法的有限元結構流程中，但邊界元實際上融入這個流程還有相當的困難，按照最終的方程來說，我們可以把P(r)前的系數陣當成剛度陣，然后也可以采用迭代法來求非平衡力，正常來說也是一次平衡，但邊界元基本不這么做，我們理解困難點在于全局剛度陣的組裝，有限元中由于節點只與跟它相連的單元節點影響，可以先求出單元剛度陣得到該單元內節點之間關系，然后只需要對結構單元一次遍歷就能把所有單元間節點關系累加得到全局剛度陣，也就是節點和所有其它節點之間的關系，但對邊界元來說，因為每個單元節點未知量都和其它所有單元相關，所以事實上沒有單元剛度陣一說，需要對單元兩次遍歷才能累加所有的節點間關系，也就無法用有限元稀疏陣的組裝方式。這看起來是個很小的改動，但當一個大型程序的流程固定后，看似一點小改動在代碼上也很難共存，當然，相信Nastran和Abaqus等商軟絕對有能力來改造，只是最終沒實現而已，不清楚具體原因，也許是應用領域比較狹窄和編程成本的平衡考慮吧，一種簡單方式是類似Nastran的Glue或者Abaqus的Tie一樣單獨開辟一個節點間關系的流程。邊界元合作團隊在實現邊界元過程中采用的是另一個單獨的exe來完成邊界元求解，求解輸入輸出數據文件按iSolver的以往的規范，這樣可以同步采用iSolver的前后處理，當然，就算融入到以結構為主的前后處理軟件中，邊界元單元、屬性、材料等的設置能否和結構一體化也是大家需要考慮的問題。不管如何，讓用戶感覺一體化是最重要的，類似Abaqus的standard和explicit求解器融入到Abaqus/CAE的前后處理一樣。

1.3 直接邊界元脈動球外場輻射模型校核

無論什么計算模塊集成到iSolver中，我們要求和自己的求解器iSolver一樣，都首先必須保證精度。我們僅以一個簡單算例考核一下集成到iSolver的聲學邊界元求解器精度。

1.3.1 模型介紹

以一個半徑為0.5m的脈動球為例來分別驗證邊界元計算結果。脈動球網格劃分網格如下所示，網格數為900，網格類型均為四邊形網格；聲學介質為水，密度為1000 Kg/m³，體積模量為2.18e+09 Pa。

1.3.2 邊界元分析步設置

iSolver中設置為邊界元分析步，且外場計算，頻率點設置為100-900Hz，取9個頻率點。

創建外場中圓柱場點，圓柱殼半徑5.0，長度10.0。

設置球表面法向振速為2m/s。

1.3.3 自主邊界元結果

采用邊界元求解后導入iSolver后處理查看云圖，選擇Magnitude幅值，iSolver中得到900Hz的聲壓幅值云圖如下，在脈動球表面聲壓幅值=2.92e6Pa：

1.3.4 和VirtualLab對比

將iSolver模型導出為bdf，然后導入聲學商業軟件Virtual.Lab.Acoustics中，并施加同樣的振動速度，Virtual.Lab中脈動球900Hz表面聲壓結果為2.93e6和2.97e6之間。

1.3.5 和理論對比

脈動球外場聲壓解析表達式為：

當r=a時，可得在脈動球表面在900Hz聲壓理論幅值= 2.92e6；

1.4 總結

本章我們首先推到了聲學邊界元的數值方程，然后合作團隊通過iSolver上的二次開發將邊界元流程融入到了以結構有限元為主的iSolver前后處理和求解器中，通過典型算例可得邊界元數值結果和邊界元商軟Virtual.Lab及理論基本一致。其實在iSolver中還實現了結構振動和聲學邊界元的耦合的操作流程，并且結構有限元結果和邊界元載荷的后臺數據實現了自動流轉，這才是真正實際工程的應用場景。最后，感謝我們的合作團隊的辛苦努力，自主CAE軟件平臺不好用，二次開發更是問題多多，理解我們的有限元的數據結構和接口是絕對需要花很多時間的，當然，數據結構都是我們自己編的，所以可以比商軟解釋的更清楚，而且起碼沒有原則上改不了的地方，這也是我們自主CAE的優勢之一。