脈沖爆轟發(fā)動機爆轟噪聲場區(qū)域劃分及其數(shù)值模擬

黃孝龍，許桂陽，翁春生，李寧

(南京理工大學 瞬態(tài)物理國家重點實驗室，江蘇 南京210094)

0 引言

脈沖爆轟發(fā)動機(PDE)是一種新型發(fā)動機，與傳統(tǒng)發(fā)動機相比，具有燃燒效率高，結(jié)構(gòu)簡單等優(yōu)點。但是，與其他發(fā)動機一樣，PDE 在工作過程中會產(chǎn)生較大噪聲，對相關(guān)工作人員和設(shè)備的正常工作會產(chǎn)生較大影響。因此，對其噪聲的研究，具有重要意義。

Boesch 等［1］測量了PDE 管外軸線上不同距離處的爆轟噪聲，并通過理想沖擊波理論，較為準確地預測出了爆轟波的爆轟噪聲峰值和頻譜變化。Glaser 等［2］對PDE 出口不同距離、不同方位角處的噪聲進行了測量，發(fā)現(xiàn):隨著出口距離的增加，聲壓峰值逐漸減小;并且沖擊波傳播速度在出爆轟管口很短的距離內(nèi)，迅速衰減為聲速;聲壓輻射曲線，呈現(xiàn)出明顯的指向性;并通過定義參考半徑，將爆轟發(fā)動機聲場分為強沖擊波區(qū)和弱沖擊波區(qū)，分別給出了PDE 爆轟噪聲峰值在參考半徑內(nèi)外的變化規(guī)律，與實驗測量所得噪聲峰值變化規(guī)律一致。He 等［3］對PDE 的爆轟噪聲做了初步的數(shù)值模擬研究，計算了爆轟管內(nèi)和管外不同位置的聲壓級，并通過傅里葉轉(zhuǎn)換初步得到了爆轟噪聲聲壓級的估計值，峰值聲壓常出現(xiàn)在頻率為330 Hz 范圍。Eerden 等［4］在求解點爆炸產(chǎn)生的沖擊波由近場強非線性區(qū)向遠場線性區(qū)轉(zhuǎn)變的過渡過程中，提出劃分計算區(qū)域的思想，采用不同的方法求解不同區(qū)域的控制方程，并與采用一種計算方法得到的結(jié)果對比，得出區(qū)域劃分可以更好地模擬點爆炸聲場傳播特性的結(jié)論。

PDE 爆轟噪聲與點爆炸相似，都會產(chǎn)生高強度的沖擊波。在此前研究的基礎(chǔ)上，本文嘗試將劃分計算區(qū)域思想應(yīng)用到PDE 爆轟噪聲的研究中。根據(jù)沖擊波強弱，將PDE 爆轟噪聲場劃分為3 個區(qū)域，并采用相應(yīng)的計算方法對PDE 爆轟噪聲場的傳播過程進行數(shù)值模擬，得到爆轟噪聲的分布規(guī)律。

1 PDE 爆轟噪聲場的分區(qū)及計算方法

PDE 爆轟噪聲主要由沖擊波噪聲和噴流噪聲組成，即爆轟波從爆轟管口逸出并退化為沖擊波引起的沖擊噪聲以及緊隨其后的高溫高壓燃氣射流引起的射流噪聲。為了便于計算與分析，爆轟噪聲由近場到遠場的傳播過程劃分為3 個區(qū)域，如圖1所示:第1 區(qū)域為近場的強非線性區(qū)，主要由沖擊波和高速射流引起的極強非線性噪聲組成;第2 個區(qū)域為近場與遠場之間的過渡區(qū)，即聲波由強非線性聲波向線性聲波過渡的區(qū)域;第3 個區(qū)域為遠場線性聲波區(qū)。

