馬瑞軒 王益民 張樹海 武從海 王勛年
1) (中國空氣動力研究與發展中心, 空氣動力學國家重點實驗室, 綿陽 621000)
2) (中國空氣動力研究與發展中心, 氣動噪聲控制重點實驗室, 綿陽 621000)
以聲波為主要表現形式的膨脹過程和以旋渦為主要表現形式的剪切過程之間的非線性耦合問題一直以來都是流體力學的研究熱點.尤其是旋渦對聲波的散射問題, 具有重要的科學意義與工程應用背景.本文通過線性緊致格式直接數值求解二維歐拉方程, 獲得了平面聲波穿過均熵Taylor渦的散射特性.與之前經典文獻中的標準算例比較, 結果極其吻合, 直接驗證了研究所采用的高精度高分辨率空間差分和時間推進格式以及遠場無反射邊界條件(緩沖區)的計算方法在時域同時解析動力學量和聲學量(量級遠遠小于動力學量)的有效性.通過引入散射截面, 將全區域的散射分為長波近似區、共振散射區和幾何聲學區.針對每個子區域,重點分析了無量綱尺度量旋渦強度和長度尺度比對散射聲場的影響, 給出了散射聲場關于上述兩個關鍵無量綱參數的尺度律關系, 并且得到了極低馬赫數極大波長時散射聲場的分布函數.在此基礎上給出了關于旋渦聲散射物理機制的一種解釋.
旋渦流動中的聲傳播問題作為氣動聲學的經典問題, 不僅可以直接應用于聲目標識別與探測,旋渦主動降噪等實際工程問題[1?5], 而且對于發展氣動聲學計算方法[6]以及認識復雜流動(如剪切層、湍流)與聲波相互作用及其發聲機理有著重要意義[7?12], 極具應用與科研價值.
聲波穿過旋渦流動時會產生強烈的非線性散射, 頻率、幅值和相位會發生顯著的變化, 而且聲波的傳播也會影響流場本身隨時間發展的運動規律.旋渦流動對聲波的散射問題有著悠久的研究歷史, 其中聲波與湍流的相互作用問題的研究歷史甚至超過氣動聲學本身.早在1941年, 前蘇聯學者Obukhov[13]就開展了湍流對聲散射的研究.Kraichnan[10]和Lighthill[11]針對湍流聲散射問題采用Lighthill聲比擬方法發展了最早的理論預估模型.Howe[14]發展了用來描述包含多尺度連續旋渦對聲波散射的動力學方程.Clifford和Brown[15]基于含源項的亥姆霍茲方程發展了描述散射聲場強度的理論模型, 并且發現流場的平均輸運速度會帶來散射聲波的多普勒頻移現象.采用拋物型的近似方程, Ostashev等[16]開展了理論研究, Dallois等[17]開展了數值研究.
目前絕大多數針對湍流聲散射的研究主要采用半經驗的理論預估模型和數值計算方法, 缺乏描述湍流的精確物理模型.而且研究中發現對聲散射起決定作用的流動結構是一系列具有特定尺度的旋渦結構, 因此從單個旋渦出發研究聲散射具有重要的意義.在早期, Fetter[18]從線性歐拉方程,Oshea[19]采用聲比擬理論都獲得了經典點渦對聲波的散射特性.Colonius[20,21]等選取多個基本旋渦流動結構, 分別采用直接數值模擬和理論研究的方法開展了聲散射研究.其中針對Rankine渦, 基于聲比擬理論分析了低馬赫數下聲散射特性, 采用聲線法獲取了高頻近似解.Ford等[22?24]和Hattori等[25]采用匹配漸近展開的方法研究了球對稱旋渦聲散射的特性, 并且采用直接數值計算的方法驗證了其理論結果.Howe[26], Kopiev和Belyaev[27]分別采用渦聲理論和諧波展開的方法研究了Rankine渦的聲散射特性.
現有的數學理論方法直接嚴格求解旋渦聲散射的控制方程時會遇到一系列技術上的難題, 發展理論模型時必須要進行數學上的近似, 比如漸近展開時要求旋渦流動是低馬赫數, 玻恩近似時要求聲波波長遠遠大于渦核半徑, 幾何聲學理論中的射線追蹤法只適用于高頻聲散射.但是在實際推導過程中, 大多時候這些近似所需要的數學條件無法同時滿足, 而且近似條件也大大地限制了結果的適用范圍.目前還不存在普遍適用且準確性高的理論模型, 絕大多數理論模型缺乏對散射聲場細節的解析, 比如相位信息的缺失.近年來, 隨著計算機軟硬件的提升, 開展精細化的數值模擬成為了一種常用的研究旋渦聲散射的手段.其中, Candel[28]通過數值求解近似的拋物型方程獲得了聲波穿過Rankine渦的畸變特性, Colonius等[20,21]通過直接數值求解Navier-Stokes方程獲得了典型緊致渦與非緊致渦的聲散射特性, Clair和Gabara[29]通過數值求解線性歐拉方程獲得了旋渦運動對散射聲場的影響特性.但是這些研究中, 大多是將數值計算作為驗證理論結果的手段, 沒有對計算結果系統分析進而給出散射聲場的定量描述.本文針對均熵Taylor渦, 通過高精度高分辨率方法直接數值求解二維歐拉方程, 獲得了散射聲場隨時間演化的全部信息.通過分析散射截面與旋渦強度和長度尺度比的尺度律關系, 將旋渦聲散射進行了分類.針對每一類型的散射, 重點分析了散射聲場關于上述兩個無量綱量的尺度效應.
首先, 給出本文所用旋渦模型的詳細數學表示.考慮二維圓對稱均熵流動, 在柱坐標 (r,θ) 下其法向和徑向速度表示如下:

