可滲透面對流FW-H 方程偽聲產生機制和抑制方法

何嘉華，王垿桁，劉秋洪，錢振昊

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

0 引言

聲比擬方法是工程中應用最廣泛的氣動噪聲數值預測方法[1-2]，先對有限區域內流動進行高保真數值計算以精確模擬噪聲的產生，然后假設聲波從近場到遠場線性運動，通過求解基于Lighthill 聲比擬思想的標準波動方程或標準對流波動方程，比如FW-H 方程[3]和對流FW-H 方程[4,5]，來快速預測噪聲的傳播。在實際應用中，FW-H 方程和對流FW-H 方程多采用非物理的封閉可滲透面作為積分面，完整解由面積分和體積分組成。如果可滲透面足夠大，聲波的產生、反射和折射等過程都包圍其中，即可忽略體積分對聲場的貢獻[6]，從而積分面外的聲場僅利用可滲透面上的厚度源和載荷源信息計算，比如解析公式F1A[7]和F1C[4]。

對繞流和射流等實際噪聲問題，為了保證可滲透面上氣動聲源具有足夠精度，選擇的可滲透面往往穿過一定非線性區域。近場流動的非線性擾動同時包含了聲波和渦波，渦波衰減緩慢且沿平均流向下游傳播，當通過可滲透面時，渦波擾動不可避免地被收集為積分面的輸入而產生偽聲[8]。學術界對偽聲產生機制的理解存在可滲透面急劇截斷積分區域導致譜泄漏和聲比擬方程求解方法破壞源項濾波功能的分歧[9-12]。Morfey 等[9]研究發現，考慮可滲透面外體積源項的貢獻即可消除偽聲污染。Rahier 等[10]通過理論分析和數值試驗明確指出，偽聲污染就是忽略四極子體積源項貢獻引起的。Obrist 等[11]的研究顯示，可滲透面急劇截斷積分區域會導致面積分源項頻譜泄漏，從而產生污染聲場的偽聲。Mao 等[12]細致分析了對流FW-H 方程的右端面源，提出了與Obrist 等不同的觀點，認為載荷源和厚度源能自動過濾自身包含的渦波偽源，但數學演繹過程即便對靜止可滲透面也不嚴格成立。

遠場噪聲的準確預測需要采用有效方法抑制偽聲干擾。最簡單的方法是直接舍棄可滲透面中部分允許渦波通過的片區[13-14]，通過對非封閉可滲透面積分來預測噪聲。然而FW-H 方程的推導建立在封閉可滲透面基礎上，因此使用非封閉積分面的方法在數學上是不正確的。Spalart 等[15]采用多個相鄰可滲透面來提取聲源，最終聲場是這些可滲透面積分結果的平均值。Mendez 等[16]給出了估計可滲透面數量和位置的表達式。Obrist 等[11]建議使用平滑變化的空間窗函數來減少截斷誤差。Wright 等[17]提出了一個廣義FW-H 方程，利用不同權重的可滲透面數據或有限厚度的體積源來降低偽聲強度。利用多可滲透面數據或有限厚度體積源貢獻抑制偽聲方法會導致聲源存儲量和聲學計算量的增加。Ikeda 等[18]將無窮薄厚度體積源貢獻轉化為面積分，提出了適用于聲學遠場的附加面源偽聲抑制方法。Zhong 等[8]從線化Euler方程出發，發展了一種考慮可滲透面局部非均勻流影響的聲遠場外推方法，其核心思想是利用擾動速度的散度 ?·u′和物質導數 D∞u′/Dτ自動過濾渦波，但增加了積分面源項的計算困難。Mao 等[12]建議采用Kirchhoff 積分方程預測噪聲傳播以抑制偽聲污染，但Kirchhoff 積分面僅適用于線性區域。

明晰渦波偽聲產生機制可為偽聲干擾抑制方法的發展指明方向。本文以運動介質中的聲傳播為例，通過耦合分析對流FW-H 方程的載荷源和厚度源，闡明渦波偽聲的產生機制。在此基礎上提出一種偽聲干擾抑制方法，并利用算例驗證其正確性。

