孔腔脈動壓力及其波數—頻率譜的大渦模擬研究

鄧玉清，張 楠

（中國船舶科學研究中心 水動力重點實驗室，江蘇 無錫214082）

孔腔脈動壓力及其波數—頻率譜的大渦模擬研究

鄧玉清，張 楠

（中國船舶科學研究中心 水動力重點實驗室，江蘇 無錫214082）

孔腔流動中包含著流動分離和失穩以及渦旋相互干擾等復雜的流動現象。孔腔渦旋流動引起的流體振蕩能夠引起脈動壓力的顯著增加從而產生強烈的噪聲，在工程實際中備受關注。湍流脈動壓力是流激噪聲的重要來源，也是湍流研究中的基礎性問題，對其進行數值計算研究是流聲耦合領域的重要內容，而湍流脈動壓力波數—頻率譜的構建更是該領域的技術難點。文章采用大渦模擬方法（LES）對孔腔脈動壓力進行了數值模擬，考察了四套網格和四種亞格子應力模型對計算結果的影響，并與試驗結果進行比較，驗證數值計算方法的可靠性。首先采用大渦模擬方法計算了孔腔的脈動壓力，并與中國船舶科學研究中心的空泡水筒試驗結果進行對比分析。接著詳細地分析孔腔脈動壓力，研究亞格子應力模型和網格數量對計算結果的影響。最后，對數值計算得到的脈動壓力多元陣列結果進行時間/空間Fourier變換，構建了三維脈動壓力波數-頻率譜。該文工作對今后流激結構振動噪聲的預報和流動控制研究奠定了基礎。

大渦模擬；孔腔；壁面脈動壓力；波數—頻率譜

0 引 言

飛機、導彈以及艦船、潛艇、魚雷等無論是在空氣中還是水中運動的航行體，在它們的外表面都會形成邊界層。邊界層前緣為層流邊界層，經轉捩后發展成為湍流邊界層。湍流和層流兩者的主要區別就在于湍流具有速度和壓力的脈動性。邊界層內隨機的脈動速度擾動就產生物面上隨機的脈動壓力[1]。這種非定常流動一方面會直接產生輻射噪聲，稱之為“流激噪聲”或“流致噪聲”；另一方面脈動壓力作為聲源激勵引起物面結構振動，從而產生二次輻射噪聲，稱之為“流激結構振動噪聲”。流體流經孔口時，在導邊脫出自由剪切層撞擊腔口隨邊產生壓力反饋，形成剪切層自持振蕩現象，脈動壓力急劇增大。劇烈的自持振蕩不但會使艇體結構振動，疲勞甚至斷裂，也會引發強烈的噪聲[2-5]。脈動壓力是分析預報和控制這兩類噪聲的重要參量，也可以說是控制噪聲的源頭。所以湍流脈動壓力一直是湍流中的研究重點，近年來國內外對孔腔流動開展了試驗與數值計算研究。

Roshko（1955）[6]在國際上首次進行孔腔流動試驗測量與分析，拉開了孔腔試驗研究的序幕。楊國晶（2009）[7]采用LES方法研究不同雷諾數、不同腔型、不同攻角、有無質量交換條件下陷入式孔腔的水動力特性、剪切層振蕩頻率特性及腔中渦體的運動頻率特性。與物理試驗的結果比較，數值計算結果較準確。張楠等（2010）[8]采用LES結合FW-H聲學類比方法，預報了五種不同尺寸的方形孔腔在水中的流動特征及流激噪聲，孔腔中的流譜、孔腔與載體上的渦量分布以及五個孔腔的輻射噪聲頻譜均與試驗結果接近。Ji等（2012）[9]采用LES方法計算三維前臺階和后臺階的流動。臺階高度取四個不同值：53%、13%、3.3%和0.83%邊界層厚度，用以研究不同高度下表面壓力的特性和壓力源的形成機理。Li等（2013）[10]采用隱式LES方法計算三維孔腔的超音速流動來探索孔腔自持振蕩的本質。算例有層流和湍流兩種來流入口狀態以及變化湍流邊界層厚度形成的不同狀態。結果驗證了反饋機制是孔腔自持振蕩的最主要機理，表明LES是研究孔腔流動物理本質的有效手段之一。張楠等（2015）[11]基于LES與Powell渦聲理論對孔腔流激噪聲進行數值模擬研究，并在中國船舶科學研究中心循環水槽進行試驗，測試了兩種帶格柵的方腔在不同水流速度時的流激噪聲。數值計算結果與試驗結果對比分析后，表明計算方法切實可行，計算結果可靠有效。