這就是經典Taylor渦 (Taylor, 1918) 在無粘情況下的表述.其中, 速度分量和空間坐標分別采用無窮遠處的聲速c∞和渦核半徑R進行無量綱化處理.旋渦的馬赫數(旋渦強度)定義為

這里Uv是旋渦的最大速度.再由均熵關系:

和聲速關系:

可以得到密度和壓強的表達式:

其中密度ρ和壓強p分別采用無窮遠處密度ρ∞和無量綱化;γ表示絕熱指數, 空氣中一般為1.4.
選取無窮遠處沿x正方向傳播的平面波作為入射波, 具體表達式如下:

這里聲學量的密度ρ′;x方向速度分量u′;y方向速度分量v′和壓強p′采用和 旋 渦動力學量相同參考量進行無量綱化;ω為平面波的角頻率;ε=10?5?1, 以保證聲學量遠遠小于動力學量.
與之前的研究類似[25], 我們認為聲傳播的過程可以忽略粘性和熱傳導效應.因此在本文的研究中, 通過數值求解二維歐拉方程獲取平面聲波穿過旋渦后的散射特性.其守恒形式如下:

其中U是守恒變量;F和G分別是x,y方向的流通矢量, 具體形式為

ρ,u,v,p分別是密度、x和y方向的速度、壓強;總能量密度
為了更好地從背景流場(動力學量)中解析出聲學量, 空間離散采用經典的6階中心緊致格式(邊界處為3階精度), 并且采用4階龍格庫塔格式進行時間推進.為了保持計算的穩定性, 采用8階中心緊致濾波技術抑制計算中產生的非物理高頻振蕩.
計算在直角坐標系(x,y)下進行, 矩形計算域如圖1所示.上節給出的均熵Taylor渦是非定常歐拉方程的精確解, 在計算中保持定常狀態, 可作為理想的背景流動.旋渦中心位于計算域的幾何中心.本文只考慮順時針旋轉的旋渦, 逆時針對應的結果只需沿x軸對稱變換即可.左邊界為入口邊界, 按照(8)式給定沿x軸正方向傳播的平面聲波.右邊界為出口邊界, 為了減少邊界反射對計算的影響, 在右邊界的外設置緩沖區(buffer zone)[30,31].在上下邊界處, 我們將物理量分成定常部分和非定常部分, 定常部分按照(1)式, (2)式, (6)式和(7)式給定, 非定常部分采取高階單邊插值的方法處理.

圖1 計算域示意圖Fig.1.Schematic diagram of computation configuration.
計算中采用等間距的笛卡爾網格, 網格間距滿足能同時解析旋渦動力學量和聲學小擾動量的要求, 取 ?x=?y=min{λ/8,R/16}, 其中λ為入射聲波的波長, 且有λ=2πc∞/ω.計算域設置為緩沖區的厚度為 2λ.柯朗數c∞?t/?x=0.4 , 總計算時長為20(R+λ)/c∞.
跟隨文獻[21]中的思想, 散射聲場按如下方式計算:

