基于關(guān)聯(lián)速度的FW-H積分四極子聲源修正模型

周志騰, 王洪平, 王士召,*, 何國威

(1.中國科學(xué)院 力學(xué)研究所 非線性力學(xué)國家重點(diǎn)實(shí)驗(yàn)室, 北京 100190; 2.中國科學(xué)院大學(xué) 工程科學(xué)學(xué)院, 北京 100049)

0 引 言

流動(dòng)噪聲既是流體力學(xué)基礎(chǔ)研究的重要難題，又是航空工程中的重要問題。Ffowcs Williams and Hawkings(FW-H)方程[1]是目前計(jì)算流動(dòng)噪聲的重要工具。隨著計(jì)算流體力學(xué)的發(fā)展，F(xiàn)W-H方程已被成功地應(yīng)用于直升機(jī)旋翼[2-3]、增升裝置[4-5]、和飛機(jī)起落架[6]等產(chǎn)生的流動(dòng)噪聲研究。特別是近年來，數(shù)值風(fēng)洞[7]的發(fā)展更是為流動(dòng)噪聲的預(yù)測提供了有力的支撐。

FW-H方程通過引入廣義函數(shù)將自由空間的Lighthill聲比擬方程推廣至有壁面的流動(dòng)問題，其中的壁面(下文中稱為FW-H積分面)既可以是真實(shí)的固體邊界，也可以是流動(dòng)中的虛擬邊界。FW-H方程是一個(gè)關(guān)于密度(或壓力)脈動(dòng)的非齊次波動(dòng)方程。遠(yuǎn)場的噪聲可在形式上表示為流場變量和格林函數(shù)的積分(下文中稱為FW-H積分)。根據(jù)FW-H積分，遠(yuǎn)場噪聲為單極子聲源、偶極子聲源和四極子聲源三部分貢獻(xiàn)的疊加。其中，與單極子聲源和偶極子聲源相關(guān)的積分分別涉及積分面上的質(zhì)量通量和非定常載荷，是相關(guān)量在FW-H積分面上的面積分。與四極子聲源相關(guān)的積分是涉及Lighthill應(yīng)力張量的體積分，其積分域應(yīng)該包含Lighthill應(yīng)力張量不為零的所有區(qū)域。計(jì)算四極子聲源相關(guān)的積分所需要的計(jì)算量和存儲(chǔ)量遠(yuǎn)大于偶極子聲源和單極子聲源相關(guān)的積分。并且，常用的實(shí)驗(yàn)測量和數(shù)值模擬數(shù)據(jù)難以包含Lighthill應(yīng)力張量不為零的全部區(qū)域。Lighthill應(yīng)力張量從FW-H積分面穿出，常會(huì)引起四極子聲源對遠(yuǎn)場噪聲的虛假貢獻(xiàn)[8-9]。

在低馬赫數(shù)流動(dòng)中，相對于非定常載荷對應(yīng)的偶極子聲源產(chǎn)生的遠(yuǎn)場噪聲，Lighthill應(yīng)力張量對應(yīng)的四極子聲源產(chǎn)生的噪聲所占比例相對較小，通常做法是忽略四極子聲源體積分項(xiàng)。但是，隨著來流馬赫數(shù)的增加或流動(dòng)邊界的復(fù)雜化，四極子聲源的影響逐漸不可忽略[10-14]。Yao等[13]研究了來流馬赫數(shù)為0.3的半球擾流，發(fā)現(xiàn)回流區(qū)尾端的四極子聲源是遠(yuǎn)場噪聲的主要來源，Tomoaki Ikeda等[10]計(jì)算了來流馬赫數(shù)為0.2、雷諾數(shù)為1×104、攻角為2°的二維NACA0012翼型繞流，結(jié)果顯示，在部分觀測點(diǎn)四極子聲源穿出FW-H面造成的虛假聲源與偶極子聲源的貢獻(xiàn)處于同一量級。