從以上可以看出，大渦模擬方法是研究孔腔流動的有效手段之一，本文將采用大渦模擬方法。

1 計算方法與模型

1.1 大渦模擬方法

大渦模擬的基本思想是通過濾波過程將湍流分為可解尺度運動和不可解尺度運動，進而得到大渦模擬的控制方程，直接計算求解大尺度脈動，而對小尺度運動利用亞格子應力模型計算。實現大渦模擬計算的第一步就是將大尺度渦和小尺度渦進行分離，通過某種方式將小尺度渦過濾掉，稱之為濾波（filtering operation）。

流動變量φ的濾波運算表達為卷積型積分式，定義如下：

對于不可壓縮流體，忽略質量力，假定濾波運算和求導運算可以交換，上述方程經過濾波的控制方程如下：

方程（3）和雷諾平均方程（RANS）有類似的形式，右端出現不封閉項：

亞格子應力表達式可寫成下式：

其中：Lij稱為Leonard應力，也稱作外射項（outscatter term），它由兩個可解尺度脈動間的相互作用產生，是大尺度渦之間的相互作用，以產生湍流小尺度渦；Cij稱為交叉應力，也稱作交叉項（cross term），它由可解尺度脈動和不可解尺度脈動的相互作用產生，表示大尺度渦和小尺度渦的能量傳遞作用，能量既可以由大尺度渦向小尺度渦傳遞，也可由小尺度渦向大尺度渦傳遞，總體平均作用結果為能量由大尺度渦傳向小尺度渦；Rij稱為Reynolds應力，也稱作反散射項或逆散射項（backscatter term），它由不可解尺度脈動的相互作用產生，表示小尺度渦間的相互作用，產生大尺度渦，并促使能量從小尺度渦反饋到大尺度渦。本文采用四種常用亞格子應力模型（SL、DSL、WALE、KET）進行計算，其中SL為原初Smagorinsky-Lilly模型，DSL為動態SL模型，WALE為壁面自適應的渦粘性模型，KET為動能輸運模型。

1.2 計算模型

如圖1所示，計算模型為三維孔腔。其上部為長方形槽道，與模型試驗保持一致，總長L=2 160 mm，寬度B=225 mm，高度H=225 mm。孔腔開口跨長為l=360 mm，寬度b=180 mm，高度h=120 mm。孔腔導邊距入口La=800 mm，孔腔隨邊距出口出Lb=1 000 mm，計算域沿流向左右對稱。在孔腔的壁面上設置6個監測點（p1點、p2點、p3點、p4點、p5點、p6點），用來監測脈動壓力等物理量，位置分布如圖2所示，其中p1點位于孔腔前緣，p2點、p3點位于孔口外部側壁面，p4點位于孔腔底部中心，p5點、p6點位于孔口后緣，與模型試驗保持一致。

邊界條件設置為速度入口（設定來流速度大小和方向），壓力出口（設定相對參考點的流體靜壓值）以及壁面無滑移粘附條件。

時間項采用二階隱式格式離散，動量方程采用二階迎風格式差分，壓力速度耦合采用SIMPLE算法。計算中時間步長Δt=5×10-5s。壁面y+≤5。采用FFT結合Hanning窗處理非定常信號的時間序列。

圖1 孔腔幾何模型以及計算域（單位：mm）Fig.1 The cavity model and computational domain(mm)

圖2 孔腔壁面6個監測點位置示意圖（單位：mm）Fig.2 Six monitoring points on wall of the cavity(mm)

2 計算方法驗證

2.1 網格數量收斂性研究

