張繹典，黃江濤，高正紅,*

1. 西北工業(yè)大學(xué) 航空學(xué)院，西安 710072 2. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

近年來，隨著全球化經(jīng)濟(jì)的不斷進(jìn)步，人們對乘坐超聲速客機進(jìn)行長途旅行的渴望愈發(fā)強烈，各國也相繼開展了不少有關(guān)超聲速客機的研究項目。音爆噪聲是制約早期超聲速客機(如“協(xié)和”和圖-144)持續(xù)運營的重要因素之一[1]，各國研究者對其關(guān)注的程度也日益增加。

音爆是超聲速飛行器在近場產(chǎn)生的一系列以激波為主的擾動波，經(jīng)過在大氣層內(nèi)的傳播，最終在遠(yuǎn)場演化成的“N”形聲波[2]。遠(yuǎn)場波形類似英文字母“N”的原因是聲波在高溫高壓處的傳播速度高于低溫低壓處，這種非線性傳播效應(yīng)造成了弱激波的堆積，最終在信號的兩端形成兩道較強的激波。準(zhǔn)確的音爆預(yù)測是超聲速飛機特別是超聲速客機設(shè)計面臨的重要技術(shù)問題。音爆模擬的關(guān)鍵是需要精確模擬音爆信號的過壓值和上升時間[3]，它們和音爆強度直接相關(guān)。通過減小過壓值，增加上升時間可以減弱音爆的強度。

目前的音爆預(yù)測方法主要分4類[4]：①基于Calson法的半經(jīng)驗預(yù)測方法[5];②基于超聲速線化理論和聲學(xué)的預(yù)測方法[6];③將計算流體力學(xué)(CFD)和聲學(xué)結(jié)合的預(yù)測方法[7];④全場直接數(shù)值求解[8]。

方法①是半經(jīng)驗法，適用于工程估測。方法②和方法③都是先計算飛行器的近場聲壓，然后通過聲學(xué)方法將預(yù)測聲波傳遞至遠(yuǎn)場。二者的不同在于前者采用超聲速線化理論計算近場而后者使用CFD的方法。由于線化小擾動理論僅適用于細(xì)長體，故方法③的適用面更廣。方法④看似最簡單直接，但因其巨大的計算量加上需要配套發(fā)展低耗散的高精度數(shù)值方法，國際上鮮有科學(xué)家進(jìn)行嘗試。2016年，東京大學(xué)的Yamashita和Suzuki[8]通過使用超算集群首次進(jìn)行了音爆全場直接求解。但其模擬的海拔僅為500 m，遠(yuǎn)低于超聲速客機的巡航高度(17 km左右)，同時一個模型的計算總時間達(dá)到了31.5 h。綜上所述，統(tǒng)籌考慮計算量和計算精度，基于CFD和聲學(xué)的預(yù)測方法是比較適合的方法，也是當(dāng)代研究者主要關(guān)注的方法。

在得到近場過壓分布后，通過聲學(xué)的方法模擬遠(yuǎn)場聲波信號也分以下3種方式來實現(xiàn)：①基于弱激波理論的波形參數(shù)法[9]；②基于時域、頻域分析的Pestorius-Anderson法[10]；③基于非線性聲學(xué)的增廣Burgers方程法[7]。其中波形參數(shù)法只能模擬音爆的超壓值，無法直接預(yù)測上升時間，也無法對結(jié)果進(jìn)行快速傅里葉變換(FFT)得到感覺噪聲級(Perceived Noise Level,PLdB)。過去，人們只能通過“3/p”定理[6]或穩(wěn)態(tài)解理論(tanh方法)[11]等經(jīng)驗公式來估計上升時間的大小，這會影響音爆強度的預(yù)測精度，可能會導(dǎo)致飛機設(shè)計師做出錯誤的判斷[7]。Pestorius-Anderson法需要在計算的過程中不斷進(jìn)行時頻域變換，這樣會造成誤差的堆積，同時計算開銷也隨之增加?；诜蔷€性聲學(xué)的增廣Burgers方程法直接在時域內(nèi)對音爆傳播過程進(jìn)行模擬求解，排除了時、頻域轉(zhuǎn)換產(chǎn)生的誤差，也減小了計算代價。同時，該方法可以考慮溫度、濕度等環(huán)境條件對音爆傳播的影響。因此增廣Burgers方程類方法成為了學(xué)者較為關(guān)注的一類音爆遠(yuǎn)場預(yù)測方法。

