馬麗璇，李恩義，孫書霞，楊慶祥

（安陽工學院飛行學院，河南 安陽 455000）

在航空應(yīng)用領(lǐng)域中，空腔流動是空氣動力學的研究熱點之一，例如飛機起落架艙和飛行器武器艙等。當高速氣流流經(jīng)空腔時，空腔外部剪切流與內(nèi)部流動相互作用，經(jīng)常會出現(xiàn)自持振蕩現(xiàn)象，會產(chǎn)生劇烈的壓力、速度脈動，進而引發(fā)強烈的噪聲。該問題涉及到流體動力學中的非定常分離流、自由剪切層的不穩(wěn)定性、渦動力學、湍流剪切層中相干結(jié)構(gòu)（擬序結(jié)構(gòu)）、聲與流動的相互作用等前沿問題。對于武器艙來說，高于170 dB 的壓力脈動可能會損害飛行器結(jié)構(gòu)和導(dǎo)彈及防礙導(dǎo)彈的正常分離。

對于空腔流動，相關(guān)學者已進行了大量理論分析、實驗測量和數(shù)值計算。Rossiter[1]引入邊緣聲調(diào)分析理論并推導(dǎo)出用于估算空腔流場自持振蕩頻率的半經(jīng)驗公式，即著名的Rossiter 半經(jīng)驗公式。隨著技術(shù)發(fā)展，大量先進的非接觸式高精度測量方法[2-3]應(yīng)用于空腔流動中，如粒子圖像速度場儀（PIV）、激光多普勒測速儀（LDV）。由于流激噪聲在流場和聲場的長度、尺度存在巨大不一致性，準確計算聲場特性具有一定難度。在低雷諾數(shù)瞬態(tài)粘性流場中采用直接數(shù)值模擬（DNS）[4-5]是可行的，但對于高雷諾數(shù)流動則不可行，工程實踐中多采用雷諾平均N-S 方法、大渦模擬方法（LES）、混合RANS/LES 方法。Temmerman 等[6]采用URANS 方法分析了均流和單音模態(tài)及時間步長對噪聲預(yù)測的影響。Aybay 等[7]采用LES 方法研究了湍流流場空腔近場噪聲，并將計算結(jié)果同實驗值和高精度格式計算值進行了對比分析。Peng[8]采用分離渦模擬（DES）方法研究了空腔內(nèi)壓力脈動和流經(jīng)空腔上方的混合層及腔內(nèi)渦運動對后壁附近混合層的影響。王顯圣等[9]對空腔可壓縮流致噪聲問題的研究進展進行綜述，分析空腔流激振蕩及其產(chǎn)生噪聲的機制、參數(shù)影響規(guī)律與噪聲控制技術(shù)發(fā)展趨勢。劉瑜等[10]使用延遲分離渦模擬方法研究了來流馬赫數(shù)為0.85 的矩形開放腔體氣動聲學環(huán)境。余培汛等[11]基于耦合湍流速度生成模型與線化歐拉方程，建立了氣動噪聲混合預(yù)測方法，對空腔標模M219 的氣動噪聲進行了預(yù)測。