其中psc,p,pvor和pinc分別表示散射聲壓、總壓強、旋渦的動力學壓強和入射波的聲壓.等到計算穩定后, 開始統計散射聲壓計算得到散射聲壓的均方根prms.
在時域直接數值模擬聲傳播的方法主要有三種: 求解Navier-Stokes方程(文獻[21])、求解線性歐拉方程(文獻[29])、還有本文所采用直接數值求解歐拉方程[32].
為了考察本文方法的可行性和計算的準確性,選取文獻[21, 29]中的標準算例計算并進行對比.算例中入射聲波波長λ=4R, 均熵Taylor渦的馬赫數Mv=0.125.為了更好的與文獻結果對比, 對所得到的散射聲壓均方根prms(r= 8R)使用入射聲壓的幅值pi進行歸一化處理, 從而得到無量綱量prms/pi, 其中pi=ε=10?5.下文中如果不作特殊說明, 在指向性分布計算中, 散射聲壓均方根prms均取r= 8R處的值.
由圖2可知, 采用本文所建方法得到的結果與Colonius等在文獻[21]和Clair等在文獻[29]中的計算結果吻合度很高, 尤其與文獻[29]中結果極度吻合, 充分證明了本文所建方法的有效性.相比直接求解Navier-Stokes方程, 歐拉方程不考慮粘性和熱傳導效應, 保證了均熵Taylor渦在計算過程中保持定常狀態(不考慮數值耗散), 有利于有效地提取出散射聲場, 因此歐拉方程更加適合具有解析解的旋渦對聲波的散射問題.而且相比線性歐拉方程的人為固定背景流動, 歐拉方程的計算是將動力學量和聲學量耦合在一起, 考慮了兩者之間的非線性相互作用, 更加接近真實物理情況.但是,采用歐拉方程計算的結果中同時也包含了聲學量和動力學量, 而將這兩部分物理量分離一直以來都是氣動聲學的難題.因此, 之前大多數聲散射研究采用將完整的歐拉方程近似為線性歐拉方程的方法達到分離動力學量和聲學量的目的.而本研究采用歐拉方程進行計算, 并將動力學量與聲學量成功分離是建立在等熵Taylor渦隨距離指數衰減的特性.在渦核區域外, 流體動力學量快速衰減, 背景流動接近靜止空氣介質.那么在遠場存在的擾動成分就只有聲波.

圖2 驗證算例Fig.2.Comparison with previous studies.
由第2節中旋渦和聲波的數學表達式中容易看出: 關于旋渦對聲波的散射問題中, 有四個重要的物理量—旋渦的渦核半徑R、旋渦最大速度Uv、聲波的波長λ和聲波的傳播速度c∞.前兩個物理量決定了渦的運動, 后兩個決了聲傳播的特性.為了衡量兩者的相互作用, 我們互相進行無量綱化, 就可以得到兩個重要的無量綱參數: 速度尺度比(這里我們稱之為旋渦強度, 即旋渦馬赫數)Mv和長度尺度比(聲波波長與旋渦渦核半徑之比)λ/R.在本文, 針對旋渦強度Mv從0.015625到0.25和長度尺度比λ/R從0.125到16范圍內旋渦聲散射特性進行計算并分析研究.這里旋渦強度和長度尺度比的具體選擇采用二分法, 從需要模擬的最大值開始, 每次縮小1/2, 直至需要模擬的最小值.前期經過大量數值模擬結果摸索出了這種取法, 可以用盡可能少的算例, 得到比較完整的規律(包含了聲散射的三個子區域).
為了更加直觀地評價兩個無量綱尺度參數Mv和λ/R對旋渦聲散射特性的影響(尺度效應), 引進散射截面Σ表征旋渦對聲波散射的強弱, 其中:

