劉秋洪，王垿桁，薛絲丹，何嘉華

(西北工業(yè)大學(xué) 翼型葉柵空氣動力學(xué)重點實驗室，西安 710072)

0 引言

經(jīng)典聲比擬理論建立了聲場與不同類型聲源之間的響應(yīng)關(guān)系，成為工程中應(yīng)用最廣泛的氣動噪聲數(shù)值預(yù)測方法。聲比擬方法先采用高精度數(shù)值方法計算有限域內(nèi)氣動聲源信息，然后利用聲比擬積分方程模擬噪聲傳播。FW-H 方程[1]及其擴(kuò)展形式[2,3]采用封閉可滲透面作為積分面，完整解由面積分和體積分兩部分組成。當(dāng)主要氣動聲源包圍在可滲透面內(nèi)時，可忽略體積分對聲場的貢獻(xiàn)[4]，從而僅利用可滲透面上厚度源和載荷源信息外推積分面外聲場，比如解析公式F1A[5]和F1C[3]。然而，經(jīng)典聲比擬理論僅關(guān)注聲壓等標(biāo)量場的預(yù)測，沒有對聲學(xué)速度矢量（聲波擾動引起的介質(zhì)粒子振動速度矢量，與聲速是兩個不同的概念）的外推做出規(guī)定，因而是不完整的[6]。聲學(xué)標(biāo)量僅表征聲波的傳播狀態(tài)，不能清晰闡述聲波的能量輸運。聲強矢量是聲壓和聲學(xué)速度的函數(shù)，表示單位時間內(nèi)通過單位面積的噪聲能量通量，其方向就是噪聲能量的輸運方向。聲強場的可視化可以詳細(xì)顯示噪聲能量的傳播路徑，有助于揭示聲能量輻射機制，Mao 等[7]利用聲強可視化說明了旋轉(zhuǎn)聲源輸出聲能的三種模式。噪聲傳播路徑上的非緊致結(jié)構(gòu)會導(dǎo)致聲散射，從而改變聲場大小、波形和指向性。盡管聲散射概念在實際應(yīng)用中有許多變化[8-9]，但都需要預(yù)測散射邊界處的聲學(xué)速度作為邊界條件。噪聲傳播路徑分析還有助于噪聲抑制方案的研究，Lee 等[10]揭示了散射面對聲能再分布的影響，Mao 等[11]研究了阻抗邊界的主要吸聲位置。

根據(jù)線化歐拉方程，采用聲壓梯度可間接表示聲學(xué)速度。對靜止的聲傳播介質(zhì)，F(xiàn)arassat 等[12]提出了一個半解析公式來計算聲壓梯度。Lee 等[13]導(dǎo)出了用于壓力梯度計算的時域公式G1 和G1A，結(jié)合等效源法評估了旋翼噪聲的聲散射[14]。為考慮均勻流的對流效應(yīng)，Ghorbaniasl 等[15]建立了一個頻域聲壓梯度解析公式，Bi 等[16]則發(fā)展了時域聲壓梯度計算公式G1-M 和G1A-M。聲壓梯度積分解析公式的提出推動了聲比擬方法向矢量場領(lǐng)域拓展。但這些公式在數(shù)學(xué)表達(dá)上都很復(fù)雜，且計算量大。

提高聲場預(yù)測效率對解決氣動聲學(xué)實際問題非常重要。對運動介質(zhì)中的聲傳播，聲壓梯度不能直接表示聲學(xué)速度，因此需要發(fā)展更簡明、高效和可矢量化的聲學(xué)公式。Ghorbaniasl 等[17]基于FW-H 方程提出了聲學(xué)速度時域計算公式V1 和V1A。Mao 等[18]將公式V1A 推廣到頻域，得到了公式FV1A 和FV2A，用于計算角速度不變的旋轉(zhuǎn)單極子和偶極子聲源誘導(dǎo)的聲學(xué)速度。Mao 等進(jìn)一步利用這些聲學(xué)速度公式分析了簡諧點源的聲強[7]和聲散射[19]。恰如Lee 等[20]指出的那樣，聲學(xué)速度公式V1 和V1A 可以直接從聲壓梯度公式G1 和G1A 導(dǎo)出。基于聲學(xué)介質(zhì)靜止假設(shè)，Mao 等[21]提出了以聲學(xué)速度為變量的矢量波動方程，引入了另一種推導(dǎo)公式V1 和V1A 的方法。質(zhì)量和動量守恒方程在形式上類似真空電動力學(xué)的控制方程—麥克斯韋方程，可通過定義“流體電場”和“流體磁場”轉(zhuǎn)換為“流體麥克斯韋方程”[22]。Dunn[6]將電磁類比概念應(yīng)用于聲比擬理論，建立了同時包含聲壓標(biāo)量和聲學(xué)速度矢量（三個分量）的四維聲比擬表述，但聲學(xué)速度的可滲透面源項處理不正確，且沒有考慮運動介質(zhì)的對流效應(yīng)。

