王一丁，陳濱琦，郭 亮，鐘范俊，童明波

空腔噪聲非線性數(shù)值模擬*

王一丁1，陳濱琦1，郭 亮2，鐘范俊2，童明波1

(1.南京航空航天大學(xué) 航空宇航學(xué)院， 江蘇 南京 210016； 2.成都飛機(jī)設(shè)計(jì)研究所， 四川 成都 610091)

將雷諾平均N-S方程與非線性噪聲求解方法相結(jié)合，對(duì)M219空腔在Ma=0.6，Ma=0.85,Ma=1.35條件下進(jìn)行了氣動(dòng)噪聲分析。通過(guò)雷諾平均N-S方程求解空腔流場(chǎng)，得到包含空腔平均流場(chǎng)基本特征以及強(qiáng)制設(shè)定的湍流脈動(dòng)統(tǒng)計(jì)描述的初始湍流統(tǒng)計(jì)平均解，采用非線性噪聲求解方法重構(gòu)噪聲源并高精度模擬壓力脈動(dòng)的傳播。通過(guò)與試驗(yàn)結(jié)果對(duì)比表明非線性噪聲求解方法能夠較好地捕捉空腔流動(dòng)中的壓強(qiáng)脈動(dòng)及噪聲水平。與分離渦模擬方法相比，非線性噪聲求解方法在保持計(jì)算精度的同時(shí)大大減少計(jì)算網(wǎng)格，對(duì)內(nèi)埋彈艙快速設(shè)計(jì)具有一定的參考意義。

空腔；非線性；噪聲源；湍流;內(nèi)埋彈艙

(1.CollegeofAerospaceEngineering,NanjingUniversityofAeronauticsandAstronautics,Nanjing210016,China;

2.ChengduAircraftDesignandResearchInstitute,Chengdu610091,China)

空腔結(jié)構(gòu)廣泛存在于飛行器中，新一代戰(zhàn)斗飛行器使用內(nèi)埋彈艙能夠使飛行器的阻力最大降低30%，雷達(dá)橫截面最低可達(dá)0.07～0.12m2[1]。內(nèi)埋彈艙給飛行器帶來(lái)了隱身、超聲速巡航等益處，但同時(shí)也產(chǎn)生了一系列復(fù)雜的空氣動(dòng)力學(xué)問(wèn)題，包括剪切層極不穩(wěn)定、邊界層分離、激波邊界相互干擾等[2-3]，這加劇了空腔內(nèi)的非定常效應(yīng)，由此產(chǎn)生的氣動(dòng)噪聲將對(duì)飛行器的性能與安全產(chǎn)生影響，嚴(yán)重時(shí)將導(dǎo)致飛行器產(chǎn)生災(zāi)難性的后果。因此對(duì)于空腔氣動(dòng)噪聲的研究具有重大而緊迫的現(xiàn)實(shí)意義。

Le等[4]采用大渦模擬方法(LargeEddySimulation，LES)、Allen等[5]采用分離渦模擬(DetachedEddySimulation，DES)以及Peng等[6]采用混合雷諾平均大渦模擬(Reynolds-AveragedNavier-Stokes/LES，RANS/LES)對(duì)空腔噪聲進(jìn)行了數(shù)值模擬研究。但對(duì)于工程實(shí)際問(wèn)題，這些方法存在計(jì)算網(wǎng)格數(shù)量大、計(jì)算時(shí)間長(zhǎng)等缺點(diǎn)。目前噪聲領(lǐng)域常用的DES方法經(jīng)常與有限傳輸算法結(jié)合使用，增加了亞格子尺度模型的耗散，有可能導(dǎo)致有效粘度過(guò)大，同時(shí)統(tǒng)計(jì)學(xué)湍流能量的傳輸也存在很大困難，這極大地限制了DES方法適用的流動(dòng)及網(wǎng)格類(lèi)型。

