莊 超， 劉漢光， 康凱旋， 蘇俊收

(1.徐工集團 江蘇徐州工程機械研究院，江蘇 徐州 221004；2.高端工程機械智能制造國家重點實驗室，江蘇 徐州 221004)

聲學(xué)邊界元法[1-2]，作為聲學(xué)仿真的一種數(shù)值計算方法，相比于聲學(xué)有限元法的技術(shù)特點在于該方法只需計算域的表面網(wǎng)格，這一特點對于計算大規(guī)模聲學(xué)問題相對于聲學(xué)有限元法具有無可比擬的優(yōu)勢，比如計算船體的輻射噪聲，有限元法須要劃分計算域的體網(wǎng)格，其計算量非常大，而邊界元法只需獲取結(jié)構(gòu)表面網(wǎng)格即可實現(xiàn)輻射噪聲計算，計算量大幅度下降，同時快速多級邊界元法可以提高計算求解速度。針對工程廣泛應(yīng)用的圓柱殼結(jié)構(gòu)，張俊杰等[3]探討了邊界元軟件 Sysnoise兩種加載方式(即基于單元(Element)加載和基于節(jié)點(Node)加載)對聲輻射影響，認(rèn)為基于單元加載方式得到的結(jié)果要好于基于節(jié)點加載方式，并且基于單元加載方式網(wǎng)格加密時收斂到解析解，而基于節(jié)點加載方式則不一定。目前，傳統(tǒng)的聲學(xué)邊界元仿真過程包含建立分析對象的CAD模型，劃分聲學(xué)表面網(wǎng)格，施加邊界條件并實施聲學(xué)邊界元數(shù)值計算等。然而，對于以下兩種情況，傳統(tǒng)聲學(xué)邊界元仿真呈現(xiàn)明顯地局限性：①分析對象無現(xiàn)成CAD模型，須要測繪或逆向獲取，無論是測繪建模還是通過逆向獲取點云數(shù)據(jù)再蒙皮建立CAD模型，都將花費大量人力和時間，尤其當(dāng)分析對象包含大量復(fù)雜曲面結(jié)構(gòu)；②分析對象尺寸規(guī)模較大并包含復(fù)雜幾何結(jié)構(gòu)，聲學(xué)表面網(wǎng)格劃分將帶來較大的工作量和技術(shù)難度。

逆向工程[4]，借助三維掃描和數(shù)據(jù)采集處理技術(shù)，將實物點云數(shù)據(jù)，通過三維模型重構(gòu)獲取三維模型，廣泛應(yīng)用于汽車、航空、機械工業(yè)產(chǎn)品開發(fā)設(shè)計，以及生物醫(yī)學(xué)領(lǐng)域的器官三維模型重建等各領(lǐng)域。尤其針對具有復(fù)雜曲面的實物模型，相比于傳統(tǒng)的測繪方法，具有天然的效率和優(yōu)勢?；邳c云數(shù)據(jù)，技術(shù)人員將逆向工程與數(shù)值分析相結(jié)合，前者獲取的三維網(wǎng)格模型作為后者的計算網(wǎng)格數(shù)據(jù)，實現(xiàn)兩者的無縫銜接。例如，在生物醫(yī)學(xué)應(yīng)用領(lǐng)域，趙巖等[5-6]，將CT掃描圖像在Mimics等逆向工程軟件進(jìn)行三維模型重建，并應(yīng)用Remesh優(yōu)化三維面網(wǎng)格以獲取規(guī)則化的三角面片，然后導(dǎo)入有限元軟件中劃分體網(wǎng)格，賦以材料屬性并添加邊界條件，實現(xiàn)有限元數(shù)值分析。

考慮到傳統(tǒng)聲學(xué)邊界元的上述缺陷，以及逆向工程與數(shù)值分析相結(jié)合的可行性，本文首次將逆向工程和聲學(xué)邊界元法相結(jié)合，將前者獲取的點云數(shù)據(jù)轉(zhuǎn)換為可用于計算的聲學(xué)邊界元表面網(wǎng)格，實現(xiàn)兩者的無縫銜接；同時，給出了小型掃路機廂體聲傳播問題的具體算例，并與傳統(tǒng)聲學(xué)邊界元仿真結(jié)果對比，得到了較高的吻合度，證明了本文方法的有效性與可行性。

1 點云及三角網(wǎng)格獲取

