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

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

(西北工業(yè)大學(xué) 翼型葉柵空氣動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，西安 710072)

0 引言

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

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

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

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

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

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

假設(shè)無(wú)窮遠(yuǎn)均勻來(lái)流的密度、壓力和速度矢量分別為 ρ∞、p∞和u∞，將當(dāng)?shù)亓鲃?dòng)分解為均勻來(lái)流和非穩(wěn)態(tài)擾動(dòng)兩部分：

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

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

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

1.1 對(duì)流矢量波動(dòng)方程

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

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

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

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

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

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

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

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

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

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

a1～a4用于簡(jiǎn)化積分公式書(shū)寫(xiě)，具體表述如下：

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

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

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

2.1 對(duì)流矢量波動(dòng)方程源項(xiàng)分析

方程（9）忽略了方程（8）右端最后一項(xiàng)。令

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

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

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

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

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

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

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

3 點(diǎn)源算例數(shù)值驗(yàn)證

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

3.1 均勻流中的靜止單極子

均勻流中的靜止單極子聲場(chǎng)可以用一個(gè)簡(jiǎn)單的諧波速度勢(shì)函數(shù)來(lái)描述：

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

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

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

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

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

3.2 均勻流中的靜止偶極子

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

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

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

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

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

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

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

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

運(yùn)動(dòng)單極子點(diǎn)源的速度勢(shì)函數(shù)為：

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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