三維聲腔模態的等幾何計算分析

薛亞強，靳國永，葉天貴，師康康

(哈爾濱工程大學 動力與能源工程學院，黑龍江 哈爾濱 150001)

三維聲腔的模態分析對壓縮機和消聲器等機械設備的聲學設計有重要的指導意義，對于比較規則的三維聲腔，如矩形、圓柱形和球形等，其特征解可以用解析法[1]求得。有限元法(finite element method，FEM)[2]被廣泛應用于復雜結構內部聲學的計算與分析。FEM中的網絡離散過程舍棄了幾何模型的精確信息，導致其精度下降，這種誤差主要來源于簡單的多項式形函數。Huttunen等[3]進行了單位分解有限元(partition of unity FEM，PUFEM)與超弱變分(ultra-weak variational formulation，UWVF)2種方法的對比研究，通過分析波在二維管道中的傳播，指出UWVF在高頻范圍具有更好的計算精度。廣義有限元法[4]通過引入特定的局部逼近形函數減弱了“污染效應”，其思想是在節點處引入廣義自由度，隨著自由度的增加，計算量逐漸增大。Grosu等[5]采用間斷富集法(discontinuous enrichment method, DEM)研究了二維區域的內部聲學問題，DEM將非協調富集函數與FEM中的多項式形函數之和作為近似解，與PUFEM和UVWF相比具有更高的精度與計算效率[6]。Desmet[7]以Trefftz法為基礎提出了波函數法(wave based method，WBM)，并進行了結構-聲耦合系統的穩態動態分析。WBM在整個求解區域內采用滿足控制方程的波函數描述場變量，與FEM相比，WBM不需要劃分單元，具有更高的收斂率，但保證其收斂的條件是計算域為凸形域[8]，這在一定程度上限制了WBM的應用范圍。He等[9]采用alpha有限元(α-FEM)研究了一維及二維聲學問題的無色散分析；α-FEM通過改變參數α的值來控制數值模型的剛度，由于波數、網格等多種因素的影響，難以選取最佳的α值。Missaoui等[10]提出了積分模態法，并研究了二維不規則聲腔的聲學特性；該方法將整個聲腔分割成多個子聲腔，通過增加子聲腔的數量來提高結果的準確性。此外，一些學者采用級數展開法，如Taylor級數[11]、Fourier級數[12]、Chebyshev級數[13]對聲壓進行計算，研究了封閉空間的聲學特性。

在結構設計中，計算機輔助設計(computer aided design，CAD)與計算機輔助工程(computer aided engineering，CAE)已經密不可分，但兩者采用不同的數學語言描述幾何模型與計算模型，這種不一致性割裂了產品的設計與分析，還會引起幾何誤差。為了克服這些障礙，Hughes等[14]提出了等幾何分析(isogeometric analysis，IGA)，將非均勻有理B樣條(non-uniform rational B-spline，NURBS)基函數作為FEM中的形函數，保證了幾何模型的精確性，從源頭消除了幾何誤差，實現了CAD和CAE數學描述的統一。近十幾年來，等幾何分析已成功應用于板殼振動[15]、流體力學[16]等多個領域。

結構形狀對域內聲場分布有明顯影響，有限元法缺乏精確的幾何模型，會導致求解失真，且網格劃分過程十分冗繁，網絡修改難度很大。NURBS基函數能靈活建立復雜結構的精確模型，IGA通過節點插入與基函數升階來增加幾何模型控制點的個數，避免了繁瑣的網格劃分過程，具有較高的潛在應用價值。目前等幾何分析在聲學領域[17-18]的應用較少，等幾何分析在聲學領域的研究有待深入。

本文介紹了NURBS的基本概念和三維內部聲場的等幾何分析，以截面為橢圓或橢圓環聲腔為例，計算了其波數和模態振型。將數值計算結果與解析解對比，驗證了等幾何分析法的正確性及收斂性，分析了幾何參數對橢圓聲腔固有特性的影響。

1 非均勻有理B樣條