為了解決四極子積分項(xiàng)引起的虛假聲源問題，Shur等[15]通過在接近控制面處添加人工黏性的方法耗散Lighthill應(yīng)力張量。但是，Nitzkorski等[9]認(rèn)為這種方法使得被打開面的法向的輻射噪聲被低估。為了避免引入人工黏性，Wang等[11]曾在Curle積分的框架下提出了一種修正四極子聲源計(jì)算的模型。該模型利用凍結(jié)流假設(shè)將所選聲源區(qū)外的Lighthill應(yīng)力張量的體積分對遠(yuǎn)場噪聲的貢獻(xiàn)轉(zhuǎn)換為計(jì)算域邊界上四極子聲源流量的面積分，有效地恢復(fù)了低馬赫數(shù)翼型繞流問題中產(chǎn)生的四極子聲源對遠(yuǎn)場噪聲的貢獻(xiàn)。Wang等[11]的工作中使用縮放的均勻來流速度作為Lighthill應(yīng)力張量穿過計(jì)算域邊界的對流速度，需要將聲源區(qū)的邊界取至下游足夠遠(yuǎn)處，以消除近場尾跡中的Lighthill應(yīng)力張量對流速度的非均勻性對修正模型的影響。隨后，Nitzkorski等[9]將Wang等[11]的思路推廣至FW-H積分，并通過利用多出口面之間的交叉關(guān)聯(lián)函數(shù)構(gòu)建Lighthill應(yīng)力張量的通量在FW-H積分面上的對流速度，避免了依賴于經(jīng)驗(yàn)的縮放均勻來流作為對流速度。

上述工作在時(shí)域方法框架下修正計(jì)算四極子聲源對遠(yuǎn)場噪聲的貢獻(xiàn)，有效地抑制了Lighthill應(yīng)力張量穿過控制面引起的虛假聲源問題。但由于FW-H積分面的運(yùn)動(dòng)和對流速度的計(jì)算方式等問題，這些方法沒有被廣泛地應(yīng)用在頻域方法中。Lockard等[16]提出了頻域方法框架下的四極子聲源積分項(xiàng)的計(jì)算模型。該模型通過連續(xù)使用分部積分，將四極子聲源的體積分近似表達(dá)為一系列面積分的和，避免了在整個(gè)計(jì)算域內(nèi)存儲(chǔ)和計(jì)算Lighthill應(yīng)力張量，并且抑制了四極子聲源體積分計(jì)算中對遠(yuǎn)場噪聲的虛假貢獻(xiàn)。在這一過程中，Lockard等[16]根據(jù)Wang等[11]的工作以對流方程的形式描述Lighthill應(yīng)力張量通過積分面的過程，并假設(shè)對流速度取為均勻來流的速度。該假設(shè)沒有考慮積分面上的Lighthill應(yīng)力張量對流速度的非均勻性，一般難以將積分面取在近尾跡區(qū)域。

針對這一問題，本文發(fā)展一種可以在FW-H積分面上計(jì)算非均勻?qū)α魉俣鹊哪Ｐ?，用于改進(jìn)Lockard等[16]提出的四極子修正方法中的均勻?qū)α魉俣燃僭O(shè)。該模型受時(shí)域方法中Nitzkorski等[9]的工作的啟發(fā)，基于FW-H積分面上的關(guān)聯(lián)函數(shù)計(jì)算對流速度。與Nitzkorski等[9]的工作不同的是，我們根據(jù)Lighthill應(yīng)力張量的關(guān)聯(lián)函數(shù)，而不是Lighthill應(yīng)力張量的通量的關(guān)聯(lián)函數(shù)，來計(jì)算對流速度，避免了依賴于經(jīng)驗(yàn)選用的均勻來流作為對流速度，并考慮了對流速度空間分布的非均勻性。

本文第1部分簡要地給出FW-H方程在頻域的解，第2部分給出FW-H積分四極子聲源項(xiàng)修正模型以及四極子聲源項(xiàng)的對流速度的具體計(jì)算方法，第3部分用二維對流渦以及低馬赫數(shù)圓柱繞流兩個(gè)標(biāo)準(zhǔn)算例檢驗(yàn)所發(fā)展的模型，第4部分為結(jié)論。

1 FW-H方程及其解

常用的FW-H方程的微分形式：

(1)

其中，

(2)

(3)

其中，yi=zi+Uit是隨FW-H積分面運(yùn)動(dòng)的坐標(biāo)下的空間位置變量，Ui為無窮遠(yuǎn)來流在yi方向的速度。不影響一般性，在本文中按照Wang等[11]和 Lockard等[16]的處理方式，取無窮遠(yuǎn)來流僅沿y1方向，即U=(U1,0,0)。則聲源項(xiàng)變?yōu)?/p>

Tij=ρuiuj

(4)

式中的速度ui和uj皆為減去無窮遠(yuǎn)來流后的速度擾動(dòng)。為了將對流波動(dòng)方程變化到頻域，我們根據(jù)Lockard[17]的工作定義傅里葉變換對如下：

(5)

利用式(5)，頻域中的FW-H方程可寫為：

(6)

其中k=ω/co為波數(shù)，Mi=Ui/co為馬赫數(shù)。

式(6)的解可寫為[17-18]：

(7)

