開式方腔可壓縮流動及噪聲LES計算分析

馮 峰，王 強

(中國航天空氣動力技術研究院，空氣動力科學中心， 北京 100074)

1 引言

戰斗機和轟炸機內埋武器艙可壓縮空腔流動噪聲可超過160 dB，能夠致使機載電子設備失靈、艙內結構疲勞損傷，甚至誘發武器在分離過程中非常規運動，導致投放精度降低[1]。此外，空腔流強噪聲輻射還會影響商用航空運輸的噪聲適航性及軍用戰斗機、轟炸機飛行隱蔽性[2]。20世紀50年代，Krishnamurty[3]最早開始研究空腔流動現象，并指出空腔流中存在聲-流耦合共振的聲反饋環機制。當高速氣流掠過空腔上方表面后，形成高度有序的剪切層渦列結構，該渦列與空腔后壁撞擊產生反射聲波擾動空腔前緣初始剪切層，當擾動波與反射壓力波頻率相匹配時，空腔內出現聲場-流場相互耦合的聲反饋環流動振蕩現象。之后，Rossiter[4]發展了半經驗公式對聲-流共振現象提供了更完整的描述，并可預測空腔流振蕩頻率，在空腔流研究中得到廣泛應用[5]。

針對超聲速飛行條件下戰斗機內埋武器分離投放可靠性問題，高保真數值技術在探索超聲速空腔流機制方面展現出良好的前景。其中LES(large eddy simulation，大渦模擬)方法可以精確求解空腔內各類渦結構及其產生的噪聲特性，具有很強的預測分析空腔流噪聲機制的能力。Sinha等[6]和Smith等[7]較早采用三維LES方法分別計算了Mach 2.0和Mach 1.5超聲速開式空腔流動噪聲。其中，Smith等[7]的LES結果顯示壓力頻譜峰值與實驗結果一致。Rizzetta和Visbal[8]采用了基于動態Smagorinsky亞格子模型的LES技術求解了長高比L/D=5，Mach數1.19，Reynolds數 200 000的超聲速方腔流動，并應用了高精度緊致格式和大量計算網格，獲得的方腔底部多個位置上壓力譜與實驗結果在幅值和特征頻率上趨勢一致。但上述研究缺乏細致的流場動力學分析，近期在深化空腔流機理認知以改進噪聲控制的需求牽引下，人們再次開始關注展向受限條件下空腔流動基本特性。Brès和Colonius[9]使用DNS(direct numerical simulation，直接數值模擬)方法和整體穩定性理論分析了低Reynolds數1500、亞聲速0.325空腔流動展向穩定性，并發現三維空腔模態具有一個與方腔寬度相當的展向波長，其頻率比二維Rossiter模態低一個量級。Sun等[10]完成了一個有限寬度超聲速開式方腔流動LES計算分析，研究了側壁效應、自由來流Mach數、Reynolds數等對空腔流動的影響。整體而言，目前對可壓縮空腔流動及噪聲特性及機理細節認識仍不足。本文采用基于計算氣動聲學技術的LES方法計算分析三維Mach 1.4和Mach 0.6超/亞聲速來流條件下可壓縮開式方腔流動。針對空腔流動中剪切層演化、多尺度渦結構、聲反饋、可壓縮效應等基本動力學過程開展分析，以深化可壓縮空腔流動及噪聲現象認識和機制理解。

2 數值方法

2.1 控制方程及SGS模型

為實現LES方法，采用無量綱化三維可壓縮Favre濾波Navier-Stokes方程

(1)

(2)

(3)

(4)

(5)

(6)

2.2 計算格式

控制方程采用基于計算氣動聲學的高精度數值格式進行離散，空間格式使用Tam等[11-12]提出的低色散低耗散誤差7點4階DRP格式，時間推進采用Berland等[13]提出的優化低存儲6步4階Runge-Kutta格式。由于DRP數值格式不能精確離散高波數短波，為避免數值高頻波導致Gibbs振蕩，引入Tam[11-12]構造的人工選擇性阻尼項對非物理數值短波進行濾波。選擇性阻尼項、黏性項及SGS項在計算中均在Runge-Kutta最后一個子步引入。對于含強湍流脈動和激波間斷的超聲速空腔流場，采用Bogey等[14]提出的自適應空間濾波法對守恒變量濾波以捕捉激波并維持計算穩健性。