根據(jù)文獻［5］中超壓在6.0×103～1.96×104Pa 可作為沖擊波與沖擊聲波過渡邊界的結(jié)論，結(jié)合PDE爆轟噪聲實驗測量結(jié)果與數(shù)值計算結(jié)果，劃定PDE爆轟噪聲由近場向遠場傳播的3 個區(qū)域界線，如圖1所示。超壓≥9 100 Pa，為強非線性區(qū)。1 000 Pa≤超壓＜9 100 Pa，為弱非線性區(qū)。超壓＜1 000 Pa 作為線性區(qū)。當PDE 噪聲傳播到線性區(qū)時，控制方程為線性聲波方程，采用拋物線方程(PE)法對其進行求解。由于線性方程的求解比較簡單，本文對線性區(qū)不再介紹，著重對前兩部分進行研究。

圖1 爆轟產(chǎn)生的沖擊波由近場向遠場傳播過程數(shù)值計算的區(qū)域劃分Fig.1 Schematic diagram of the propagation of a shock wave from a source of detonation to the far field

1.1 強非線性區(qū)控制方程及計算方法

不考慮流體微團間的熱傳導和熱輻射等耗散效應(yīng)，考慮粘性影響下的PDE 外流場軸對稱控制方程［6-7］為

強非線性區(qū)域內(nèi)PDE 噪聲具有極強的非線性，采用時空守恒元/求解元(CE/SE)方法求解軸對稱納維-斯托克斯方程。CE/SE 方法是近年來求解強間斷的一種新的數(shù)值方法。它最初由美國國家航空航天局科學家Chang［8］提出，該方法在構(gòu)造思想上與其它傳統(tǒng)的數(shù)值方法有所不同，它將空間與時間統(tǒng)一對待，從守恒律積分方程出發(fā)，設(shè)立守恒元和求解元，使其計算格式在局部和全部計算區(qū)域內(nèi)嚴格保證物理意義上的守恒。CE/SE 方法可以準確捕捉非線性區(qū)的重要特征:強間斷性。

PDE 外流場計算區(qū)域如圖2所示。軸向和徑向計算總網(wǎng)格數(shù)為1 200 ×1 200，再增加網(wǎng)格數(shù)量對計算結(jié)果影響極小。ABFG 所圍區(qū)域為PDE 管，BCDEFG 所圍成的區(qū)域為外流場計算區(qū)域。FG 為PDE 管管壁，BG 為PDE 管管口，AB、BC 段采用軸對稱邊界條件，CD、DE、EF 段采用非反射自由邊界，AF、FG 段采用固壁反射邊界條件，BG 采用入口邊界條件，將脈沖爆轟發(fā)動機管內(nèi)流場中爆轟管出口參數(shù)作為外流場計算的入口參數(shù)。

圖2 強非線性區(qū)計算區(qū)域示意圖Fig.2 Schematic diagram of computational domain

1.2 弱非線性區(qū)控制方程的推導及計算方法

第2 個區(qū)域為弱非線性區(qū)，是由強沖擊波向線性聲波的轉(zhuǎn)變過程，這個區(qū)域的聲壓變化采用非線性行波方程(NPE)進行描述，NPE 方程由McDonald和Kuperman 首先提出，該模型被最初廣泛應(yīng)用于預測高強度沖擊波在水中向遠場線性區(qū)的傳播過程，近幾年也被成功的應(yīng)用于預測弱非線性聲波在空氣中的傳播，Piacsek 等［9］將NPE 方程用于計算飛行器在高空突破音速時產(chǎn)生的音爆向地面?zhèn)鞑サ倪^程。NPE 方程是基于弱非線性聲波擾動的單一方程，同時考慮了波傳播過程中的折射和非線性等因素，且其計算區(qū)域是以當?shù)芈曀傧騲 正方向移動的，可以很好捕捉計算區(qū)域有限振幅的變化。

軸對稱柱坐標系下的NPE 方程［10］為

將p'關(guān)于ρ'做2 階展開可得

定義無量綱壓強的Q =p'/ρ0c20，由(3)式可知，Q=R+O(ρ'2)，若取1 階近似，則二者相等，因此，在理想流體中，可以將NPE 方程演變?yōu)殛P(guān)于Q 的方程。

本文采用時間算子分裂格式對NPE 方程進行求解，即對(4)式右邊兩項采取不同格式進行求解。其中第1 項為非線性項，采用2 階迎風通量校正格式(FCT);對于右邊第2 項，即衍射項，采用Crank-Nicholson 方法進行離散求解，對于其中的積分項，采用梯形積分法進行求解。