為考慮均勻流對流效應(yīng)，Mao 等[23]發(fā)展了一個對流矢量波動方程，并提出了相應(yīng)的聲學(xué)速度半解析公式，本文將其命名為公式V1-M。當(dāng)均勻流速度為零時，公式V1-M 退化為公式V1。然而對均勻流中的偶極子點源聲輻射，半解析公式V1-M 并不能得到正確的預(yù)測結(jié)果。他們將數(shù)值結(jié)果的錯誤歸結(jié)為偶極子會誘導(dǎo)渦波擾動，使得線化歐拉方程不能描述聲學(xué)速度。事實上，均勻流中偶極子點源的聲輻射是純聲學(xué)問題，不會誘發(fā)渦波擾動。

本文研究目的是通過分析對流矢量波動方程的右端源項，闡明半解析公式V1-M 的理論不足，提出修正的聲學(xué)速度積分解析公式CV1A-M，并利用均勻流中單極子和偶極子點源聲輻射的理論解驗證公式CV1A-M 的正確性。

1 聲學(xué)速度解析公式V1A-M

假設(shè)無窮遠(yuǎn)均勻來流的密度、壓力和速度矢量分別為 ρ∞、p∞和u∞，將當(dāng)?shù)亓鲃臃纸鉃榫鶆騺砹骱头欠€(wěn)態(tài)擾動兩部分：

可滲透面f=0 的運動速度為v，且可滲透面單位外法線矢量n=?f。那么，可滲透面外f＞0區(qū)域的連續(xù)方程和動量方程可寫為：

其 中，H(·)和 δ(·)分別表示Heaviside 函數(shù)和Dirac delta 函數(shù)，c∞是均勻運動介質(zhì)中的聲速，物質(zhì)導(dǎo)數(shù)符號 D∞/Dt、Lighthill 應(yīng)力張量Tij以及參數(shù)Li和Q分別定義為：

式（5）和式（6）中的 σij和 δij分別表示黏性應(yīng)力張量和Kronecker delta 函數(shù)。

1.1 對流矢量波動方程

從方程（2）和方程（3）出發(fā)，文獻(xiàn)[23]得到了一個聲學(xué)矢量波動方程：

其中t*為積分變量。對線性小振幅波動，聲學(xué)密度ρ′遠(yuǎn) 小于均勻平均來流密度 ρ∞，因此在可滲透面外存在近似關(guān)系 ρu′≈ρ∞u′。Mao 等[23]認(rèn)為方程（8）右端最后一項的積分核函數(shù)為有旋矢量，對聲學(xué)速度場沒有實際貢獻(xiàn)，因而將其忽略，從而將聲學(xué)控制方程最終表述為：

方程（9）的左端是以聲學(xué)速度為變量的標(biāo)準(zhǔn)對流波動算子，右端七項為聲源項，其中第一項為單極子厚度源，緊接著的三項為偶極子載荷源，最后三項為四極子體積源。

1.2 聲學(xué)速度公式V1A-M

方程（9）的求解可采用對流格林函數(shù)，以便考慮均勻流的對流作用。假設(shè)均勻流亞聲速運動，時域?qū)α鞲窳趾瘮?shù)表達(dá)式為[2-3]：

其中，g為延遲時間函數(shù)，x和t表示觀察點的位置與時間，y和 τ表示聲源的位置與時間，聲學(xué)半徑 ?和R定義為：

式中，M∞=u∞/c∞為均勻來流馬赫數(shù)，α為均勻來流的對流放大系數(shù)，r表示聲源y與觀察點x之間的距離，即：