表1 四種網格數量Tab.1 Grid number of different meshes

圖3 孔腔的計算網格Fig.3 Computational meshes of cavity

孔腔脈動壓力的測量是在中國船舶科學研究中心的小型多功能高速空泡水筒中完成的，受試驗條件限制，僅能得到100 Hz-8 kHz范圍的可靠的測量結果。對上述模型的脈動壓力進行數值計算，入口水速為5.83 m/s，與試驗保持一致，水中參考聲壓為1 μPa。

本節詳細地分析網格數量對計算結果準確度的影響，圖4給出了6個監測點在不同網格數量下的脈動壓力聲壓級的1/3 Oct頻譜計算結果，并與試驗結果進行比較。

總體定性來看，四種網格的計算結果在譜型上均表現出孔腔脈動壓力頻譜的典型特性，即能量主要分布在低頻段，隨著頻率的增加，譜級減少，高頻段的聲壓級低于低頻段的。再定量地分析，網格1和網格2由于網格數量較少，網格比較粗糙，計算得到的脈動壓力聲壓級頻譜圖與試驗結果相比，差別較大。隨著網格加密，網格數量增加，越來越多的更小尺度的渦被辨識出來，網格3和網格4的預報結果更加精確。

表2給出了四種不同網格數量的孔腔各個監測點脈動壓力總聲級計算偏差，負號表示計算結果小于試驗結果。從表中我們可以看出，網格1的平均偏差幅值為10.4 dB，網格2的平均偏差幅值為8.9 dB，網格3的平均偏差幅值為4.9 dB，網格4的平均偏差幅值為4.6 dB，說明網格3和網格4總體精度比網格1和網格2都要好，隨著網格數量的增加，偏差減小。網格4的計算偏差較穩定，計算結果是令人滿意的。p1處由于還未發生流動分離，流動形式與平板流動相似，其聲壓級比其他各點均低。綜上所述，6個監測點的脈動壓力聲壓級譜型和幅值與試驗結果較吻合。因此，下面進行亞格子應力模型計算時，采用網格4。

圖4 不同網格數量下脈動壓力頻譜計算結果與試驗對比Fig.4 Comparison between the calculated frequency spectra of wall pressure fluctuations with different meshes and measurement

表2 不同網格數量孔腔脈動壓力總聲級計算與試驗對比（單位：dB，頻段：100 Hz-8 kHz）Tab.2 Comparison between calculated OASPL and measured OASPL with different meshes(dB，100 Hz-8 kHz)

2.2 亞格子應力模型影響性研究

本節詳細地分析亞格子應力模型對計算結果的影響。圖5是不同亞格子應力模型脈動壓力頻譜計算結果，并與試驗進行對比。

圖5給出了6個監測點的脈動壓力聲壓級計算結果，四種亞格子應力模型計算結果得到的譜型與試驗結果較相符，脈動壓力隨著頻率的增加而衰減，能量集中在低頻段。高頻段的計算偏差要比低頻段的計算偏差大，這是由于脈動壓力高頻段的頻譜特征主要受小尺度渦的影響。小尺度渦的產生、輸運、演化并且對周圍流場產生強烈相互作用，是一個很復雜的過程。不同性質的小渦用同一種亞格子渦模型進行模擬，會引起一定的偏差。

圖5 不同亞格子應力模型脈動壓力頻譜的計算結果與試驗對比Fig.5 Comparison between the calculated frequency spectra of wall pressure fluctuations with sub-grid models and measurement

亞格子渦模型均為湍流假設下的半理論半經驗模型，模型中參數的選取往往借助一些較簡單的基礎性試驗，采用不同亞格子渦模型計算得到的孔腔脈動壓力結果也產生明顯的差別。表3給出了不同亞格子應力模型計算得到的各個監測點孔腔脈動壓力總聲級的偏差，負號表示計算結果小于試驗結果。從表中可以看出，SL模型的平均偏差幅值為4.6 dB，最大不超過9 dB；DSL模型的平均偏差幅值為15.2 dB，KET模型的平均偏差幅值為9.4 dB，WALE模型的平均偏差幅值為11.1 dB。對于本文的計算頻段，四種亞格子應力模型中SL模型的計算結果是最令人滿意的。