可以看出散射截面依賴于積分圓的半徑r.考慮到散射波以球面波的形式向外輻射, 在遠場(r?R)基本上以衰 減[21,29].因 此 以 πr進 行 無 量 綱 化,既使得在遠場時散射截面保持與半徑r無關, 又可以保證入射波與散射波強度相當時, 散射截面的大小與1接近.不失一般性, 計算散射截面時, 當時,r取8R; 當λ>πR時,r取16R.
圖3給出了當長度尺度比λ/R為0.125, 0.25,0.5, 1, 2, 4, 8, 16時, 散射截面 Σ 與旋渦強度Mv之間的關系.對于所計算范圍的所有λ/R, Σ 都會隨著Mv的增加而增加, 直至達到1的量級時, 不再會有明顯的變化.一旦散射截面超過100, 其與旋渦強度的平方關系將不再成立.散射截面與旋渦強度的定量關系依賴于長度尺度比的大小, 當λ/R>1 時,Σ與Mv的平方呈正比; 當λ/R<1 時,只有對旋渦強度較小時成立.顯然, 散射波的強弱與旋渦的強度呈正相關, 且最大只能達到與入射波相同的量級.

圖3 不同長度尺度比下散射截面與旋渦強度的關系(對數坐標系)Fig.3.Scattering cross-section Σ potted against M v at different λ /R (Logarithmic coordinate system).
圖4 給出了當旋渦強度Mv為0.015625, 0.3125,0.0625, 0.125, 0.25時, 散射截面Σ與長度尺度比λ/R之間的關系.對于所計算范圍的所有Mv, 散射截面Σ都會隨著λ/R的增加而減小, 且減小的速度隨著λ/R的增加而越來越快, 在長度尺度比較大時, 衰減速度的接近 (λ/R)?4.經典的理論分析證明在玻恩近似(即長波近似,λ/R→+∞時)下,衰減的速度可達到最大值 (λ/R)?4.按照文獻[22]針對希爾球渦聲散射的理論分析方法, 對于Taylor渦的聲散射結果, 在低馬赫數大波長的理想情況下, 同樣可以得到:


圖4 不同旋渦強度下散射截面與長度尺度比的關系(對數坐標系)Fig.4.Scattering cross-section Σ plotted against λ /R at different M v (Logarithmic coordinate system).
下面考慮旋渦強度Mv和長度尺度比λ/R兩者共同對散射截面 Σ 的影響.將整個旋渦聲散射定義在關于Mv和λ/R的平面上, 并且將該平面劃分為:長波近似區, 共振散射區和幾何聲學區.具體劃分如圖5所示, 以為特征變量.在長波近似區近似有這是由經典的長波近似理論, 即玻恩近似所得到的.在共振散射區有而λ/R的–4次方律不再成立.在該區域里入射聲波波長與渦核尺度相當, 散射聲場強度與入射聲場也相當, 產生了類似“共振”的散射效應; 在幾何聲學區,Mv的2次方律和λ/R的–4次方律都不再成立, 該區域的聲散射現象可以通過幾何聲學理論進行分析.其中, 長波近似區和共振散射區的分界線大致在λ/R=2π , 正好對應于亥姆霍茲數kR=2πR/λ=1 , 其中k為波數.而共振散射區和幾何聲學區通常以S= 1為分界線.下面分別對這三個區域的散射特性的尺度效應做詳細分析.

圖5 旋渦強度和長度尺度比共同對散射截面的影響(對數坐標系)Fig.5.Scattering cross-section as a function of M v and λ/R(Logarithmic coordinate system).
圖6 給出了無量綱時間t= 280時的無量綱化的散射聲壓psc/pi.在長波近似區, 旋渦對聲波的散射很弱, 散射聲場的強度遠遠小于入射聲場.散射聲場的指向性類似于四極子聲場, 出現了四道波束(beam), 分別位于θ=45?,?45?,?135?,135?附近, 依次記作第一波束(first beam), 第二波束(second beam), 第三波束(third beam), 第四波束(fourth beam).靠近x正方向的兩道波束強于靠近x軸負方向的兩道波束, 且任意相鄰的兩道波束之間相位相反.