其中IT、IL分別為厚度聲源與載荷聲源，分別對應(yīng)于單極子聲源和偶極子聲源，表達(dá)式如下：

(8)

其中x、y分別表示觀測點(diǎn)與聲源所在位置。

四極子聲源積分項(xiàng)為：

(9)

2 FW-H積分四極子聲源項(xiàng)的修正

準(zhǔn)確求解四極子聲源對遠(yuǎn)場噪聲的貢獻(xiàn)需要保存FW-H積分面以外的所有Lighthill應(yīng)力張量不為零的區(qū)域流場，如式(9)所示。當(dāng)四極子聲源的積分域不能包含所有的Lighthill應(yīng)力張量不為零的區(qū)域時(shí)，Lighthill應(yīng)力張量會(huì)穿過FW-H積分面，導(dǎo)致四極子聲源對遠(yuǎn)場噪聲的虛假貢獻(xiàn)。Lockard等[16]通過分部積分的方法將四極子體積分項(xiàng)(式(9))近似地表達(dá)為面積分的和：

IQ(x;ω)≈

(10)

在計(jì)算中一般取前三項(xiàng)。Uc為Lighthill應(yīng)力張量的輸運(yùn)速度，表示四極子聲源穿出FW-H積分面的速度，注意，式(10)在推導(dǎo)過程中假設(shè)對流速度沿流向與時(shí)間的變化可以忽略不計(jì)。Lighthill應(yīng)力張量穿出FW-H積分面的過程用如下方程模化：

(11)

按照式(10)計(jì)算四極子聲源項(xiàng)對遠(yuǎn)場噪聲的貢獻(xiàn)的關(guān)鍵是計(jì)算對流速度Uc。Lockard等[16]選取對流速度為均勻來流速度。以均勻來流為對流速度沒有考慮FW-H積分面上的對流速度的非均勻性，會(huì)導(dǎo)致遠(yuǎn)場噪聲依賴于FW-H的位置。不同于Lockard的選用均勻來流為對流速度，我們采用多個(gè)出口面Lighthill應(yīng)力張量的關(guān)聯(lián)函數(shù)構(gòu)建對流速度，其基本思想是利用關(guān)聯(lián)函數(shù)得到四極子聲源在出口面之間的輸運(yùn)時(shí)間，繼而根據(jù)出口面之間的距離得到穿出出口面的速度。與直接選用均勻來流為對流速度相比，此方法不僅提高了遠(yuǎn)場噪聲指向性預(yù)測的準(zhǔn)確性，而且減小了預(yù)測結(jié)果對FW-H面選取位置的依賴性。

構(gòu)建對流速度Uc的 Lighthill應(yīng)力張量的關(guān)聯(lián)函數(shù)為：

Rxy(τ)=

(12)

其中，Tij(y1,y2,t)與Tij(y1+Δy1,y2,t+τ)分別為上下游兩個(gè)面對應(yīng)位置的信號，Δy1為所選取的兩出口面之間的距離。T為用來計(jì)算Uc的采樣時(shí)間，在文中，我們?nèi)c僅隨空間位置變化。如圖1所示，基于泰勒凍結(jié)假設(shè)，渦在輸運(yùn)過程中不發(fā)生變形。因此，Lighthill應(yīng)力張量的輸運(yùn)速度Uc與渦的運(yùn)動(dòng)速度保持一致，并且本文考慮的渦沿方向穿過出口面，對應(yīng)單元的位置僅在流向變化。需要注意的是，由于使用了泰勒凍結(jié)流假設(shè)，下游出口面的位置應(yīng)取在速度梯度較小的區(qū)域[14]，并且用于求解出口速度的兩個(gè)出口面之間的距離不能太大。在3.2節(jié)的圓柱繞流算例中，第一個(gè)出口面到圓柱中心的距離大于或等于2.5D，兩個(gè)出口面的間距為 0.1D。為了便于計(jì)算，本文用于計(jì)算出口速度的兩個(gè)出口面都垂直于流向。式(12)的離散形式寫為：

圖1 四極子聲源修正模型的求解示意圖

Tij(y1+Δy1,y2,iΔ′+jΔ′)

(13)

Uc=Δy1/τr

(14)

3 結(jié)果與討論

3.1 二維對流渦

二維對流渦為可壓縮歐拉方程的一個(gè)精確解，是檢驗(yàn)FW-H積分修正方法的正確性的典型算例[16,20-21]。在圖2和圖3所示的對流渦算例中，渦以勻速直線運(yùn)動(dòng)沿y1方向向下游運(yùn)動(dòng)，嚴(yán)格滿足凍結(jié)流假設(shè)。壓力與速度場的表達(dá)式為:

