冷 巖，張勁柏，錢戰(zhàn)森，*

(1.中國航空工業(yè)空氣動力研究院，沈陽 110034；2.高超聲速氣動力/熱技術(shù)重點實驗室，沈陽 110034；3.北京航空航天大學 航空科學與工程學院，北京 100191)

0 引言

自“協(xié)和號”首飛以來，聲爆問題很大程度上決定了超聲速客機能否進入商業(yè)運營及取得商業(yè)成功，因而成為了研制超聲速民用飛機必須解決的關(guān)鍵難題[1-3]。特別地，針對起飛階段不可避免的超聲速加速過程和快速機動飛行時產(chǎn)生的聚焦聲爆現(xiàn)象（如圖1 所示[4]）的研究成為了研制新一代環(huán)保型超聲速飛行器的重要課題。

圖1 不同飛行狀態(tài)下四種聚焦聲爆示意圖[4]Fig.1 Illustration of four types of caustics under supersonic flight operations[4]

只有在超聲速巡航狀態(tài)下地面聲爆特征才是標準N 波。超聲爆現(xiàn)象一般產(chǎn)生于飛機起飛后的超聲速加速階段或變馬赫數(shù)狀態(tài)下發(fā)生轉(zhuǎn)彎、拉升等機動動作時，一般呈非對稱U 型，其過壓峰值通常是超聲速巡航狀態(tài)下產(chǎn)生的N 波的2～5 倍[5-6]。因此，量化焦散面上的聚焦聲爆對超聲速飛行器研制具有重要意義。經(jīng)典理論在幾何聲學框架下開展聲爆數(shù)值預測研究[7-8]，其中相位函數(shù)由射線路徑?jīng)Q定，信號幅值由射線管面積決定，沿每條射線的非線性效應解釋了壓力信號從近場復雜激波/膨脹波系到地面簡單N 波的變化。以焦散面為中心的超聲爆區(qū)域射線管面積為零，導致幾何聲學方法在焦散面附近因奇異性而失效。

國外超聲爆現(xiàn)象的研究歷史悠久，試驗測量主要是在20 世紀六七十年代。Haglund 和Kane 針對加速引起的聚焦聲爆開展了詳細研究[9-10]。其他研究團隊針對超聲速飛行器加速和轉(zhuǎn)彎引起的聚焦聲爆開展飛行測量研究，并成功捕捉到超聲爆現(xiàn)象[11]。Sanai[12]等證明利用彈道靶實現(xiàn)超聲爆測量的可能性。在開展試驗研究的同時，也從數(shù)值建模和分析的角度對超聲爆進行了研究。數(shù)值建模研究始于對具有解析解[12]的線性Tricomi 方程和非線性Tricomi方程[13-14]（nonlinear Tricomi equation,NTE）的推導，之后眾多學者[15-19]用NTE 方程對超聲爆進行建模，進而預測得到焦散線附近的聲爆信號。在聲爆預測及遠場傳播過程中，必須考慮包含損耗機制在內(nèi)的真實大氣效應影響，熱黏吸收效應和分子弛豫效應是兩個主要耗散機制。通過在NTE 方程中添加吸收/色散效應，Salamone[20]等提出廣義非線性Tricomi 方程（lossy nonlinear Tricomi equation,LNTE）并數(shù)值求解，與飛行試驗數(shù)據(jù)對比表明，在聚焦聲爆的研究中添加耗散機制是非常必要的。非線性Tricomi 方程為混合型雙曲/橢圓方程，Sescu[21]等通過添加虛擬時間變量以保證守恒形式的方程為雙曲型，并采用非線性限制器來控制求解域的振蕩。在此基礎(chǔ)上，Kanamori[22]等開發(fā)了超聲爆評估工具FFnoise。

國內(nèi)在近場聲爆CFD 計算與基于非線性Burgers方程的遠場聲爆預測方面取得了較大進展。西北工業(yè)大學[23-26]、中國空氣動力研究與發(fā)展中心[27]、中國航空工業(yè)空氣動力研究院[28-30]、中國航空研究院[31]、北京航空航天大學[32]等多家單位都開展了相關(guān)研究，發(fā)展了一系列聲爆的近、遠場預測方法。但是，針對超聲速飛行器機動狀態(tài)下的超聲爆現(xiàn)象的研究，尚未見公開發(fā)表的文章。