對于NPE 算法計算區(qū)域，同樣采用軸對稱方程建立模型及矩形網(wǎng)格劃分方法。由于NPE 方程建立時進行了小角度的近似，因此從聲源點開始，計算區(qū)域與對稱軸所成的夾角不能大于10°. 如圖1中標注區(qū)域所示，NPE 計算的區(qū)域采用動網(wǎng)格法，計算區(qū)域的移動速度和當?shù)芈曀僖恢?，隨時間的積累向前移動。NPE 算法程序的計算網(wǎng)格的長度的選取中，x 代表對稱軸方向，r 代表徑向方向，軸向上網(wǎng)格長度dx =0.25 cm，徑向上dr =0.125 cm. 選取340 m/s 為當?shù)匾羲?。對于計算區(qū)域的上和右邊界，均采用p =0，對稱軸上采用軸對稱邊界條件:

2 計算結(jié)果及分析

2.1 強非線性區(qū)計算結(jié)果及分析

圖3 PDE 管外的壓力云圖Fig.3 Pressure nephogram outside PDE

圖3為爆轟發(fā)動機出口處，不同時刻，CS/SE 方法計算得到的爆轟管外的壓力云圖。其中圖3(a)、圖3(b)和圖3(c)分別為0.138 ms、1.324 ms、4.014 ms時刻的爆轟管外的壓力云圖。從圖3(a)中可以看出，爆轟波離開出口以后，迅速退化為無化學反應(yīng)的“弓形”激波，在激波掠過區(qū)域，形成高壓區(qū)?！肮巍奔げň哂忻黠@的指向性，在軸向位置上壓力最高，隨著角度的增大，壓力下降。同時，管內(nèi)噴射的高壓爆轟產(chǎn)物在出口附近形成球狀膨脹波向四周膨脹，由于過渡膨脹，在膨脹波后產(chǎn)生負壓區(qū)。由圖3(b)可以看出，隨著時間的推移，激波沿軸向和徑向方向上繼續(xù)傳播，最先形成的激波逐漸衰減，激波后面的膨脹波繼續(xù)向外傳播，并由于燃燒產(chǎn)物是以噴流方式向外噴出，與管外靜止空氣邊界層形成的強烈的湍流脈動，導致渦系的產(chǎn)生，在x =0.35 m 上方，存在一明顯渦系。隨著時間的推移，渦系沿軸向傳播。當t=4.014 ms 時，退化的激波傳播到1.8 m，渦系傳播到0.8 m 處。

2.2 弱非線性區(qū)計算結(jié)果及分析

圖4為在4.014 ms 時，管外對稱軸上的聲壓曲線。此時，沖擊波的波陣面距離爆轟管口1.8 m，峰值為9 100 Pa，聲壓級達到173.16 dB，根據(jù)分區(qū)判斷標準，爆轟噪聲開始進入弱非線性區(qū)。如圖1中所示，強非線性區(qū)與弱非線性區(qū)的計算區(qū)域存在交匯區(qū)，利用交匯區(qū)內(nèi)第1 區(qū)的計算結(jié)果，作為第2 區(qū)計算的初始值。由于除了第1 個完整波形之外的曲線上升與下降，均為渦系所導致，而渦系會隨傳播距離的增加而迅速衰減，因此計算中只取第1 個完整波形。取強非線性區(qū)計算結(jié)果中，x 方向上1.2 ～2.0 m，r 方向上0 ～0.8 m 區(qū)域作為NPE 計算區(qū)域的初始值，如圖5所示。因為NPE 方程是關(guān)于聲壓的單一方程，因此，將強非線性區(qū)計算結(jié)果減去大氣壓的值作為下一步計算的初值。

圖4 軸線上的聲壓變化曲線Fig.4 Curve of on-axis sound pressure