如果誘發(fā)聲輻射的主要氣動聲源位于可滲透面以內(nèi)，忽略體積源貢獻(xiàn)不會引起明顯數(shù)值誤差，從而可滲透面外的聲場可僅根據(jù)積分面上的厚度源和載荷源來預(yù)測，以提高數(shù)值計算效率。根據(jù)Farassat 方法[3,5]，文獻(xiàn)[23]提出了一個聲學(xué)速度半解析公式，其中來自厚度源貢獻(xiàn)的厚度噪聲u′T為：

MR代表聲源向觀察者方向的運動馬赫數(shù)：

來自載荷源貢獻(xiàn)的載荷噪聲u′L為：

a1～a4用于簡化積分公式書寫，具體表述如下：

方程（14）和方程（18）需要對面積分結(jié)果求觀察點的時間導(dǎo)數(shù)，將其轉(zhuǎn)化為積分符號內(nèi)對聲源的時間導(dǎo)數(shù)可以提高計算精度。此外，對可滲透面固定不動的特殊情況，如CAA/CFD 計算或風(fēng)洞聲學(xué)問題，參數(shù) ?和R、可滲透面外法向矢量n和 聲源運動速度v（等于零）是恒定的，不依賴延遲時間，且延遲時間可以采用顯式方法得到。對靜止可滲透面，本文將方程（14）和方程（18）改寫為：

2 聲學(xué)速度解析公式CV1A-M

文獻(xiàn)[23]利用公式V1-M 對均勻流中的單極子和偶極子聲輻射進(jìn)行研究，結(jié)果顯示單極子聲學(xué)速度預(yù)測結(jié)果滿足線化歐拉方程，而偶極子聲學(xué)速度的數(shù)值解不滿足線化歐拉方程。鑒于方程（14）和方程（18）的推導(dǎo)是正確的，那么對流波動方程（9）的右端載荷源必然會丟失。利用對流格林函數(shù)對丟失的載荷面源積分，即可實現(xiàn)對公式V1A-M 的修正。

2.1 對流矢量波動方程源項分析

方程（9）忽略了方程（8）右端最后一項。令

那么根據(jù)矢量恒等關(guān)系，有：

其中矢量F=n×(ρu′)。顯然，對小振幅擾動有近似關(guān)系?×(ρu′)≈ρ∞?×u′，因此對純聲學(xué)問題，式（24）右端第一項近似為零，可被忽略；但第二項為對聲場有實際貢獻(xiàn)的載荷面源，當(dāng)平均流速度不為零時不能舍棄。利用式（24），將對流矢量波動方程（9）修正為：

2.2 聲學(xué)速度公式CV1A-M

SA中 載荷面源的貢獻(xiàn)用表示，積分結(jié)果可寫為：

其中Ei j=εijmFm，εi jm為置換算子。利用關(guān)系式：

方程（26）的積分結(jié)果為：

對靜止可滲透面，方程（30）可進(jìn)一步表達(dá)為：

3 點源算例數(shù)值驗證

點源聲輻射算例常用于氣動聲學(xué)理論公式的數(shù)值驗證。考慮自激頻率 ω=340 rad/s的簡諧單極子和偶極子點源聲輻射，密度 ρ∞=1.2 kg/m3的均勻流沿x1軸 正向運動，馬赫數(shù)為0.5，聲速c∞=340 m/s。為保持對聲學(xué)速度公式的關(guān)注，并避免與流動模擬準(zhǔn)確性有關(guān)的任何偏差，可滲透面上的瞬態(tài)流動參數(shù)，包括壓力、密度和速度，都由點源產(chǎn)生的流場精確解產(chǎn)生。聲學(xué)速度公式中的時間導(dǎo)數(shù)采用二階有限差分方法近似，觀察點聲學(xué)信號采用行進(jìn)時間算法[24]計算。根據(jù)文獻(xiàn)[17,23]，取t0=0。

3.1 均勻流中的靜止單極子

均勻流中的靜止單極子聲場可以用一個簡單的諧波速度勢函數(shù)來描述：

其中參數(shù) ?和R分別由式（11）和式（12）定義。點源誘導(dǎo)的聲學(xué)速度、壓力和密度場為：