Batten等于2002年提出了一種非線性噪聲求解(NonlinearAcousticSolver,NLAS)方法，該方法通過(guò)對(duì)湍流物理量進(jìn)行重構(gòu)兼顧了亞格子尺度聲源的影響，在保持計(jì)算精度的同時(shí)降低了網(wǎng)格需求[7-9]。王一丁等將NLAS方法引入內(nèi)埋彈艙噪聲預(yù)測(cè)中，計(jì)算了英國(guó)QinetiQ公司M219空腔[10]，空腔計(jì)算條件為Ma=0.6，Ma=0.85，Ma=1.35，將NLAS仿真結(jié)果與試驗(yàn)數(shù)據(jù)以及Allen[5]使用DES方法計(jì)算結(jié)果進(jìn)行了對(duì)比，驗(yàn)證了NLAS方法用于空腔噪聲預(yù)測(cè)的有效性與準(zhǔn)確性。NLAS方法計(jì)算網(wǎng)格數(shù)相對(duì)DES方法大幅減少，降低了計(jì)算成本，具有一定的工程應(yīng)用價(jià)值。

1 數(shù)值方法

(1)

式中，

(2)

(3)

忽略密度脈動(dòng)項(xiàng)，對(duì)以上方程取時(shí)間平均可得：

(5)

(6)

式中，Ri是標(biāo)準(zhǔn)雷諾應(yīng)力張量和湍流熱通量相關(guān)項(xiàng)。求解噪聲的關(guān)鍵是通過(guò)RANS計(jì)算求得到這些未知項(xiàng)，不能求解的小尺度量則通過(guò)RANS計(jì)算得到的湍流統(tǒng)計(jì)結(jié)果重構(gòu)出來(lái)，以此生成亞格子源項(xiàng)。Batten提出的湍流重構(gòu)方法為：

(7)

(8)

2 計(jì)算模型

2.1 腔體構(gòu)型

M219模型為典型的開(kāi)式空腔，在QinetiQ風(fēng)洞進(jìn)行了一系列風(fēng)洞試驗(yàn)，腔體長(zhǎng)深比L/D=5，寬深比W/D=1.0。圖1為腔體在DERA風(fēng)洞試驗(yàn)的照片，圖2為M219空腔風(fēng)洞試驗(yàn)件構(gòu)型圖。

圖1 M219空腔在DERA風(fēng)洞噪聲試驗(yàn)Fig.1 Noise test of M219 cavity in wind tunnel

圖2 M219空腔風(fēng)洞試驗(yàn)件構(gòu)型圖Fig.2 Sketch of M219 wind tunnel test pieces

2.2 計(jì)算網(wǎng)格

圖3為RANS計(jì)算網(wǎng)格區(qū)域，整個(gè)計(jì)算域由4個(gè)結(jié)構(gòu)塊組成，來(lái)流方向(x)從-8D到12D，展向從-2D到2D，壁面法向方向從-1D到4D。為了更好地模擬湍流脈動(dòng)的統(tǒng)計(jì)平均結(jié)果，RANS計(jì)算應(yīng)采用非線性的各向異性湍流模型，選取cubick-epsilon模式,該模式通過(guò)矩陣近似各個(gè)位置對(duì)渦粘系數(shù)的影響，更加符合物理本質(zhì)。流場(chǎng)物面第一層網(wǎng)格尺度為5×10-2mm，網(wǎng)格數(shù)量為260萬(wàn)，基于空腔深度的雷諾數(shù)Re=7×106。RANS求解的超聲速時(shí)來(lái)流為固定超聲速來(lái)流入口條件，遠(yuǎn)場(chǎng)為特征線條件，出口邊界不指定，由內(nèi)層網(wǎng)格物理量推得。亞聲速時(shí)來(lái)流入口、遠(yuǎn)場(chǎng)及出口位置均采用特征線邊界條件。

圖3 RANS計(jì)算區(qū)域Fig.3 Computational domain of RANS