1 可滲透面渦波偽聲產生機制

假設無窮遠均勻來流的密度、壓力和速度分別為ρ∞、p∞和u∞i（對應馬赫數為Ma∞i），將當地流動分解為均勻來流和非穩態擾動兩部分：

引入運動速度為vi（對應運動馬赫數為Mai）的可滲透面f=0，可滲透面內部和外部區域分別定義為f＜0和f＞0，且可滲透面單位外法線矢量滿足n=?f。對流FW-H 方程可表示為[4-5]：

方程左端c∞和H(·)分別表示均勻運動介質中的聲速和Heaviside 函數；方程右端 δ(·)表示Dirac delta 函數；物質導數 D∞/Dt定義為：

右端3 項分別代表四極子體積源、偶極子載荷源和單極子厚度源，張量Tij、Li和Q定義為：

式中，δij為 Kronecker delta 函數，σij為黏性應力張量，下標i,j=1,2,3表示直角坐標方向（對二維問題i,j=1,2）。

1.1 對流FW-H 方程源項的渦波自動過濾功能

可滲透面的選擇應遵循一定規則，即在瞬態流動參數精確解析前提下，可滲透面應覆蓋聲波的產生、反射和折射等主要過程，并使得可滲透面處的平均流近似均勻。如果渦在流動中快速發展或衰減，或彼此之間有強相互作用，就會產生不可忽視的聲輻射，可滲透面應包圍這些區域。也就是在一個恰當的可滲透面處，渦波衰減緩慢，且不會產生明顯的壓力和密度擾動。由于主要的氣動聲源位于可滲透面內，忽略體積源貢獻不會引起明顯的數值誤差，因此公式F1A[7]和F1C[4]僅根據積分面上的厚度源和載荷源評估可滲透面外聲場，以提高計算效率。

對高雷諾數流動，黏性應力張量因對遠場的貢獻小而被忽略，從而載荷源源參數Li可近似為：

非加熱流體的運動速度可以分解為無旋分量和有旋分量兩部分，其中擾動速度的無旋分量和有旋分量分別對應于聲波和渦波分量[19]，也就是：

其中下標a和v分別表示聲波和渦波分量。利用式（1）和式（8）可將式（6）等價表述為：

其中，右端前兩項分別代表一階的聲波擾動和渦波擾動，后兩項分別表示二階的聲-聲干擾和渦-渦干擾。同理，式（7）可等價寫為：

其中，右端前兩項分別代表一階的聲波擾動和渦波擾動，緊接著的3 項分別表示二階的聲-聲干擾、渦-渦干擾和渦-聲干擾，最后兩項分別為三階的聲-聲干擾和渦-渦干擾。

理論上講，式（9）和式（10）中非線性的渦-渦干擾和渦-聲干擾能夠誘導聲輻射，應該視為聲源的一部分。但在通常情況下，可滲透面處的聲波和渦波擾動速度遠小于平均流速度，因而可滲透積分面上的厚度源和載荷源強度主要取決于一階擾動項。已有的研究表明，均勻流中的一階渦波擾動項僅表征對流效應[20]，不能產生聲波擾動，即不會誘發聲輻射，屬于偽聲源。定義：

當可滲透面外法線方向與均勻平均流方向并非完全一致時，D∞[?H(f)]/Dt非零，因此方程（16）并不嚴格成立。

方程（14）減去方程（15）即可得到對流FW-H 方程（2）右端的偽聲源之和為：

方程（17）并沒有限定可滲透面靜止，可見對流FWH 方程（2）的面源項確實能夠自動過濾渦波偽聲源。

忽略黏性和非等熵效應，應力張量Tij可展開為：

由于非線性渦-渦干擾和渦-聲干擾屬于聲源，且應力張量Tij不存在一階渦波擾動，故體積源項不含有偽聲源。忽略體積源對聲場的貢獻會導致一定數值誤差，但不會產生偽聲輻射[20]。因此，可滲透面急劇截斷空間積分區域不是偽聲產生的根源。

1.2 對流FW-H 方程求解過程的偽聲引入