鑒于此，受Salamone 等[20,33]工作的啟發(fā)，本文基于廣義Tricomi 方程建立表征焦散線附近壓力波演化過程的空間模型，采用偽時間推進和分裂算法，在時域和頻域分開求解以達到數(shù)值模擬超聲爆的目的，并開發(fā)相應計算機代碼ARI_Superboom，線性解析解分析和飛行試驗數(shù)據(jù)驗證表明了模型及算法的準確性。

1 超聲爆及其模型方程

1.1 超聲爆

當飛行器超聲速飛行時，由近場復雜激波/膨脹波系引起的壓力脈動沿聲射線傳播。聲射線垂直于馬赫錐，因此在固定參考系下，巡航狀態(tài)生成的聲射線間彼此平行，互無交叉。當飛行器加速飛行時，隨著馬赫數(shù)的增加，馬赫角逐漸減小，聲射線間的交叉導致焦散的形成，如圖2 所示。進一步，在衍射和非線性效應的共同作用下會在焦散附近產(chǎn)生最大聚焦區(qū)域。當聲射線進入聚焦區(qū)域，沿射線傳播的波形就會發(fā)生畸變，進而產(chǎn)生聚焦聲爆或非對稱U 波。這種變形發(fā)生在聚焦區(qū)域的任何位置，包括焦散面和地面之間的交叉點周圍，因此，如圖2 所示，可在地面觀測到由加速度引起的超聲爆現(xiàn)象。

圖2 超聲速飛行器加速狀態(tài)聚焦聲爆形成示意圖[20]Fig.2 Illustration of the focus boom formation during the acceleration phase of a supersonic flight vehicle[20]

以勻加速為例（圖2），Ma從0.96 加速到1.35，馬赫數(shù)隨時間變化率為0.004/s，飛行高度為13.716 km。加速過程中，當馬赫數(shù)大于截止馬赫數(shù)后出現(xiàn)聚焦區(qū)域。標準大氣無風狀態(tài)下截止馬赫數(shù)為1.153 2，根據(jù)真實飛行狀態(tài)，初始計算狀態(tài)設(shè)為Ma=1.156，每隔1 s 計算一次聲射線，得到如圖3 所示聲射線射線和焦散線。

圖3 ARI_Boom 預測所得聲射線和焦散線Fig.3 Rays and caustic predicted by ARI_Boom

1.2 模型方程

在三維空間中，焦散是一個面；在二維空間中，焦散是一條線，如圖3 和圖4 中紅色線所示，其中x軸與焦散線相切于O點，z軸與焦散線垂直于O點。焦散線上方為強干擾區(qū)（明亮區(qū)），反映了入射壓力波在焦散線附近的強非線性相互作用，在此區(qū)域內(nèi)的每個位置都有入射和出射兩束射線通過。焦散線下方為消散波區(qū)（陰影區(qū)），反映由于衍射效應投射的波，一般強度較弱。

圖4 射線在焦散附近變化示意圖[33]Fig.4 Illustration of rays in the vicinity of the caustic[33]

廣義Tricomi 方程表征的是焦散附近壓力波的演化圖像，控制焦散線附近聲爆信號變化的無量綱形式的廣義Tricomi 方程如下[20]：

方程（1）左端前兩項表示衍射效應，屬于線性項；第三項和第四項表示自然來流影響，也屬于線性項，Max=u0x/c0、Maz=u0z/c0分別表示自然來流馬赫數(shù)在x向和z向的分量，u0為飛行器所在高度的來流速度；第五項為非線性項；第六項為耗散項，包含了熱黏性和分子弛豫效應。無量綱聲學壓力=(p?p0)/pac，聲學馬赫數(shù)pac為輸入波形特征參考壓力（最大過壓峰值）；ρ0、p0、c0分別為大氣環(huán)境密度、壓力和聲速，ε為與焦散線相切于O點的聲射線的曲率半徑。

