張樹海, 李 虎, 王益民
(中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點實驗室, 四川 綿陽 621000)
激波和旋渦是可壓縮流動的基本結(jié)構(gòu)。在許多流動中,比如超聲速混合層、超聲速射流和燃燒等,激波與旋渦共存,它們之間會發(fā)生相互作用,如激波與激波之間、激波與旋渦之間以及旋渦與旋渦之間的相互作用,這種作用,不僅會引起激波和旋渦的變形,還會產(chǎn)生噪聲。因此,吸引了很多學(xué)者開展理論、實驗和數(shù)值模擬研究。
1955年,Hollingsworth 和Richards[1]首先研究了激波與旋渦的相互作用,他們利用激波管產(chǎn)生平面激波,該激波掃過一個具有一定攻角的機(jī)翼后,產(chǎn)生一個旋渦,激波碰到激波管端部的固壁后反射回來,與旋渦發(fā)生干擾。該實驗發(fā)現(xiàn),激波與旋渦干擾后,伴隨著激波與旋渦變形會產(chǎn)生聲波。Dosanjh 和Week[2]對此問題進(jìn)行了較細(xì)致的實驗,測量了第一道聲波的強(qiáng)度,成為后來實驗研究和數(shù)值模擬研究的一個標(biāo)準(zhǔn)。Barbosa 和Skews[3]設(shè)計了一個分叉型激波管,激波管所產(chǎn)生的激波被分成兩支,一支通過激波管內(nèi)的拐角產(chǎn)生旋渦,另一支直接與旋渦發(fā)生干擾,研究發(fā)現(xiàn),激波與旋渦干擾后,在旋渦中心附近一個很小的區(qū)域,出現(xiàn)壓力脈沖區(qū),其脈沖壓力是周圍壓力的兩倍以上。
為了解釋聲波的產(chǎn)生機(jī)理,早在20世紀(jì)60年代,Ribner[4]就發(fā)展了線性分析理論。他利用傅立葉變換將旋渦分解成平面正弦剪切波,這種單正弦波與激波相互作用后產(chǎn)生聲波,經(jīng)過合成后得到的激波與旋渦干擾產(chǎn)生的聲波在垂直于激波的方向反對稱分
布,其理論結(jié)果與Dosanjh和Week[2]的測試結(jié)果的反對稱部分一致。后來,Week和Dosanjh[5]發(fā)展了Lightill[6]的聲比擬理論,將周向壓力分布展開成四極、兩極和單極聲源之和。單極的作用使得聲波偏離了單純反對稱態(tài),他們所得到的聲波與其實驗一致。
到20世紀(jì)90年代,伴隨著高階精度捕捉激波格式的成熟,采用數(shù)值模擬研究激波和旋渦干擾并捕捉聲波已成為研究聲波產(chǎn)生機(jī)理的重要手段。Ellzey等[7-9]采用四階差分格式,通過求解非定常可壓縮Euler方程,系統(tǒng)地模擬了激波與旋渦的干擾,得到的第一道聲波強(qiáng)度與實驗測得的聲波強(qiáng)度基本一致。Inoue等[10-12]利用六階緊致格式直接模擬了不同強(qiáng)度激波與旋渦相互作用及聲波的產(chǎn)生過程,給出了弱激波與旋渦相互作用聲波的產(chǎn)生及傳播機(jī)理。他們通過精細(xì)的模擬發(fā)現(xiàn),當(dāng)激波碰到旋渦的瞬間,旋渦邊緣發(fā)生局部壓縮和膨脹,形成雙極的第一道聲波,隨著激波與旋渦干擾的繼續(xù),出現(xiàn)第二組壓縮和膨脹區(qū),此壓縮和膨脹區(qū)與前面出現(xiàn)的壓縮和膨脹區(qū)域交替出現(xiàn),形成第一道完整的聲波,第一道聲波演變成四極波。激波通過旋渦后,在旋渦兩側(cè),形成兩道反射激波,反射激波后伴隨第二道和第三道聲波。Grasso和Pirrozli[13-16]也系統(tǒng)地模擬了較大范圍激波與旋渦的干擾問題,并根據(jù)激波變形的形狀,將激波與旋渦的干擾分成三類:弱干擾,具有規(guī)則反射的強(qiáng)干擾,具有馬赫反射的強(qiáng)干擾。
除了激波與旋渦之間的相互作用外,在包含激波和旋渦的復(fù)雜流場中,還存在旋渦與旋渦之間的相互作用,兩個或多個旋渦之間的相互作用也是噪聲產(chǎn)生的重要機(jī)制。比如,兩個同向旋轉(zhuǎn)的旋渦的合并過程,是剪切層噪聲的重要部分。
這種包含激波、旋渦的復(fù)雜流場的非定常運動產(chǎn)生的噪聲是激波噪聲。目前,雖然氣動聲學(xué)理論已有相當(dāng)?shù)陌l(fā)展,但是關(guān)于激波噪聲的產(chǎn)生機(jī)制認(rèn)識還很不夠。伴隨高階精度數(shù)值格式和大型并行計算機(jī)的快速發(fā)展,數(shù)值模擬已成為研究激波噪聲產(chǎn)生機(jī)制的重要手段,但是激波噪聲的計算仍然是計算氣動聲學(xué)所面臨的相當(dāng)嚴(yán)峻的挑戰(zhàn),具體表現(xiàn)在:
(1) 流場的強(qiáng)非定常特性。激波噪聲是激波與其它流場結(jié)構(gòu)相互作用產(chǎn)生的,其中存在激波變形和產(chǎn)生過程,與激波相互作用的流場結(jié)構(gòu)也隨之發(fā)生變化。比如,在激波與旋渦的相互作用過程中[17-20],入射激波會變形并產(chǎn)生反射激波,在旋渦中心附近還產(chǎn)生小激波串,同時旋渦被壓扁,而且出現(xiàn)扭轉(zhuǎn)運動,當(dāng)激波強(qiáng)度達(dá)到一定程度之后,會引起旋渦的破碎。在激波與剪切層相互作用過程中,存在有激波的局部震蕩[21]。因此,傳統(tǒng)的定常流體動力學(xué)問題計算方法或是大尺度問題的非定常方法已不再適用。盡管已發(fā)展了多種高階精度Runge-Kutta方法[22],但大多Runge-Kutta方法是顯式格式,由于穩(wěn)定性條件的限制,允許的時間步長非常小,對包含邊界層的流動計算效率非常低,因此仍需要發(fā)展時間離散精度和效率更高的計算方法。
(2) 流場的多尺度特性。激波噪聲是流場的一部分,但是在很多情況下,它的擾動尺度比流體動力學(xué)尺度小很多[23]。比如在馬赫數(shù)為1.5的噴流中,距噴流直徑40倍處測得的噪聲強(qiáng)度是124 dB,聲場的脈動速度與噴流的速度之比是1.5×10-4。特別地,當(dāng)流場中存在強(qiáng)激波時,聲波引起的物理量脈動比激波引起的變化尺度差距更大。這就要求不僅能精確計算宏觀尺度的流體動力學(xué)量,能夠捕捉激波,還能精確計算微弱尺度的聲波,對計算方法和計算網(wǎng)格的要求遠(yuǎn)比動力學(xué)計算高。
(3) 模擬激波噪聲所需的計算區(qū)域遠(yuǎn)大于模擬氣動力所需的計算區(qū)域。研究激波噪聲既關(guān)心近場流體動力學(xué)的變化,挖掘噪聲產(chǎn)生的機(jī)制,又關(guān)心遠(yuǎn)場激波噪聲的傳播特性,計算區(qū)域邊界與聲源的距離非常遠(yuǎn),這是與傳統(tǒng)計算流體動力學(xué)的鮮明差異。
(4) 激波與聲波的計算對格式需求存在矛盾的因素。捕捉激波需要格式具有一定的耗散,而聲波是近乎無色散、無耗散的等熵運動,聲波的模擬需要高分辨率和無耗散格式。因此,現(xiàn)在還不存在一種較為理想的高階精度、高分辨率、低耗散的激波捕捉格式。目前,聲波的計算普遍采用線性緊致格式,精度和分辨率非常高,耗散性小,適合于聲波的計算。但是線性格式不能捕捉流場中的強(qiáng)激波,而捕捉激波普遍采用WENO格式[24]。盡管WENO格式的魯棒性非常好,可以模擬很強(qiáng)的激波,但同所有的捕捉激波格式一樣,在激波下游會產(chǎn)生微弱的非物理波動。雖然這種波動不影響流場的整體特性,但其強(qiáng)度常常會比流場中的噪聲高。另外,雖然可以設(shè)計很高精度的WENO格式[25],但格式的耗散性比較大,對聲波的分辨率不滿意。盡管已有很多種簡單的混合格式[26]、優(yōu)化WENO格式[27]、非線性加權(quán)緊致格式[28-29]以及DRP格式[30],試圖提高捕捉激波類格式的分辨率,但實際應(yīng)用比WENO格式并沒有明顯的優(yōu)勢,因此目前還不存在一種既能捕捉激波又適用于聲波計算的較理想數(shù)值方法。
針對這些問題,我們近年開展了較為系統(tǒng)性的研究,解決了計算激波噪聲的一些方法上存在的問題,對包含有激波、旋渦和聲波的一系列問題進(jìn)行了直接數(shù)值模擬,揭示了激波噪聲產(chǎn)生機(jī)制。本文對這些研究做簡要回顧。
WENO格式是一種高階精度捕捉激波的格式[24],被廣泛用于工程實際和科學(xué)計算。但是像所有其它的激波捕捉格式一樣,用WENO格式計算定常激波時,在激波下游會出現(xiàn)微弱的非物理波動,導(dǎo)致殘差不能收斂到滿意程度。雖然這種波動不會影響整體的氣動特性,但是至少會造成兩方面的影響,第一,殘差是結(jié)束計算的重要判據(jù),殘差無法收斂到很小的值會使計算無法判斷何時停止,只能靠經(jīng)驗結(jié)束進(jìn)程;第二,雖然激波后的非物理波動是微弱的,不影響氣動力等宏觀流動特性,但是對于激波噪聲,這種微弱的非物理波動則非常大,有時比激波噪聲的幅度高幾個數(shù)量級,因此計算氣動噪聲問題,應(yīng)該消除這種微弱的非物理波動。經(jīng)過系統(tǒng)的數(shù)學(xué)物理分析和大量的數(shù)值實驗,我們提出了兩種方法消除這種微弱的非物理波動[31-32]。
1) 針對五階WENO格式,我們提出了一種新的光滑測試因子如下:
(1)
具體形式為:
(2)
采用這種光滑測試因子后,消除了激波下游的微弱非物理波動,對典型問題計算可以收斂到機(jī)器零。
2) 前面的方法只適合于五階精度的格式,對更高階精度格式不適用,而且會引起格式在高階極值點附近降階。經(jīng)過大量的數(shù)值實驗,我們發(fā)現(xiàn)WENO格式在特征投影過程中,采用Roe平均方法計算半點上的物理量確定特征值和特征矩陣,是影響WENO格式不收斂的重要原因。為此,我們提出采用迎風(fēng)插值方法代替原來的Roe平均方法。即:
具體的迎風(fēng)插值公式如下:
一階插值:U(1)=UiU(2)=Ui+1
(3)
(4)
最優(yōu)加權(quán)插值:
(5)
(6)
最優(yōu)插值精度與WENO重構(gòu)過程精度一致,如五階精度的WENO格式采用五階插值,七階精度的WENO格式采用七階插值,九階精度的WENO格式采用九階插值。圖1和圖2是采用兩種方法對馬赫數(shù)2的一維定常激波計算得到的密度分布和收斂歷程,圖1中ResA為殘差的平均值??梢钥吹?,改進(jìn)的WENO格式可以收斂到機(jī)器零。值得指出的是,采用最優(yōu)加權(quán)插值方法求解半點上的值,使WENO格式具有自相似特性,并一致高階精度,詳見文獻(xiàn)[31-32]。
WENO格式是一種高階精度捕捉激波格式,雖然格式的精度可以很高,但是對聲波的分辨率并不理想,而且耗散性也比較大,采用WENO格式模擬湍流結(jié)構(gòu)、激波噪聲等復(fù)雜多尺度流動并不理想。為了提高格式的分辨率,降低耗散,并保持捕捉激波的魯棒性,我們將WENO格式的設(shè)計思想應(yīng)用到緊致格式中,設(shè)計了一種類譜分辨率的高階精度緊致格式。這一格式是基于Lele的半點型線性緊致格式[33]設(shè)計的,Lele的半點型線性緊致格式形式是:
(7)
該格式左端是網(wǎng)格點上的函數(shù)導(dǎo)數(shù)值,右端是半點的函數(shù)值,它最大的優(yōu)點是分辨率非常高,接近譜方法分辨率。但是右端半點的函數(shù)值是未知的,Lele[33]設(shè)計了一種線性緊致插值的方法求得半點的值。圖3是Lele的半點型線性緊致格式的設(shè)計模板,可以看出,該格式有兩個缺點:(1)設(shè)計模板包括了整點和半點的函數(shù)值,但Lele僅用了部分信息,格式?jīng)]有達(dá)到該模板上的最優(yōu)精度;(2)半點上的值需要插值方法求得,任何插值方法都存在傳遞誤差,導(dǎo)致格式對高波數(shù)的分辨率大幅度下降。
為克服Lele格式的缺點,我們設(shè)計了一種類譜分辨率的高階精度緊致格式。在網(wǎng)格節(jié)點上,我們采用如下格式計算其一階導(dǎo)數(shù):
(8)
采用同樣的格式計算半點上的一階導(dǎo)數(shù)值:
(9)
求得網(wǎng)格點和半點處的一階導(dǎo)數(shù)之后,通過求解相應(yīng)的控制方程,即可求得下一時刻的網(wǎng)格點和半點處函數(shù)值。
與Lele的格式相比,式(8)、式(9)給出的數(shù)值格式有兩個優(yōu)點:1) 格式的精度有一定的提高,Lele格式最高達(dá)到十階精度,我們設(shè)計的格式可以達(dá)到十四階;2) 網(wǎng)格點和半點上的函數(shù)值都是通過同一格式計算得到,消除了原格式緊致插值引起的傳遞誤差,保證了計算過程的精度和分辨率。圖4給出了六階和八階格式修正波數(shù)及其與 Lele 格式的比較,相比之下我們改進(jìn)的格式分辨率與精確值相差無幾,詳見文獻(xiàn)[34]。
線性格式不能捕捉流場中的強(qiáng)激波,為了使格式具有捕捉強(qiáng)激波的功能,我們將WENO重構(gòu)思想應(yīng)用到插值技術(shù),發(fā)展了高階WENO插值方法,并將該方法應(yīng)用于式(8)、式(9)右端,初步發(fā)展了一種非線性類譜分辨率的緊致格式。具體的WENO插值方法是:在格式的設(shè)計模板S2r-1={xi-r+1,...,xi+r-1}上,采用插值函數(shù):
(10)
可以得到半點上的通量值:
(11)
與WENO格式重構(gòu)類似,將模板S2r-1分成三個子模板,在每個子模板上,存在r階插值:
(12)
采用加權(quán)思路,可以得到半點上的加權(quán)插值:
(13)
在式(8)的右端,采用如下混合插值方法替代原來的精確值fi+1/2,即:
(14)
而在式(9)中,采用類似的混合插值方法。得到的非線性格式,能夠保證格式的精度,具有較好的捕捉激波能力,同時格式的分辨率比現(xiàn)有的WENO明顯提高,耗散性也大幅度降低,詳見文獻(xiàn)[35]。
Gaussian渦是一種典型的旋渦,在一定的條件下,兩個同向旋轉(zhuǎn)的Gaussian渦會合并[36],對于合并條件,已有很多研究[37]。1995年,Mitchell等[36]采用直接數(shù)值模擬方法,研究了兩個同向旋轉(zhuǎn)的
(a)r=λ/2
(b)r=2λ
圖7兩個同向旋轉(zhuǎn)的Gaussian渦合并過程中
產(chǎn)生的聲波及其比較
Fig.7Far-fieldpressuretracesatr=λ/2andr=2λandthecomparisonbetweenourdirectnumericalsimulation(DNS)andthatbytheM?hring’sequation
值得指出的是,我們的結(jié)果與Mitchell的結(jié)果并不一致,其原因可能是由于初始參數(shù)的微弱不一致性,導(dǎo)致兩個旋渦合并時間的差異。為了驗證我們結(jié)果的正確性,我們進(jìn)行了網(wǎng)格收斂性研究,初始條件的影響以及采用不同方法對此問題進(jìn)行了模擬,結(jié)果一致,詳見文獻(xiàn)[20]。
Taylor渦也是一種典型的旋渦,與Gaussian渦結(jié)構(gòu)不同,Taylor渦具有雙層結(jié)構(gòu)。為了深入了解渦聲的產(chǎn)生機(jī)理,我們系統(tǒng)研究了兩個Taylor旋渦的相互作用,并根據(jù)旋渦特性,將兩個Taylor旋渦的相互作用分成四類:(1) 兩個反向旋轉(zhuǎn)的Taylor渦的相互作用;(2) 兩個同向旋轉(zhuǎn)Taylor渦的相互作用;(3) 兩個同向旋轉(zhuǎn)Taylor渦的合并;(4) 兩個尺度或強(qiáng)度相差較大的旋渦的相互作用。
2.2.1 兩個反向旋轉(zhuǎn)的旋渦的相互作用
Taylor渦具有雙層結(jié)構(gòu),內(nèi)層外層渦量符號相反。如果兩個反向旋轉(zhuǎn)的旋渦距離足夠近,它們之間會發(fā)生耦合,結(jié)果導(dǎo)致內(nèi)層外層分離,兩個渦核逐漸靠攏,形成一個渦偶極子,兩個Taylor渦的外層逐漸靠攏,形成一個反向運動的渦偶極子,如圖8所示,是強(qiáng)度為Mv=0.5、初始距離為4倍旋渦半徑的兩個反向渦相互作用的渦量演化歷程。詳見文獻(xiàn)[20]。
2.2.2 兩個同向旋轉(zhuǎn)Taylor旋渦的相互作用
圖9是兩個強(qiáng)度為Mv=0.5、初始旋渦間距是4倍旋渦半徑的兩個同向旋轉(zhuǎn)的Taylor旋渦相互作用過程的渦量演化歷程。這種相互作用存在兩種特性,1) 兩個同向旋轉(zhuǎn)的旋渦將演化成兩個不對稱的旋渦區(qū),每一個渦區(qū)存在三個渦核結(jié)構(gòu);2) 與兩個反向旋轉(zhuǎn)的旋渦相互作用具有本質(zhì)區(qū)別的是,每一個渦區(qū)的兩個渦核來自同一個初始旋渦,強(qiáng)的渦核來自初始旋渦的內(nèi)層,而初始旋渦的外層逐漸從初始旋渦中分離形成圍繞渦核的兩個獨立的渦核。
2.2.3 兩個同向旋轉(zhuǎn)Taylor旋渦的合并過程
與兩個同向旋轉(zhuǎn)的Gaussian旋渦相似,在一定條件下兩個同向旋轉(zhuǎn)的Taylor旋渦也會合并。如圖10所示,是兩個強(qiáng)度為Mv=0.5的同向旋轉(zhuǎn)的Taylor旋渦的合并過程,兩個旋渦的初始距離是旋渦半徑的2倍。兩個渦核合并形成主渦核,而在渦核的外部,形成兩個旋轉(zhuǎn)臂結(jié)構(gòu),兩個渦的外層逐漸集中于旋轉(zhuǎn)臂附近,形成一個三核渦結(jié)構(gòu)。
2.2.4 兩個強(qiáng)度或尺度相差較大的旋渦的相互作用
兩個強(qiáng)度相差很大的旋渦的耦合過程是最令人感興趣的。1) 兩個強(qiáng)度或尺度相差較大的旋渦的相互作用會產(chǎn)生多核渦結(jié)構(gòu);2) 強(qiáng)度或尺度的較大差異,導(dǎo)致它們之間相互作用所產(chǎn)生的多核渦結(jié)構(gòu),具有非對稱特性。這種相互作用包括有兩種典型狀態(tài),一是強(qiáng)度相差較大的兩個Taylor渦的相互作用,另一是尺度相差較大的Taylor渦的相互作用。
圖11是兩個強(qiáng)度相差較大的Taylor渦的相互作用過程,其中強(qiáng)旋渦的強(qiáng)度為Mvu=-0.8,順時針旋轉(zhuǎn),弱旋渦的強(qiáng)度為Mvd=0.25,逆時針旋轉(zhuǎn),它們之間的初始距離為半徑的4倍。可以看到,與兩個等強(qiáng)度的反向旋渦相互作用類似,兩個強(qiáng)度相差較大的Taylor渦相互作用會產(chǎn)生兩個渦偶極子。兩個旋渦的渦核逐漸靠攏,形成一個強(qiáng)的渦偶極子,兩個旋渦的外層逐漸從原旋渦中分離,形成一個弱的渦偶極子。由于其強(qiáng)度相差較大,兩個渦偶極子都呈現(xiàn)非對稱特性,導(dǎo)致弱的渦核圍繞強(qiáng)渦核旋轉(zhuǎn),較弱的渦偶極子圍繞強(qiáng)渦偶極子旋轉(zhuǎn)。
圖12是兩個尺度相差較大的反向旋渦相互作用的過程,其中兩個旋渦的強(qiáng)度都是0.5,一個逆時針旋轉(zhuǎn),一個順時針旋轉(zhuǎn)。兩個旋渦的尺寸相差較大,上面的旋渦半徑為1, 而下面的旋渦半徑為0.2, 兩個旋渦的初始距離為2.2??梢钥吹剑⌒郎u的渦核穿過大旋渦的外層,進(jìn)入到大旋渦的渦核外圍,并圍繞大旋渦的渦核不斷旋轉(zhuǎn),最終形成一個橢圓形的渦核,而旋渦的外層形成兩個渦核結(jié)構(gòu),圍繞大的渦核旋轉(zhuǎn),形成一個三核渦結(jié)構(gòu),并持續(xù)做蛙跳運動。
2.2.5 兩個Taylor渦相互作用過程中聲波產(chǎn)生機(jī)理分析
兩個Taylor渦的相互作用包括豐富的動力學(xué)特性,這種過程會產(chǎn)生噪聲。圖13和圖14 是兩個典型狀態(tài)兩個Taylor渦相互作用的聲壓Δp=(p-p0)/p0等值線和其徑向分布。旋渦的強(qiáng)度均是0.5,兩個旋渦初始距離都是半徑的4倍。其中圖13是兩個反向Taylor渦相互作用所產(chǎn)生的噪聲,圖14是兩個同向Taylor渦相互作用所產(chǎn)生的噪聲??梢钥吹剑瑑煞N旋渦相互作用產(chǎn)生的噪聲具有很大區(qū)別,兩個反向旋渦相互作用僅產(chǎn)生幾個聲脈沖,而兩個同向旋轉(zhuǎn)的Taylor渦相互作用,與兩個同向旋轉(zhuǎn)的Gaussian渦相互作用很類似,會產(chǎn)生持續(xù)的強(qiáng)噪聲。
以兩個尺度相差較大的旋渦相互作用為例,分析噪聲產(chǎn)生機(jī)理。我們計算了Lamb量·(ρω×u),它代表著噪聲產(chǎn)生的源。如圖15所示,是Lamb量的等值線演化過程??梢钥吹?,旋渦相互作用具有三個明顯的運動:1) 初始圓形旋渦被扭曲,弱旋渦被撕扯成兩個渦帶,他們與強(qiáng)旋渦繼續(xù)相互作用;2) 盡管強(qiáng)旋渦的渦心位置幾乎不變,但整個渦核被壓成橢圓形;3) 兩個渦帶與強(qiáng)旋渦的渦核之間相互作用產(chǎn)生一種蛙跳式運動。圖16是監(jiān)測點(x,y)=(100,0)和(0,100)的聲壓隨時間變化曲線,壓力峰值P1、P2、P3、P4、P5、P6和P7位于時間軸103、121、139、163、190、223和262。考慮到聲從渦心附近產(chǎn)生聲波到達(dá)監(jiān)測點所需要的時間,它們對應(yīng)于長軸在y軸的方向,而聲波谷M1、M2、M3、M4、M5、M6和M7位于時間軸114、130、151、176、206、241和285,它們對應(yīng)于長軸位于x的方向。這意味著,噪聲主要由旋渦的蛙跳運動產(chǎn)生。
激波與旋渦相互作用是激波噪聲的一個簡化模型,通過研究激波與旋渦相互作用,可以深入了解激波噪聲產(chǎn)生機(jī)制。通過對激波與旋渦相互作用系統(tǒng)的模擬,我們發(fā)現(xiàn)激波與強(qiáng)旋渦相互作用存在多級干擾。第一級干擾是入射激波和初始旋渦的相互作用,第二級干擾是反射激波和變形旋渦的相互作用,第三級干擾是渦心附近產(chǎn)生的小激波串和變形旋渦的干擾,如圖17所示。與此同時,旋渦的變形也存在多級特征。在第一級干擾中,初始圓形旋渦被壓成橢圓形;在第二級干擾中,橢圓形旋渦被壓回成圓形旋渦;在第三級干擾中,圓形旋渦又被壓縮成橢圓形。此外,這種多級干擾會引起旋渦在水平線附近做振蕩運動。根據(jù)激波與強(qiáng)旋渦相互作用過程中的多級干擾過程,當(dāng)時我們預(yù)測,這種多級干擾過程會產(chǎn)生更為復(fù)雜的聲波。這一預(yù)測被Chatterjee和Vijayaraj[39]的研究所證實,如圖18所示。此外,我們還系統(tǒng)研究了激波與渦列相互作用過程流場結(jié)構(gòu)和聲波產(chǎn)生機(jī)制[17-19]。
激波與剪切層相互作用是一個典型的問題,也可以看做是研究超聲速噴流激波噪聲的一個簡化模型。為了了解超聲速噴流激波噪聲的產(chǎn)生機(jī)制,我們設(shè)計了一個激波與剪切層相互作用的模型。如圖19所示,一道斜激波打到一個剪切層上,下邊界是一個反射固壁,透射激波打到固壁上后,形成反射激波,反射激波再次與剪切層相互作用發(fā)生反射,形成類似于超聲速噴流的激波串結(jié)構(gòu)。通過較為系統(tǒng)的數(shù)值模擬,我們揭示了兩種噪聲產(chǎn)生機(jī)制,一是激波與旋渦相互作用,發(fā)生在圖20所示的區(qū)域1, 另一是激波泄露機(jī)制,發(fā)生在圖20所示的區(qū)域2。
圖21是第一個區(qū)域入射激波與剪切層中旋渦相互作用過程的一個演化周期中幾個典型時刻的脹量圖,其中紅色虛線標(biāo)記的是剪切層下層流體中被追蹤的兩道聲波。我們看到由于流體以超聲速向下游傳播,聲波也隨著流體向下游方向?qū)α?,并且聲波的輻射半徑在逐漸增大。當(dāng)t=91.69時,產(chǎn)生了第二道聲波,可以明確地看到聲波的產(chǎn)生源在渦與小激波相互作用的位置,這和張樹海等[17-19]的激波與旋渦相互作用的研究結(jié)果相同。剪切層的渦列穿過激波,每個渦與激波發(fā)生干擾時,根據(jù)這種激波噪聲產(chǎn)生和傳播的機(jī)制,在噪聲源處向外輻射出一道道的聲波,同時我們可以看到當(dāng)聲波傳播到下壁面時聲波反射回內(nèi)場,從而聲波再次與激波和剪切層相互干擾。
圖22是第二個區(qū)域反射激波與剪切層相互作用過程中一個演化周期中幾個典型時刻的脹量圖,其中紅色虛線標(biāo)記的是剪切層下層流體中的一道聲波。我們看到由于流體以超聲速向下游傳播,聲波也隨著流體向下游方向?qū)α?,并且聲波的輻射半徑在逐漸增大。當(dāng)t=91.69時,噪聲源處產(chǎn)生了第二道聲波,可以明確地看到聲波的產(chǎn)生源在旋渦之間的鞍點處,即辮子區(qū)的位置,這和Manning等[41-43]的激波泄漏機(jī)制完全相同。在剪切層的渦列穿過激波時,渦對與激波發(fā)生干擾,在渦對之間的鞍點處激波透過剪切層,同時在鞍點處激波被辮子區(qū)的旋臂所阻礙,激波與辮子區(qū)的相互干擾使得激波強(qiáng)度變?nèi)?,部分激波穿過剪切層,還有一部分激波沒有穿過剪切層,而是在鞍點處以聲波的形式發(fā)生泄漏。這說明了強(qiáng)激波與剪切層相互作用中也存在激波泄漏機(jī)制,在鞍點處激波泄漏并向外輻射聲波。根據(jù)這種激波噪聲產(chǎn)生和傳播的機(jī)制,當(dāng)剪切層穿過激波時,在噪聲源處向外輻射出一道道的聲波,同時我們可以看到當(dāng)聲波傳播到下壁面時聲波反射回內(nèi)場,從而聲波再次與激波和剪切層相互干擾。
超聲速噴流噪聲是一種典型的激波噪聲,與很多工程密切相關(guān)。采用五階WENO格式直接求解軸對稱Navier-Stokes方程,數(shù)值模擬了射流馬赫數(shù)(完全膨脹)為Mj=1.19(嘯聲呈現(xiàn)軸對稱模態(tài),三維效應(yīng)不明顯)的軸對稱欠膨脹射流,其對應(yīng)的射流聲速馬赫數(shù)為Ma=1.05,雷諾數(shù)為Re=6.216×105,噴管直徑為D=25.4 mm,噴管厚度為0.4D,環(huán)境溫度為288.15 K。
圖23給出了噴管唇口壁面[0.0,0.642D]和噴管出口平面[0.0,2.0D]上壓力監(jiān)測點處噪聲信號的頻譜信息,包括頻率和聲壓級。結(jié)果顯示,在大于5000 Hz的頻率區(qū)間內(nèi),有4個尖銳的聲壓級突起:前兩個聲壓級突起的頻率為6567 Hz(128 dB)和8637 Hz(123 dB),分別是嘯聲A1模態(tài)和嘯聲A2模態(tài);另外兩個聲壓級突起的頻率為12087 Hz(125 dB)和14310 Hz(123 dB),分別是嘯聲B模態(tài)的諧頻和嘯聲A0模態(tài)。
(a) 噴管唇口壁面
(b) 噴管出口平面
圖23監(jiān)測點處壓力信號的譜分析
Fig.23Soundpressurelevel
綜合流場結(jié)構(gòu)和聲場信息可以得到的嘯聲模態(tài)的圖像及其產(chǎn)生位置,其中流場結(jié)構(gòu)用數(shù)值紋影(密度梯度的模)表示,聲場由脹量場表示。根據(jù)嘯聲頻率可以得到嘯聲的波長,利用波長可以在脹量場中清晰地分辨出嘯聲A1模態(tài)和A2模態(tài),相應(yīng)的波長分別為λ=2.04D和λ=1.51D。從圖24中可以清晰地分辨出向上游傳播的嘯聲,并且嘯聲產(chǎn)生位置都是在第三個激波柵格和第五個激波柵格之間的區(qū)域。
聲波產(chǎn)生在近流場區(qū)域,流體介質(zhì)的非定常運動產(chǎn)生的噪聲經(jīng)過復(fù)雜流場結(jié)構(gòu)傳播達(dá)到遠(yuǎn)場。當(dāng)噪聲穿過復(fù)雜流場結(jié)構(gòu)過程中,會與流場相互作用,發(fā)生畸變。其中,聲波穿過激波是一個典型問題,我們采用直接數(shù)值模擬方法,研究了聲波穿過激波的畸變特性。
初始條件為:
(15)
(16)
(17)
計算區(qū)域為x∈[-10,10],馬赫數(shù)M=2.0,聲波幅值ε=1.0d-5,頻率ω=0.6π,比熱比取γ=1.4,計算結(jié)果如圖25所示,由圖可以看出,聲波穿過馬赫數(shù)M=2.0的激波,振幅放大了約3.64倍多。
當(dāng)聲波穿過旋渦時,聲波的幅值和頻率都會發(fā)生變化。采用我們發(fā)展的線性緊致格式,通過求解Navier-Stokes方程,我們研究了平面聲波穿過等熵渦的散射特性。
圖26是計算模型,平面聲波與一個二維Taylor渦相互作用。計算區(qū)域的左端入射聲波為:
(18)
二維Taylor渦置于計算區(qū)域的中心。圖27是入射聲波頻率ω=0.6π與強(qiáng)度為Mv=0.25的Taylor渦相互作用的瞬態(tài)脈動壓力場云圖。可以看到,旋渦對聲波的影響非常明顯。在旋渦下游,形成兩條強(qiáng)噪聲干擾條帶,在強(qiáng)條帶外側(cè),有二次條帶,由于旋渦以逆時針方向旋轉(zhuǎn),干擾條帶關(guān)于y=0不對稱。圖28是散射聲壓的均方根(RMS)等值線圖,圖29是離渦心r=10的圓周上散射聲壓的指向性分布。
近年來,針對激波噪聲計算的需求,我們開展了數(shù)值方法研究,通過對典型問題的直接數(shù)值模擬,揭示了激波噪聲產(chǎn)生機(jī)理。
(1) 在數(shù)值方法方面,針對高階精度捕捉激波的WENO格式不收斂問題,提出了一種新的光滑測試因子,發(fā)展了一種迎風(fēng)插值技術(shù)計算半點上物理量,用以確定WENO重構(gòu)過程中特征投影所需物理量。所改進(jìn)的WENO格式,消除了激波下游的微弱非物理波動,對典型問題計算可以收斂到機(jī)器零。發(fā)展了一種類譜方法的緊致格式,這一格式包括線性和非線性兩部分,其中線性格式具有高階精度、高分辨率和低耗散特性,是模擬低速氣動聲學(xué)問題的理想算法。非線性格式捕捉激波能力與WENO格式相當(dāng),分辨率明顯提高,耗散性大幅度降低。
(2) 采用直接數(shù)值模擬方法,系統(tǒng)研究了兩個旋渦相互作用、激波與旋渦/渦列相互作用、激波與剪切層相互作用以及超聲速噴流激波噪聲。發(fā)現(xiàn)激波與強(qiáng)旋渦相互作用具有多級干擾特性,提出每級干擾都有自己的聲波產(chǎn)生。
當(dāng)然,對激波噪聲的認(rèn)識還相當(dāng)有限,有許多問題有待進(jìn)一步解決。
參 考 文 獻(xiàn):
[1]Hollingsworth M A, Richards E J. On the sound generated by the interaction of a vortex and a shock wave[R]. British Aeronautical Research Council, Fluid Motion Sub-Committee, 18257, FM2371, 1956.
[2]Dosanjh D S, Weeks T M. Interaction of a starting vortex as well as a vortex street with a travelling shock wave[J]. AIAA J, 1965, (3): 216-223.
[3]Barbosa F J, Skews B W. Shock wave interaction with a spiral vortex[J]. Phys Fluids, 2001, 13: 3049-3060.
[4]Ribner H S. The sound generated by interaction of a single vortex with a shock wave[R]. UTIA Report No. 61, 1959.
[5]Weeks T M, Dosanjh D S. Sound generation by shock vortex interaction[J]. AIAA J, 1967, (5): 660-669.
[6]Lighthill M J. On sound generated aerodynamically I: general theory[J]. Proceedings of the Royal Society of London, 1952, 211: 564-587.
[7]Ellzey J L, Henneke M R, Picone J M, et al. The interaction of a shock with a vortex: shock distortion and the production of acoustic waves[J]. Phys Fluids, 1995, 7: 172-184.
[8]Ellzey J L, Henneke M R. The shock-vortex interaction: the origins of the acoustic wave[J]. Fluid Dyn Res, 1997, 21: 171-184.
[9]Ellzey J L, Henneke M R. The acoustic wave from a shock-vortex interaction: comparison between theory and computation[J]. Fluid Dyn Res, 2000, 27: 53-64.
[10]Inoue O, Hattori Y. Sound generation by shock-vortex interaction[J]. J Fluid Mech, 1999, 380: 81-116.
[11]Inoue O. Propagation of sound generated by weak-shock interaction[J]. Phys Fluids, 2000, 12: 1258-1261.
[12]Inoue O, Takahashi T, Hatakeyama N. Separation of reflected shock waves due to secondary interaction with vortices: Another mechanism of sound generation[J]. Phys Fluids, 2002, 14: 3733-3744.
[13]Grasso F, Pirozzoli S. Shock-wave thermal inhomogeneity week shock-vortex interaction in a mixing zone[J]. AIAA J, 1999, 33: 1797-1802.
[14]Grasso F, Pirozzoli S. Shock-wave-vortex interactions: shock and vortex deformations, and sound production[J]. Theroet Comput Fluid Dyn, 2000, 13: 421-456.
[15]Grasso F, Pirozzoli S. Simulations and analysis of the coupling process of compressible vortex pairs: Free evolution and shock induced coupling[J]. Phys Fluids, 2001, 13: 1343-1366.
[16]Pirozzoli S, Grasso F, D’Andrea A. Interaction of a shock wave with two counter-rotating vortices: Shock dynamics and sound production[J]. Phys Fluids, 2001, 13: 3460-3481.
[17]Zhang S, Zhang Y T, Shu C W. Multistage interaction of a shock wave and a strong vortex[J]. Phys Fluids, 2005, 17: 116101.
[18]Zhang S, Zhang Y T, Shu C W. Interaction of an oblique shock wave with a pair of parallelvortices: Shock dynamics and mechanism of sound generation[J]. Phys Fluids, 2006, 18(12): 126101.
[19]Zhang S, Jiang S, Yong-Tao Zhang, et al. The mechanism of sound generation in the interaction between a shock wave and two counter-rotating vortices[J]. Phys Fluids, 2009, 21: 076101.
[20]Zhang S, Li H, Liu X, et al. Classification and sound generation of two-dimensional interaction of two Taylor vortices[J]. Physics of Fluids, 2013, 25: 056103.
[21]Lui C V. C. M. A numerical investigation of shock-associated noise[D]. Stanford University, 2003.
[22]Gottlib S, Shu C W. Total variation diminishing Runge-Kuttaschemes[J]. Mathematics of Computation, 1998, 67: 73-85.
[23]Tim Colonius, SanjivaK Lele. Computational aeroacoustics: progress on nonlinear problems of sound generation[J]. Progress in Aerospace Sciences, 2004, 40: 345-416.
[24]Jiang G S, Shu C W. Efficient implementation of weighted ENO schemes[J]. Journal of Computational Physics, 1996, 126: 202-228.
[25]Balsara D S, Shu C W. Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order accuracy[J]. Journal of Computational Physics, 2000, 160: 405-452.
[26]Pirozzoli S. Conservative hybrid compact-WENO schemes for shock-turbulence interaction[J]. J Comput Phys, 2002, 178: 81-117.
[27]Wang Z J, Chen R F. Optimized weighted essentially non-oscillatory schemes for linear waves with discontinuity[J]. J Comput Phys, 2001, 174: 381-404.
[28]Deng X, Zhang H. Developing high-order weighted compact nonlinear schemes[J]. J Comput Phys, 2000, 165: 22-44.
[29]Zhang S, Jiang S, Shu C W. Development of nonlinear weighted compact schemes with increasingly high order accuracy[J]. J Comput Phys, 2008, 227: 7294-7321.
[30]Tam C K W, Webb J C. Dispersion-Relation-Preserving finite difference schemes for computational acoustics[J]. J Comput Phys, 1993, 107: 262-281.
[31]Zhang S, Shu C W. A new smoothness indicator for the WENO schemes and its effect on the convergence to steady state solutions[J]. Journal of Scientific Computing, 2007, 31: 273-305.
[32]Zhang S, Jiang S, Shu C W. Improvement of convergence to steady state solutions of Euler equations with the WENO schemes[J]. Journal of Scientific Computing, 2011, 47: 216-238.
[33]Lele S K. Compact finite difference schemes with spectral-like resolution[J]. Journal of Computational Physics, 1992, 103: 16-42.
[34]Liu X, Zhang S, Zhang H, et al. A new class of central compact schemes with spectral-like resolution I: Linear schemes[J]. Journal of Computational Phys, 2013, 248: 235-256.
[35]Liu X, Zhang S, Zhang H, et al. A new class of central compact schemes with spectral-like resolution II: Hybrid weighted nonlinear schemes[J]. Journal of Computational Phys, 2015, 284: 133-154.
[36]Mitchell B E, Lele S K, Moin P. Direct computation of the sound from a compressible co-rotating vortex pair[J]. Journal of Fluid Mechanics, 1995, 285: 181-202.
[37]Cerretelli and Williamson C H K. The physical mechanism for vortex merging[J]. J Fluid Mech, 2003, 475: 41-77.
[38]Mohring W. On vortex sound at low Mach number[J]. J Fluid Mech, 1978, 85: 685-691.
[39]Chatterjee A, Vijayaraj S. Multiple sound generation in interaction of shock wave with strong vortex[J]. AIAA Journal, 2008, 46: 2558-2567.
[40]Liu X, Zhang S. Direct numerical simulation of the interaction of 2D shock wave and shear layer[J]. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45: 61-75. (in Chinese)劉旭亮, 張樹海. 二維激波與剪切層相互作用的直接數(shù)值模擬研究[J]. 力學(xué)學(xué)報, 2013, 45: 61-75.
[41]Manning T A. A numerical investigation of sound generation in supersonic jet screech[D]. Stanford University, 1999.
[42]Suzuki T, Lele S K. Shock leakage through an unsteady vortex-laden mixing layer: application to jet screech[J]. Journal of Fluid Mechanics, 2003, 490: 139-167[43]Lui C C M. A numerical investigation of shock-associated noise[D]. Stanford University, 2003.
[44]Li H. Numerical study of shock-associated noise in axisymmetric supersonic jet[D]. China Aerodynamics Research and Development Center, 2016. (in Chinese)李虎. 軸對稱超聲速射流激波噪聲數(shù)值模擬研究[D]. 中國空氣動力研究與發(fā)展中心, 2016.