三維聲腔模態(tài)的等幾何計算分析

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

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

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

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

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

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

1 非均勻有理B樣條

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

(1)

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

圖1 二次B樣條基函數(shù)

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

(2)

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

(3)

(4)

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

(5)

(6)

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

2 三維聲場等幾何分析

2.1 控制方程

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

(7)

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

(8)

2.2 基于NURBS的離散方程

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

(9)

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

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

(10)

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

(11)

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

(K-k2M)p=0

(12)

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

(13)

(14)

3 數(shù)值結(jié)果與分析

3.1 橢圓臺形聲場特性分析

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

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

圖2 橢圓臺形聲場

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

圖3 橢圓臺形聲場前8階模態(tài)振型

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

3.2 同心橢圓環(huán)柱形聲場特性分析

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

圖5 橢圓環(huán)柱形聲場前8階模態(tài)振型

表3 不同網(wǎng)格下橢圓環(huán)柱形聲場的前8階波數(shù)

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

圖6 離心率e1對橢圓環(huán)柱形聲腔波數(shù)的影響

3.3 橢圓偏心環(huán)柱形聲場特性分析

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

表4 不同圓心位置下橢圓偏心環(huán)形聲場的前8階波數(shù)

圖7 橢圓偏心環(huán)柱型聲場的截面

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

圖8 不同圓心位置下橢圓偏心環(huán)柱形聲場的第1階與第3階模態(tài)振型

4 結(jié)論

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

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

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

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