噪聲計(jì)算采用單獨(dú)的計(jì)算網(wǎng)格，物面邊界為黏性無(wú)滑移絕熱壁，采用壁面函數(shù)法求解保證物面區(qū)域求解精度。RANS計(jì)算得到當(dāng)?shù)乩字Z應(yīng)力張量和熱通量的統(tǒng)計(jì)平均值，將它插值到噪聲計(jì)算網(wǎng)格，根據(jù)這一統(tǒng)計(jì)平均結(jié)果對(duì)湍流進(jìn)行人工重構(gòu)。NLAS方法的優(yōu)勢(shì)在于噪聲求解器可以在各項(xiàng)同性更好的網(wǎng)格單元上進(jìn)行計(jì)算，特別是在近壁面區(qū)域。在進(jìn)行噪聲計(jì)算時(shí)，計(jì)算域選取為包含噪聲源周邊的區(qū)域。圖4為RANS計(jì)算得到的最大湍流動(dòng)能kmax的10%等值面，這部分區(qū)域是湍流脈動(dòng)最為劇烈區(qū)域、也是主要的噪聲源區(qū)域、噪聲網(wǎng)格重點(diǎn)關(guān)注區(qū)域。新的邊界被設(shè)置為吸收層邊界，它的遠(yuǎn)場(chǎng)及衰減層數(shù)據(jù)由之前RANS計(jì)算提供。

圖4 RANS計(jì)算10% kmax等值面(背景為來(lái)流速度)Fig.4 Iso-surface of 10% kmax calculated by RANS(shaded with streamwise velocity)

由于近壁面網(wǎng)格要求放寬以及計(jì)算域的縮小，噪聲計(jì)算網(wǎng)格數(shù)量為120萬(wàn)，較RANS計(jì)算的260萬(wàn)網(wǎng)格有了顯著減少。RANS與NLAS計(jì)算的網(wǎng)格對(duì)比如圖5所示。

圖5 RANS網(wǎng)格與噪聲計(jì)算網(wǎng)格對(duì)比Fig. 5 Comparison of RANS and acoustics meshes

RANS計(jì)算控制方程采用有限體積法求解，無(wú)黏項(xiàng)采用二階精度TVD格式離散，黏性項(xiàng)采用中心差分格式離散，時(shí)間推進(jìn)采用隱式方法。NLAS計(jì)算空間和時(shí)間離散格式與RANS計(jì)算相同，時(shí)間步長(zhǎng)Δt=2×10-5s，共計(jì)算20 000步。

來(lái)流馬赫數(shù)與M219風(fēng)洞試驗(yàn)一致，取Ma=0.6，Ma=0.85，Ma=1.35，覆蓋了亞、跨、超聲速以充分驗(yàn)證NLAS方法在各種來(lái)流條件下模擬空腔噪聲的有效性與準(zhǔn)確性。在空腔底面中心線處設(shè)置10個(gè)點(diǎn)記錄壓力的變化，分別表示為K20～K29，具體位置如圖6所示。

圖6 脈動(dòng)壓力監(jiān)測(cè)點(diǎn)位置Fig.6 Monitoring locations of oscillating pressure

計(jì)算使用4個(gè)計(jì)算機(jī)節(jié)點(diǎn)，每個(gè)計(jì)算機(jī)節(jié)點(diǎn)包含1個(gè)8核2.6GHz處理器和24G內(nèi)存。

3 計(jì)算結(jié)果對(duì)比分析

空腔底部監(jiān)測(cè)點(diǎn)壓強(qiáng)均方根是測(cè)量脈動(dòng)壓力的常用指標(biāo)。圖7～9為Ma=0.6，Ma=0.85，Ma=1.35條件下由NLAS方法計(jì)算得到的壓強(qiáng)均方根值與QinetiQ風(fēng)洞試驗(yàn)值以及Allen等[5]采用DES方法計(jì)算的結(jié)果對(duì)比。

圖9 Ma=1.35壓強(qiáng)均方根試驗(yàn)、DES及NLAS對(duì)比圖Fig.9 Comparison of prms between experiment, DES and NLAS at Ma=1.35

