同向旋轉渦對發聲的數值研究

包振忠， 秦國良， 和文強， 王亞洲， 穆毅偉

(西安交通大學 流體機械及工程系， 西安 710049)

Lighthill[1]提出了著名的聲比擬理論，其后的60多年里，聲比擬理論得到了快速發展，并且在氣動聲學問題中得到了廣泛應用。而Ribner[2-3]和Meecham[4]則從偽聲壓的角度對Lighthill方程作了修改，得到了一個新的方程，命名為膨脹理論(Dilatation Theory)，其等價于后來Hardin等[5]提出的黏聲分離技術。與Lighthill四極子源項不同，膨脹理論以偽聲壓的時間二階導數作為方程源項，簡化了源項求解時的復雜性；并且，與聲比擬理論相比，膨脹理論能夠更好的解釋發聲機理。Flemming等[6]結合大渦模擬和膨脹理論研究了湍流火焰的燃燒噪聲，數值結果與試驗結果吻合良好，驗證了膨脹理論在求解氣動聲學問題時的適用性。Hiramoto等[7]應用流動可視化技術和靜壓測量技術，并結合膨脹理論研究了湍流剪切層的發聲情況，結果表明膨脹理論的源項與渦結構以及渦結構的運動有極強的關聯性。Escobar等[8]用有限元法計算了同向旋轉渦對的發聲問題，并比較了應用Lighthill聲比擬理論和膨脹理論時的求解結果，結果表明膨脹理論求解過程更加簡捷，且能夠得到更好的數值結果。

由于所面臨問題的特性不同，傳統的CFD方法無法滿足計算氣動聲學問題的研究需要[9-10]。因此，需要發展一種具有高分辨率、低色散和低耗散的時空離散方法，用于對氣動聲學問題進行高精度數值模擬。譜元法[11](Spectral Element Method, SEM)具有有限元方法和偽譜法的思想，兼備有限元可以模擬復雜介質模型的韌性和偽譜法的精度，能夠較好地捕捉和傳播聲波。譜元法最初是Patera在1984年提出，后由 Canuto等[12]多位學者發展并應用于求解流動問題，獲得了較為理想的結果。 Seriani等[13]應用譜元方法求解了聲波動方程并對其進行了色散、耗散分析。 Hasbestan等[14]應用最小二乘譜元法求解了層流壓縮流動問題。Canuto等利用 DG譜元法求解了非線性雙曲問題，為了捕捉到間斷點加入了人工黏性，成功捕捉到了激波的發生位置。耿艷輝等[15-18]將譜元法應用于對氣動聲學問題、波動問題的數值模擬中，結果體現出了該方法的優越性。

本文將膨脹理論擴展到有背景流場存在的情形下，并應用 Chebyshev譜元法對控制方程進行空間離散，時間離散采用隱式 Newmark法,邊界采用吸收邊界條件，對同向旋轉渦對在均勻流場中的發聲問題進行了數值模擬，同時研究了不同背景流場對渦對發聲的影響，為深入理解流致噪聲的產生機理聲提供了一定的理論基礎。

1 控制方程

均勻流動狀態下的Lighthill方程為

(1)

在只考慮水平方向來流的情況下

(2)

式中：c0和U分別是聲速和均勻來流速度。p′為壓力波動，Tij為Lighthill應力張量，在低馬赫數、高雷諾數的情況下可表示為

Tij?ρ0uiuj

(3)

取馬赫數Ma=U/c0,則式(1)可寫成如下形式

(4)

根據Lighthill聲比擬理論，聲壓波動將滿足上述非齊次對流波動方程。但是Ribner和Meecham分別提出在不可壓縮流動情況下，壓力波動p′可以分解成聲壓pa和偽聲壓pinc(又被稱為不可壓縮壓力波動或水動力學壓力波動)

p′=pinc+pa

(5)