為改善空腔流場內(nèi)的聲學環(huán)境，可采用一些措施對流場進行控制，常用方法包括：主動控制和被動控制，如幾何修型、底部泄壓管、前緣立齒、前緣射流。Bastrzyk 等[12]通過在空腔前緣位置處放置圓桿來削減剪切層，達到降低噪聲的目的，并考慮了不同圓桿位置的降噪效果。Lada 等[13]通過分析3 種不同微射流速度工況，得出射流改變了剪切層的穩(wěn)定性和空腔內(nèi)的反饋回路，可實現(xiàn)降噪及減小腔內(nèi)阻力。Soemarwoto 等[14]分析兩種不同幾何形狀工況，得出修形后可使諧振頻率增大但振幅減小。Wilcox 等[15]通過在空腔底板開泄壓腔，使空腔后部的高壓區(qū)向前部的低壓區(qū)流動，從而達到平衡壓力的目的，進而有效地改善空腔流場特性。Ukeiley 等[16]采用在空腔前緣添加立齒的方法對空腔流動進行控制，結(jié)果表明僅將空腔前緣邊界層抬高并不能很好地降低腔內(nèi)噪聲。吳繼飛等[17]在高速風洞中對空腔流場氣動聲學特性進行實驗研究，對空腔后壁進行倒角，以降低氣流在該處的撞擊強度，從而達到抑制空腔流場氣動噪聲的目的。周方奇等[18]通過開展高速風洞實驗研究跨超聲速來流條件下前緣直板裝置對空腔（長深比為6）流動和噪聲的控制機理，通過對比多種前緣直板控制條件下的腔內(nèi)噪聲聲壓級分布，確定直板控制參數(shù)的優(yōu)化選擇方法及最優(yōu)參數(shù)；利用靜態(tài)／動態(tài)壓力傳感器和油流實驗采集腔內(nèi)靜壓、脈動壓力和壁面流譜，著重分析前緣直板對腔內(nèi)流動結(jié)構(gòu)、聲壓級和聲壓頻譜的影響規(guī)律。張群峰等[19]利用基于SST k-ω 湍流模型的改進延遲分離渦模擬方法，對開式空腔氣動聲學特性進行數(shù)值模擬研究，在腔體尾緣采用后壁開孔板加耗能腔的方法進行流動控制。

由于聲波頻率范圍很寬，高頻和低頻的聲波波長尺度相差很大，為能捕捉到最大頻率波，要求聲場計算網(wǎng)格要很密，時間步長取值要很小，聲場的計算時間較長。上述研究中，為取得良好的噪聲數(shù)值模擬結(jié)果，大多都采用DES 或LES 湍流模型，甚至于DNS 模型，所需計算網(wǎng)格分辨率要求高，數(shù)值計算代價大，不利于工程實踐應(yīng)用。為此，提出一種基于CFD/CAA 方法來求解聲學信息，即Cubic k-ω 湍流模型結(jié)合非線性聲學求解器方法（NLAS, nonliner acoustic solver）。首先，通過與實驗對比驗證了該方法的有效性。然后，為了降低噪聲，對比分析了4 種不同構(gòu)型，得到最優(yōu)結(jié)構(gòu)，該方法減少了網(wǎng)格需求量且精度較高，有利于工程實踐應(yīng)用。

1 數(shù)值計算方法

1.1 湍流方程

采用Cubic k-ω 湍流模型[20]，其湍動能k 和耗散率ε 的輸運方程為

其中：ρ 為流體密度；t 為時間；ui為速度在i 方向上的分量；xi為笛卡爾坐標系下i 方向的位移；Tt為時間尺度；μ 為粘度，μt為渦粘度；Pk為速度梯度產(chǎn)生的湍動能；E 為內(nèi)能；常量Cε1=1.44，Cε2=1.92，σk=1.0，σε=1.3。

1.2 非線性聲學求解器

小尺度湍流運動產(chǎn)生的噪聲源在一定程度上可較好地模擬。然而，對于影響相干擬序結(jié)構(gòu)或輻射噪聲的大尺度結(jié)構(gòu)模擬不太理想，對此Batten 等[21]提出了非線性聲學求解器的基本原理。

NLAS 計算噪聲的產(chǎn)生和傳播基于初始的統(tǒng)計穩(wěn)態(tài)湍流，數(shù)據(jù)可由傳統(tǒng)的RANS 模型提供。初始步提供了平均流場基本特征和強制設(shè)定的湍流脈動統(tǒng)計描述。NLAS 通過對統(tǒng)計結(jié)果的重構(gòu)獲得噪聲源來模擬壓力脈動傳播。

其中

其中

式中，cp為定壓比熱容。