圖6、圖7和圖8分別為擾動傳播到距離爆轟發(fā)動機管口2.3 m、5.8 m 和10.0 m 位置處，由NPE方程計算得到的聲壓云圖。由于聲波的衍射作用，聲波在軸向和徑向上傳播。隨著能量的擴散，壓力擾動的區(qū)域越來越大。徑向上的壓力擾動由最開始時的0.8 m 處傳播到3.4 m 處，軸向上的壓力擾動由1.8 m 處傳播到10.0 m 處，徑向上的壓力擾動傳播速度明顯較慢。這是因為PDE 爆轟噪聲具有明顯的指向性，能量主要集中在軸向上，因此在軸向上的傳播速度快于徑向上的傳播速度。由于傳播過程中，聲波的能量擴散衰減，聲壓峰值也從9 100 Pa(173.16 dB)，衰減為6950 Pa(170.82 dB)、2 645 Pa(162.43 dB)和1 515 Pa(157.58 dB)，且在能量較高區(qū)域，衰減速度較快。衰減速度最快時可達3 756 Pa/m(33 dB/m). 然后，衰減速度逐漸減慢，在接近線性區(qū)時，衰減速度最小為198.18 Pa/m(1.02 dB/m).

3 計算結(jié)果與實驗結(jié)果的對比及分析

圖5 弱非線性區(qū)的初值壓力云圖和軸線上的聲壓曲線圖Fig.5 Pressure nephogram and curve of initial sound pressure in weak nonlinear region

圖6 波振面到達2.3 m 時的壓力云圖Fig.6 The pressure nephogram at 2.3 m

為了驗證PDE 爆轟噪聲場計算結(jié)果的準確性，對相應(yīng)距離處PDE 管外的爆轟噪聲聲壓進行了測量。所用PDE 試驗系統(tǒng)主要包括PDE、供氣系統(tǒng)、供油系統(tǒng)、點火控制系統(tǒng)和測試系統(tǒng)等。以93 號汽油為燃料，空氣為氧化劑，脈沖爆轟發(fā)動機爆轟管直徑為80 mm，工作頻率為10 Hz，傳感器布置于距離管口0.4 m、0.8 m、1.2 m 、1.6 m、2.4 m、3.0 m、5.0 m 和10.0 m 處，測試得到的噪聲電信號通過信號調(diào)理器處理后由同步數(shù)據(jù)采集系統(tǒng)采集，采樣率為500 kS/s.

聲壓級的表達式為

圖7 波振面到達5.8 m 時的壓力云圖Fig.7 The pressure nephogram at 5.8 m

圖8 波振面到達10.0 m 時的壓力云圖Fig.8 The pressure nephogram at 10.0 m

式中:SPL 為聲壓級(dB);p 為聲壓(Pa);pref=2 ×10-5Pa，為基準聲壓。

分別取強非線性區(qū)和弱非線性區(qū)內(nèi)對稱軸上不同位置處的觀測點計算得到的壓力峰值，連同實驗測量所得數(shù)據(jù)，繪制峰值壓力隨軸向距離x 變化的曲線，如圖9所示。在強非線性區(qū)壓力峰值迅速下降，衰減速度較快。這是由于在此區(qū)域，沖擊噪聲相比射流噪聲要大很多，但由于沖擊噪聲在大幅值時衰減較快，導致該區(qū)域的壓力峰值迅速衰減。到1.6 m 處，隨著沖擊噪聲的迅速衰減，噴流噪聲占據(jù)主導地位，其衰減速度相對較慢。但是，當脫離燃氣噴流核心區(qū)時，由于湍流脈動效應(yīng)的減小，噴流噪聲較快衰減。此時小振幅沖擊噪聲的峰值衰減速度相對較慢。在2.4 m 處，沖擊噪聲占據(jù)主要地位，進入弱非線性區(qū)，由于小振幅沖擊噪聲的衰減速度減慢，因此之后的噪聲峰值的衰減速度也相對較慢。

PDE 近場爆轟噪聲峰值ppeak與軸向距離x 的倒數(shù)滿足一定的變化規(guī)律，計算公式［2］為

圖9 PDE 爆轟噪聲聲壓衰減曲線Fig.9 Curves of PDE detonation noise amplitudes on axis