2.3 模型設置

計算模型為一個有限寬度的開式方腔，如圖1所示，方腔長深比L/D=6，寬深比W/D=2，D、L、W分別為方腔深度、長度和寬度。以方腔深度為特征長度的Reynolds數Re=10 000，來流邊界層厚度為δ=0.167。

圖1 方腔模型及計算網格示意圖

亞/超聲速方腔流動使用相同Cartesian網格，900萬個網格點求解空腔流動及噪聲。空腔內設置199×103×91個網格點，計算域范圍為x1=0～6，x2=-1～0，x3=-1～1。空腔外計算域339×141×183個網格點，計算域范圍為x1=-3～13，x2=-1～9，x3=-2.5～2.5。壁面附近對網格進行加密，各方向最小網格尺度為Δx1=0.002，Δx2=0.001，Δx3=0.002。超聲速方腔算例進行500 000次迭代，時間步長Δt=0.000 5，流場統計在第100 000步后啟動。亞聲速方腔算例進行250 000次迭代，時間步長Δt=0.000 5。為加速計算，LES程序采用MPI方法，使用84個進程，每網格點每次迭代約花費0.5 μs CPU計算時間，一次計算花費約480 h。

2.4 Rossiter空腔頻率預測模型

Rossiter[4]于1964年提出了預測空腔流噪聲特征頻率的半經驗公式，或稱Rossiter模態預測，其形式如下

(7)

式中:fm為第m個模態的頻率，u、Ma為自由來流速度和Mach數，γ為比熱比，α為相位移。κν為大尺度渦結構對流速度與自由來流速度的比值，是依賴于空腔幾何形狀和測試條件常數。本文LES計算中在方腔底部(3，-1，0)和中心(3，-0.5，0)設置兩個測壓點，并將計算壓力脈動值與Rossiter預測結果對比，相互驗證并確認空腔噪聲特征頻率。

3 驗證與校核

將計算統計得到的空腔流時均流向速度和湍流脈動速度與Dudley和Ukeiley[15]的實驗結果和DES數值結果進行對比，以驗證LES結果的準確性。圖2展示了相同來流條件下本文LES計算和Dudley和Ukeiley的實驗測試的方腔x3=0截面上平均流向速度，兩圖等值云圖范圍均在u/U∞=-0.25～1之間，增量均為0.25。圖中LES和實驗獲得的超聲速空腔流流向速度沿空間分布完全一致，均展現了明顯的剪切流動特征，即方腔上方剪切層沿下游逐漸增厚，受方腔后壁遮擋而浸入腔內的動力學現象，空腔前部及底部平均流向速度相對較低。

圖2 時均流向速度場云圖

圖3展示了方腔x3=0截面上流向脈動速度均方根，LES結果與文獻中DES結果[15]采用等值云圖范圍一致，均為〈u1rms〉/U∞=0～0.25，增量0.05。計算與實驗結果均顯示了湍流脈動速度沿剪切層向空腔后部逐漸增強，在方腔后部剪切層區域湍流脈動達到最大值的。受方腔后緣誘導，后緣表面也形成一個小范圍的強湍流區。值得注意的是，在靠近空腔后底部小的湍流區域，文獻DES獲得的湍流脈動速度強度高于本文LES計算統計結果。這是由于DES計算的底壁面邊界層較薄，易于在腔內剪切層于后壁回流誘導的逆壓梯度下，形成較大范圍湍流分流區，致使后底壁區呈現較LES結果高的湍流脈動強度。