近年來，音爆的相關(guān)研究也逐漸引起了國內(nèi)學(xué)者的重視。朱自強和蘭世隆[1]對音爆預(yù)測、優(yōu)化方法做了詳細(xì)梳理，指出了高精度流場模擬和增廣Burgers方程等方法的優(yōu)勢。王剛等[12]對典型標(biāo)模的音爆預(yù)測方法作了分析，指出采用適當(dāng)幾何修形和采用熵相容格式可以保證近場音爆預(yù)測精度。馮曉強等[13-14]發(fā)展了一整套基于CFD和波形參數(shù)法的音爆優(yōu)化軟件，并對低音爆機理進(jìn)行了探究。但是，在音爆傳播的模擬方法的選擇上，國內(nèi)學(xué)者幾乎都采用了波形參數(shù)法等基于弱激波理論的方法，尚未看到通過增廣Burgers方程來模擬音爆傳播過程的研究。

本文主要就音爆遠(yuǎn)場計算展開研究。通過推導(dǎo)增廣Burgers方程的無量綱形式，明確了其中各項的物理意義。采用算子分裂法，建立了增廣Burgers方程的數(shù)值求解方法。接著，提出了一些方法來提高過壓值、上升時間等音爆關(guān)鍵參數(shù)的計算精度。在此基礎(chǔ)上，用自研的基于增廣Burgers方程的遠(yuǎn)場傳播程序計算了第二屆國際音爆預(yù)測研討會[15](SBPW-2)的兩個標(biāo)模(LM1021、Axibody)以驗證程序的準(zhǔn)確性。最后基于以上方法，研究了不同大氣環(huán)境對地面音爆波形的影響，可為飛越不同氣候環(huán)境的遠(yuǎn)程超聲速客機的音爆預(yù)測提供參考。

1 Burgers方程理論

20世紀(jì)60年代，Blackstock[16]最早使用Burgers方程來模擬波在有損耗介質(zhì)(Lossy Media)中的傳播，提出了經(jīng)典Burgers方程。在此基礎(chǔ)上，他又和合作者Carlton[17]在70年代將經(jīng)典Burgers方程擴展為廣義Burgers方程(Generalized Burgers Equation)，使其能夠模擬幾何擴散(Geometrical Spreading)和非均勻介質(zhì)對波傳播的影響。1981年，Pierce[18]又將分子馳豫效應(yīng)引入廣義Burgers方程，并將其稱作增廣Burgers方程(Augmented Burgers Equation)。至此，該方程用以模擬音爆信號傳播的形式基本固定。

1.1 經(jīng)典Burgers方程

從Navier-Stokes方程可以得到忽略二階以上小量的非線性波動方程——Westervelt方程[19]，其一維形式為

(1)

式中：p′為擾動聲壓；z為傳播距離;t為傳播時間;b和β分別為經(jīng)典吸收系數(shù)和非線性系數(shù);ρ0和c0分別為環(huán)境密度和環(huán)境波速，隨海拔高度而變化。

為簡化式(1)，引入坐標(biāo)變換：

(2)

式中:t′為延遲時間，可以看做是波傳播過程中某點的當(dāng)?shù)貢r間(即以波傳播到當(dāng)?shù)氐臅r間為時間原點)。

由鏈?zhǔn)椒▌t可以得到式(1)中的各階導(dǎo)數(shù)為

(3)

