方腔渦流運動對壓力脈動噪聲影響的數值分析

郝維，劉正先，陳麗英

(1.天津大學 機械工程學院，天津300072;2.天津大學 化工學院，天津300072)

流體流過方腔引起的脈動噪聲是實際工程中經常遇到的一類問題，如風管道、渦輪機中的空腔以及水下航行體的空隙和排水孔等.對腔體內部流動機制的研究得到許多研究者的重視[1-4]，近年來，隨著計算機技術的飛速進步帶動計算流體力學(CFD)快速發展，基于計算流體力學與氣動聲學交叉學科的計算氣動聲學得到迅速發展.采用數值模擬的方法研究流體在固體腔中由流動形成的噪聲，以及非定常湍流的渦流形成機理，為研究渦流噪聲進而削弱和消除噪聲提供了基礎的數值依據.Powell等[5]認為低馬赫條件下的等熵絕熱流體，其流體動力場和輻射聲場的基本且唯一的源是渦，可以通過研究流場中渦的運動分析其聲壓變化規律.

本文通過對方腔內流場的數值模擬，分析渦旋產生、發展的規律及與其相應關聯位置處壓力脈動變化之間的耦合性，采用快速傅里葉變換得到頻域下的脈動壓力級值，并與實驗數據進行對比和驗證.

1 數值分析方法

1.1 湍流模型

以方腔為研究對象建立計算模型，通過數值模擬求解湍流N-S方程得到腔體內的湍流流場和渦流運動信息.分別利用標準k-ε模型和大渦模擬(LES)實現對湍流脈動量的模化，標準k-ε模型用于確定流場的定常平均參數，為進一步應用LES方法模擬非定常流動條件下的渦旋產生、發展和耗散提供基礎信息.

1.1.1 標準 k-ε模型

標準k-ε模型是目前應用較廣泛的兩方程湍流模型[6]，它對運動方程中的所有脈動量實施模化，湍流粘度通過求解湍動能k方程和湍動耗散率ε方程得到，k和ε的求解方程如下:

式中:t為時間，s;Xi為空間坐標(i=1，2，3);ρ為流體密度，kg/m3;μ 為流體的動力粘度，Pa·s;μt為湍動能粘度，Pa·s;Gk為由平均速度梯度引起的湍動能產生項;Gb為由浮力影響引起的湍動能產生項;YM為湍流脈動膨脹對耗散率的影響.其中，Gk和YM對不可壓縮流動為 0;C1=1.44，C2=1.92，C3=0.09;湍流普朗特系數分別為:σk=1.0，σε=1.3.

1.1.2 大渦模擬

大渦模擬方法(LES)是介于直接數值模擬(DNS)與雷諾應力模型(RSM)之間的一種湍流數值模擬方法.LES假設湍流中的動量、質量、能量及其他物理量的輸運主要受大尺度渦影響，其湍流運動通過N-S方程直接求解，而小尺度渦對湍流的影響則通過模化在N-S方程中體現出來.

LES的運動方程通過對N-S方程在波數空間進行濾波得到.過濾原則是去掉比過濾寬度小的渦旋，大渦控制方程如下:式中:ui表示坐標維度下平均速度(i=1，2，3);τij為亞格子尺度應力，

方程式(3)～(4)與雷諾平均方程相似，不同的是方程的變量是過濾后的物理量，而非時間平均量.式中的湍流應力采用Smagorinsky[7]的基本亞格子應力(SGS)模型，該模型在LES方法中應用較廣，且取得了有效的結果[8-9].

亞格子應力具有下列形式:

1.2 壓力脈動分析方法

1.2.1 壓力級的傅里葉變換

數值模擬得到的關鍵位置點的噪聲分析數據來源于確定時間段內的壓力脈動信息，該脈動信息是時域下的流動信號，而湍流噪聲是頻域信息，故采用快速傅里葉變換(FFT)算法[10]實現LES結果從時域到頻域空間的轉換.

對有限長序列x(n)進行離散傅里葉變換(DFT)表示為