表3 不同亞格子應力模型孔腔脈動壓力總聲級計算與試驗對比（單位：dB，頻段：100 Hz-8 kHz）Tab.3 Comparison between calculated OASPL and measured OASPL with different sub-grid models(dB,100 Hz-8 kHz)

綜合以上分析可以看出：

（1）對孔腔壁面脈動壓力進行預報和定量分析可知，隨著網格加密，計算偏差減小，500萬和1 415萬兩套數量的精細網格的預報準確度明顯提高。

（2）采用不同亞格子應力模型對模擬孔腔流動的計算結果也產生明顯的差別。通過實際計算可知SL亞格子應力模型的效果較好。

3 孔腔脈動壓力波數—頻率譜的計算

前面已對脈動壓力計算方法進行了驗證，確定了計算網格數量與亞格子應力模型。本節在前面研究的基礎上利用數值模擬方法構建孔腔脈動壓力的波數—頻率譜。

湍流脈動壓力波數頻率譜是計算流激結構振動噪聲的重要輸入，其構建技術是流聲耦合基礎研究領域的技術難點。傳統的方法是利用模型試驗進行脈動壓力的傳感器陣列式聯合測量得到波數頻率譜，試驗花費大且非常耗時，利用數值計算方法獲得波數—頻率譜的研究是對傳統試驗方法的拓展，目前還非常稀少。張曉龍（2014）[16]利用大渦模擬方法對平板壁面脈動壓力波數—頻率譜進行計算，得到的譜級、遷移速度和含能渦分布等主要參數，與試驗結果相吻合。Abraham和Keith(1998)[17]將脈動壓力波數頻率譜定義為時間/空間相關函數在時間-空間的Fourier變換。在實際中，首先對時間-空間離散的脈動壓力進行離散Fourier變換（Discrete Fourier Transform，DFT），然后再求解幅值平方的系綜平均值得到波數-頻率譜。數學表達式如下：

其中：kx為縱向波數，ω為圓頻率，〈〉表示通過系綜平均得到的期望值，p（ xm,tn）表示第m個傳感器測在tn時刻的脈動壓力；N為時間結點數，M為傳感器數量；Δx為傳感器間距，Δt為時間步長；

圖6 傳感器陣列布置示意圖Fig.6 Distribution of sensors

本文在孔腔后緣中心線布置48元等間距線陣傳感器計算壁面脈動壓力，具體布置位置如圖6所示。Δl為第l個傳感器距孔腔后緣的距離，相鄰兩個傳感器的間距為4.22 mm。

圖7給出了水速為5.83 m/s時，在不同位置處（線陣a：Δl1=10 mm、線陣b：Δl2=30 mm）孔腔壁面脈動壓力的波數—頻率譜。左邊圖為計算得到的波數-頻率譜的三維視圖，譜級最高的部分稱為波數—頻率譜的“遷移脊”。定義遷移速度為Uc=dω/dk，即遷移脊的斜率表征遷移速度，表示湍流主要能量的平均遷移速度。遷移脊集中了大部分的湍流能量（含能渦），遷移速度反映能量輸運和傳遞的快慢，主要受外部流速影響。

圖7 不同位置的孔腔壁面脈動壓力波數—頻率譜Fig.7 Wavenumber-frequency spectra of cavity wall pressure fluctuations at different locations