圖2 初始時(shí)刻二維對流渦壓力分布云圖

圖3 無量綱時(shí)刻t=25時(shí)二維對流渦壓力分布云圖

u=M-a0a1y2exp((1-r2)/2)

v=a0a1(y1-Mt)exp((1-r2)/2)

(15)

其中，po為無窮遠(yuǎn)處的壓力；γ為空氣的熱容比,取1.4；y1、y2表示二維笛卡爾坐標(biāo)系；M為對流馬赫數(shù)；t為時(shí)間。

(16)

式(15)式(16)所有的量皆經(jīng)過無量綱化，參考物理量為單位長度L、聲速co以及密度ρ。Lighthill應(yīng)力張量滿足：

(17)

不失一般性，我們在計(jì)算中將FW-H積分面取為矩形，其在y1、y2方向的取值范圍分別為[-5L,5L]、[-5L,5L]，如圖3所示。當(dāng)來流馬赫數(shù)為0.2時(shí)，在無量綱時(shí)刻t=25，僅有一部分四極子聲源在FW-H積分面內(nèi)，穿出下游出口面的四極子聲源使常規(guī)的FW-H積分的計(jì)算結(jié)果產(chǎn)生誤差。我們采用式(10)和式(14)修正了Lighthill應(yīng)力張量穿出FW-H積分面所引起的誤差，如圖4所示。需要注意的是，圖中的黑色實(shí)線表示FW-H方程所計(jì)算出的(IT+IL)的相反數(shù)。由于觀察點(diǎn)的真實(shí)壓力脈動(dòng)趨于0，則不進(jìn)行修正的FW-H方程的解實(shí)際上遠(yuǎn)遠(yuǎn)高估了壓力脈動(dòng)的大小，因此誤差近似為此解的相反數(shù)，四極子聲源修正模型所得結(jié)果應(yīng)逼近于此誤差。在計(jì)算式(10)的過程中，我們?nèi)×瞬煌膎值，當(dāng)n=3時(shí)(與Lockard等[16]所取的n值相同)，我們的方法在二維渦的上下游遠(yuǎn)場觀測點(diǎn)均能夠修正Lighthill應(yīng)力張量穿出FW-H積分面引起的誤差。這一結(jié)論證實(shí)了我們的模型在二維對流渦運(yùn)動(dòng)中的有效性。

(a)Observer(-100,0)

3.2 二維圓柱繞流遠(yuǎn)場噪聲

低雷諾數(shù)圓柱繞流的遠(yuǎn)場輻射噪聲是一個(gè)廣泛應(yīng)用的檢驗(yàn)流動(dòng)噪聲計(jì)算模型和方法的算例。直接數(shù)值模擬[22]、FW-H方程[9]以及精確格林函數(shù)[23]等方法皆被應(yīng)用于此問題當(dāng)中，并且相應(yīng)的結(jié)果顯示出一致性。其中，Inoue等[22]較早地利用直接數(shù)值模擬的手段系統(tǒng)研究了雷諾數(shù)150下，二維圓柱繞流輻射噪聲的指向性、遠(yuǎn)場衰減特性以及遠(yuǎn)場壓力脈動(dòng)關(guān)于馬赫數(shù)的標(biāo)度率，并常被用作比對算例。本文的結(jié)果將與Inoue 等[22]的直接數(shù)值模擬結(jié)果進(jìn)行比較，以檢驗(yàn)所發(fā)展的模型的有效性。

我們通過數(shù)值求解不可壓流動(dòng)的Navier-Stokes方程作為遠(yuǎn)場噪聲的聲源。聲源數(shù)值模擬的設(shè)置如圖5所示，圓柱中心在原點(diǎn)，直徑為D，來流速度為U0，基于來流速度與圓柱直徑定義的雷諾數(shù)為150，坐標(biāo)軸y1、y2分別對應(yīng)流向與垂向。流場計(jì)算采用浸入邊界方法，計(jì)算區(qū)域?yàn)閇-10D,15D]，[-12D,12D]，均勻網(wǎng)格的尺度為0.01D，流場計(jì)算的時(shí)間步長為0.001D/U0。FW-H積分面選為長方形，上游以及兩側(cè)的面位置分別為y1=-D與y2=±2D，計(jì)算聲場的采樣步長為0.027D/U0。為了檢驗(yàn)結(jié)果對出口面位置的敏感度，我們選取了不同位置的出口面，分別為y1=4.5D、y1=2.5D，計(jì)算了馬赫數(shù)為0.2下，位于(0D,100D)處觀測點(diǎn)的壓力脈動(dòng)，并與Inoue 等[22]的直接數(shù)值模擬結(jié)果進(jìn)行了對比。結(jié)果如圖6所示，其中Tp為所選觀察點(diǎn)處壓力脈動(dòng)的周期，對于不同位置的出口面，我們的模型均能夠消除四極子聲源穿出FW-H積分面引起的誤差，與直接數(shù)值模擬的結(jié)果相符。此結(jié)果表明我們的模型不依賴于出口面位置的選取。