為了避免求解多維空間問題，Tricomi 方程對焦散線附近的現(xiàn)象進行了簡化，建立了一維空間模型。通過合理地選擇計算區(qū)域，可以借助定義與2 個維度空間位置相關(guān)的無量綱時間參數(shù)來表征壓力波從入射到出射的演化過程；這一無量綱參數(shù)既可以看成是沿著焦散線的空間發(fā)展，也可以看作是壓力自身的時間發(fā)展。式（1）中無量綱時間變量的定義為：

其中Rcau為焦散線曲率半徑。

方程（1）中 ε為小參數(shù)，反映特征波長λac=c0/fac與衍射層厚度 δ的比率：

定義 μ為與衍射相關(guān)的非線性效應的衡量值，對于超聲爆情形其量值在0.1 左右。根據(jù)式（5）非線性項系數(shù)可改寫為：

其中，下標 υ可取為1 和2，分別代表氧氣和氮氣分子弛豫效應；τυ為分子弛豫時間；c∞,υ為凍結(jié)聲速，c∞,υ=c0+?c∞,υ，?c∞,υ為聲速增量，其近似表達式及參數(shù)可見參考文獻[23，36-37]。

2 廣義Tricomi 方程的數(shù)值解法

2.1 數(shù)值求解

廣義Tricomi 方程是關(guān)于z和的雙曲/橢圓混合型方程，通常引入偽時間變量后采用時間推進求解，當偽時間收斂時，就可以得到壓力波信號在焦散線附近的分布。數(shù)學上的求解區(qū)域及圖像如圖5 所示，其中，即計算區(qū)域在方向為2 倍衍射邊界層厚度。

圖5 Tricomi 方程求解區(qū)域示意圖[22]Fig.5 Illustration of the computational domain of the Tricomi equation[22]

方程（1）的數(shù)值求解過程中，引入虛擬時間變量σ，則方程變?yōu)椋?/p>

基于傅里葉變換，數(shù)值求解通過兩步分裂算法分別在時域和頻域中實現(xiàn)。針對每一個無量綱角頻率ωn，在頻域中求解衍射和吸收/色散效應項，非線性項在時域內(nèi)求解。在分裂算法的每一步中，設(shè)置偽時間增量的大小一致并限制其量值以避免非線性求解中出現(xiàn)激波過沖現(xiàn)象[29]。

第一步：在頻域內(nèi)求解衍射項、耗散項、x向大氣風項和z向大氣風項；

在該步的離散求解時，對于偽時間項采用一階前差，z方向的導數(shù)項取中心差分，具體計算格式為：

第二步：在時域內(nèi)采用激波捕捉法求解非線性項，如下式：

該步的時間離散承接上步，右端項采用經(jīng)典的二階Roe 格式。簡略形式如下：

2.2 邊界條件

上述數(shù)值求解過程中每一步都需要設(shè)置邊界條件。時域內(nèi)，假設(shè)輸入波形到達前和輸出波形離開后聲學壓力為0，因此，

式（14）中，F(xiàn)表示經(jīng)過pac無量綱化的輸入波形。值得注意的是，上述輻射邊界條件可允許在不知道輸出波形的情況下使用。頻域內(nèi)，對應的輻射邊界條件經(jīng)傅里葉變換后形式為[34]：

2.3 焦散線附近壓力波演化過程模擬方法

圖6 給出超聲速飛行器勻加速狀態(tài)聚焦聲爆預測過程及Tricomi 方程求解區(qū)域。整個過程包括以下三步：第一步，基于航空工業(yè)氣動院自主研發(fā)的ARI_Overset 高精度數(shù)值模擬軟件，采用CFD 手段數(shù)值求解三維Navier-Stokes（N-S）方程得到近場流場結(jié)構(gòu)并提取模型下方設(shè)定位置處空間壓力分布，壓力分布提取位置一般為1～3 倍飛行器特征長度處；第二步，基于自主研發(fā)的ARI_Boom 聲爆預測程序開展射線追蹤和聲爆遠場傳播計算，獲得焦散線及所需聲射線信息和距離地面Zδ處聲爆信號，此時所得聲爆信號類似于標準N 波；第三步，將第二步所得聲爆信號作為聚焦計算初始輸入波形，通過自主研發(fā)的基于廣義Tricomi 方程的ARI_Superboom 超聲爆預測程序獲得地面或指定位置聚焦聲爆信息，由于聲衍射等效應的聯(lián)合作用，聲爆信號從幅值、形狀等方面發(fā)生了明顯的變化。第一、第二步涉及的近場CFD 預測和遠場聲爆在文獻[28-29]中均有驗證，本文的重點主要體現(xiàn)在第三步焦散線附近壓力波演化過程模擬方法中，為了驗證其準確性，聚焦區(qū)域輸入波形采用文獻數(shù)據(jù)。

