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

郝維，劉正先，陳麗英

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

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

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

1 數(shù)值分析方法

1.1 湍流模型

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

1.1.1 標(biāo)準(zhǔn) k-ε模型

標(biāo)準(zhǔn)k-ε模型是目前應(yīng)用較廣泛的兩方程湍流模型[6]，它對運動方程中的所有脈動量實施?；牧髡扯韧ㄟ^求解湍動能k方程和湍動耗散率ε方程得到，k和ε的求解方程如下:

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

1.1.2 大渦模擬

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

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

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

亞格子應(yīng)力具有下列形式:

1.2 壓力脈動分析方法

1.2.1 壓力級的傅里葉變換

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

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

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

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

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

數(shù)值模擬模型除與水槽實驗?zāi)Ｐ途哂邢嗤膸缀纬叨?、流場邊界條件外，還應(yīng)保持相同的數(shù)據(jù)處理方法.實驗采用1/3倍頻程數(shù)據(jù)處理方法，總壓力級等效為

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

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

2 計算模型

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

圖1 方腔模型結(jié)構(gòu)及邊界條件Fig.1 The cavity model and the boundary conditions

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

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

圖2 腔體網(wǎng)格分布Fig.2 The grid distribution for the cavity hole

3 數(shù)值模擬結(jié)果分析

圖3 監(jiān)測點位置Fig.3 The location of the reference points

3.1 流場分析

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

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

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

3.2 壓力脈動數(shù)值分析

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

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

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

圖 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 傅里葉壓力級

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

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

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

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

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

4.2 壓力級的數(shù)值與實驗數(shù)據(jù)比較

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

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

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

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

5 結(jié)論

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

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

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

[1]耿冬寒，劉正先.大渦模擬-Lighthill等效聲源法的空腔水動噪聲預(yù)測[J].哈爾濱工程大學(xué)學(xué)報，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]朱習(xí)劍，何祚鏞.水洞中突出矩形腔的流激駐波振蕩研究[J].哈爾濱工程大學(xué)學(xué)報，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，F(xiàn)REUND 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]王福軍.計算流體動力學(xué)分析-CFD軟件原理與應(yīng)用[M].北京:清華大學(xué)出版社，2004:116-123.WANG Fujun.Computational fluid dynamics analysis-principle and application of CFD software[M].Beijing:Tsinghua University Publisher，2004:116-123.

[7]張兆順，崔桂香，許春曉.湍流大渦數(shù)值模擬的理論和應(yīng)用[M].北京:清華大學(xué)出版社，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].國外電子測量技術(shù)，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算法實現(xiàn)[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].測控技術(shù)，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.