對流FW-H 方程的求解可以采用對流格林函數，以便將均勻流對流作用包含在格林函數中。假設均勻流亞聲速運動，三維時域格林函數可表示為[4-5]：

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

式中，R 表示歐式空間，G為對流格林函數。根據方程（17）可知，對源項直接求物質導數和散度可自動過濾渦波，因此依據方程（22）積分得到的結果能夠消除渦波干擾，但需要數值計算右端項的空間導數，其過程復雜且易導致不必要的數值誤差[21]。令Qa=Q?Qv和La=L?Lv，利用分部積分得到：

根據Farassat 提出的積分方法[7]，方程（23）右端前兩項的積分結果均為零，因而被忽略。利用格林函數性質D∞G/Dτ=?D∞G/Dt與?G/?yi=??G/?xi，并將聲場分解為厚度噪聲和載荷噪聲兩部分，有：

其中下標T和L則分別表示厚度噪聲和載荷噪聲分量。方程（24）和方程（25）中，利用了函數 δ(f)的性質將體積分降階為f=0處的面積分。厚度噪聲和載荷噪聲積分公式的詳細推導過程見參考文獻[4-5]。

由于Qa和La與一階渦波擾動無關，方程（24）和方程（25）的積分結果不含渦波偽聲。然而，對數值模擬得到的擾動速度場u′，一般難以實現渦波速度分量的分離，因此聲壓積分公式（比如公式F1C[4]）直接采用變量Q和L計算厚度噪聲和載荷噪聲：

當有渦波通過可滲透面時，渦波擾動不可避免地被收集為積分面的輸入，從而導致偽聲輻射。偽聲的厚度分量和載荷分量可表示為：

大量的研究結果[4-5]顯示，方程（26）和方程（27）對聲學問題是成立的，即忽略方程（23）右端前兩項積分貢獻是正確的?？蓾B透面處存在渦波擾動時，如果聲傳播仍然采用方程（26）和方程（27）計算，則下述積分項貢獻也被忽略了：

根據方程（17），可滲透面處的渦波偽聲源是可以自動過濾的，那么有：

也就是說，被方程（26）和方程（27）忽略的積分項貢獻可以抵消渦波偽聲。

通過上述分析可以知道，盡管載荷源和厚度源能夠自動過濾渦波，但按照Farassat 方法，波動方程求解采用了分部積分，最終的聲傳播積分方程忽略了部分積分項的貢獻，因難以實現渦波速度分量的分離，使得積分結果引入了渦波偽聲輻射。

2 可滲透面偽聲抑制方法

Farassat 積分方法對僅含聲波分量的源項積分有效，而對含有渦波分量的源項積分失效。將可滲透面源項拆分為僅含聲波分量項和可含渦波分量項兩部分。對流FW-H 方程（2）可等價變形為：

2.1 時域積分公式

可滲透面外聲場的積分解可表示為：

將三維對流格林函數式（19）代入上述兩個方程，得到厚度噪聲和載荷噪聲的半解析積分公式為：

對可滲透面固定不動的特殊問題，如CAA/CFD計算或風洞聲學問題，參數 ?與R、可滲透面外法向矢量n和聲源運動速度v（等于零）不依賴于延遲時間，且延遲時間可采用顯式方法獲取。方程（38）和方程（39）可進一步寫為：

2.2 頻域積分公式

若可滲透面周期性運動，如勻速平移、恒速旋轉和靜止不動，噪聲預測可在頻域下實施。對時域信號ξ(t)和頻域信號(ω)，若Fourier 變換對為：

其中，ω為角頻率，j為虛數單位。那么三維頻域對流格林函數為：

其中k=ω/c∞為聲學波數。對方程（38）和方程（39）兩端進行Fourier 變換，得到聲壓頻域積分解為：

當可滲透面靜止時，可進一步簡化為：

將頻域聲壓信號執行Fourier 逆變換，即可得到聲壓的時間歷程。如果聲學問題是二維的，則采用二維頻域對流格林函數[22]進行積分，此處不再詳述。