圖5 二維圓柱繞流渦量分布云圖與遠(yuǎn)場噪聲計(jì)算所用FW-H面的位置示意圖

圖6 馬赫數(shù)0.2下，不同位置的出口面得到的(0D,100D)處觀測點(diǎn)的壓力脈動(dòng)，并與Inoue的直接數(shù)值模擬結(jié)果對比

圖7比較了利用各種方法得到的遠(yuǎn)場壓力脈動(dòng)的指向性，|Δp|為壓力脈動(dòng)幅值的均方根。參考Inoue 等[22]對觀測點(diǎn)位置的選取，圖7中的觀測點(diǎn)到圓柱中心的距離定義為75[1+Mcos(θ)]D，其中M為來流馬赫數(shù)，此處為M=0.2，方向角θ為觀測點(diǎn)位置對應(yīng)向量與y1軸正向的夾角。圖中的黑色實(shí)線為Inoue 等[22]的直接數(shù)值模擬的結(jié)果。紅色點(diǎn)劃線為利用本文所發(fā)展的模型的計(jì)算結(jié)果，藍(lán)色雙點(diǎn)劃線為四極子聲源對流速度取為來流速度的結(jié)果，紫色虛線為未經(jīng)FW-H積分四極子聲源修正的結(jié)果。在本文所選用的FW-H積分面下，如果在用FW-H方程計(jì)算遠(yuǎn)場噪聲中不考慮四極子聲源修正，將會(huì)引起指向性和壓力脈動(dòng)幅值嚴(yán)重偏離DNS值。如表1所示，在本文坐標(biāo)系下，Inoue等[22]給出的最大壓力脈動(dòng)在101.5°，未經(jīng)修正的結(jié)果在120.0°，經(jīng)均勻來流為對流速度的四極子聲源修正后，最大壓力脈動(dòng)在92.0°處，可經(jīng)本文所發(fā)展的方法修正后，最大壓力脈動(dòng)在100.0°，指向性的最大誤差降至2°以下。在壓力脈動(dòng)幅值誤差明顯的120°方向角處，未經(jīng)修正的幅值的誤差為31%，經(jīng)均勻來流為對流速度的四極子聲源修正后，誤差為12%，可經(jīng)本文所發(fā)展的方法修正后，誤差則降為4%以下。這一結(jié)果表明了本文所提出的四極子聲源修正模型的有效性，并且本文所提出的模型改進(jìn)了Lockard等[16]所提出的頻域FW-H積分修正方法。

圖7 二維圓柱繞流中，利用各種方法計(jì)算出的遠(yuǎn)場壓力脈動(dòng)指向性

表1 最大壓力脈動(dòng)方向

4 結(jié) 論

針對FW-H積分中的四極子聲源項(xiàng)常在遠(yuǎn)場噪聲的計(jì)算中引起虛假聲源的問題，本文在頻域方法的框架下改進(jìn)了四極子聲源項(xiàng)的修正模型，用于修正Lighthill應(yīng)力張量穿過積分域邊界引起的誤差。該工作基于泰勒凍結(jié)流假設(shè)?；疞ighthill應(yīng)力張量穿過積分域邊界的過程，特色是基于關(guān)聯(lián)函數(shù)計(jì)算Lighthill應(yīng)力張量的對流速度，而不是將對流速度取為均勻來流。所提出的模型考慮了對流速度的空間非均勻性，改善了非均勻流動(dòng)區(qū)域FW-H積分面對遠(yuǎn)場噪聲的影響。利用二維對流渦算例以及低雷諾數(shù)圓柱繞流算例檢驗(yàn)了所提出的模型，結(jié)果表明模型準(zhǔn)確修正了四極子聲源穿出FW-H積分面所導(dǎo)致的遠(yuǎn)場噪聲計(jì)算誤差。作為在頻域方法中考慮FW-H積分面上的非均勻?qū)α魉俣鹊牡谝徊焦ぷ?，我們在二維流動(dòng)的例子中驗(yàn)證了基于關(guān)聯(lián)函數(shù)構(gòu)建Lighthill應(yīng)力張量的可行性，并將進(jìn)一步把該模型用于三維流動(dòng)的四極子聲源的修正。