忽略二階以上小量，經(jīng)過化簡后得到經(jīng)典Burgers方程[18]為

(4)

這里值得注意的是，經(jīng)過坐標(biāo)變換后，耗散項為聲壓對延遲時間的二階導(dǎo)，這和一般形式的Burgers方程中耗散項為聲壓對空間的二階導(dǎo)有所區(qū)別。

1.2 增廣Burgers方程的無量綱形式

在經(jīng)典Burgers方程中加入非均勻介質(zhì)、幾何擴散以及分子馳豫效應(yīng)的影響，即可得到增廣Burgers方程[11]為

(5)

式中:S為聲管面積[20];(Δc)ν為分子馳豫效應(yīng)造成的聲速變化量；τν為馳豫時間，下標(biāo)ν表示不同的大氣組成成分(如氧氣、氮氣)的馳豫過程。式(5)等號右邊5項分別對應(yīng)：非線性效應(yīng)、經(jīng)典吸收、大氣分層、幾何擴散以及分子馳豫效應(yīng)對音爆的影響。

從增廣Burgers方程的組成可以看出，該方程可以模擬聲波在分層、有損耗的真實大氣中的傳播。

為了便于求解，首先對式(5)進(jìn)行無量綱化處理：

(6)

增廣Burgers方程的數(shù)值解法在附錄A中給出。

2 算例驗證

基于附錄A中的的數(shù)值方法，開發(fā)了自研軟件來進(jìn)行遠(yuǎn)場預(yù)測，軟件的框架如圖 1所示。針對研發(fā)的過程中遇到的一些問題，提出了一些方法來提高過壓值和上升時間的計算精度。使用SBPW-2的兩個標(biāo)準(zhǔn)算例LM1021和Axibody來對軟件進(jìn)行校核。

2.1 驗證算例簡介

本文采用了音爆預(yù)測研討會會議提供的兩個標(biāo)準(zhǔn)算例來校核軟件。該會議是美國AIAA學(xué)會組織的專門研討音爆預(yù)測方法以及可信度校核的會議，至今召開兩屆(2014[21]年、2017[15,22]年)。該算例提供了完整的飛行數(shù)據(jù)、近場數(shù)據(jù)以及大氣數(shù)據(jù)，可以準(zhǔn)確設(shè)置增廣Burgers方程的各項輸入?yún)?shù)。由于該會議的影響力和權(quán)威性，國際上的學(xué)者廣泛采用會議提供的標(biāo)準(zhǔn)算例來檢驗自己的音爆預(yù)測工具。

2.2 “緩沖信號”對上升時間的影響

以SBPW-2中的LM1021算例為例，選取周向角φ為0°的音爆信號為近場輸入，周向角的定義如圖 2所示。官方給出的近場輸入信號如圖 3中的黑色點劃線所示。首先，直接使用官方給出的風(fēng)洞試驗數(shù)據(jù)作為近場輸入進(jìn)行了計算，并和NASA的sBOOM[23]軟件的計算結(jié)果作了比較。如圖 4所示，發(fā)現(xiàn)計算得到的遠(yuǎn)場信號的首段上升時間會偏小。導(dǎo)致這一現(xiàn)象的原因主要如下：增廣Burgers方程中存在強非線性，初始波形中幅值不同的位置傳播速度不同，幅值越大的部分傳播速度越快。這一速度差就會導(dǎo)致時域中幅值為正的部分波形向左偏移，幅值為負(fù)的部分波形向右偏移。如果僅僅將計算域限制在有擾動的區(qū)域，最前端的正值信號就沒有向左偏移的空間，這就導(dǎo)致了信號的失真。極端情況下會導(dǎo)致遠(yuǎn)場信號的首個上升段和縱軸重合從而無法正確預(yù)測上升時間。