上述渦波偽聲抑制方法的核心思想與文獻[8]是一致的，即利用了擾動速度的散度 ?·u′和物質導數D∞u′/Dτ能自動過濾渦波的特性。需要進一步指出的是，前文分析過程默認可滲透面處渦波不會產生壓力和密度擾動。恰如文獻[8]所強調的，若湍流中存在對流的渦波壓力或密度擾動，當其通過可滲透積分面時，預計仍會存在偽聲輻射。更具體地說，因參數Q′和L′中的壓力和密度擾動含有非聲學的成分，必將對偽聲的產生做出貢獻，此時應重新選擇可滲透面，盡量避免或弱化非聲學的壓力或密度擾動。

3 均勻流算例數值驗證

利用具有理論解的流動算例驗證渦波偽聲的產生機制和抑制方法。假設密度 ρ∞=1.216 kg/m3的均勻流沿x1軸正向運動，馬赫數Ma∞=0.5。首先考慮自激角頻率 ω=680 rad/s的單極子聲輻射，證明算法對純聲學問題的有效性；然后利用均勻流中的對流渦算例檢驗算法的渦波抑制效果；最后將單極子與對流渦疊加在一起，數值驗證渦波偽聲的產生機制。所有算例中，取聲速c∞=340 m/s。

3.1 均勻流中的單極子

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

其中，A為速度勢函數的幅值，參數 ?和R分別由式（20）和式（21）定義。點源誘導的聲學速度、壓力和密度場為：

單極子位于坐標原點，選擇邊長a=4.0 m的正方體表面作為可滲透面包圍點源（如圖1 所示），采用均勻的結構化網格離散可滲透面，網格總數38 400?？蓾B透面處瞬態流動參數（包括壓力、密度和速度）由流場精確解產生。對線性聲學問題，忽略源項中的非線性作用。36 個觀察點均勻布置在x3=0平面內一個半徑為r=34 m的圓上，將基于x1軸測量到的觀察點與坐標原點間的幾何角度定義為觀察角 θ。取A=0.1 m2/s，單極子誘導的瞬時聲學壓力和聲學速度場在x3=0平面內的空間分布分別如圖2 和圖3 所示。因平均流對流作用，點源上游的聲波被壓縮，而下游的聲波則被拉伸。

圖1 單極子、可滲透面和觀察點示意圖Fig.1 Schematic of monopole,permeable surface,and observers

圖2 靜止單極子誘導的聲學壓力場Fig.2 Pressure field induced by a stationary monopole

圖3 靜止單極子誘導的聲學速度場Fig.3 Velocity field induced by a stationary monopole

時域聲壓預測公式中的時間導數和空間導數采用二階有限差分方法近似，在每個聲源周期T內布置128 個時間步，空間步長則取為5 mm（與可滲透面的網格間距不一致），以保證差分算法的數值精度。觀察點聲學信號采用超前時間算法[23]計算。將含有偽聲抑制算法的聲場數值預測結果稱為CFW-H 解。圖4（a）和圖4（b）分別為點源下游 θ=30?和上游θ=150?處的聲壓時間歷程預測結果，CFW-H 數值解與理論解、FW-H 積分方程解都吻合。對僅含聲模態的噪聲問題，應盡量采用經典方法去外推聲傳播，比如對流FW-H 積分方程，以避免空間導數的計算。

圖4 靜止單極子聲學壓力時間歷程Fig.4 Pressure time histories induced by a stationary monopole

3.2 均勻流中的對流渦

均勻流中的對流渦算例常用于檢驗渦波偽聲抑制方法。對流渦隨平均流運動，假設在t=0時刻通過(0,0)點，其誘導的速度場為：

仍采用圖1 所示的正方體6 個表面（S1～ S6）作為可滲透面。實際物理問題中，渦波隨平均流向下游運動，僅在 S4面上存在渦波干擾。因此本算例忽略S4之外的5 個面上的渦波速度場，而且不考慮對流渦誘導的壓力和密度擾動。對小振幅渦波擾動，不計入偽聲源非線性項的貢獻，分別采用FW-H 積分方程和偽聲抑制方法預測遠場噪聲。在點源下游 θ=30?和上游 θ=150?處，對流渦誘導的偽聲時間歷程數值預測結果分別如圖6（a）和圖6（b）所示。均勻流中的對流渦不會誘發能輻射到遠場的聲波。FW-H 數值解在對流渦通過 S4面前后會產生偽聲輻射，噪聲預測結果與理論值差異顯著，而偽聲抑制方法能有效過濾渦波偽聲源，CFW-H 數值解與理論值保持一致。

