伴流聲場計算的有限元方法及其應用

范一良，季振林

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

聲波是由于介質的彈性效應引起的，相對于介質的質點在運動[1]。在消聲器和管路系統中都存在氣體的流動，氣體的流動直接影響聲波的傳播，進而影響消聲器的聲學特性。對于等截面或漸變截面管道，假定管道內流體在垂直于軸線的橫截面上以均勻速度流動，此時管道內的聲傳播可用解析方法求解[1-2]。但是多數消聲器內部的流場是非均勻的，需要使用三維數值方法計算消聲器內的流場并考慮流場對聲場影響。由于問題的復雜性，通常假設傳輸聲波的是無旋、無粘的運動介質，于是速度可以由一個標量函數的梯度來表示，這個標量函數被稱為速度勢。速度勢可以進一步分解為不隨時間發生變化的平均量和小擾動的聲學量，結合流體的連續性方程、動量方程和狀態方程得到伴流下的聲波波動方程。

Danda等[3-4]在研究渦輪發動機風扇通道的聲輻射時給出了勢流中的波動方程，討論了與之相關的聲學邊界條件，并在邊界上采用近似方法用FEM編程計算出聲場的分布，但是省略了詳細的推導過程，后續相關研究[2]也是直接引用其結論。Myers[5]建立了掠過流下無粘流體中邊界面上法向聲質點振速與聲質點法向位移、法向位移的梯度、平均流和圓頻率之間的關系，法向位移又與聲壓通過聲阻抗建立起聯系，但表達式比較復雜，很難在FEM中編程實現。Eversman[6]在Myers的基礎上研究聲阻抗邊界時，受到Moehring[7]的啟發，得到了更加簡化、準確、容易編程實現的表達式，隨后該表達式被Eversman應用在可壓縮勢流的逆流定理與聲學互易性[8]和管道內可壓縮勢流聲學互易的數值實驗[9]研究中，但Eversman只探討了阻抗邊界條件，對于聲激勵邊界、穿孔阻抗邊界和無反射邊界則沒有涉及。徐貝貝等將有限元法應用于計算穿孔管消聲器的傳遞損失[10]，分析了穿孔率對消聲器聲學性能的影響，但僅限于無流情況。

Eversman省略了推導勢流中的聲波方程的過程，并且聲邊界條件的探討也僅限于阻抗邊界。本文補充勢流中聲波方程推導過程中的相關細節，運用伽遼金加權余量法建立相應的有限元弱形式，討論和分析勢流狀態下管道聲學計算的各種聲邊界條件，通過裝配和組裝得到弱形式的矩陣表達式，并給出有限元編程計算的相關細節。為驗證理論和計算程序的正確性，使用自行編寫的有限元程序計算H-Q管的傳遞損失，并與一維理論計算結果進行比較，之后將有限元法應用于分析流速和結構變化對H-Q管消聲性能的影響。

1 伴流聲學方程

假設消聲器內部是無旋、無粘、定常的氣體流動，對于任意空間位置固定的無窮小微團，微分形式的質量守恒方程為：

(1)

由于無粘，微團表面力只有壓力，忽略體積力，則動量方程可表示為：

(2)

聲音傳播過程是等熵過程，對于等熵流體，狀態方程可表示為：

(3)

式中，γ為比熱比。

聲速可以表示為：

(4)

(5)

把式(5)代入到式(2)中，得：

(6)

因為聲波傳遞是等熵過程，所以密度只是壓力的函數，式(6)可以寫成：

(7)

式(7)只與時間有關，與空間無關，所以有：

(8)

式(8)為正壓流的伯努利方程，F(t)為力勢。對狀態方程式(3)求導，并代入到式(8)的第三項中，可得：

(9)

把式(9)代入到式(8)得到：

(10)

用c和v分別表示遠離反射體和吸聲體等均勻無擾動流體中的常速聲速和流體速度。在未被聲擾動處，速度勢函數不隨時間發生變化，所以有：

(11)

將式(11)改寫成：

(12)