從圖7中我們可以清晰地看出能量大部分集中在200Hz的低頻段范圍內，這與前面計算脈動壓力的結果相符。圖7（a）Δl1=10 mm，其歸一化譜峰值為-19.2 dB，顯示存在兩條脊，一條集中在波數為0的區域，另一條則是一般認為的遷移脊，其形式是傾斜的。計算得到遷移速度為4.23 m/s，無量綱化即為Uc/U0≈0.73；圖7（b）Δl2=30 mm，其歸一化譜峰值為-8.7 dB。可以看到僅有一條脊，計算得到遷移速度為4.23 m/s，無量綱化即為Uc/U0≈0.73。對比（a）、（b）圖，隨流動往下發展，遷移速度幾乎未發生變化。圖7（a）因為孔腔的存在，對流動產生阻礙作用，譜峰值比圖7（b）小，這是合理的。對比文獻[18]的計算結果，孔腔后緣波數-頻率譜的譜型和遷移脊特征與平板相似，說明隨著水流經孔口往下游發展，孔腔對流動的擾動作用減弱，流動形式與經典的平板流動相似。可以初步推測圖7（a）零波數脊的存在是由于渦旋結構與孔腔后緣撞擊，引起渦的振蕩和脫落而引起的。上述計算結果表明利用大渦模擬方法構建孔腔壁面脈動壓力波數—頻率譜是可靠的。

4 結 論

本文采用大渦模擬方法，對孔腔的壁面湍流脈動壓力及其波數—頻率譜進行計算。首先計算了孔腔流動脈動壓力的頻譜并進行了詳細的分析，從而研究網格數量和四種亞格子渦模型對計算準確度定性和定量的影響，并構建孔腔脈動壓力的波數-頻率譜。得到結論主要如下：

（1）隨著網格數量增加，分辨率隨之提高，脈動壓力計算的準確度也相應提高。

（2）壁面湍流脈動壓力的預報精度與亞格子渦模型直接相關，建立可靠的數值方法是進行數值計算的關鍵。對于孔腔，在本文的計算頻段內，SL亞格子應力模型計算結果較優，計算準確度較高。

（3）構建的孔腔脈動壓力的波數-頻率譜的譜級、遷移脊特征、遷移速度和能量分布域等結果與經典理論和已有研究結果相比是合理的，說明用此數值方法建立波數—頻率譜是切實可行的。

綜上所述，本文所建立的基于大渦模擬方法的數值研究用于孔腔流動脈動壓力計算切實可行，為下一步流激結構振動噪聲預報和流動控制的研究奠定了基礎。

[1]Kim J,Hussain F.Propagation velocity of perturbations in turbulent channel[J].Physics of Fluids,1993,A5(3):695-698.

[2]Zhuang N,Alvi F S,Alkislar M B,Shih C.Supersonic cavity flows and their control[J].AIAA,2006,44(9):2118-2128.

[3]Rowley C W,Williams D R.Dynamics and control of high-Reynolds-number flow over open cavities[J].Ann.Rev.Fluid Mech,2006,38:251-276.

[4]Lawson S,Barakos G N.Review of numerical simulations for high-speed,turbulent cavity flows[J].Prog.Aerosp.Sci,2011,47(3):186-216.

[5]Cattafesta L N,Song Q,Williams D R,Rowley C W,Alvi F S.Active control of flow-induced cavity oscillations[J].Prog.Aerosp.Sci.,2008,44:479-502.

[6]Roshko A,Some measurements of flows in a rectangular cutout[R].NACA TN 3488,1955.

[7]楊國晶.陷落式腔體水動力特性研究[D].哈爾濱:哈爾濱工程大學,2009.Yang Guojing.Research on the hydrodynamic characters of the cave-in cavity[D].Harbin:Harbin Engineering University,2009.

[8]張 楠,沈泓萃,朱錫清,姚惠之,謝 華.三維孔腔流激噪聲的大渦模擬與聲學類比預報與驗證研究[J].船舶力學,2010,14(1-2):181-190.Zhang nan,Shen Hongcui,Zhu Xiqing,Yao Huizhi,Xie hua.Validation and prediction of flow induced noise of 3-dimensional cavity with large eddy simulation and acoustic analogy[J].Journal of Ship Mechanics,2010,14(1-2):181-190.

[9]Ji M,Wang M.Surface pressure fluctuations on steps immersed in turbulent boundary layers[J].Fluid Mech,2012,7(12):471-504.

[10]Li W,Nonomura T,Fujii K.On the feedback mechanism in supersonic cavity flows[J].Physics of Fluids,2013.