圖6 對流渦偽聲時間歷程Fig.6 Time histories of spurious noise induced by a convective vortex

3.3 均勻流中的單極子與對流渦

考慮均勻流中同時存在靜止單極子和恒速運動對流渦的情況。將單極子誘導的壓力和密度，以及單極子和對流渦誘導的速度擾動和，作為可滲透面的流動參數（假設僅在 S4面存在渦波擾動），將其輸入程序代碼預測噪聲傳播。因對流渦誘導的壓力場不會傳播到遠場，觀察點的聲場理論上僅含單極子貢獻。不考慮聲波與渦波間相互作用，即忽略厚度源和載荷源中的非線性項，僅考慮一階線性擾動項對聲學遠場的貢獻。圖7（a）和圖7（b）分別為 θ=30?和 θ=150?處聲場預測結果。FW-H 數值解在渦波通過 S4面時產生了顯著的偽聲輻射，且受均勻流的對流效應影響，點源下游的偽聲污染比上游更加強烈。CFW-H 數值解與理論解很好地吻合，證明了渦波偽聲產生機制的正確性和抑制方法的有效性。

圖7 靜止單極子和對流渦誘導的聲壓時間歷程Fig.7 Time histories of acoustic pressure induced by a stationary monopole and a convective vortex

需要指出的是，如果可滲透面處存在強烈的渦-聲和渦-渦非線性相互作用，CFW-H 數值解仍將產生明顯的數值誤差。不妨考慮 S4面上因對流渦運動誘導的壓力擾動和密度擾動，解析表達式為：

圖8 計入對流渦壓力和密度擾動影響的聲壓時間歷程Fig.8 Time histories of acoustic pressure time accounting for the pressure and density fluctuations induced by the convective vortex

4 結論

可滲透面對流FW-H 方程右端源項包含了單極子厚度源、偶極子載荷源和四極子體積源。雖然緩慢衰減的渦波穿過可滲透面后沿平均流向下游傳播，但由于體積源項不含有偽聲源，只要主要聲源被可滲透面包圍，忽略體積源的遠場貢獻僅會產生數值誤差，而不會誘發渦波偽聲輻射，因此可滲透面急劇截斷空間積分區域不是產生渦波偽聲的根源。盡管厚度源的物質導數和載荷源的散度操作不能分別過濾偽厚度源和偽載荷源，但兩者的耦合作用可自動過濾偽聲源。偽聲產生的真正原因在于對流FW-H 方程的求解采用了Farassat 提出的分部積分公式，部分積分項被忽略，從而導致偽聲源自動過濾功能失效。如果直接利用流動數值模擬所得擾動速度計算積分面厚度源參數Q和載荷源參數L，對流FW-H 方程積分公式就引入了渦波偽聲輻射。

將對流FW-H 方程進行恒等變形，從厚度源和載荷源中抽取含有渦波擾動的項，在對流波動方程求解過程中保留相應的物質導數和散度操作，即可過濾渦波偽聲源，從而實現渦波偽聲傳播的有效抑制。根據這一思想，發展了修正的時域和頻域聲學積分方程。對均勻流中的單極子和對流渦算例，利用時域積分解討論了渦波偽聲抑制方法的可靠性。數值結果顯示，FW-H 積分方程數值解在渦波通過可滲透面時間產生了顯著的偽聲輻射，且受均勻流對流效應影響，點源下游的偽聲污染比上游更強烈；修正后的數值解與理論解相吻合，驗證了偽聲產生機制分析的正確性。

本文關注渦波偽聲的產生機制，但僅研究了具有理論解的簡單算例，后續將利用三維復雜流動算例開展進一步的時域與頻域驗證工作。此外，本文偽聲抑制方法也適用于近場預測，只是文中沒有給出相應的驗證結果。