聲學理論是研究流場變量的小擾動行為。流體變量被分解為不隨時間變化的平均量(用下標0表示)和小擾動的聲學量(用下標a表示)。因此速度勢、速度、密度和壓強可以表示為：

(13a)

(13b)

(13c)

(13d)

式中:v0=φ0,va=φa。

將式(13(a))和(13(b))代入式(12)中，忽略高階聲學量，得：

(14)

由式(14)可知，沒有聲擾動時，流動流體質中的聲速為：

(15)

(16)

把式(16)代入式(15)中，可建立起靜止流體和運動流體中的聲速關系：

(17)

把式(15)代入到式(14)中，可得：

(18)

將式(13(c))代入到式(4)中，并將式(4)在ρ0處進行泰勒展開，保留一階項，得：

(19)

綜合式(18)和(19)可得聲密度與聲速度勢之間的關系：

(20)

把式(13(c))和(13(d))代入式(3)中，并將式(3)在ρ0處進行泰勒展開，僅保留一階項，得到：

(21)

把式(20)代入到式(21)中，得到聲壓與聲速度勢之間的關系：

(22)

把式(13(a))和(13(d))代入質量守恒方程中，即式(1)中，去掉守恒項和高階項，得：

(23)

把式(20)代入到式(23)中，可得：

(24)

對于簡諧聲場，則有：

φa(x,y,z,t)=φ(x,y,z)ejωt

(25)

把式(25)代入到式(24)中，可得：

(26)

式(26)即為勢流中的聲波方程，當平均流速為零時，則簡化為亥姆霍茲方程。

2 聲學有限元方程

對式(26)應用伽遼金加權余量法，并對第一項使用格林第一公式，得：

(27)

式中:N為權函數，n為面單元的單位外法向量，上式即為弱形式的變分公式。

2.1 聲學邊界條件

消聲器常用的邊界條件如圖1所示，包括入口處的激勵面Si，聲導納面SR，穿孔面Sp1和Sp2，剛性壁面Sw，以及出口端So，將消聲器分為兩個域，分別用Ⅰ和Ⅱ標識。流場的入口面在聲學上被當成剛性壁面來處理，這樣會使聲音只向下游傳播，不向上游傳播。邊界將分為背景流與邊界面平行和垂直兩種情況討論，當邊界面法向方向與平均流速相切時，則有：

圖1 消聲器邊界條件

φ0·n=0

(28)

把式(28)代入式(27)的等式右端中，則有：

(29)

式中，un為法向聲質點振速。

Myers[5]給出了掠過流下無粘流體中簡諧聲激勵下邊界面上法向聲質點振速un與法向位移ζn，平均流v0以及壁面單位法向量之間n的關系為：

un=jωζn+v0·ζn-n·(n·v0)ζn

(30)

把式(30)代入到式(29)中，可得：

n·(n·v0)ζn]dS

(31)

當邊界面與平均流相切時，Eversman簡化處理了式(31)，即：

(32)

式中，ζn為法向位移。

對于非滑移壁面，壁面上的背景流速為0，此時式(32)為：

(33)

對于如圖1中所示的邊界面，當與背景流相切時，可在式(33)的基礎上進一步推導；當與背景流垂直時，需重新處理式(27)的等號右端部分。

(1)聲激勵邊界

在消聲器進口管的壁面上，如圖1中Si所在區域，需要給定一個聲激勵計算整個聲學域內的聲學響應。邊界面與背景流相切，當激勵為法向位移時，有：

(34)

當激勵為加速度時，a=-ω2ζ，則有：

(35)

(2)聲阻抗邊界

在面SR上，聲阻抗Z0已知，壁面的法向位移與聲壓的關系為[11]：

jωZ0ζn=pa

(36)

把式(22)代入式(36)中，壁面上背景流速為0，得到位移與聲速度勢之間的關系，即：

(37)

把式(37)代入式(33)中，可得：

(38)

(3)穿孔聲阻抗邊界