對復數序列x(n)中的一個k值，按上式計算X(k)需要N次復數乘法和(N-1)次復數加法.對所有N個k值，共需要N2次復數乘法和N(N-1)次復數加法.利用旋轉因子的對稱性)和周期性，將長序列大點數的DFT分解為小點數的DFT，利用多個小點數DFT的計算代替整體的DFT計算，達到降低運算量，明顯提高DFT運算速度，縮短運算時間的目的[11].

1.2.2 壓力級的倍頻程處理方法

與本數值模擬結果進行對比分析的實驗數據是由國內某大型水槽實驗室通過對相同方腔模型的水動噪聲測試提供，實驗針對方腔模型的關鍵位置點進行噪聲監測，并對監測數據進行1/3倍頻程處理.

數值模擬模型除與水槽實驗模型具有相同的幾何尺度、流場邊界條件外，還應保持相同的數據處理方法.實驗采用1/3倍頻程數據處理方法，總壓力級等效為

式中:i為倍頻程帶寬內3個1/3倍頻程帶寬信號，Lpi為倍頻程帶寬內第i個1/3倍頻程壓力級測量值.

完成一個標準的1/3倍頻程分析，需要分別進行32次測量[12]，因此實驗測量得到32個頻率及對應的脈動壓力級.為便于數值模擬結果與實驗數據的比較，對數值模擬數據經傅里葉變換后再進行相同方法的1/3倍頻程頻處理，使其轉化為與實驗測量頻率同值下的噪聲壓力級數據.

2 計算模型

方腔計算模型如圖1所示，計算區域總尺度為4.3 ×0.9 ×0.8(長 × 深 × 寬，單位:m);其中方腔尺度為0.3 ×0.3 ×0.2;流動介質為水;計算中入口給定均勻流速為7 m/s;對應雷諾數為2.1×106;出口按照不可壓縮流動條件滿足質量守恒;下邊界為不滲透固體壁面，滿足無滑移流動條件，上邊界設速度與壓力梯度均為0.

圖1 方腔模型結構及邊界條件Fig.1 The cavity model and the boundary conditions

圖2為腔體在X-Y與Y-Z流面上的網格分布.在網格生成過程中重點對渦旋生成、發展區域進行加密，而對主流區域采用逐漸過渡方法減少網格總數.根據模型對壁面處網格的要求，壁面第1層網格的無量綱垂直距離y+保持在10左右，整體計算域的網格總數為186萬，其中方腔內部網格數為94.5萬.非定常流場的計算時間步長為0.000 1 s，共計算10 000個時間步長，即1.0 s間流場的運動過程.

圖3為方腔關鍵點位置，前緣點p1用于考察左側來流經過腔體前的流動參數變化，腔底p2點用于監測流體在腔體內部流動對底壁面的影響，后緣點p3用于分析來流撞擊臺階時流動參數變化.首先分析3點處的脈動壓力值，再進行聲壓分析.

圖2 腔體網格分布Fig.2 The grid distribution for the cavity hole

3 數值模擬結果分析

圖3 監測點位置Fig.3 The location of the reference points

3.1 流場分析

圖4為方腔中間面在t=1.0 s時刻速度矢量圖.可看到，圖中一個明顯的順時針大尺度渦，同時在腔體下側兩邊角處分別有一個逆時針角渦.腔體左側前緣臺階的存在，使流動在此處出現分離，并在下游形成渦流.不同時刻的渦流運動分析發現(圖5)，該渦流在向后運動過程中逐漸發展和合并，渦的尺寸不斷增大，最終在腔體右側后緣臺階處破裂，一部分能量沿水平壁面向右傳遞，并在后緣誘發產生新的渦旋;另一部分能量則沿垂直壁面向下傳遞，成為腔內中心大渦形成的主要因素.在流場中，圖4的中心渦和所有角渦保持相對穩定的位置和尺度，而前緣和后緣處的渦則進行生成、運動、發展和耗散變化，呈明顯的周期變化規律.

圖4 腔體t=1.0 s時刻的速度矢量Fig.4 The velocity vector at the end of t=1.0 s

圖5為3個不同時刻方腔渦流隨時間的形成和發展過程，其中不等間隔時間t1、t2、t3分別取為方腔前緣臺階渦生成、后緣臺階渦到達和破裂時刻.在圖5(a)中可明顯觀察到渦A后方連續的渦旋存在，且渦的尺度不斷增大;渦C在t2時刻接近后緣處并具備了分裂趨勢，在t3時刻可明顯看到分裂后的2個渦分別向水平和垂直方向移動;渦B則為后緣臺階上新的誘導渦.

