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

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

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

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

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

空腔；非線性；噪聲源；湍流;內(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)埋彈艙給飛行器帶來了隱身、超聲速巡航等益處，但同時也產(chǎn)生了一系列復(fù)雜的空氣動力學(xué)問題，包括剪切層極不穩(wěn)定、邊界層分離、激波邊界相互干擾等[2-3]，這加劇了空腔內(nèi)的非定常效應(yīng)，由此產(chǎn)生的氣動噪聲將對飛行器的性能與安全產(chǎn)生影響，嚴(yán)重時將導(dǎo)致飛行器產(chǎn)生災(zāi)難性的后果。因此對于空腔氣動噪聲的研究具有重大而緊迫的現(xiàn)實(shí)意義。

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

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

1 數(shù)值方法

(1)

式中，

(2)

(3)

忽略密度脈動項，對以上方程取時間平均可得：

(5)

(6)

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

(7)

(8)

2 計算模型

2.1 腔體構(gòu)型

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

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

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

2.2 計算網(wǎng)格

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

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

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

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

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

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

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

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

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

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

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

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

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

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

x方向

y方向

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

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

對于長深比L/H=5的開式空腔，來流氣體流經(jīng)前緣時，因為腔體深度較大，氣流未觸及空腔底部，剪切層跨越空腔中部與后壁發(fā)生碰撞，空腔前部和中部受剪切層影響較小，壓力不會發(fā)生大的變化，空腔后部壓力上升，對于超音速流動，會誘發(fā)激波產(chǎn)生?？涨坏牧鲃犹匦砸约扒粌?nèi)復(fù)雜的流動環(huán)境會導(dǎo)致空腔后部發(fā)生顫振，從而產(chǎn)生噪聲。噪聲通過腔內(nèi)循環(huán)氣流傳播到空腔前緣，導(dǎo)致剪切層分離，當(dāng)滿足一定相位條件時，形成聲波反饋循環(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三個監(jiān)測點(diǎn)NLAS仿真得到的聲壓頻譜曲線與試驗值對比圖。由于壓力脈動值在計算起始的一段時間內(nèi)不具有周期性，而在該段時間之后壓力脈動呈現(xiàn)出一定的周期特性并一直保持下去，而這正是所需要的壓力脈動數(shù)據(jù)，也是噪聲傳播對應(yīng)的數(shù)據(jù)。為了防止起始時間的不規(guī)則壓力脈動影響噪聲求解，將起始段的壓力脈動截掉，使用0.1s～0.4s的脈動數(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)測點(diǎn)聲壓頻譜特性計算與試驗對比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)測點(diǎn)聲壓頻譜特性計算與試驗對比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)測點(diǎn)聲壓頻譜特性計算與試驗對比Fig.13 Comparison of spectrum between calculation and test at Ma=1.35

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

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

表1 最高聲壓級比較

表2 主頻位置比較

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

4 結(jié)論

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

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

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

References)

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

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

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

[3] 吳繼飛，羅新福，范召林．內(nèi)埋式彈艙流場特性及武器分離特性改進(jìn)措施[J]．航空學(xué)報，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—)，男，四川樂山人，博士研究生，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