圖4展示了方腔內7個流向站位上LES與實驗測量獲得的時均流向速度型和流向脈動速度均方根。如圖所示，二者很好地吻合，表明LES計算準確地捕捉到了剪切層和方腔流動對流場平均速度及湍流脈動的影響。在空腔下游第六和第七個站位處，數值和實驗的時均流向速度和流向脈動速度均方根在靠近空腔底壁面時出現最大約10%的相對誤差。對比本文LES及文獻DES結果[15]，發現在空腔后底壁面區域，LES計算湍流脈動速度強度略低于實驗測試結果，而DES結果則略高于實驗結果。可見未來仍需不斷深入探索以完善更精確的空腔流動數值分析方法。

圖3 時均流向脈動速度均方根云圖

圖4 空腔內不同流向站位時均流向速度和 流向脈動速度均方根示意圖

4 LES結果與分析

4.1 瞬態流場

LES方法可精確計算空腔流動渦、聲、激波等全部流場細節，便于對復雜振蕩空腔流動的演化特征和機制開展分析探索。如圖5，使用Q準則等值面(Q=3)和壓力云圖展示了Mach 1.4超聲速開式空腔瞬態流場。來流進入計算域后，在壁面誘導的氣流壓縮效應下計算域前緣形成一道激波。氣流經過空腔前緣后，與腔內氣體混合使得邊界層向剪切層轉化，且由于流動不穩定性剪切層卷起形成展向渦結構。渦結構在向下游運動過程中沿展向失穩，產生大量Λ渦及其他更小尺度的結構，在高速來流裹挾下繼續向下游對流。在空腔后緣，剪切層被撞擊產生大量小尺度湍流渦結構，一部分侵入空腔后部與腔內流體急劇混合，另一部分溢出空腔沿壁面繼續向下游流出計算域。剪切層撞擊空腔后緣還形成向上游傳播的聲波，聲波輻射出空腔后被超聲速來流的對流作用下向下游遠場傳播，在腔內聲波穿過亞聲速流動輻射至空腔前壁面。

圖5 瞬態Q準則等值面(Q=3)及壓力云圖

圖6展示了空腔中心x1x2截面(x3=0)上4個連續時刻密度梯度幅值云圖，由圖可見復雜超聲速空腔流場動力學演化。首先，圖中清晰地展示了剪切層大尺度結構的運動發展過程。在t=139時刻大尺度結構形成并向下游對流，t=139.5時刻該大尺度結構向方腔后緣靠近，同時新的渦結構卷起，t=140時刻原大尺度結構繼續逼近方腔后緣且頭部結構開始變形，新的渦結構顯著增強，t=140.5原大尺度結構與方腔后緣撞擊，新大尺度渦結構形成并向下游對流，進而完成了一次大尺度渦結構演化周期。同時，如圖6(c)所示，可觀察到剪切層附近的亞聲速流動區域內形成向上游輻射的聲波，其波長與大尺度渦結構相當。這證實了空腔振蕩流動中剪切層與空腔后壁面撞擊產生向空腔前方傳播的聲波，并擾動空腔前緣剪切層，形成自持振蕩的聲反饋環現象。此外，超聲速來流在開式空腔誘導下在空腔前緣也形成一道激波，如圖6(a)所示，由于與剪切層相互作用，其始終處于振蕩狀態，并反向作用影響剪切層的發展演化。整體而言，三維超聲速開式空腔流動的聲反饋環特性較為顯著。

圖6 密度梯度幅值(紋影圖)演化云圖

4.2 空腔壓力脈動

在空腔底部壁面中心P1點(3，-1，0)和空腔空間中心P2點(3，-0.5，0)兩個典型位置設置測壓點以觀測腔內壓力變化。如圖7(a)所示，P2點壓力脈動幅值略高于P1點，這是由于湍流在空腔空間中的脈動強度強于空腔底壁所致。使用了Fourier變換計算兩點上壓力頻譜，由于壓力信號計算及存儲時間有限，對原始壓力脈動信號做Hamming窗濾波處理以方便識別主頻信息。