式中:k 為衰減系數(shù)，A 為非線性系數(shù)。當在非線性區(qū)時，1 ＜A≤3. 非線性程度越強，A 的取值越接近3;非線性程度弱，A 的值接近1. 利用數(shù)值計算和實驗測量所得聲壓級繪制隨軸向距離x 的變化規(guī)律，如圖10 所示。其中x 采用對數(shù)坐標，A =1.0 直線為以第2 區(qū)最后1 個點為參考點、A 的值取1.0 繪制的曲線?？梢钥闯觯谧铋_始的時刻，衰減速度很快，此時A=2.1. 隨著軸向距離的增加，速度減慢。但在第2 區(qū)中，衰減速度逐漸變緩，其中到x =10 m處最后一個點時，A =1.02，此時聲壓級的衰減規(guī)律已經(jīng)很接近線性區(qū)的衰減規(guī)律。

圖10 衰減規(guī)律擬合曲線Fig.10 Fitting curves of attenuation law

4 結(jié)論

將PDE 爆轟噪聲場進行分區(qū)求解，可以有效的模擬出脈沖爆轟發(fā)動機爆轟噪聲的傳播特性，避免了單一方法計算的不準確性，對PDE 爆轟噪聲場的傳播特性的研究具有重要的參考價值。并得出以下結(jié)論:

1)以超壓9 100 Pa 作為強非線性區(qū)和弱非線性區(qū)的判斷標準是合理的，可以有效避免CE/SE 方法計算誤差所產(chǎn)生的衰減速度較快現(xiàn)象。

2)爆轟噪聲峰值聲壓級，在出口處迅速衰減。隨著軸向距離的增加，衰減速度逐漸放緩。在強非線性區(qū)，衰減速度最大值可達33 dB/m. 在弱非線性區(qū)，最小衰減速度為1.02 dB/m. 非常接近線性區(qū)的衰減速度。

References)

［1］Boesch H E，Reiff C G，Benwell B T.Modification of the acoustic spectrum of detonation tube shock waves by timed multiple-pulse addition，ARL-TR-2203［R］. America:Army Research Laboratory Report，2000.

［2］Glaser A J，Caldwell N，Gutmark E.A fundamental study on the acoustic behavior of pulse detonation engines［C］∥45th AIAA Aerospace Sciences Meeting and Exhibit. Reno，Nevada，US:AIAA，2007.

［3］He X，Karagozian A R . Performance and noise characteristics of pulse detonation engines［C］∥42nd AIAA Aerospace Sciences Meeting and Exhibit. Reno，Nevada，US:AIAA，2004.

［4］Eerden F V，Védy E. Propagation of shock waves from source to receiver［J］.Noise Control Engineering，2005，53(3):87 -93.

［5］王秉義.槍炮噪聲與爆炸聲的特性和防治［M］. 北京:國防工業(yè)出版社，2001:31 -32.WANG Bing-yi. Gun muzzle noise ＆ explosive sound-Character，Protection and control［M］. Beijing:National Defence Industry Press，2001:31 -32.(in Chinese)

［6］王杰，翁春生.脈沖爆轟發(fā)動機沖擊波外流場數(shù)值模擬［J］.火炮發(fā)射與控制學報，2009，30(3):53 -56.WANG Jie，WENG Chun-sheng. Numerical simulation shock wave exteral flow of pulse detonation engine［J］. Journal of Gun Launch＆ Control，2009，30(3):53 -56.(in Chinese)

［7］王研艷，翁春生. 尾噴管構(gòu)型對多循環(huán)兩相脈沖爆轟發(fā)動機流場及性能影響［J］.航空動力學報，2013，28(10):66 -69.WANG Yan-yan，WENG Chun-sheng. Effects of nozzle on flow field and performance of multi-cycle two-phase pulse detonation engines［J］. Journal of Aerospace Power，2013，28(10):66 -69. (in Chinese)

［8］Chang S C. The method of space-time conser vation element and solution element-a new approach for solving the navier-stokes and euler equations［J］. Journal of Computational Physics，1995，119(2):295 -324.

［9］Piacsek A A，Locey L L，Sparrow V W. Time-domain modeling of atmospheric turbulence effects on sonic boom propagation［C］∥14th AIAA/CEAS Aeroacoustics Conference. Vancouver，British Columbia，Canada:AIAA，2008.

［10］Piacsek A A，Plotkin K J. Application of nonlinear progressivewave equation to sonic boom transition focus［C］∥51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition. Grapevine，Texas，US:AIAA，2013.