給定一個單調不減的節點矢量E= {ξ1,ξ2,…,ξi, …,ξn+p+1}, 其中ξi為第i個節點，n和p分別為基函數的個數與階次，B樣條基函數定義為：

(1)

圖1給出了定義在E={0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1}上的二次B樣條基函數。一般地，B樣條基函數僅在節點0與1處為插值函數。

圖1 二次B樣條基函數

對每個B樣條基函數乘以對應的權因子ωi,i= 1, 2, …,n, 并除以權函數W(ξ)，可以得到B樣條基函數的有理形式，即NURBS基函數：

(2)

張量積形式的二維及三維NURBS基函數為：

(3)

(4)

式中：ωi,j、ωi,j,k是二維和三維權因子；Mj,q(η)、Lk,r(ζ)分別為定義在節點矢量H=(η1,η2, …,ηj, …,ηm+q+1)和Z=(ζ1,ζ2, …,ζk, …,ζl+r+1)上的B樣條基函數。NURBS曲面和實體分別定義為：

(5)

(6)

式中Bi,j和Bi,j,k是曲面和實體的控制點。

2 三維聲場等幾何分析

2.1 控制方程

考慮邊界為剛性壁面的三維內部聲場，其控制方程及邊界條件為：

(7)

?p/?n=-jρ0ωun，邊界Г上

(8)

2.2 基于NURBS的離散方程

等幾何分析的核心思想是將精確描述幾何的NURBS基函數作為有限元中場變量的形函數。因此，式(7)中的聲壓可表示為：

(9)

式中：m、n及l分別代表ξ、η和ζ方向的基函數個數；pA為控制點A處的聲壓；R=[R1R2…Rm×n×l]T為基函數向量；p=[p1p2…pm×n×l]T為聲壓向量。

采用伽遼金法，將NURBS基函數乘以式(7)并在區域Ω上積分，得到：

(10)

應用格林公式，式(10)可以寫為：

(11)

將式(8)、(9)代入式(11)，可得聲學問題的特征方程為：

(K-k2M)p=0

(12)

式中K、M為整體剛度與質量矩陣，表示為：

(13)

(14)

3 數值結果與分析

3.1 橢圓臺形聲場特性分析

考慮截面連續變化的橢圓臺形聲場，如圖2所示，兩端橢圓的長短半徑分別為a1、b1、a2、b2，橢圓臺高度為h。采用NURBS建立其精確的幾何模型，初始節點矢量為E=H=(0, 0, 0, 1, 1, 1)，Z=(0, 0, 1, 1)，初始控制點數為3×3×2。當a1=a2，b1=b2時，聲腔為橢圓柱，其特征值問題存在解析解。取尺寸參數a1=a2=0.15 m，b1=b2=0.09 m，h=0.4 m，聲速c0=1 m/s，前8階波數的計算結果如表1所示，隨著基函數階次的提高和控制點數的增加，等幾何分析的計算結果收斂并與解析解[19]吻合?；瘮档碾A次越高，前8階波數收斂的速率越快。采用10×10×10(自由度為1 000)和12×12×12(自由度為1 728)的計算結果優于自由度為1 361和3 672時有限元[19]得到的結果，因此等幾何分析與傳統有限元相比具有更高的準確性和計算效率。

表1 橢圓柱形聲場前8階波數的收斂情況

圖2 橢圓臺形聲場

接下來分析離心率對橢圓臺聲腔固有特性的影響，截面的連續變化增加了解析法求解的難度。取a1=0.10 m，a2=0.15 m，b2=0.09 m，h=0.4 m，基函數階次p=q=r=3，控制點10×10×10，不同離心率下橢圓臺聲腔的前8階波數如表2所示。可以看出，隨著離心率逐漸增大，前8階波數呈上升趨勢，這是因為離心率的增大使b1減小，橢圓臺形聲場的幾何變化導致其波數增大。此外，離心率在0～0.4波數的增幅小于離心率在0.4～0.8波數的增幅，原因是離心率較小時，離心率對短半徑值的影響不大。當e1=0.4時，橢圓臺形聲場的模態振型如圖3所示。