為了驗證上述原因，在官方給出的近場輸入前加入一段無幅值的“緩沖信號”(如圖 3所示)，就可以有效提升上升時間的模擬精度。圖 4同時也給出了加入了“緩沖信號”后的遠(yuǎn)場輸出。可以發(fā)現(xiàn)在加入“緩沖信號”后，遠(yuǎn)場輸出的首個上升段基本和sBOOM的結(jié)果重合。由于非線性扭轉(zhuǎn)的作用，前端信號又向前位移了一小段距離，約為0.005 s。換個角度來說，如果不加入“緩沖信號”，就相當(dāng)于限制了音爆信號的總長度，在傳播的過程中信號實際長度始終不變，從物理上來說也是不合理的。這種不合理的現(xiàn)象在近幾年的一些期刊論文[7]中較為常見。

2.3 網(wǎng)格收斂性研究

同樣以LM1021算例為例，選取收斂歷程較有代表性的周向角為30°的音爆信號作為近場輸入進(jìn)行網(wǎng)格收斂性研究。

增廣Burgers方程在時間和空間上各有一個維度，在這兩個維度上的網(wǎng)格都是均勻的。下面依次對這兩個維度進(jìn)行網(wǎng)格收斂性分析：首先固定空間網(wǎng)格密度，對時間網(wǎng)格進(jìn)行收斂性分析，然后固定時間網(wǎng)格密度，對空間網(wǎng)格進(jìn)行收斂性分析。表 1和表 2分別給出了時間網(wǎng)格收斂性以及空間網(wǎng)格收斂性研究中的網(wǎng)格量和計算時間。其中：numz為空間網(wǎng)格的網(wǎng)格點數(shù)目：numt為時間網(wǎng)格的網(wǎng)格點數(shù)目。通常，人們更喜歡使用時間采樣率(kHz)來表示時間方向的網(wǎng)格密度。

圖5和圖6分別給出了遠(yuǎn)場信號在時間和空間上的收斂歷程。從中可以發(fā)現(xiàn)，遠(yuǎn)場波形對時間網(wǎng)格的密度更加敏感。加密時間網(wǎng)格后，遠(yuǎn)場波形的形態(tài)發(fā)生了較大的變化。需要將網(wǎng)格密度加密到一定程度后，得到的遠(yuǎn)場聲壓波形才逐漸收斂。相反，在空間方向加密網(wǎng)格后，地面波形并沒有出現(xiàn)很大的變化。在SBPW-2[15]會議中，P4和P10這兩位參會者可能因為時間方向的網(wǎng)格加密不充分，提交的結(jié)果出現(xiàn)了和低采樣率的波形形態(tài)相類似的失真現(xiàn)象，如圖 5所示。

表1 時間網(wǎng)格收斂性分析算例Table 1 Cases for temporal mesh convergence study

表2 空間網(wǎng)格收斂性分析算例Table 2 Cases for spatial mesh convergence study

綜上所述，在運用上述的數(shù)值方法求解增廣Burgers方程時，建議適當(dāng)加密計算網(wǎng)格，尤其是時間方向的網(wǎng)格，建議的信號采樣率為85 kHz以上。

2.4 LM1021算例驗證

LM1021為SBPW-2的第一個標(biāo)準(zhǔn)算例，這是洛克西德·馬丁公司設(shè)計的一款未來超聲速客機，其幾何外形如圖 7所示，風(fēng)洞測得的近場過壓值分布(φ=0°,30° )如圖 8所示。采用的大氣參數(shù)為SBPW官方給出的標(biāo)準(zhǔn)大氣，空氣相對濕度為70%。

圖9分別給出了自研程序和sBOOM在不同周向角下計算得到的地面波形。可以發(fā)現(xiàn)，二者預(yù)測的地面波形外形相似。當(dāng)周向角為0°時(即飛機的正下方)，二者的計算結(jié)果基本重合。當(dāng)周向角為30°時，二者地面波形的最后一個上升段斜率稍有偏差。

2.5 Axibody算例驗證