利用三維激光掃描儀獲取小型掃路機廂體部件(見圖1)的點云數(shù)據(jù)，在Geomagic逆向處理軟件中，同實物比對以獲取連續(xù)的、無離散數(shù)據(jù)的點云為原則；通過直接觀察法、弦高差法，采用自動修復(fù)、去除特征、清除、填充孔等操作，識別并去除點云的噪聲點和冗余數(shù)據(jù)[7-8]。

圖1 小型掃路機Fig.1 Small sweeper

考慮計算量基于均勻采樣或者弦偏差采樣實現(xiàn)點云精簡，采用曲線、曲面插值補充法對缺失的幾何特征完成數(shù)據(jù)補齊；基于曲率均勻法借助快速光順對其進(jìn)行濾波光順處理[9]，獲得最終的點云數(shù)據(jù)三角網(wǎng)格，為后續(xù)的數(shù)值計算提供最基本的網(wǎng)格數(shù)據(jù)，該三角網(wǎng)格是連續(xù)的且能夠表達(dá)實物幾何特征，如圖2所示。圖2中點云三角化模型為了保留實體幾何特征，曲面網(wǎng)格數(shù)量較大，整個模型范圍內(nèi)網(wǎng)格單元分布不均與、網(wǎng)格密度相差大，同時三角面片大小與形狀不均勻，存在大量狹長型、內(nèi)角多為小銳角和大鈍角的三角面片，以及階數(shù)大于6的頂點。

圖2 三角網(wǎng)格模型Fig.2 Triangular mesh model

2 三角網(wǎng)格的重網(wǎng)格化

由圖2可知，三角網(wǎng)格的質(zhì)量并不滿足插值與積分計算的需求，無法直接用于邊界元數(shù)值計算，須要進(jìn)行網(wǎng)格精簡和單元質(zhì)量優(yōu)化處理，以獲取既符合幾何表達(dá)又滿足數(shù)值計算要求的三角網(wǎng)格單元模型。沈建國等[10]提出面向有限元分析的三角網(wǎng)格迭代優(yōu)化，通過網(wǎng)格細(xì)分、簡化和規(guī)則化處理，在給定的誤差范圍內(nèi)有效地改善三角形形態(tài)，提高了網(wǎng)格質(zhì)量，為有限元單元的生成提供合適的網(wǎng)格模型。

為實現(xiàn)三角網(wǎng)格的重網(wǎng)格化[11]；首先，根據(jù)數(shù)值計算的類型、精度、模型整體規(guī)模等，給定目標(biāo)數(shù)量或單元尺寸簡化模型；然后，通過重網(wǎng)格化實現(xiàn)單元質(zhì)量優(yōu)化，包括幾何和拓?fù)鋬?yōu)化兩個方面，前者以三角面片形狀和尺寸均勻為目標(biāo)，后者以網(wǎng)格頂點具有相同的階數(shù)為目的。

三角網(wǎng)格的重網(wǎng)格化的基本思想是對于給定的初始網(wǎng)格，基于采樣均勻性、規(guī)則化、尺寸和單元形狀等準(zhǔn)則，重新生成網(wǎng)格在較好近似原網(wǎng)格的基礎(chǔ)上獲取較高的網(wǎng)格單元質(zhì)量。本文采用半邊數(shù)據(jù)結(jié)構(gòu)(Halfedge Data)以及相應(yīng)的庫函數(shù)組成的表示和管理多邊形網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)OpenMesh[12-13]，來存儲、管理以及遍歷模型的網(wǎng)格數(shù)據(jù)，并根據(jù)Botsch等[14]提出的各向同性重網(wǎng)格化實現(xiàn)單元質(zhì)量優(yōu)化。其算法實現(xiàn)過程為：①分裂所有長度大于目標(biāo)邊長的邊；②消除所有長度短于目標(biāo)邊長的邊；③翻轉(zhuǎn)邊以獲取最優(yōu)的平均頂點階數(shù)；④在切平面中移動頂點至重心位置。

將上述步驟迭代數(shù)次，獲得目標(biāo)邊長、頂點階數(shù)為6的三角面片網(wǎng)格。然后，采用基于面積的切向光滑處理，將頂點移動至重心加權(quán)的形心位置，加權(quán)系數(shù)計算如式(1)，最終得到理想的三角網(wǎng)格模型，具有網(wǎng)格單元密度與尺寸均勻、形狀均接近等邊三角形、絕大多數(shù)頂點階數(shù)為6等特點。原廂體三角面片網(wǎng)格的重網(wǎng)格化后的網(wǎng)格模型，作為聲學(xué)邊界元分析的單元模型，如圖3所示。