NLAS 可以運行由原始RANS 計算結(jié)果截斷得到減小的子域，截斷的邊界被定義為遠場吸收邊界。該邊界不僅提供湍流統(tǒng)計，還可以提供變化的平均場數(shù)據(jù)，且外邊界更接近能表現(xiàn)流場重要特征的區(qū)域。此外，還可以進一步放寬對網(wǎng)格的要求。這是因為穩(wěn)態(tài)RANS 解可以插值到用于后續(xù)瞬態(tài)模擬中較粗糙的網(wǎng)格。所以，非線性聲學計算求解器在求解時對網(wǎng)格的需求較小，與DNS、LES、Hybird RANS/LES 的需求量對比如圖1所示。

1.3 湍流人工重構(gòu)

圖1 DNS、LES、Hybird RANS/LES、NLAS 所需近壁網(wǎng)格分辨率Fig.1 Required near-wall mesh resolutions with DNS,LES,Hybrid RANS/LES and NLAS

在NLAS 方法中，大尺度脈動可被直接描述；對于小尺度壓力脈動采用人工重構(gòu)方法。Smirnov 等[22]基于雷諾應(yīng)力張量的相似變換張量嵌入尺度概念，使其成功的應(yīng)用于各向異性湍流中。傅里葉模式方程[23]如下

其中：αik為應(yīng)力張量的Cholesky 分解；n 為模態(tài)階數(shù)；N 為階數(shù)上限；pni和qni分別為余弦和正弦的系數(shù)；ω 為湍流頻率；di為波矢。

2 NLAS 方法驗證

2.1 幾何模型

驗證算例采用文獻[24]中的M219 矩形開式空腔，其尺寸長：寬：深（L：W：D）=5：1：1，空腔深為0.101 6 m，在空腔底部設(shè)有10 個壓力傳感器（K20～K29），如圖2所示。計算域：進口到空腔前緣的距離是7.75 D，空腔尾緣到出口的距離是5.25 D，空腔開口到頂部邊界的距離是17 D，兩側(cè)邊到空腔的距離均為1 D，如圖3所示,其中箭頭為自由來流方向。

圖2 空腔結(jié)構(gòu)及傳感器位置Fig.2 Cavity structure and sensor location

圖3 空腔流動計算域Fig.3 Calculation domain of cavity flow

2.2 邊界條件和求解方法

自由來流流動條件為：馬赫數(shù)Ma∞=0.85，靜壓P∞=62 096 Pa，靜溫T∞=266.53 K，基于空腔長度的雷諾數(shù)Re=6.8×106。對于流動邊界，進口施加1 個穩(wěn)態(tài)湍流邊界層，在空腔前緣的邊界層厚度為0.37×Re-0.2，即10.16 mm；空腔壁面采用絕熱無滑移壁面；其余為基于流動特征的邊界并設(shè)置3 層遠場吸收層。對于聲學邊界，外部來流設(shè)為NLAS 外部邊界并設(shè)置3 層遠場吸收層，以抑制聲波在遠場邊界反射；空腔壁面采用絕熱無滑移壁面。

基于有限體積法離散控制方程，空間離散采用具有二階精度的TVD 格式，TVD 限制器為Min-mod，時間積分采用雙時間步長隱式方法。同時采用多重網(wǎng)格法加速收斂和并行計算減小計算代價。時間步長Δt=2×10-5s，共計算20 000 步，考慮到壓力的時間序列和時間均化，為獲得平均流動特性將前5 000 步數(shù)據(jù)舍棄。

2.3 M219 空腔聲學預(yù)測

空腔噪聲的主要特征：空腔開口處混合層的不穩(wěn)定性和空腔內(nèi)部瞬態(tài)渦運動。不穩(wěn)定性是由于剪切層與空腔后部壁面碰撞造成自持振動。在文獻[24]實驗中，測量了空腔底部10 個位置（K20～K29）的壓力脈動，其被均勻地分布在空腔底部X/L=0.05 ～0.95 位置上。

功率譜密度（PSD）無法直接測得，可通過Burg 方法將采樣壓力脈動轉(zhuǎn)化為功率譜密度，而聲壓級(SPL)可通過功率譜密度計算得到，即