[11]張 楠,李 亞,王志鵬,王 星,張曉龍.基于LES和Powell渦聲理論的孔腔流激噪聲數值模擬研究[J].船舶力學,2015,19(11):1393-1408.Zhang nan,Li Ya,Wang Zhipeng,Wang Xing,Zhang Xiaolong.Numerical simulation on the flow induced noise of cavity by LES and Powell vortex sound theory[J].Journal of Ship Mechanics,2015,19(11):1393-1408.

[12]Smagorinsky J.General circulation experiments with primitive equations[J].Monthly Weather Review,1963,91(3):99-164.

[13]Germano M,Piomelli U,Cabot W H.A dynamic subgrid-scale eddy viscosity model[J].Phys.Physics of Fluids,1991,A(3):1760-1765.

[14]Lilly D K.A Proposed modification of the Germano Subgrid scale closure method[J].Physics of Fluids,1992,A(4):633-635.

[15]張 楠.孔腔流動和流激噪聲機理及耦合計算方法研究[D].無錫:中國船舶科學研究中心,2010.Zhang Nan.Research on mechanism and hybrid computation approach for cavity flow and flow induced noise[D].Wuxi:China Ship Scientific Research Center,2010.

[16]張曉龍.壁面湍流脈動壓力及其波數—頻率譜大渦模擬計算研究[D].無錫:中國船舶科學研究中心,2014.Zhang Xiaolong.Investigation of turbulent wall pressure fluctuation and its wavenumber-frequency spectrum using large eddy simulation[D].Wuxi:China Ship Scientific Research Center,2014.

[17]Abraham B M,Keith W L.Direct measurements of turbulent boundary layer wall pressure wavenumber-frequency spectra[J].Journal of Fluids Engineering,1998,120(3):29-39.

[18]張曉龍,張 楠,吳寶山.平板壁面湍流脈動壓力及其波數頻率譜的大渦模擬計算分析研究[J].船舶力學,2014,18(10):1151-1164.Zhang Xiaolong,Zhang Nan,Wu Baoshan.Computation of turbulent wall pressure fluctuation and its wavenumber-frequency spectrum using large eddy simulation[J].Journal of Ship Mechanics,2014,18(10):1151-1164.

Computation of wall pressure fluctuations and wavenumberfrequency spectrum of cavity using large eddy simulation

DENG Yü-qing,ZHANG Nan

(Key Laboratory of Hydrodynamics,China Ship Scientific Research Center,Wuxi 214082,China)

Cavity flow always contains phenomenon of flow separation,flow instability and vorties interaction,which can cause fluid oscillations.The fluid oscillations increase pressure fluctuations intensely and induce obvious noises.It is a serious concern in naval application.Wall pressure fluctuations beneath turbulent boundary layers are important sources of flow induced noise.The research of wall pressure fluctuation is one of the groundworks in turbulent research.Thus,numerical simulation of pressure fluctuations has become an important issue in the flow-acoustic coupling field and construction of wavenumber-frequency spectra of turbulent pressure fluctuations is the technological difficulty in turbulence research.In this paper,wall pressure fluctuations of cavity are calculated by using large eddy simulation(LES)and the effects of different grid numbers and four different sub-grid scale models on calculated results are discussed.The computed results are compared with experimental data to verify the reliability of the numerical method.First,wall pressure fluctuations of cavity are computed by LES.The computed results are compared with the experimental data in cavitation tunnel of CSSRC.Then,wall pressure fluctuations are analyzed particularly.The effects of different grid numbers and different sub-grid scale models are also discussed in detail.Finally,wavenumber-frequency spectra of wall pressure fluctuations of sensors array are calculated by timespace Fourier Transformation.It lays a foundation for further research in the prediction of flow induced vi-bration noise and flow control approach.

large eddy simulation(LES);cavity;wall pressure fluctuations;wavenumber-frequency spectrum

O35

A

10.3969/j.issn.1007-7294.2017.10.003

1007-7294（2017）10-1199-11

2017-06-19

鄧玉清（1992-），女，碩士研究生，E-mail:dengyuqing@hust.edu.cn； 張 楠（1977-），男，博士，研究員。