通過(guò)圖7～9可以看到NLAS計(jì)算得到的均方根值略大于試驗(yàn)值，與DES方法計(jì)算值精度基本相當(dāng)。而NLAS方法所用網(wǎng)格數(shù)量?jī)H為120萬(wàn)，DES方法所用網(wǎng)格數(shù)量為260萬(wàn)，在保證計(jì)算精度的同時(shí)，NLAS方法大大減少了網(wǎng)格需求，縮短了計(jì)算時(shí)間。

x方向

y方向

z方向圖10 NLAS方法得到的渦量分量瞬時(shí)等值面Fig.10 Instantaneous streamwise vorticity iso-surfaces calculated by NLAS

圖10為空腔在Ma=0.85下x，y，z方向渦量瞬時(shí)等值面，渦量相同均為6×103，背景代表流向速度。

對(duì)于長(zhǎng)深比L/H=5的開(kāi)式空腔，來(lái)流氣體流經(jīng)前緣時(shí)，因?yàn)榍惑w深度較大，氣流未觸及空腔底部，剪切層跨越空腔中部與后壁發(fā)生碰撞，空腔前部和中部受剪切層影響較小，壓力不會(huì)發(fā)生大的變化，空腔后部壓力上升，對(duì)于超音速流動(dòng)，會(huì)誘發(fā)激波產(chǎn)生。空腔的流動(dòng)特性以及腔內(nèi)復(fù)雜的流動(dòng)環(huán)境會(huì)導(dǎo)致空腔后部發(fā)生顫振，從而產(chǎn)生噪聲。噪聲通過(guò)腔內(nèi)循環(huán)氣流傳播到空腔前緣，導(dǎo)致剪切層分離，當(dāng)滿足一定相位條件時(shí)，形成聲波反饋循環(huán)，發(fā)生腔內(nèi)自持振蕩。

圖11～13分別為Ma=0.6，Ma=0.85，Ma=1.35下，位于空腔底部中心線X/L=0.25，X/L=0.55，X/L=0.95三個(gè)監(jiān)測(cè)點(diǎn)NLAS仿真得到的聲壓頻譜曲線與試驗(yàn)值對(duì)比圖。由于壓力脈動(dòng)值在計(jì)算起始的一段時(shí)間內(nèi)不具有周期性，而在該段時(shí)間之后壓力脈動(dòng)呈現(xiàn)出一定的周期特性并一直保持下去，而這正是所需要的壓力脈動(dòng)數(shù)據(jù)，也是噪聲傳播對(duì)應(yīng)的數(shù)據(jù)。為了防止起始時(shí)間的不規(guī)則壓力脈動(dòng)影響噪聲求解，將起始段的壓力脈動(dòng)截掉，使用0.1s～0.4s的脈動(dòng)數(shù)據(jù)。

(a) K22(X/L=0.25)

(b) K25(X/L=0.55)

(c) K29(X/L=0.95)圖11 Ma=0.6監(jiān)測(cè)點(diǎn)聲壓頻譜特性計(jì)算與試驗(yàn)對(duì)比Fig.11 Comparison of spectrum between calculation and test at Ma=0.6

(a) K22(X/L=0.25)

(b) K25(X/L=0.55)

(c) K29(X/L=0.95)圖12 Ma=0.85監(jiān)測(cè)點(diǎn)聲壓頻譜特性計(jì)算與試驗(yàn)對(duì)比Fig.12 Comparison of spectrum between calculation and test at Ma=0.85

(a) K22(X/L=0.25)

(b) K25(X/L=0.55)

(c) K29(X/L=0.95)圖13 Ma=1.35監(jiān)測(cè)點(diǎn)聲壓頻譜特性計(jì)算與試驗(yàn)對(duì)比Fig.13 Comparison of spectrum between calculation and test at Ma=1.35