其中，pref=2×10-5Pa 是最小可聽聲壓，也是基準聲壓。

圖4分別給出了位置點K20、K22、K29的計算聲壓級，對比實驗數(shù)據(jù)，NLAS 方法可與其很好地吻合，驗證了該方法在求解氣動噪聲的合理性。在圖5中，SADES[25]和NLAS 方法計算所得功率譜密度幅值大部分情況下均大于實驗值，可知兩種方法對壓力脈動的預(yù)測有些偏高。此外，SA-DES 方法求解獲得曲線較平滑，NLAS 方法求解獲得曲線的波動較大，可見NLAS方法能更好地捕捉到壓力脈動。

諧振頻率可用Rossiter 經(jīng)驗關(guān)系式評估

其中：U∞和M∞分別為自由來流速度和馬赫數(shù)；γ 為壓力波與空腔尾緣碰撞后產(chǎn)生的相位移，取0.29；L 為空腔長度，α 和κ 分別為相位延遲量和對流渦與自由流流速比，κ 取0.57?；赗ossiter 公式，Larcheveque 等[26]提出，誘導(dǎo)諧振與混合層渦在下游的移動速度和壓力波在空腔上游聲速相關(guān)。

圖4 空腔底部K20、K22、K29 位置點的聲壓級對比Fig.4 SPL comparison at locations K20，K22 and K29 on cavity floor

表1給出了K29處模態(tài)頻率的實驗值、Rossiter 公式估算值、4 階高精度格式算值[27]和NLAS 計算值。將NLAS 方法計算所得的1～4 階模態(tài)頻率與實驗值對比，可知1 階、2 階、4 階模態(tài)頻率略高于實驗值且誤差分別為2.6%、3.2%和2.7%，3 階模態(tài)頻率與實驗值吻合良好，誤差僅為0.3%。然而，4 階高精度格式的1 ～3階模態(tài)頻率計算值均偏低和4 階模態(tài)頻率偏高，誤差分別為-13.2%、-10.3%、-8.6%和2.7%?？芍琋LAS 方法可以很好地預(yù)測空腔噪聲的模態(tài)頻率。

表1 K29 處的模態(tài)頻率對比Tab.1 Modal frequency comparison at K29

3 被動控制

為抑制空腔噪聲，改善空腔流場環(huán)境，采取3 種不同結(jié)構(gòu)：①前緣壁面改為45°斜邊（Case 1）；②尾緣壁面改為45°斜邊（Case 2）；③前緣和尾緣壁面均改為45°斜邊（Case 3）。對于初始模型算例標記為Original。

圖6給出了K20、K22、K293 個位置點4 種工況下的聲壓級對比?？梢钥闯?，在Case 2、Case 3 工況下，聲壓級有所下降，實現(xiàn)降噪目的。對于Case 1 工況下，特別是空腔中部（K22），其聲壓級有所增大，未實現(xiàn)降噪目的。對于Case 3 來說，由于具有Case 1、Case 2 的構(gòu)型，因而一方面實現(xiàn)了降噪目的，另一方面又增大了空腔中部壓力脈動，使空腔中部聲壓級提高。此外還使其幾何結(jié)構(gòu)變得復(fù)雜，不利于機械加工。造成這種現(xiàn)象的原因是，由于前緣斜壁改變了氣流方向，自由剪切層向空腔底部移動，造成壓力脈動增大；尾緣斜壁使剪切層對空腔后壁的撞擊減緩，進而使壓力脈動減小，聲壓級降低。

時間平均總聲壓級計算如下

由圖7對比實驗值和NLAS 計算值（Original）可知，從K20到K21測量點總聲壓級稍微下降，然后沿著流動方向（從K21到K29）總聲壓級增大，即壓力脈動增大。實驗值與計算值趨勢一致，且計算值相比實驗值被過高預(yù)測約5 dB。