3.2 壓力脈動數值分析

圖5 腔體不同時刻渦量Fig.5 The vorticity distribution in the cavity hole on three time points

圖6為方腔3個位置點即p1、p2和p3在時間步長為4 000～10 000時段的壓力脈動強度分布.可以看出，三者的壓力值均顯示出周期或準周期波動特征.就波動幅值而言，p3點最大，p2點最小，p1點的周期性效果最不明顯;從平均值看，p1點基本在零值附近，p2和p3點均為負值，p3點的負壓值最大.分析認為各點的壓力脈動特征與其位置相關:p3點位于腔體后緣臺階處，前緣處流體產生的渦旋經生長和合并，在此處破裂，發生較劇烈的變化，同時沿后緣水平壁面再次形成渦流，由此導致壓力波動劇烈且幅值較大;p2點位于方腔底部中心，由于腔體的中心渦尺寸雖大但較穩定，因此p2點的壓力波動幅值很小;方腔進口端的p1點，受水流入腔時邊界條件變化影響，突然從穩定流態演變為劇烈的后臺階渦旋流動，并受中心渦的制約，在臺階下方的垂直壁面附近誘發出逆向小渦，由此形成了雖呈周期性波動但幅值不穩定的壓力變化規律.

以p3點的壓力脈動為例，提取其波峰值和波谷值對應時刻的渦量位置如圖7所示，可以明顯看到，圖6中的波峰a點對應圖7(a)中渦流到達后緣處對壁面形成的擠壓狀態，而波谷b點正是圖7(b)渦旋沿水平和垂直方向分裂為2個渦時的對應狀態.

圖 6 p1、p2、p3 位置點壓力脈動Fig.6 The pressure pulsation for p1，p2and p3

圖7 p3點壓力峰、谷時刻渦量Fig.7 The vorticity at peak and valley point of pressure for p3

4 渦流壓力級分析

4.1 傅里葉壓力級

數值模擬得到流場的壓力信息通過以下處理得到脈動聲壓值:以MATLAB程序為平臺，確定壓力點的采樣點數為6 000(由數值模擬數據分析，確定0.4～1.0 s時間段，壓力脈動呈現準周期且穩定的波動規律)，采樣頻率為10 000，對所得壓力信號進行傅里葉轉換并提取幅值，再進行頻率轉換得到頻譜圖，即頻率-脈動壓力級圖.

圖8分別為 p1、p2、p3位置點的頻率-脈動壓力級.可看到，隨著頻率值的增大，壓力級呈下降趨勢，同時出現若干峰值點.把大于10 Hz的第1高峰值做為主頻值，則 p1、p2、p3點的主頻值分別為42、26.67、37.14 Hz.各點的次高和其余主頻值見表 1.由壓力級峰值對應的頻率可以發現，p1的頻率普遍大于p2和p3;p3的脈動壓力級最大，即渦旋撞擊臺階產生的壓力脈動對聲壓值貢獻最大.

另外，還考察了p1位置處渦旋隨時間的變化過程，發現渦旋在向右發展過程中，不斷與周圍小渦旋合并、生長，這導致p3的脈動頻率減小.

圖8 監測點的頻率-脈動壓力級對比Fig.8 The frequency-SPL diagrams of the reference points

表1 參考點頻率與脈動壓力級峰值對應表Table 1 The correspondence between frequency and SPL peak value for the key points

4.2 壓力級的數值與實驗數據比較

實驗測量脈動壓力數據和數值模擬計算數據均運用前節的噪聲處理方法得到關鍵位置的脈動壓力級.圖9分別為p1、p2、p3這3個位置點壓力級的數值與實驗數據比較，其中橫坐標為各頻率對應的序號Ni，頻率對應值見表2.

從圖9可以看出，數值與實驗結果的總體趨勢一致，均隨參考點頻率值的增大，脈動壓力級呈下降趨勢.3個壓力點的數值與實驗數據在40 Hz之后呈現很好的符合度，在40 Hz之前則一致地表現出模擬值大于實驗值.

