同向旋轉(zhuǎn)渦對發(fā)聲的數(shù)值研究

包振忠， 秦國良， 和文強(qiáng)， 王亞洲， 穆毅偉

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

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

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

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

1 控制方程

均勻流動(dòng)狀態(tài)下的Lighthill方程為

(1)

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

(2)

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

Tij?ρ0uiuj

(3)

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

(4)

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

p′=pinc+pa

(5)

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

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

由不可壓縮連續(xù)性方程和動(dòng)量方程得偽聲壓滿足下式

(6)

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

(7)

初始條件為

(8)

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

(9)

式中:n為邊界外法線方向;上游邊界Ma前為負(fù)號(hào);下游邊界Ma前為正號(hào);垂直流場方向Ma為零。

2 數(shù)值方法

空間方向采用高精度Chebyshev譜元法進(jìn)行離散，以Chebyshev多項(xiàng)式的極值點(diǎn)作為插值點(diǎn)，利用Chebyshev多項(xiàng)式逼近待求函數(shù)。結(jié)合吸收邊界條件，式(7)相應(yīng)的Galerkin變分形式為

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

(10)

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

(11)

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

3 數(shù)值算例

3.1 數(shù)值模型

湍流產(chǎn)生的噪聲大都起源于流動(dòng)中的渦結(jié)構(gòu)，渦不僅被稱為流動(dòng)的肌腱[21]，還被稱為是流動(dòng)的聲音[22]。因此，在計(jì)算氣動(dòng)聲學(xué)中常使用渦源作為人工聲源來分析氣動(dòng)噪聲產(chǎn)生的機(jī)理。其中同向旋轉(zhuǎn)渦對[23-27]的發(fā)聲問題已經(jīng)被廣泛用于測試各種計(jì)算氣動(dòng)聲學(xué)的數(shù)值格式。

圖2 同向旋轉(zhuǎn)渦對配置圖

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

(12)

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

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

(13)

(14)

式中：Re表示對復(fù)數(shù)求實(shí)部。

同向旋轉(zhuǎn)渦對流場固有的不穩(wěn)定性是其發(fā)聲的根本原因，其聲壓波動(dòng)的解析解可通過匹配漸近展開法獲得

(15)

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

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

3.2 數(shù)值結(jié)果

圖3 點(diǎn)(0, 58.2)處聲壓隨時(shí)間的變化

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

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

(a) Mr=0.13

(b) Mr=0.08

(c) Mr=0.03

數(shù)值解與解析解在基頻處均吻合良好，但在其它頻率成分處，數(shù)值解均略低于解析解，這也是由于對源項(xiàng)進(jìn)行數(shù)值截?cái)嗨鸬摹Ａ硗鈹?shù)值解中也出現(xiàn)了諧波成分，而解析解中并沒有出現(xiàn)，這是由于八極子和更高階的諧波項(xiàng)在推導(dǎo)解析解的時(shí)候被忽略的緣故。

同時(shí)為了比較不同旋轉(zhuǎn)馬赫數(shù)時(shí)渦對的發(fā)聲情況，圖5給出了取上述Mr時(shí)100 s時(shí)刻的聲場數(shù)值模擬結(jié)果。可以看出隨著Mr的增大，聲源強(qiáng)度隨之增大，聲場中各個(gè)位置的聲壓也隨之增大，而聲波的波長減小，也就需要更精細(xì)的網(wǎng)格來對其進(jìn)行計(jì)算。高旋轉(zhuǎn)馬赫數(shù)的聲壓云圖中出現(xiàn)了不光滑的現(xiàn)象，也表明了在對高波數(shù)問題進(jìn)行數(shù)值模擬時(shí)需要提供更精細(xì)的網(wǎng)格。

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

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

為了觀察渦對在剪切流場中的發(fā)聲情況，對渦對在剪切流場存在的情況下進(jìn)行了數(shù)值模擬。剪切速度表達(dá)式如下，其分布如圖8所示。

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

(16)

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

(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時(shí)x軸方向聲壓分布

圖8 剪切速度沿y軸分布

(a) ΔU=0.1c0

(b) ΔU=0.3c0

(c) ΔU=0.5c0

4 結(jié) 論

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