圖6 t = 280時的散射聲壓, 其中 M v=0.125 ,λ/R=10Fig.6.Snapshot of scattered pressure withMv=0.125 and λ /R=10 at t = 280.
圖7 給出了長波近似區五個典型狀態的聲壓均方根 (r= 16R) 的指向性特征, 采用無量綱量作為特征量.和瞬時散射聲壓類似, 均方根值也出現了四道波束, 且每道波束的形狀類似一個半正弦函數.散射聲壓均方根prms與旋渦強度Mv呈正比.但是隨著Mv增加, 散射聲場的指向性分布表現出更多的不對稱性, 第一波束的強度逐漸增強, 相應地, 第二波束減弱.這與文獻[29]中的完全對稱的描述不符.這種差異主要來源于計算所使用的的控制方程不同.本文所采用歐拉方程模型是比文獻[29]中忽略高階小項的線性歐拉方程更加貼合真實的物理模型, 而且在旋渦的馬赫數較高時, 這種非線性的高階項愈發重要, 因此本文所得到的結果更為可信.在極低馬赫數時,指向性分布關于聲波入射方向完全對稱, 這與基于低馬赫數大波長的理論分析結果(玻恩近似[18,19]和匹配漸近展開[22?25])類似.同時, 隨著Mv減小,在θ=0?處的散射逐漸趨于0.散射聲場關于長度尺度比λ/R的“–2次方律”基本符合.同樣隨著λ/R增加, 不同波束的強度逐步接近, 散射聲場表現了良好的對稱性.

圖7 長波近似區的散射聲壓均方根指向性分布Fig.7.Directivity for root-mean-square pressure of the scattered fields in long-wavelength domain.
綜合以上結果, 可以得到在長波近似區散射聲場強度關于旋渦強度和長度尺度比的尺度律關系:

圖8給出了低馬赫數 (Mv=0.015625) 、大波長 (λ/R8) 時的散射聲場指向性分布(r= 16R).從圖8可以看出, 隨著波長的增加, 指向性分布逐步接近四個等強度的波束(參考曲線).我們可以得到在極低馬赫數、極大波長時散射聲場的分布函數有如下形式:

其中AT是一個常數, 針對等熵Taylor渦,AT=36.2038671967512.

圖8 低馬赫數大波長下的散射聲壓均方根指向性分布Fig.8.Directivity for root-mean-square pressure of the scattered fields with low Mach number and long wavelength.
隨著長度尺度比λ/R的減小, 旋渦對聲波的散射逐漸過渡到共振散射區.與長波散射相比, 共振散射區的散射聲場中, 靠近x軸負方向的兩道波束持續衰減, 直至消失.如圖9所示, 采用無量綱量作為特征量, 會發現靠近x軸正方向的兩道波束其強度仍然滿足λ/R的“–2次方律”,并逐漸向聲波入射方向靠攏.這是長波近似區向共振散射區過渡的顯著特征.從圖10可以看出, 當λ/R=4時, 散射聲場已經完全表現共振散射的特性, 只存在靠近x軸方向的兩道反相的波束.

圖9 共振散射區較大波長散射聲場指向性分布Fig.9.Directivity for root-mean-square pressure of the scattered fields in resonance domain with relatively long wavelength.

圖10 t = 140時的散射聲壓, 其中 M v=0.125 ,λ/R=4Fig.10.Snapshot of scattered pressure withMv=0.125 and λ /R=4 at t = 140.

圖11 旋渦強度對共振散射區散射聲場指向性的影響Fig.11.Directivity for root-mean-square pressure of the scattered fields in resonance domain at different vortex strength.
下面采用無量綱量prms/(piMv) 作為特征量考察散射聲場與旋渦強度的關系.從圖11可以看出,在共振散射區, 散射聲壓均方根值prms與旋渦強度Mv呈正比.和長波散射類似, 隨著Mv的增加, 散射聲場逐漸表現出不對稱性.結合以上結果, 可以得到在共振散射區散射聲場的強度關于旋渦強度有如下關系式:

針對共振散射區的短波散射, 如圖12所示,兩道主要的波束逐漸向聲波入射方向靠攏,x軸負方向的散射徹底消失.從圖13可以看出隨著旋渦強度的增加, 兩道波束發生干涉, 標志著散射由共振散射區轉向幾何聲學區.
在共振散射區域, 散射聲波的強度相比于長波近似區有明顯的增強, 接近入射聲波的強度.這個區域的聲散射特性對于采用聲波探測旋渦結構、基于旋渦陣列主動降噪等工程實用的聲學技術具有很強的現實意義.