圖7 空腔內P1、P2兩測點的壓力脈動和頻譜曲線

從圖7(b)可見，壓力脈動頻譜具有多個峰值，其中前3個最優峰值頻率被標記為f1=0.88、f2=1.88和f3=2.84(圓頻率)。兩觀測點上前2個主導頻率一致，但第3個主導頻率在空腔中心位置較為顯著，展示了腔內復雜的非線性渦系運動特性。結合第4.1節流場演化特征，分析可知主頻f2與聲反饋導致的剪切層渦運動相關，而頻率f1與空腔后部運動較緩慢的大尺度渦運動相關，據信f3應是2種運動交互作用衍生頻率。值得說明的是，相似的脈動主頻在Dudley和Ukeiley[15]和Sheta等[16]的工作中也被發現過，表明超聲速空腔流存在固有動力學特性。

此外，如圖7(b)所示，虛線為Rossiter公式[4]預測的空腔Rossiter模態，其中設定常參數κν=0.66，相位移由經驗公式α=0.062×L/D計算得到，本文方腔模型α=0.397。LES計算得到的3個壓力脈動主導頻率與Rossiter預測頻率非常吻合，證明Rossiter公式對超聲速空腔流動共振模態也具有預測能力。

4.3 壓縮性效應

為明確壓縮效應對空腔流動的影響，對Mach 0.6亞聲速空腔流進行LES計算及對比分析。值得說明的是，除來流速度外，該算例其余來流參數以及計算網格均與Mach 1.4超聲速空腔流一致。圖8展示了亞聲速空腔流動密度梯度幅值云圖隨時間的演化。由圖8(a)可見，該流動仍然受剪切層渦結構運動主導，在空腔前緣位置卷起一列規則小渦。約在x1=2～5站位，細小渦積聚形成大尺度渦團結構。如圖8(b)，包含大量小渦的大尺度渦團開始向空腔后緣撞擊，破碎成大量的小湍流渦，同時形成向上游傳播的聲波。聲波在空腔中沿腔內向上游傳播，則與初始剪切層耦合形成聲-流共振狀態。聲波向腔外輻射，形成上游遠場占優的強聲場，這與超聲速空腔流中聲波無法向上游傳播不同。另外，由于沒有激波出現，缺少了激波與剪切層相互干擾現象。

圖8 亞聲速(Ma=0.6)空腔流動密度梯度幅值云圖

圖9 亞/超聲速空腔流壓力脈動及其頻譜曲線

5 結論

1) 超聲速(Mach 1.4)空腔流動LES數值結果與文獻中實驗及數值結果對比，表明基于低色散低耗散計算氣動聲學技術的LES方法對空腔流計算的精確性、魯棒性，解決了常規CFD方法或商用軟件計算精度低不易捕捉噪聲及邊界非物理聲反射等問題。

2) 超聲速來流條件下，剪切層大尺度渦結構與空腔后壁發生撞擊并產生經腔內向前傳播的反射聲波，該聲波與空腔前緣初始剪切層相互作用，導致空腔內形成聲-流耦合自持振蕩的聲反饋環現象。基于該認識，可設置專門裝置在超聲速戰斗機武器艙等工程應用中破壞來流剪切層內大尺度渦結構的一致性，以削弱聲-流耦合現象，降低腔內噪聲水平。

3) 超聲速空腔流展現出復雜非線性渦系運動導致多個特征頻率。觀測發現前3個頻率較為顯著，與流場演化歷程結合分析可知第2個主頻f2最優，其與聲反饋導致的剪切層渦運動相關，而第1個主頻f1與空腔后部較緩慢的大尺度渦運動相關，第3個主頻f3是前2種運動相互作用衍生。

4) 當來流壓縮效應減弱時，亞聲速(Mach 0.6)空腔流動剪切層渦團與后緣撞擊形成的聲-流共振現象與超聲速情形基本一致，但亞聲速空腔流動的壓力脈動強度明顯減弱，且流動振蕩主導頻率減小，表明腔內渦運動緩慢。