單極子位于坐標(biāo)系原點，速度勢函數(shù)的幅值為A=1 m2/s。采用中心在坐標(biāo)系原點且半徑為0.5 m的球面作為可滲透積分面，并利用數(shù)量為7 840 的四邊形網(wǎng)格離散積分面，確保空間解析精度。在每個聲源周期T內(nèi)布置128 個時間步，保證差分算法的數(shù)值精度。36 個觀察點均勻布置在x1-x2平面內(nèi)半徑為10 m 的圓上，將基于x1軸測量到的觀察點與坐標(biāo)原點間的幾何角度定義為觀察角θ。

所有觀察點的時間步長與聲源時間步長保持一致，輸出一個聲源周期（128 個時間步）的聲學(xué)信號。將觀察角 θ=120°處聲學(xué)速度分量的時間歷程預(yù)測值與理論解進(jìn)行對比，結(jié)果如圖1 所示。公式V1A-M 的數(shù)值解接近理論值，但在波峰與波谷附近仍存在一定差異；公式CV1A-M 的數(shù)值解則與理論值一致。

圖1 θ=120°靜止單極子聲學(xué)速度時間歷程Fig.1 Acoustic velocity time histories of a stationary monopole at θ=120°

圖2 不同觀察點處靜止單極子聲學(xué)速度均方根值Fig.2 RMS value of the acoustic velocity of a stationary monopole at different observation points

3.2 均勻流中的靜止偶極子

考慮均勻流中靜止偶極子的聲輻射，假設(shè)偶極子的軸線與x2軸對齊，速度勢函數(shù)可表示為：

偶極子的位置、自激頻率、速度勢函數(shù)幅值、可滲透面網(wǎng)格，以及觀察點的設(shè)置，均與單極子點源算例相同。誘導(dǎo)的聲學(xué)速度、壓力和密度場理論解分別采用式（34）、式（35）和式（36）計算。

圖3 為1 20°觀察點聲學(xué)速度數(shù)值解隨時間的變化歷程及其與理論解的對比。由于可滲透面有效聲源部分丟失，公式V1A-M 預(yù)測的聲學(xué)速度分量和與理論解存在顯著差異；修正后的公式CV1A-M 能獲得與理論解一致的數(shù)值結(jié)果。改變平均流的速度（包括大小與方向）、點源的自激頻率或可滲透面的大小，公式CV1A-M 均能得到與理論解完全吻合的結(jié)果（這里不再給出）。圖4 為不同觀察點處偶極子聲學(xué)速度分量均方根值數(shù)值解與理論解的對比。公式CV1A-M 數(shù)值解再一次與理論值保持一致，而公式V1A-M 數(shù)值解在點源上游與理論解存在巨大差異。公式V1A-M 丟失的源項為載荷面源，因而偶極子點源的數(shù)值誤差遠(yuǎn)大于單極子算例的數(shù)值誤差。

圖3 θ=120°靜止偶極子聲學(xué)速度時間歷程Fig.3 Acoustic velocity time histories of a stationary dipole at θ=120°

圖4 不同觀察點處靜止偶極子聲學(xué)速度均方根值Fig.4 RMS value of the acoustic velocity of a stationary dipole at different observation points

3.3 均勻流中的旋轉(zhuǎn)單極子

對運動點源算例，假設(shè)初始方位角為 0°的點源以恒定角速度 Ω=ω/4繞x3軸逆時針旋轉(zhuǎn)，旋轉(zhuǎn)半徑d=0.25 m，如圖5 所示。點源自激頻率、速度勢函數(shù)幅值、可滲透面網(wǎng)格和觀察點位置均與靜止點源算例相同，輸出512 個時間步（對應(yīng)4 個聲源周期）的聲學(xué)速度信號。

圖5 均勻平均流中旋轉(zhuǎn)點源示意圖Fig.5 The schematic of a rotating point source in a uniform mean flow

運動單極子點源的速度勢函數(shù)為：

其中MR為延遲時刻點源向觀察點方向運動的馬赫數(shù)，由式（17）定義。圖6 為1 50°觀察點聲學(xué)速度時間歷程的數(shù)值解與理論解對比。公式V1A-M 預(yù)測的聲學(xué)速度分量與 理論值基本吻合，而與理論值存在大小和相位上的明顯差異。公式CV1A-M 的數(shù)值解與理論解的一致性進(jìn)一步驗證了其正確性。