圖12 t = 60時的散射聲壓, 其中 M v=0.125 ,λ/R=1Fig.12.Snapshot of scattered pressure withMv=0.125 and λ /R=1 at t = 60.

圖13 旋渦強度對共振散射區小波長散射聲場指向性的影響Fig.13.Directivity for root-mean-square pressure of the scattered fields in resonance domain with relatively small wavelength at different vortex strength.
隨著旋渦強度Mv增加和長度尺度比λ/R的減小, 聲散射逐漸進入到幾何聲學區, 其散射聲場強度高于入射聲場, 具體參考圖14和圖15.由于主波束之間的相互干涉, 產生許多強度較弱的次級波束, 且隨著λ/R減小和Mv的增加, 次級波束的數量逐漸增多, 散射聲場的指向性更加不規則.在幾何聲學區, 散射聲場的強度不再與旋渦強度呈正比.與長波近似區和共振散射區不同, 對于低馬赫數的旋渦, 幾何聲學區的散射聲場也表現出了不對稱性.

圖14 散射聲壓, 其中旋渦強度 M v=0.25 (a) λ /R=1 ; (b) λ /R=0.5 ; (c) λ /R=0.25 ; (d)λ/R=0.125Fig.14.Snapshot of scattered pressure with M v=0.25 : (a) λ /R=1 ; (b) λ /R=0.5 ; (c) λ /R=0.25 ; (d) λ /R=0.125.

圖15 不同旋渦強度對幾何聲學區散射聲場指向性的影響Fig.15.Directivity for root-mean-square pressure of the scattered field in geometrical acoustics domain at different vortex strength.
旋渦對聲波的散射包含了兩種聲散射的機制.首先是非線性散射的效應, 即旋渦與聲波非線性耦在渦核內.另外一種是長程折射效應, 主要發生在渦核外部.表現為背景流動對渦核內所產生的散射聲波傳播路徑長距離的偏轉作用, 是造成非對稱性的原因, 稱之為旋轉模態.
采用以上分析結果可以很好地解釋: 為什么隨著波長的減小, 反向散射越來越不明顯.原本由渦核所產生的散射波在前向和反向的強度是一致的,但是與入射波疊加后各個方向出現了差異.而這種差異隨著波長的減小越發明顯.使得反向的折射效應強于正向, 原本反向的散射波便被偏轉為接近入射波的方向.合發聲機制.由于入射聲波的存在, 旋渦不再保持原來定常狀態, 而是受到聲波的激發向外輻射出散射聲波以抵消入射聲波的影響.散射聲波的強度與旋渦強度成正比, 此外還與聲波和旋渦的尺度比值密切相關, 當兩者接近時, 非線性發聲達到最大值即通常意義上的共振.旋渦受激輻射聲波有良好的對稱性, 可以稱之為對稱模態, 且該過程主要發生
本文通過數值求解二維歐拉方程, 研究了無量綱尺度參數(旋渦強度Mv和長度尺度比λ/R)對均熵Taylor渦聲散射特性的影響.隨著旋渦強度的增加與長度尺度比的減小, 旋渦聲散射現象依次經歷長波近似區、共振散射區以及幾何聲學區, 并且散射聲場強度逐漸增加, 指向性表現出更多的非對稱性.
1) 在長波近似區, 散射聲場很微弱, 由四道波束組成, 相鄰兩道波束反相.散射強度與旋渦強度呈正比, 與長度尺度比的平方呈反比; 旋渦強度較弱時, 散射聲場關于入射方向表現出良好的對稱性;
2) 在共振散射區, 散射主要發生在靠近入射波的方向, 包含兩道反相的波束, 其強度與旋渦強度呈正比;
3) 在幾何聲學區, 散射聲場強度和入射聲場相當, 主要集中在旋渦后方區域, 沒有表現出明確的指向性;
4) 旋渦對聲波的散射包含了兩種不同的機制,即非線性的散射效應與線性的長程折射效應.
需要注意的是本文只考慮了相對靜止的旋渦對聲波的散射.因此聲波穿過旋渦后, 頻率并沒有改變.下一步研究將考慮運動旋渦對散射聲場的影響, 尤其是不同形式的運動(平動、振動等)給散射聲場帶來的多普勒頻移效應.