(1)

圖3 三角網(wǎng)格的重網(wǎng)格化Fig.3 Remeshing of triangular mesh

3 編寫數(shù)值計算文件

經(jīng)重網(wǎng)格化獲取的廂體三角網(wǎng)格單元模型，仍然是一種表達(dá)幾何的模型，不能直接用于計算。表達(dá)幾何的三維網(wǎng)格數(shù)據(jù)包含網(wǎng)格頂點的三維坐標(biāo)、編號以及三角面片的拓?fù)溥B接關(guān)系等。

本文通過編程，將上述數(shù)據(jù)寫成通用的、包含數(shù)值計算信息的INP格式數(shù)據(jù)文件，包含節(jié)點、單元、單元類型等基本數(shù)據(jù)。圖4給出了表達(dá)幾何及數(shù)值計算的數(shù)據(jù)文件格式。通過編寫數(shù)值計算數(shù)據(jù)文件，我們獲取了可以導(dǎo)入聲學(xué)邊界元計算程序的單元模型數(shù)據(jù)。

圖4 幾何與計算數(shù)據(jù)文件格式Fig.4 File formats for geometric and computational data

4 聲學(xué)邊界元理論

聲學(xué)問題控制方程Helmholtz微分方程

2p+k2p=0

(2)

式中:k為波數(shù)，k=ω/c?？刂品匠贪珼irichlet,Neumann以及Robin等三類邊界條件。對于穩(wěn)態(tài)聲場外問題，還須要滿足遠(yuǎn)場Sommerfeld輻射邊界條件。

(3)

借助Helmholtz方程基本解u*與格林公式[15]，

(4)

將式(2)和式(4)代入式(5)，得到間接邊界元網(wǎng)格上的單層勢σ與雙層勢μ表示的Helmholtz間接邊界積分方程

(6)

式中:單層勢σ代表間接邊界元兩側(cè)的速度差即

(7)

雙層勢μ代表間接邊界元兩個的聲壓差即

μ=p+-p-

(8)

將式(6)中的源點P移到求解域的邊界，獲得邊界上的聲壓及其外法向質(zhì)點振速之間的邊界積分方程。

(9)

(10)

為了數(shù)值求解邊界積分方程，須要將聲學(xué)域的邊界劃分成有限個單元，這里用二維三角形面單元。將邊界離散成Ne個單元，則邊界積分式(9)變成

(11)

式中：Γe為單元e的邊界。對于邊界上所有N個節(jié)點，可以得到N個方程，用矩陣形式表示為

(12)

式中：B,C和D為系數(shù)矩陣;σ為單層勢列向量與μ為雙層勢列向量;f和g為外界激勵向量。求解上述方程得到邊界單元各節(jié)點的單層勢σ與雙層勢μ，然后聲場區(qū)域內(nèi)任一點聲壓為

p={Aσ}T{σ}+{Aμ}T{μ}

本文用于離散計算域的二維三角形面單元，即是點云數(shù)據(jù)經(jīng)重網(wǎng)格化獲取的分布均勻、形狀規(guī)則的二維三角形網(wǎng)格單元。

5 算 例

為了驗證本文所述方法的可行性，將其用于模擬小掃路機廂體內(nèi)放置點聲源的聲傳播問題，并與傳統(tǒng)聲學(xué)邊界元方法的結(jié)果進(jìn)行對比。計算流體域為空氣，邊界條件為位于廂體內(nèi)某一位置的點聲源，其頻譜特性如圖5所示。頻率范圍16～1 800 Hz，在廂體左、右側(cè)7.5 m某位置處放置測點，拾取聲壓響應(yīng)。

首先，建立基于點云數(shù)據(jù)的聲學(xué)邊界元計算模型，將上述重網(wǎng)格化后的小掃路機廂體三角形網(wǎng)格模型經(jīng)數(shù)據(jù)格式轉(zhuǎn)換得到的INP數(shù)據(jù)模型，并用于聲學(xué)邊界元數(shù)值計算。

圖5 點聲源頻譜曲線Fig.5 Spectrum curve of point sound source