圖9 壓力級數值與實驗數據比較Fig.9 Comparison between simulation and experiment

表2 頻率與Ni對照表Table 2 The correspondence between frequency and Ni

5 結論

1)通過對方腔內部渦流流場的非定常數值模擬和分析，確定腔體不同位置點處渦旋產生、發展和破裂運動具有脈動周期性，壓力脈動的峰、谷特征與渦流狀態具有良好的對應性.

2)采用快速傅里葉變換將數值計算的時域信息轉換為壓力脈動頻域數據，通過1/3倍頻程處理方法得到各頻率下的壓力級值.方腔不同的位置點具有不同的壓力脈動頻率和壓力級特征，與其渦流運動密切相關.

3)與實驗測量結果的對比表明，除40 Hz以下低頻區外，兩者具有很好的符合度.表明采用數值方法模擬方腔內的渦流流場，進而結合快速傅里葉變換和倍頻程處理方法實施對方腔內渦流引起壓力脈動的研究是可行的，可以為方腔類的水動力噪聲和消聲控制提供參考.

[1]耿冬寒，劉正先.大渦模擬-Lighthill等效聲源法的空腔水動噪聲預測[J].哈爾濱工程大學學報，2010，31(2):182-187.GENG Donghan，LIU Zhengxian.Predicting cavity hydrodynamic noise using a hybrid large eddy simulation-Lighthill＇s equivalent acoustic source method [J].Journal of Harbin Engineering University，2010，31(2):182-187.

[2]朱習劍，何祚鏞.水洞中突出矩形腔的流激駐波振蕩研究[J].哈爾濱工程大學學報，1993，14(3):41-52.ZHU Xijian，HE Zuoyong.Study of flow-induced standing wave resonance of rectangular cavity in water tunnel[J].Journal of Harbin Engineering University，1993，14(3):41-52.

[3]GHAR I M，ROSHKO A.The effect of flow oscillations on cavity drag[J].Journal of Fluid Mechanics，1987，177(10):510-530.

[4]WANG M，FREUND J B，LELE S K.Computational prediction of flow generated sound[J].Annual Review of Fluid Mechanics，2006，38(1):483-512.

[5]WANG Chunxu，ZHANG Tao，HOU Guoxiang.Noise prediction of submerged free jet based on theory of vortex sound[J].Journal of Ship Mechanics，2010，14(6):670-677.

[6]王福軍.計算流體動力學分析-CFD軟件原理與應用[M].北京:清華大學出版社，2004:116-123.WANG Fujun.Computational fluid dynamics analysis-principle and application of CFD software[M].Beijing:Tsinghua University Publisher，2004:116-123.

[7]張兆順，崔桂香，許春曉.湍流大渦數值模擬的理論和應用[M].北京:清華大學出版社，2008:57-59.ZHANG Zhaoshun，CUI Guixiang，XU Chunxiao.The theory and application of large eddy simulation of turbulent flows[M].Beijing:Tsinghua University Publisher，2008:57-59.

[8]GERMANO M，P IOMELL I U，MO IN P，et al.Dynamic sub-grid-scale eddy viscosity model[J].Physics of Fluid A，1991，3(7):1760-1765.

[9]LILLY D K.A proposed modification of the Germano subgrid-scale closure model[J].Physics of Fluid，1992，4(3):633-635.

[10]張廣超，宋文愛.基于遞歸法的FFT計算機仿真[J].國外電子測量技術，2008，27(6):9-11.ZHANG Guangchao，SONG Wenai.Computer simulating of FFT based on recursion[J].Foreign Electronic Measurement Technology，2008，27(6):9-11.

[11]譚子尤，張雅彬.離散傅里葉變換快速算法的研究與MATLAB算法實現[J].中國科技信息，2006(22):316-317，321.TAN Ziyou，ZHANG Yabin.Research of fast Fourier transformation and realization of MATLAB[J].China Science and Technology Information，2006(22):317-321.

[12]楊昌棋，秦樹人，何輝.基于FFT的虛擬實時噪聲倍頻程分析儀[J].測控技術，2000，19(9):25-27.YANG Changqi，QIN Shuren，HE Hui.Virtual real-time noise octave analyzer based on FFT[J].Measurement＆Control Technology，2000，19(9):25-27.