計(jì)算結(jié)果表明，不同馬赫數(shù)下不同位置的模態(tài)一致，這與文獻(xiàn)中給出的典型頻譜符合很好。從頻譜模態(tài)中可以看出，采用本方法進(jìn)行氣動(dòng)噪聲計(jì)算，前4階頻譜模態(tài)均可以被捕捉到，且除個(gè)別模態(tài)略有差別外，主頻均被精確捕捉到，說(shuō)明本方法具有較高精度。Ma=1.35時(shí)，噪聲測(cè)點(diǎn)聲壓級(jí)分布集中在130～170dB之間，上游聲壓級(jí)略低，而下游聲壓級(jí)略高。對(duì)于類(lèi)似結(jié)構(gòu)的內(nèi)埋彈艙而言，如此高的聲壓級(jí)會(huì)對(duì)艙體以及艙內(nèi)武器造成疲勞損傷，且這種分布會(huì)使艙內(nèi)武器生成一定的抬頭力矩，這主要是由于上游剪切層與下游壁面邊界層互相作用產(chǎn)生不穩(wěn)定壓力波，該不穩(wěn)定壓力波主要集中于下游區(qū)域，并且從下游沿壁面向上游傳播至上游前緣，再與剪切層互相作用使之與壁面分離從而形成聲學(xué)反饋。在空腔中主要的不穩(wěn)定區(qū)域集中于下游，使得下游噪聲聲壓級(jí)明顯高于上游。

氣動(dòng)噪聲計(jì)算中最重要是主頻位置及其對(duì)應(yīng)的最高聲壓級(jí)的預(yù)測(cè)，重點(diǎn)對(duì)這兩個(gè)參數(shù)的仿真與試驗(yàn)結(jié)果進(jìn)行了對(duì)比。表1為Ma=0.6，Ma=0.85，Ma=1.35下空腔底部前中后三個(gè)位置的數(shù)值仿真結(jié)果與試驗(yàn)結(jié)果中的最高聲壓級(jí)對(duì)比，表2為主頻對(duì)比。

表1 最高聲壓級(jí)比較

表2 主頻位置比較

通過(guò)對(duì)比不同馬赫數(shù)下NLAS仿真與試驗(yàn)的最高聲壓級(jí)以及主頻位置，可以發(fā)現(xiàn)NLAS方法能較為精確地模擬出亞、跨、超聲速情況下空腔噪聲主頻及其對(duì)應(yīng)的最高聲壓級(jí)。

4 結(jié)論

1)將非線性噪聲求解方法NLAS應(yīng)用于空腔噪聲預(yù)測(cè)，模擬了Ma=0.6，Ma=0.85，Ma=1.35三種來(lái)流條件下的空腔噪聲。應(yīng)用cubick-ε湍流模型，遠(yuǎn)場(chǎng)吸收邊界及壁面函數(shù)法計(jì)算得到了三種來(lái)流條件下空腔噪聲特性，數(shù)值結(jié)果與試驗(yàn)結(jié)果基本吻合，特別是準(zhǔn)確模擬主頻及其對(duì)應(yīng)的最高聲壓級(jí)，表明NLAS方法在亞、跨、超聲速條件下對(duì)空腔噪聲有較好的預(yù)測(cè)能力。

2)非線性噪聲求解方法對(duì)于近壁面網(wǎng)格要求低，聲場(chǎng)計(jì)算域比RANS小，可減少噪聲計(jì)算網(wǎng)格數(shù)量，降低計(jì)算成本。將NLAS計(jì)算結(jié)果與國(guó)外文獻(xiàn)中DES計(jì)算結(jié)果進(jìn)行對(duì)比，NLAS的計(jì)算精度與DES相當(dāng)，但是網(wǎng)格數(shù)大大降低，因此NLAS方法對(duì)于內(nèi)埋彈艙工程快速設(shè)計(jì)具有一定的意義。

3) 對(duì)于L/D=5的典型開(kāi)式空腔，通過(guò)對(duì)亞、跨、超聲速情況下空腔噪聲數(shù)值計(jì)算與試驗(yàn)對(duì)比，隨著馬赫數(shù)的增大，各監(jiān)測(cè)點(diǎn)的噪聲主頻位置，總聲壓級(jí)都有所增大。

References)

[1]MicharelJH. 戰(zhàn)術(shù)導(dǎo)彈空氣動(dòng)力學(xué)[M].洪金森,楊其德，毛國(guó)良，譯.北京：宇航出版社,1999.

