均勻流中聲學(xué)速度預(yù)測的時域解析公式

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

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

0 引言

經(jīng)典聲比擬理論建立了聲場與不同類型聲源之間的響應(yīng)關(guān)系，成為工程中應(yīng)用最廣泛的氣動噪聲數(shù)值預(yù)測方法。聲比擬方法先采用高精度數(shù)值方法計算有限域內(nèi)氣動聲源信息，然后利用聲比擬積分方程模擬噪聲傳播。FW-H 方程[1]及其擴展形式[2,3]采用封閉可滲透面作為積分面，完整解由面積分和體積分兩部分組成。當(dāng)主要氣動聲源包圍在可滲透面內(nèi)時，可忽略體積分對聲場的貢獻[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)通過單位面積的噪聲能量通量，其方向就是噪聲能量的輸運方向。聲強場的可視化可以詳細顯示噪聲能量的傳播路徑，有助于揭示聲能量輻射機制，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é)表達上都很復(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 等進一步利用這些聲學(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è)無窮遠均勻來流的密度、壓力和速度矢量分別為 ρ∞、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ā)，文獻[23]得到了一個聲學(xué)矢量波動方程：

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

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

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

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

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

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

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

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

來自載荷源貢獻的載荷噪聲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

文獻[23]利用公式V1-M 對均勻流中的單極子和偶極子聲輻射進行研究，結(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）右端第一項近似為零，可被忽略；但第二項為對聲場有實際貢獻的載荷面源，當(dāng)平均流速度不為零時不能舍棄。利用式（24），將對流矢量波動方程（9）修正為：

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

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

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

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

對靜止可滲透面，方程（30）可進一步表達為：

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

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

圖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é)論

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

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

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