假設穿孔處只有掠過流，通過流忽略不計，則式(33)依然成立。當穿孔部分掠過流速變化較大時，可以分段處理，假設每一段的平均流速相等。在圖1所示的穿孔面Sp1和Sp2上，聲壓、聲速度勢、聲質點法向位移、介質密度和形函數分別用p1、p2、φ1、φ2、ζn1、ζn2、ρ1、ρ2、N1和N2來表示。穿孔面之間的穿孔導納記為ζp，穿孔面Sp1和Sp2上聲質點法向位移與聲壓和穿孔導納之間的關系為：

(39)

將式(39)中的ζn1和ζn2用p1和p2來表示，并代入到式(33)中，得：

(40)

(4)剛性壁面

剛性壁面上，如圖1中Sw所在區域，背景流速和法向質點振速均為0。把式(22)代入式(27)的右端的積分中，有：

(41)

式中，vn=v0·n為邊界面上的法向平均流的速度。

(5)無反射端

計算消聲器的傳遞損失時，需將圖1所示的出口面So設置為無反射端。有時在出口管的壁面設置為聲阻抗已知，使反射聲波被聲阻抗壁面吸收，也是一種解決方案。如果分析的上限頻率在出口管的平面波截止頻率以下，則出口管內傳遞的只有平面波。本文從解析的角度推導無反射時出口面上聲壓和聲質點振速之間的關系。

假設出口管為一段足夠長的等截面直管，管道內流體介質沿軸線z方向以均勻馬赫數M運動，并與出口管橫截面垂直，管道內介質密度的變化忽略不計。管道內傳遞的聲波為平面波，聲速度勢φ沿x和y方向的梯度為0，因此式(26)可簡化為沿管道軸線的一維聲傳播方程，即：

(42)

式中，k0=ω/c0為波數。

式(42)的通解可寫成如下形式：

φ=(C1e-jk0z/(1+M)+C2ejk0z/(1-M))

(43)

式中，C1和C2為待定系數。

無反射時，C2=0，則獨立于時間項的聲壓p和法向聲質點振速un可表示為：

(44)

把式(44)代入到式(27)等號右端，可得：

(45)

式(45)即為管道中平面波的無反射端的表達式。

2.2 離散化與矩陣裝配

對于任意一種m節點的單元，其內部任意一點的場變量都可以用該單元上所有節點的場變量值來表示，聲速度勢則表示為：

(46)

式中：N為單元形函數組成的列向量，上標T表示轉置，φi為單元節點i上待求的聲速度勢。

平均流速度的矩陣形式表示為：

(47)

把式(46)和(47)代入到式(27)的左端，并在聲學域內Ⅰ和Ⅱ離散，用φ1和φ2分別表示區域Ⅰ和Ⅱ中的聲速度勢，則有：

(48)

式中:矩陣Lhs表示在域上的積分組裝而成的矩陣，矩陣Si、Mi和Ci是由單元矩陣組裝而成，即

(49)

(50)

(51)

其中，

(52)

將式(48)～(52)進行組裝，得到如下整體有限元矩陣表達式：

(53)

式(53)即為勢流中聲波方程的有限元矩陣表達式。由式(53)可知，伴流聲場的FEM矩陣為非對稱復矩陣。在實際的計算中，由于穿孔面的距離很窄，兩個對應面面積大小很接近，因此可以從數值上把穿孔矩陣P處理成對稱矩陣。所以其非對稱性不是源于穿孔矩陣P，而是式(48)中體積分中由于流的存在得到的非對稱矩陣Ci。相比有流狀態，無流時的聲學波動方程——亥姆霍茲方程，離散化后得到的對稱非共軛的復矩陣，因此無流時聲場的計算速度大概是有流時的兩倍。

3 聲學量計算

3.1 聲壓和質點振速

求解式(53)即可計算出每一個節點的聲速度勢φ，將聲速度勢代入到式(22)中便可計算出節點上的聲壓值，對聲速度勢求解梯度便可計算出聲質點振速。通過坐標變換可知速度勢對全局坐標的導數可由局部坐標表示，即：