圖6 θ=150°旋轉(zhuǎn)單極子聲學(xué)速度時間歷程Fig.6 Acoustic velocity time histories of a rotating monopole at θ=150°

圖7 對比了不同觀察點處旋轉(zhuǎn)單極子聲學(xué)速度分量均方根值的數(shù)值解與理論解。公式CV1A-M 的數(shù)值解與理論值是一致的，公式V1A-M 的數(shù)值解在點源上游方向，尤其是1 80°觀察角附近，與理論值差異明顯。

圖7 不同觀察點處旋轉(zhuǎn)單極子聲學(xué)速度均方根值Fig.7 RMS value of the acoustic velocity of a rotating monopole at different observation points

3.4 均勻流中的旋轉(zhuǎn)偶極子

仍然假設(shè)運動偶極子點源的軸線與x2軸對齊，速度勢函數(shù)可根據(jù)旋轉(zhuǎn)單極子的勢函數(shù)得到：

圖8 為1 50°觀察點旋轉(zhuǎn)偶極子聲學(xué)速度時間歷程的數(shù)值解與理論解對比。公式CV1A-M 的數(shù)值解與理論解一致；聲學(xué)速度分量的公式V1A-M 數(shù)值解雖然大小與理論值不吻合，但相位相同，而分量的公式V1A-M 數(shù)值解大小和相位與理論值均存在巨大差異。圖9 比較了不同觀察點處旋轉(zhuǎn)偶極子聲學(xué)速度分量均方根值的數(shù)值解與理論解。公式CV1A-M數(shù)值解與理論值一致，公式V1A-M 數(shù)值解在點源上游與理論解差異巨大。這是意料之中的結(jié)果。

圖8 θ=150°旋轉(zhuǎn)偶極子聲學(xué)速度時間歷程Fig.8 Acoustic velocity time histories of a rotating dipole at θ=150°

圖9 不同觀察點處旋轉(zhuǎn)偶極子聲學(xué)速度均方根值Fig.9 RMS value of the acoustic velocity of a rotating dipole at different observation points

均勻流中點源聲輻射誘導(dǎo)的聲壓計算公式（35）是根據(jù)線化歐拉方程和聲學(xué)速度計算公式（34）得到的，因此對純聲學(xué)問題，聲學(xué)速度和聲學(xué)壓力滿足線化歐拉方程。在上述4 個典型算例中，聲學(xué)速度公式CV1A-M 的數(shù)值解均與理論值一致，證明了該公式的理論正確性。公式V1A-M 不能實現(xiàn)均勻流中偶極子聲學(xué)速度的預(yù)測，并非線化歐拉方程不能描述偶極子聲學(xué)速度，而是可滲透載荷源部分丟失引起的。

4 結(jié)論

文獻(xiàn)[23]提出的對流矢量波動方程丟失了部分對聲場有貢獻(xiàn)的可滲透面載荷源，使得公式V1-M 和V1A-M 難以準(zhǔn)確外推聲學(xué)速度，即便對單極子算例，數(shù)值解與理論值仍存在明顯誤差；受均勻流對流作用影響，點源上游方向的公式V1A-M 數(shù)值解誤差遠(yuǎn)大于下游方向誤差。

通過考慮被丟失的可滲透載荷源貢獻(xiàn)，發(fā)展了積分解析公式CV1A-M。對均勻流中靜止或旋轉(zhuǎn)單極子和偶極子而言，公式CV1A-M 的數(shù)值解與理論值一致。均勻流中的偶極子點源聲輻射不會誘發(fā)渦波擾動，線化歐拉方程可以準(zhǔn)確描述偶極子誘導(dǎo)的聲學(xué)速度與聲學(xué)壓力間的關(guān)系。

對真實的流動誘發(fā)氣動噪聲問題，可滲透面往往穿過一定的非線性區(qū)域。盡管主要的氣動聲源被可滲透面包圍，但當(dāng)渦波通過可滲透面時，渦波擾動會被收集為積分面的輸入而產(chǎn)生偽聲。聲學(xué)速度公式CV1A-M 的實際應(yīng)用還需發(fā)展相應(yīng)的偽聲抑制方法，這是未來需要關(guān)注的研究方向。