圖6 聚焦聲爆預測過程示意圖[33]Fig.6 Illustration of the focus boom prediction[33]

3 方法驗證

3.1 線性Tricomi 方程分析

一種快速的機動飛行聲爆遠場傳播預測方法是采用線性Tricomi 方程的解析解計算聚焦區(qū)域的聲爆。線性Tricomi 方程相對于廣義Tricomi 方程僅保留了聲衍射項，方程變?yōu)椋?/p>

其中，IFFT 表示傅里葉逆變換，sgn為符號函數(shù)，Ai為艾里函數(shù)。

圖7 初始輸入N 波波形Fig.7 N-shape wave for the initial input

圖8 計算所得空間壓力云圖ⅠFig.8 Pressure contour Ⅰobtained by the numerical simulation

圖9 預測結(jié)果與解析解對比Fig.9 Comparison between the numerical and analytical solutions

3.2 SCAMP 項目飛行試驗

NASA Dryden 中心基于F-18 飛機開展的SCAMP（superboom caustic analysis and measurement project）飛行試驗提供大量聚焦聲爆數(shù)據(jù)以驗證數(shù)值預測方法和代碼。試驗過程中由81 個麥克風組成的地面線性陣列測量超聲速飛機產(chǎn)生的聚焦聲爆現(xiàn)象。通過飛行軌跡記錄、高空/地面大氣條件監(jiān)測、聲學數(shù)據(jù)推算等手段獲得射線追蹤和CFD 模擬所需的飛行條件。

本文選取“狀態(tài)C”[33]作為驗證算例進行討論。F-18 飛機的飛行馬赫數(shù)為1.23，馬赫數(shù)變化率為0.003 5/s，俯仰角速度為?0.25°/s，飛行高度為12.8 km，迎角為2.3°，飛機特征長度為46.3 m，近場特征提取位置為3 倍特征長度處，近場特征如圖10 所示[28]。

圖10 聚焦聲爆預測初始輸入波形Fig.10 Initial input wave for the focus boom prediction

圖11 計算所得空間壓力云圖ⅡFig.11 Pressure contour Ⅱ obtained by the numerical simulation

圖12 預測結(jié)果與飛行試驗數(shù)據(jù)對比Fig.12 Comparison between the numerical solution and the flight test data

圖13 非線性項對預測結(jié)果影響對比Fig.13 Comparison of the numerical results calculated with and without the nonlinear term

圖14 耗散項對預測結(jié)果影響對比Fig.14 Comparison of the numerical results calculated with and without the viscous term

4 結(jié)論

本文基于廣義Tricomi 方程建立真實大氣效應的聚焦區(qū)域壓力波演化過程模型開發(fā)自研程序ARI_Superboom，采用頻域和時域相結(jié)合的方式離散模型并數(shù)值求解，提高了計算精度和效率；鑒于飛機超聲速機動飛行條件聲爆預測方法的復雜性，基于航空工業(yè)氣動院自研軟件建立超聲速飛行器勻加速狀態(tài)下聚焦聲爆預測分析手段，進而獲得聚焦區(qū)域超聲爆信息，并通過線性Tricomi 方程分析和SCAMP 項目飛行試驗兩個算例驗證了ARI_Superboom 程序的準確性和過程分析的可行性；基于SCAMP 項目飛行試驗實測數(shù)據(jù)開展效應分析研究，對比數(shù)據(jù)表明，在超聲爆數(shù)值預測中必須考慮非線性效應和耗散機制的影響，但是影響區(qū)域集中在明亮區(qū)，對于衍射效應主導的陰影區(qū)影響較小。