聲壓的特征是它的傳播特性，而偽聲壓是由于流體流動的慣性所引起的，沒有波的傳播特性，這也是其“偽聲壓”名稱的由來。聲壓與偽聲壓的衰減特性也極為不同，如圖1所示。從圖中可以看出，在近場偽聲壓占據壓力波動主導地位，而隨著傳播距離的增加，偽聲壓迅速衰減，聲壓則可以傳播至很遠的距離。

圖1 聲壓與偽聲壓衰減曲線

由不可壓縮連續性方程和動量方程得偽聲壓滿足下式

(6)

將式(3)、(6)代入式(4)可得以聲壓為未知量的非齊次對流波動方程

(7)

初始條件為

(8)

邊界條件采用E-D-T[19]吸收邊界條件，其表達式為

(9)

式中:n為邊界外法線方向;上游邊界Ma前為負號;下游邊界Ma前為正號;垂直流場方向Ma為零。

2 數值方法

空間方向采用高精度Chebyshev譜元法進行離散，以Chebyshev多項式的極值點作為插值點，利用Chebyshev多項式逼近待求函數。結合吸收邊界條件，式(7)相應的Galerkin變分形式為

?pa,v∈H1(Ω),t∈[0,T]

(10)

式中：Γxup，Γxdown和Γy分別為上、下游和y方向邊界；v為檢驗函數。

(11)

式中：M，C，K，S分別代表質量矩陣、阻尼矩陣、剛度矩陣和源項，具體離散過程和表達式可參見文獻[15-16]。時間推進方式采用隱式Newmark法[20]。

3 數值算例

3.1 數值模型

湍流產生的噪聲大都起源于流動中的渦結構，渦不僅被稱為流動的肌腱[21]，還被稱為是流動的聲音[22]。因此，在計算氣動聲學中常使用渦源作為人工聲源來分析氣動噪聲產生的機理。其中同向旋轉渦對[23-27]的發聲問題已經被廣泛用于測試各種計算氣動聲學的數值格式。

圖2 同向旋轉渦對配置圖

在不可壓縮、無黏的條件下，渦對產生的流場可以用復位勢函數φ(z,t)來表示

(12)

式中：z=x+iy=reiθ,b=r0eiωt。

從式(12)中可以得到計算聲源時需要的流場量

(13)

(14)

式中：Re表示對復數求實部。

同向旋轉渦對流場固有的不穩定性是其發聲的根本原因，其聲壓波動的解析解可通過匹配漸近展開法獲得

(15)

式中：k=2ω/c0；J2(kr)、Y2(kr)分別為第一和第二類二階貝塞爾函數；ψ=2(ωt-θ)。

計算區域選為200 m×200 m的方形區域，將源項區域限定在中間范圍的100 m×100 m區域。取旋轉半徑r0=1 m，環量Γ=1.005 31m2/s，聲速c0=1 m/s。根據耿艷輝等對譜元法求解波動方程時的影響因素研究，本文中所有計算均取均勻單元，且x和y方向的單元數Nm=Nn=30，單元內的插值節點數Nx=Ny=2，因此在全部區域所形成的網格都為均勻網格；時間步長Δt取0.1 s.

3.2 數值結果

圖3 點(0, 58.2)處聲壓隨時間的變化

圖3給出了Ma= 0,Mr= 0.079 6時，聲壓在點(0，58.2)處數值解和解析解沿時間的變化歷程曲線。從圖中可以看出，在聲波到達該點后，該點處聲壓數值解與解析解吻合良好，色散誤差極小，只有在波峰與波谷處，數值解幅值略微小于解析解，這是由于為消除原點奇異性對源項進行數值截斷(即認為在r/r0<=1.5范圍內不存在源項)，從而導致了總體源項強度偏小的緣故，這也與其它文獻的結果一致。但本文采用的單元數遠小于上述文獻中的單元數，因此這也從側面也體現了譜元方法的高精度優勢。