圖3 橢圓臺形聲場前8階模態振型

表2 不同離心率下e1橢圓臺形聲場的前8階波數

3.2 同心橢圓環柱形聲場特性分析

同心橢圓環柱形聲場的固有特性，柱體截面如圖4所示。三維柱體的幾何參數為a1=0.10 m，b1=0.06 m，a2=0.15 m，b2=0.09 m，h=0.4 m。采用NURBS建立精確的幾何模型，初始節點矢量為E=(0, 0, 0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1, 1, 1)，H=Z=(0, 0, 1, 1)，周向、徑向和軸向的初始控制點個數分別為9、2、2，記為9×2×2。通過基函數升階和節點插入后，不同網格下的計算結果如表3所示，NURBS基函數階次為p=q=r=3?？梢钥闯鲭S著控制點數的逐漸增加，計算結果逐漸收斂，當控制點個數達到25×7×23后，繼續細化網格對結果的影響很小，所以該控制點網格將用于本節后續的計算與分析。圖5給出了橢圓環柱形聲場的前8階模態振型，可以看出第1階與第6階為軸向模態，第2階與第3階為周向模態。

圖5 橢圓環柱形聲場前8階模態振型

表3 不同網格下橢圓環柱形聲場的前8階波數

接下來分析橢圓內環的離心率對波數的影響，取a1=0.06 m，橢圓外環及柱形高度的尺寸參數同上，圖6給出了離心率e1對橢圓環柱形聲腔波數的影響。對于給定內環橢圓長半徑的柱形聲腔，內環橢圓短半徑隨著離心率e1的增大而減小，第1階與第5階波數未出現明顯變化，第2、3、7階波數逐漸增大，第4、6、8階波數逐漸減小。

圖6 離心率e1對橢圓環柱形聲腔波數的影響

3.3 橢圓偏心環柱形聲場特性分析

最后分析橢圓偏心環柱形聲場的固有特性，柱體截面如圖7所示，幾何參數為內圓半徑r=0.06 m，橢圓長半徑a=0.15 m，橢圓短半徑b=0.09 m，柱體高度h=0.4 m。采用NURBS基函數精確建立此聲腔初始幾何模型所用到的節點矢量與上節同心橢圓環柱相同。取NURBS基函數階次為p=q=r=3，控制點個數為25×7×23，表4給出了不同圓心位置下此聲腔的前8階波數。從表4可以看出，改變內環圓心O2的位置會引起聲腔波數的變化。

表4 不同圓心位置下橢圓偏心環形聲場的前8階波數

圖7 橢圓偏心環柱型聲場的截面

在圖7所示的柱體截面所在直角坐標系中分析圓心位置的影響，當圓心O2從原點(0, 0)向右沿著x軸方向移動到(0.02, 0)，第2、3、4、6、7階波數逐漸增加；然后圓心O2從(0.02, 0)向上移動到(0.02, 0.02)，再向左移動到(0.01, 0.02)，在此過程中第2、3、7階波數逐漸減小，第4階波數繼續增加，第6階與第8階波數先增大后較小。不同圓心位置下橢圓偏心環柱形聲場的第1階與第3階模態振型如圖8所示。由此可見，不同圓心位置下的第1階模態振型相似；然而對于第3階模態振型，內部偏心圓導致了幾何模型上的不對稱，故其聲場分布存在局部變化。

圖8 不同圓心位置下橢圓偏心環柱形聲場的第1階與第3階模態振型

4 結論

1)對于橢圓臺形聲腔，其體積隨著離心率增加而減小，導致前8階波數逐漸增大。

2)對于同心橢圓環柱形聲腔，隨著內環橢圓離心率增大，第1階與第5階波數變化不明顯，第2、3、7階波數逐漸增大，第4、6、8階波數逐漸減小。

3)由于幾何的不對稱，偏心圓所在位置對橢圓偏心環柱形聲場的前8階波數影響較為復雜。

4)等幾何分析具有良好的收斂性和計算精度，適用于復雜聲場固有特性的預測與評估。