Axibody為SBPW-2的第2個標(biāo)準(zhǔn)算例，其幾何外形就是一個簡單的幾何回轉(zhuǎn)體，如圖 10所示。風(fēng)洞測得的近場過壓值分布如圖 11所示，因為其對稱特性，各個周向角測得的近場過壓分布完全一致。采用的大氣參數(shù)和前一算例一致。

圖12分別給出了自研程序和sBOOM在不同周向角下計算得到的地面波形。可以發(fā)現(xiàn)，和LM1021算例的結(jié)果類似，二者的音爆預(yù)測結(jié)果達(dá)到了基本相當(dāng)?shù)木取?/p>

3 增廣Burgers方程在音爆模擬中的應(yīng)用

基于上述自研程序，首先研究了經(jīng)典吸收和分子馳豫效應(yīng)這兩種大氣聲吸收[24]效應(yīng)對地面波形的影響，接著通過改變大氣參數(shù)，分別模擬了超聲速客機在大氣層中飛行可能遇到的干燥、濕潤、高溫以及低溫的大氣環(huán)境，研究并比較了濕度和溫度對于地面波形的影響。

3.1 大氣聲吸收對地面波形的影響

增廣Burgers方程將經(jīng)典吸收和分子馳豫效應(yīng)這兩種大氣聲吸收效應(yīng)納入考慮。理論[25]和試驗[26]均證實，在低于10 MHz的頻率下，經(jīng)典吸收和分子馳豫效應(yīng)是可以疊加的，因此在數(shù)值模擬中可以對其進(jìn)行解耦。經(jīng)典吸收是指聲能因大氣的黏性和導(dǎo)熱性所造成的損失，而分子馳豫效應(yīng)是因部分聲能從聲波的集總運動轉(zhuǎn)移到分子內(nèi)自由度中(轉(zhuǎn)動和振動)而造成的[24]。這兩種效應(yīng)都導(dǎo)致了聲波能量在傳播的過程中被耗散為內(nèi)能。

圖13給出了分別使用增廣Burgers方程法和波形參數(shù)法計算上文提到的LM1021這一算例的結(jié)果，周向角取0°?？梢园l(fā)現(xiàn)，由于考慮了傳播過程中的損失，所以通過增廣Burgers方程計算得到的遠(yuǎn)場信號中的壓縮波是有厚度的，可以直接得到音爆的上升時間。而傳統(tǒng)的波形參數(shù)法并不考慮聲能在傳播過程中因耗散而造成的損失，計算得到的地面波形是自然界中并不存在的無厚度的壓縮波，需要通過一些經(jīng)驗公式來估計上升時間。這是增廣Burgers方程類方法與波形參數(shù)法等線化理論方法最大的不同。另外由于耗散的原因，增廣Burgers方程計算得到的過壓值的峰值較低。綜上所述，考慮了傳播介質(zhì)中的損耗的基于增廣Burgers方程的模擬方法能夠更真實地反映聲爆在大氣層中傳播的物理過程。

圖13還比較了經(jīng)典吸收和分子馳豫效應(yīng)各自對地面波形的影響?？梢园l(fā)現(xiàn)，如果在增廣Burgers方程中僅考慮大氣經(jīng)典吸收效應(yīng)的影響，最終得到的波形和不考慮大氣聲吸收的波形參數(shù)法的結(jié)果比較接近，波形的壓縮波厚度并不明顯。如果僅考慮分子馳豫效應(yīng)的影響，相比于僅考慮大氣經(jīng)典吸收的情況，壓縮波的厚度明顯增加，并且最終的波形與增廣Burgers方程的結(jié)果已經(jīng)非常接近?；谝陨系慕Y(jié)果可以得出結(jié)論：分子馳豫效應(yīng)對于音爆信號的耗散更為強烈。

3.2 濕度對地面波形的影響