[([N]′[Coord])-1[N]′]{φi}e

(54)

式中:[N]′為單元所有節點的形函數組成的列向量對局部坐標的導數組成的矩陣，[Coord]為單元所有節點的全局坐標組成的矩陣，{φi}e為單元的聲速度勢組成的列向量。

由式(54)可知，對于一個三維問題，必須將場變量置于體單元中才能計算出該變量的梯度。因此在計算橫截面上聲速度勢的梯度時，首先需要找到該面單元所屬的體單元，然后找到待求節點在體單元中的局部坐標，最后將局部坐標代入到式(54)中計算出聲速度勢的梯度。下面以如圖2所示的四面體單元為例加以說明。

圖2 一、二階四面體單元節點及局部坐標系

圖2中單元節點的局部坐標，如表1所示。

表1 四面體單元節點的局部坐標

找到面單元上待求節點在體單元中的節點編號，并以此在表1中找到其局部坐標，最后將點的局部坐標代入到式(54)中即可求解出該節點的聲速度勢梯度，進而計算出聲壓。

3.2 傳遞損失

傳遞損失定義為出口面為無反射端時，消聲器的入射聲功率級與透射聲功率級之差，即：

TL=10lg(Wi/Wt)

(55)

式中:Wi為入射聲功率，Wt為透射聲功率。

如圖1所示，在面S1計算入射聲功率，在面So上計算透射聲功率。由式(22)，聲壓p可表示為：

(56)

聲質點振速在法向上的分量可表示為：

un=n·φ

(57)

式中:n為面的法向量，并約定其方向與聲傳播方向相同。

假設入射波為平面波，根據有流時平面波聲壓和質點振速的解析表達式，采用聲波分解法，可以得到入口面上的入射聲壓pi和聲質點振速ui的表達式為：

(58)

當管道內介質均勻沿著軸向流動時，沿管道軸向的聲強I可表示為[12]：

(M|p|2/ρ0c0+ρ0c0|uz|2)]

(59)

式中,uz為聲質點振速沿著管道軸向z方向的分量。

在面S1上，把式(58)代入到式(59)中便可得到任意一點的聲強，并在面上積分，便可計算出入射聲功率。在出口面So上，聲學邊界為無反射端，因此由式(53)計算出的面上節點的聲速度勢即為透射聲速度勢。將其代入式(44)中，即可計算出透射聲壓和透射聲質點振速，代入式(59)計算出透射聲強，并在面So上積分，得到透射聲功率。把透射聲功率和入射聲功率代入到式(55)中計算出傳遞損失。

4 計算與分析

作為本文方法的驗證和應用，本節使用自編的有限元法計算程序計算和分析干涉式消聲器的聲學性能，此類消聲器是利用相位差產生聲波的相互干涉，從而實現消聲效果[1]。如圖3所示的H-Q管是一種典型的干涉式消聲器，由兩條支路并聯的管道組成，管道在第一個分叉點分成兩條支路，在第二個分叉點處再合成一條。聲波經過不同長度的管道產生了相位差，在第二個分叉點匯合處“相互干涉”實現了消聲效果。分支管道內的氣體流動直接影響兩列聲波間的相位差，從而影響共振頻率和消聲特性。

圖3 Herschel-Quincke管

本文使用的H-Q管為圓形管道，其幾何尺寸為d1=d2=0.05 m，l1=0.6 m。背景流場由CFD軟件計算得到，當入口馬赫數為0.1時，主管道與側支管道的平均流速分別為30.61 m/s和3.61 m/s；當入口馬赫數為0.2時，主管道與側支管道的平均流速分別為61.92 m/s和6.7 m/s。參與計算的聲學網格尺寸為4 mm，其流場信息通過CFD軟件直接映射而來。

圖4～6分別比較了無流、入口處馬赫數為0.1和0.2時，使用一維平面波理論和三維有限元法計算得到的傳遞損失。H-Q消聲器的傳遞矩陣可以在文獻[1]中找到，此處不再贅述。此外，側支管道處的端部修正對共振頻率有較大的影響，使用一維理論計算時需要予以考慮。