為了比較不同旋轉馬赫數對渦對發聲頻率的影響，分別取Mr= 0.03，0.08，0.13，并在三種Mr時對點(0，58.2)處的聲壓進行監視，得到各自的時間歷程曲線，然后對其進行快速傅里葉變換，結果如圖4所示。

(a) Mr=0.13

(b) Mr=0.08

(c) Mr=0.03

數值解與解析解在基頻處均吻合良好，但在其它頻率成分處，數值解均略低于解析解，這也是由于對源項進行數值截斷所引起的。另外數值解中也出現了諧波成分，而解析解中并沒有出現，這是由于八極子和更高階的諧波項在推導解析解的時候被忽略的緣故。

同時為了比較不同旋轉馬赫數時渦對的發聲情況，圖5給出了取上述Mr時100 s時刻的聲場數值模擬結果。可以看出隨著Mr的增大，聲源強度隨之增大，聲場中各個位置的聲壓也隨之增大，而聲波的波長減小，也就需要更精細的網格來對其進行計算。高旋轉馬赫數的聲壓云圖中出現了不光滑的現象，也表明了在對高波數問題進行數值模擬時需要提供更精細的網格。

大部分氣動聲學問題中，聲波均在均勻流場中進行傳播。圖6中計算了X方向來流Ma=0.3時，不同時刻渦對的發聲情況。隨著時間的推進，聲波逐漸向外傳播，但往上游和下游的傳播速度差距很大，且上游的幅值和波數都得到增加，而在下游區域則相應降低，由于背景流場的存在，多普勒效應得到了良好的體現。

為更好的闡釋圖6中所觀察到的現象，圖7給出了縱向位置y分別在±40 m處沿著水平方向的壓力曲線圖。

為了觀察渦對在剪切流場中的發聲情況，對渦對在剪切流場存在的情況下進行了數值模擬。剪切速度表達式如下，其分布如圖8所示。

u(y)=ΔUtanh(2y/δ)

(16)

式中:ΔU和δ分別是剪切最大速度和剪切層厚度。分別選取ΔU=0.1c0, 0.3c0和0.5c0，δ=10 m進行了數值計算。模擬得到的100 s時刻聲壓云圖如圖9所示。由于剪切速度的作用，波陣面由圓形變成了橢圓形，且隨著剪切速度的增加，橢圓越來越扁，圖9(c)剪切馬赫數雖然達到了0.5，但仍然得到了較好的數值結果。本算例直觀揭示了剪切流場對聲傳播的影響。

(a) Mr=0.03

(b) Mr=0.08

(c) Mr=0.13

(a) t=33 s

(b) t=66 s

(c) t=99 s

圖7 來流Ma=0.3時x軸方向聲壓分布

圖8 剪切速度沿y軸分布

(a) ΔU=0.1c0

(b) ΔU=0.3c0

(c) ΔU=0.5c0

4 結 論

本文實現了將高精度譜元法結合膨脹理論應用于求解不可壓縮流動引起的氣動聲學問題。將膨脹理論擴展到有背景流場存在的情形，并應用E-D-T吸收邊界條件，成功對由兩個點渦形成的同向旋轉渦對的聲場進行了數值模擬與分析。另外計算了同向旋轉渦對在不同旋轉馬赫數、均勻流場、剪切流場下的發聲情況，并且對其頻譜成分進行了分析。結果表明：高精度譜元法具有極低的色散和耗散特性，能夠對計算氣動聲學問題進行高精度的數值求解；渦對所輻射的聲場隨著旋轉馬赫數的增大，聲波波長減小，聲源強度增大，峰值頻率也隨之增大，也就需要更精細的網格來對其進行計算；在均勻流場和剪切流場影響的情況下，輻射聲場呈現出了典型的多普勒效應；另外，E-D-T吸收邊界條件的使用，有效降低了聲波在邊界處的反射，適用于存在背景流場時的聲場數值模擬。