MicharelJH.Tacticalmissileairdynamics[M].TranslatedbyHONGJinseng,YANGQide,MAOGuoliang.Beijing:AstronavigationPress, 1999. (inChinese)

[2]MurrayNE,UkeileyLS.Flowfielddynamicsinopencavityflows[C]//Proceedingsof12thAIAA/CEASAeroacousticsConference,AIAA2006-2428, 2006.

[3] 吳繼飛，羅新福，范召林．內(nèi)埋式彈艙流場(chǎng)特性及武器分離特性改進(jìn)措施[J]．航空學(xué)報(bào)，2009，30(10)：1840-1845．WUJifei,LUOXinfu,FANZhaolin.Flowcontrolmethodtoimprovecavityflowandstoreseparationcharacteristics[J].ActaAeronauticaetAstronauticaSinica，2009，30(10)：1840-1845. (inChinese)

[4]LeTH,MaryI,TerracolM.LESofpressureloadssuppressioninweaponsbayflow[C]//Proceedingsof43rdAIAAAerospaceSciencesMeetingandExhibit,AIAA2005-794, 2005.

[5]AllenR,Mendon?aF.DESvalidationsofcavityacousticsoverthesubsonictosupersonicrange[C]//Proceedingsof10thAIAA/CEASAeroacousticsConference,AIAA2004-2862, 2004.

[6]PengSH,HaaseW.AdvancesinhybridRANS-LESmodelling,notesonnumericalfluidmechanicsandmultidisciplinarydesign[M].USA:Springer, 2008 : 132-141.[7]BattenP，RibaldoneE,CasellaM,etal.Towardsageneralizednon-linearacousticssolver[C]//Proceedingsof10thAIAA/CEASAeroacousticsConference,AIAA2004-3001, 2004.

[8]BattenP，GoldbergU,ChakravarthyS.Reconstructedsub-gridmethodsforacousticspredictionsatallReynoldsnumbers[C]//Proceedingsof8thAIAA/CEASAeroacousticsConference&Exhibit,AIAA2002-2511, 2002.

[9]DaSilvaCRI,DeAlmeidaO,BattenP.Investigationofanaxi-symmetricsubsonicturbulentjetusingcomputationalaeroacousticstools[C]//Proceedingsof13thAIAA/CEASAeroacousticsConference,AIAA2007-3656, 2007.

[10]HenshawMJ.M219cavitycase:verificationandvalidationdataforcomputationalunsteadyaerodynamics[R].ResearchandTechnologyOrganization,RTO-TR-26,AC/323(AVT)TP/19, 2002.

Nonlinear numerical simulation of cavity noise

WANG Yiding1,CHEN Binqi1, GUO Liang2, ZHONG Fanjun2,TONG Mingbo1

InordertoevaluatetheM219cavitynoiseat0.6, 0.85and1.35Machnumber,nonlinearacousticsolveriscombinedwithReynolds-averagedNavier-Stokesequations.TheflowfieldofacavityiscalculatedbymeansofReynolds-averagedNavier-Stokesequations,whichcontainsbasiccharacteristicsofaverageflowfieldandturbulencestatisticalaveragesolutionofstatisticsdescriptionofturbulencefluctuation.Noisesourceisrefactoredbythenonlinearacousticsolver.Spreadofpressurefluctuationissimulatedprecisely.Acomparisonshowsthatthesimulationresultsofnonlinearacousticsolveragreewellwiththeexperimentresults.Comparedwithdetachededdysimulation，nonlinearacousticsolvercangreatlyreducetheamountofmesh.Inaddition,themethodcanprovidesomereferenceforinternalweaponsbaydesign.

cavity;nonlinearity;source;turbulence;internalweaponsbay

2014-02-10

王一丁(1985—)，男，四川樂(lè)山人，博士研究生，E-mail：wyding127@163.com；童明波(通信作者)，男，教授，博士，博士生導(dǎo)師，E-mail: tongw@nuaa.edu.cn

10.11887/j.cn.201504025

http://journal.nudt.edu.cn

V

A