未來超聲速客機力求做到能在大陸上空超聲速飛行，沿海和內(nèi)陸的大氣濕度相差很大，而濕度和分子馳豫效應(yīng)的強度相關(guān)，所以有必要研究濕度對地面波形的影響。

為了簡化研究，從近場到地面的大氣相對濕度取一常數(shù)，選取了3個相對濕度(10%、50%、90%)來進(jìn)行計算，其中10%的相對濕度模擬了沙漠、戈壁等干旱地區(qū)的濕度，90%的相對濕度則接近雨林、沿海等濕潤地區(qū)的最高濕度。在研究的過程中固定溫度、氣壓分布為標(biāo)準(zhǔn)大氣分布。

圖14給出了音爆信號在上述3種相對濕度的大氣中傳播得到的地面波形。總體來說，濕度對于波形幅值的影響不大。相對而言，干燥的空氣對于聲波的耗散更強，地面信號的過壓值峰值更低。在起耗散作用的經(jīng)典吸收和分子馳豫效應(yīng)中,濕度的降低僅僅會使分子振動馳豫效應(yīng)增強，而不會改變經(jīng)典吸收的強度，因此總的聲吸收強度會增強。

3.3 溫度對地面波形的影響

基于和3.2節(jié)相同的原因，對溫度對地面波形的影響作了研究。

以SBPW-2官方標(biāo)準(zhǔn)大氣的溫度分布為基準(zhǔn)，分別給所有海拔處的溫度增加15 ℃，使地面溫度為40 ℃以模擬高溫地區(qū)的溫度。再給所有海拔處的溫度相對標(biāo)準(zhǔn)大氣減少15 ℃，使地面溫度為-10 ℃以模擬低溫地區(qū)的溫度。在研究的過程中固定相對濕度(70%)、標(biāo)準(zhǔn)大氣分布等環(huán)境參數(shù)。

圖15給出了音爆在上述高溫、低溫以及標(biāo)準(zhǔn)的大氣中傳播得到的地面波形?？傮w來說，溫度對地面波形的影響比濕度更加明顯。相對而言，低溫大氣對聲波的耗散更強，傳播得到的最大過壓值最低。在70%的相對濕度下，考慮起耗散作用的經(jīng)典吸收和分子馳豫效應(yīng)，溫度的降低雖然會導(dǎo)致經(jīng)典吸收效應(yīng)減弱，但是低溫亦導(dǎo)致了分子振動馳豫效應(yīng)的增強，且因其增加的耗散比因經(jīng)典吸收效應(yīng)減弱而減少的耗散要多，這兩種效應(yīng)的疊加最終使總耗散增強。同時可以發(fā)現(xiàn)：低溫的大氣增強了傳播中的非線性效應(yīng)，這將導(dǎo)致上升時間的減小。

4 結(jié) 論

1) 相比于波形參數(shù)法等線化理論方法，增廣Burgers方程類方法可以更準(zhǔn)確地模擬音爆信號在真實大氣中傳播的各個物理過程，能夠有效提高上升時間的模擬精度。

2) 在近場聲壓信號前加入一段無幅值的“緩沖信號”可以有效提升地面波形上升時間的模擬精度。

3) 基于算子分解法的增廣Burgers方程的數(shù)值方法對于時間方向的網(wǎng)格數(shù)量更敏感，需要對其進(jìn)行適當(dāng)加密，建議的時間采樣率為85 kHz以上，否則可能導(dǎo)致波形的失真。

4) 在經(jīng)典吸收和分子馳豫效應(yīng)這兩種大氣環(huán)境造成的聲吸收機制中，分子馳豫效應(yīng)對聲能的損耗更大。

5) 干燥、低溫的大氣對地面音爆信號的過壓值有抑制作用，但低溫大氣會減少音爆信號的上升時間。

在本文工作的基礎(chǔ)上，未來擬開展結(jié)合音爆、氣動的多學(xué)科優(yōu)化設(shè)計[27]及考慮不確定因素的穩(wěn)健設(shè)計[28-29]。