圖4 Ma=0時有限元法和一維理論計算結果對比

圖5 Ma=0.1時有限元法和一維理論計算結果對比

圖6 Ma=0.2時有限元法和一維理論計算結果對比

由圖4～6可以看出，對于本文給定尺寸的H-Q管消聲器，當入口處的馬赫數分別為0、0.1和0.2時，一維理論與有限元法計算的傳遞損失大約分別在2 400 Hz、2 200 Hz和2 000 Hz以下吻合良好；高于這些頻率，兩種方法計算結果開始出現差異，這是因為在主管道和側支管道的連接處的三維流場和三維聲場效應所致，此時一維理論已不再適用。圖3中A1和A2橫截面上聲壓在2 200 Hz，馬赫數為0.1時的聲壓的最大相位差分別為為0.036弧度和6.24弧度。A1面上的聲壓相位差極小，可認為是數值計算誤差，因此A1上傳遞的是平面波。A2面上的聲壓相位差則是由兩管交接處的三維波效應引起。

此外，進口處不同的馬赫數所對應的平面波截止頻率不同，截止頻率隨著馬赫數的增加而減小，因而進一步縮減了平面波的適用范圍。

圖7顯示了介質流動速度對H-Q管消聲器傳遞損失的影響，這些數據來源于圖(4)～(6)中有限元計算結果。介質流動改變了部分共振頻率和通過頻率，低頻時部分無流時的共振頻率和通過頻率分別成為了有流時通過頻率和共振頻率；馬赫數越高，傳遞損失在受到影響的通過頻率和共振頻率附近的變化越大；在中高頻，介質流動對消聲性能的影響非常明顯。

圖7 馬赫數對傳遞損失的影響

接下來研究結構改變對H-Q管消聲性能的影響。H-Q1即為圖3所示的結構；H-Q2的結構如圖8所示，包含了兩個旁支管，其中d1=d2=d3=0.05 m，l1=0.6 m，l2=0.3 m；H-Q3的結構如圖9所示，主管中增加一個簡單膨脹腔，其中d1=d2=0.05 m，d3=0.15 m，l1=0.6 m，l2=0.3 m。

圖8 H-Q2管

圖9 H-Q3管

圖10給出了三種H-Q管消聲器傳遞損失的計算結果?？梢钥闯?，結構的改變明顯影響了H-Q管的消聲特性，特別是共振頻率。H-Q2因為包含了兩個旁支管，所以傳遞損失中共振峰的數量多于H-Q1，并且共振峰向低頻方向移動，提升了低頻段的聲學性能。H-Q3主管道增加了一個簡單膨脹腔，使得整個頻段上共振峰的數量增多，從而改善了多數頻段內的消聲性能。

圖10 結構改變對H-Q消聲器的聲學性能的影響

5 結 論

本文推導了勢流中的聲波方程，應用伽遼金加權余量法建立了有限元的弱形式表達式。針對管道聲學的計算，討論了所需邊界條件的處理方法，通過離散和裝配得到了有限元矩陣方程。自行編程實現了推導的有限元算法，計算和分析了介質流動速度和結構形式的改變對H-Q管消聲器聲學性能的影響。結果表明，當入口處的馬赫數分別為0、0.1和0.2時，一維理論與有限元法計算的傳遞損失大約分別在2 400 Hz、2 200 Hz和2 000 Hz以下吻合良好；高于這些頻率，兩者的計算結果開始出現較大的差異，這是由于主管和支管交界處的三維流場和三維聲場效應所致，此時一維方法不再適用。馬赫數越高，平面波截止頻率越低，一維理論適用的頻率范圍越窄。介質流動影響H-Q管消聲器的聲學特性，特別是在共振頻率和通過頻率附件，馬赫數越高，影響越顯著。具有兩級的H-Q管，共振峰數量增多，共振頻率降低。主管道包含了簡單膨脹腔的H-Q管，共振峰的數量增多，從而提升了在多數頻段內的消聲性能。