然后，建立傳統(tǒng)聲學(xué)邊界元分析模型。①將逆向獲取的點云數(shù)據(jù)導(dǎo)入CAD軟件Pro Enginner，通過表面蒙皮、光順等實現(xiàn)三維幾何重建以獲得廂體的三維CAD模型，如圖6所示。②利用網(wǎng)格劃分工具HyperWorks得到用以計算的表面網(wǎng)格，導(dǎo)出NASTRAN類型的Standard format格式文件，為了滿足計算精度需求，須要花費一定的時間，基于質(zhì)量評價準(zhǔn)則，通過單元質(zhì)量優(yōu)化來獲取整體單元質(zhì)量較好的網(wǎng)格單元模型，如圖7所示。③導(dǎo)入聲學(xué)邊界元程序LMS Virtual. Lab Acoustics，定義為聲學(xué)網(wǎng)格，定義聲學(xué)材料為空氣，將聲學(xué)材料定義到聲學(xué)網(wǎng)格上，定義場點，定義點聲源并導(dǎo)入圖5所示的頻譜曲線，先后計算聲場分布、場點響應(yīng)以及場點上的聲壓頻率響應(yīng)函數(shù)。

用于仿真的計算機處理器為Intel 3.60 GHz，內(nèi)存16 GB。表1給出兩種聲學(xué)仿真方法各步驟及時間花費，實物逆向時間相同，同時由于目標(biāo)計算網(wǎng)格數(shù)量差不多，所以分析計算時間基本一致；主要時間差異來源于從點云數(shù)據(jù)獲取計算網(wǎng)格的時間，傳統(tǒng)聲學(xué)邊界元須建立CAD模型并劃分網(wǎng)格，共須240 min，而基于點云數(shù)據(jù)的聲學(xué)邊界元方法通過程序?qū)崿F(xiàn)點云三角網(wǎng)格到計算網(wǎng)格的轉(zhuǎn)換，僅須15 min。

表1 基于點云數(shù)據(jù)的與傳統(tǒng)聲學(xué)BEM仿真時間Tab.1 Simulation time of based on point cloud data and traditional acoustic BEM min

圖6 CAD模型Fig.6 CAD model

圖7 邊界元表面網(wǎng)格Fig.7 Surface mesh for boundary element

圖8和圖9，分別給出了上述應(yīng)用基于點云數(shù)據(jù)的聲學(xué)邊界元與傳統(tǒng)聲學(xué)邊界元仿真計算結(jié)果，即廂體左、右側(cè)測點的聲壓響應(yīng)頻譜曲線。由圖可知，兩種方法仿真結(jié)果吻合度較高，這也證明了本文提出的基于點云數(shù)據(jù)的聲學(xué)邊界元分析方法的有效性。

圖8 左側(cè)點聲壓響應(yīng)頻譜Fig.8 Sound pressure response spectrum of the left side

圖9 右側(cè)點聲壓響應(yīng)頻譜Fig.9 Sound pressure response spectrum of the right side

6 結(jié) 論

本文提出的基于點云數(shù)據(jù)的聲學(xué)邊界元分析方法，將逆向工程與聲學(xué)邊界元數(shù)值計算相結(jié)合，把掃描獲取的不規(guī)則的幾何網(wǎng)格模型優(yōu)化處理，得到形狀規(guī)則的網(wǎng)格模型以滿足數(shù)值計算的需求，然后轉(zhuǎn)換為包含單元類型和拓?fù)潢P(guān)系的計算數(shù)據(jù)模型以用于聲學(xué)邊界元分析，實現(xiàn)了CAD幾何模型與CAE仿真分析的無縫銜接。

于是，對于實物模型的聲學(xué)分析，在通過逆向掃描獲取點云數(shù)據(jù)之后，不再需要根據(jù)點云建立CAD幾何模型并劃分聲學(xué)邊界元表面網(wǎng)格的傳統(tǒng)路徑，只須將點云數(shù)據(jù)進(jìn)行優(yōu)化處理并添加單元數(shù)值計算信息，即可得到適用于仿真計算的網(wǎng)格單元數(shù)據(jù)，該方法對于具有復(fù)雜曲面的實物模型尤其能體現(xiàn)出巨大優(yōu)勢。文中給出的小型掃路機廂體聲傳播問題具體算例，與傳統(tǒng)聲學(xué)邊界元仿真對比，得到了較高的吻合度，證明了該方法的可行性與有效性。