由于湍流剪切層的Kelvin-Helmholtz 不穩(wěn)定性會產(chǎn)生多種不穩(wěn)定性現(xiàn)象，如旋渦的破碎和運動渦核等[28]。為更好地顯示旋渦特征，可采用準則方法[29]。根據(jù)Hunt 理論，如果渦張量對流體變形的貢獻大于應(yīng)變率張量認為該區(qū)域有渦存在，準則定義如下

其中：Sij和Ωij分別為Δu 的對稱分量和反對稱分量。根據(jù)速度梯度張量具有廣義Galilean 不變性，可知其值與坐標系的選取無關(guān)。Q>0 時表示有渦存在，等值面包圍的區(qū)域定義了一個渦核，反之則無。圖8給出了4種結(jié)構(gòu)下空腔附近流場旋渦結(jié)構(gòu)（Q 準則）。Powell[30]和Howell[31]從渦動力學角度考慮流體發(fā)聲機理，建立渦聲理論，認為氣動噪聲主要源于流場中渦系的拉伸與破裂。后緣斜面起到渦流發(fā)生器的作用。

圖6 4 種工況下空腔底部K20、K22、K29 位置點的聲壓級對比Fig.6 SPL comparison at locations K20,K22,and K29 on cavity floor under four conditions

圖7 沿空腔底部總聲壓級Fig.7 OASPL comparison along cavity floor

由于渦結(jié)構(gòu)圖不能直觀的表示噪聲的大小，故引入VLamb的模。圖9為不同結(jié)構(gòu)下的VLamb模大小對比圖，其中VLamb=u×Ω，即速度叉乘以渦量，其可以反映聲源的大小。結(jié)合上述分析，由圖可知：①對于Original 結(jié)構(gòu)，剪切層脫落渦撞擊空腔尾緣壁面產(chǎn)生一次聲波并向前傳播至前緣，引起剪切層擾動，造成二次聲波以及渦脫落，完成1 個反饋回路，進而形成自持振蕩；②對于Case 1 結(jié)構(gòu)，結(jié)合圖5可知，由于前緣壁面變?yōu)閮A斜使氣流方向改變，造成空腔中部壓力脈動增強，噪聲增大；③對于Case 2 結(jié)構(gòu)，由于尾緣壁面變?yōu)閮A斜，降低了壁面對一次聲波的反射，進而降低了二次聲波的輻射，實現(xiàn)了降低噪聲的目的；④對于Case 3 結(jié)構(gòu)，由于前、尾緣壁面均變?yōu)閮A斜，使其一方面降低了一次聲波反射，減小了噪聲，另一方面空腔中部壓力脈動增大，噪聲變大。此外，當采用斜壁面時，增大了空腔的長度，使空腔流場特性由開式向閉式過度。

圖8 湍流結(jié)構(gòu)等值面（Q ＝7×106）Fig.8 Iso-surface of turbulence structure（Q ＝7×106）

圖9 不同構(gòu)型VLamb 的模大小Fig.9 Magnitude of VLamb with different configurations

4 結(jié)語

基于CFD/CAA 方法對三維可壓縮開式跨聲速空腔M219 進行了數(shù)值分析可知，聲壓級風洞實驗值與數(shù)值計算值可以很好地吻合；對比風洞實驗值、NLAS 和SA-DES 計算值的功率譜密度，可知NLAS 方法可更好地捕捉到壓力脈動；分析對比Rossiter 經(jīng)驗關(guān)系式、4 階高精度格式與實驗值的模態(tài)頻率的誤差率；驗證了非線性噪聲求解器方法（NLAS）可以很好地預(yù)測空腔噪聲。

采取被動控制措施幾何修型以實現(xiàn)降低空腔氣動噪聲的目的，分析和對比4 種工況下，其聲壓級、總聲壓級、Q 準則等值面、Lamb 矢量?？芍熬壭扌驮斐煽涨恢胁繅毫γ}動增大，噪聲也隨之增大；尾緣修型降低了壁面對一次聲波的反射，進而降低了二次聲波的輻射，實現(xiàn)了降低噪聲的目的，即Case 2 降噪效果最好，大約降低5 dB，可應(yīng)用在實際工程中。