顧奕然，黃江濤，陳樹生，3，*，劉德園，高正紅

1.西北工業(yè)大學 航空學院，西安 710072

2.中國空氣動力研究與發(fā)展中心，綿陽 621000

3.航空工業(yè)第一飛機設計研究院，西安 710089

受到諸多工程技術因素的影響，超聲速客機自20 世紀70 年代出現以來發(fā)展極其緩慢，而聲爆帶來的地面噪聲問題是超聲速客機發(fā)展的重要障礙之一［1］。以歐洲聯合研發(fā)的超聲速客機協和號為例，其巡航狀態(tài)下地面感覺噪聲約為105 dB［2-3］，這會對航線下方的生物和建筑造成難以忽視的影響。

飛行器在超聲速飛行下于近場產生的擾動波經過大氣傳播，會在遠場演化成形如英文字母“N”的遠場波形［4］，如圖 1 所示。聲爆在遠場演化的N 形聲波直接影響著地面，因此低聲爆設計工作的最終目的即通過對聲爆產生的N 形波進行修型，降低地面感覺噪聲。依據NASA 在SSBD（Shaped Sonic Boom Demonstrator）［5］項目中的研究成果，認為對超聲速飛行器的外形進行優(yōu)化設計能夠改變其對應的地面聲爆波形，這為聲爆的設計優(yōu)化提供了基礎和參考，意味著在低聲爆設計工作中引入氣動外形優(yōu)化的相關方法是行之有效的。

圖1 聲爆傳播過程示意圖Fig.1 Schematic diagram of sonic boom propagation

低聲爆設計方法依據其推進的方向，可以大致分為正設計［6-9］和反設計［10-13］兩類。正設計通過針對超聲速飛行器外形的優(yōu)化設計，經過多次迭代計算從而得到符合要求的遠場聲爆波形；反設計則是通過選擇一個符合要求的遠場聲爆波形，利用相應的反演技術反向推演出超聲速飛行器的氣動外形。通常情況下，人為設計一個理想波形然后反向尋找其對應的飛行器氣動外形，比使用正向算法搜索最優(yōu)的地面波形更加高效，因此反設計方法在低聲爆設計領域受到極大重視。

有關聲爆反設計的方法通?？梢苑譃榈孛婺繕瞬ㄐ纹ヅ洌?4-15］和近場目標波形匹配［16-17］兩類。地面目標波形匹配方法通過指定一個目標地面聲爆波形，并使用反演技術將地面聲爆波形反向傳播至近場，以完成反設計任務；近場目標波形匹配方法指定一個目標近場過壓分布，然后通過反設計算法來尋找其對應的氣動外形。由于飛行器的氣動外形和近場過壓分布是緊密耦合的，因此通常情況下，實現近場目標波形匹配更為容易。但倘若考慮大氣傳播過程，一個優(yōu)化后的較好的低聲爆近場波形未必能對應一個較好的遠場波形，這限制了近場目標波形匹配方法的應用。與此相比，地面目標波形匹配方法能夠更準確地反映地面聲爆波形與飛行器外形的關聯性，但其對聲波大氣傳播過程的模擬技術要求更為苛刻。

對于模擬聲波在大氣層內的傳播過程，通?？梢杂? 種方式實現：一為基于弱激波理論的波形參數法［18］；二為基于時頻域分析技術的Pestorius-Anderson 方法［19］；三為基于非線性聲學的增廣Burgers 方程方法［20］。3 種方法中，前2種方法由于需要引入經驗公式［21］和需要交替進行時域頻域變換，模擬精度往往難以達到工程應用的需求，因此增廣Burgers 方程方法使用更為廣泛［22］。

反設計工作的重心在于尋找與目標波形相對應的近場波形?；谙嗖ㄐ螀捣?，Seebass 和George［23］提出了一種相對較為粗糙的反設計方法?；谠鰪VBurgers 方程方法，Li 等［24］提出了一種基于逆向增廣Burgers 方程的方法來反向傳播聲波并以此完成反演計算。Rallabhandi［25］發(fā)展了基于反向等效面積目標的離散伴隨方法。張繹典和高正紅等［26］提出了基于本征正交分解（Gappy Proper Orthogonal Decomposition，GPOD）的聲爆反設計方法。各方法均能在一定程度上使設計外形的地面波形接近目標波形，從而降低地面聲爆強度。

本文基于一種逆向增廣Burgers 方程的數值反演方法展開研究，對求解過程中的各影響因素進行探索，并在已有研究成果的基礎上，進一步對基于此方法的中場反設計方法進行了研究，對聲爆傳播過程中高頻分量的演化及其對反設計方法的影響進行了討論。

1 增廣Burgers 方程

增廣的Burgers方程由Pierce［27］于20世紀80年代提出，在經典Burgers 方程中引入了額外的幾何擴散、大氣分層及分子弛豫項，能夠更好地模擬聲壓信號在大氣中的傳播。逆向增廣Burgers 方程可用于模擬聲爆信號的反向傳播，能夠將中場壓力分布反向傳播到近場位置。本節(jié)描述逆向增廣Burgers 方程及其求解方法。

1.1 經典正向增廣Burgers 方程

增廣Burgers 方程包括非線性、經典耗散、分子弛豫現象等效應。無量綱的Burgers 方程如式（1）所示：

式中：P為無量綱壓力；τ為無量綱時間；σ為無量綱聲管距離；Γ、Cv、θv分別為無量綱氣體耗散參數、無量綱松弛系數和無量綱松弛時間；S、c0與ρ0為聲管面積、聲速和大氣密度。

方程（1）右側5 項分別為非線性效應項、經典耗散效應項、分子弛豫效應項、幾何擴散效應項和大氣分層效應項，分別表示對應的物理效應對聲壓信號的影響。

1.2 逆向增廣Burgers 方程

經典正向增廣Burgers 方程中，非線性效應項、經典耗散效應項和分子弛豫效應項與時間有關，而幾何擴散效應項和大氣分層效應項只與聲管路徑有關。構建逆向方程可以看作令時間反向流動，因此只需要修改非線性效應項、經典耗散效應項和分子弛豫效應項。

對于波形的反向傳播，增廣Burgers 方程可以修改為

逆向增廣Burgers 方程與正向方程類似，在形式上只有非線性效應項、經典耗散項和分子弛豫項的正負發(fā)生了變化［24］。

1.3 數值求解方案

在對逆向增廣Burgers 方程進行求解之前，首先需要對時間和空間進行離散。時間維度上需要對過壓信號進行離散處理；空間維度上需要對聲管傳播路徑進行離散處理。使用均勻的網格對時間和空間進行離散，并使用參考時間1/ω0，參考長度為時間采樣頻率，β為非線性系數）和參考壓力p0進行無量綱化，無量綱化過程由式（3）～式（5）所示：

使用算子分裂法［28］對逆向增廣Burgers 方程進行求解。式（2）可由算子分裂法分裂成以下5 項：

求解過程中，依次使用式（6）～式（10）更新無量綱壓力P，推進至下一高度后繼續(xù)循環(huán)上述過程，直至高度推進至目標高度，從而完成求解。

對于式（6），其等價于：

使用以下離散方法進行求解：

式（7）與式（8）是病態(tài)方程，為穩(wěn)定方程求解過程，可以引入偽拋物線型正則化方法改寫方程，其形式如式（14）所示：

正則化參數ε是一個正的標量，使得式（14）數值可解，可以通過試錯法選取ε=10-4。使用偽拋物線型方程改寫經典耗散項和分子弛豫項，能夠較好的穩(wěn)定數值求解的過程。

使用待定系數法改寫經典耗散項，可以寫為

分子弛豫項也可使用類似的方法寫為

式中：α取0.5。

式（15）與式（17）均為三對角方程，可以使用TDMA 方法（Tridiagonal Matrix Algorithm）高效率求解。

式（9）與式（10）本質是一個守恒方程，可以分別使用式（20）與式（21）進行求解。

1.4 反向等效面積求解

通過使用前文中討論的反向傳播方法，中場位置的壓力信號分布可以反向傳播到近場位置，近場過壓分布可以逆變換為等效面積分布，稱為反向等效面積Ae，r。使用以下等式［24］：

式中：Ma為巡航馬赫數；R為近場位置到飛行器的距離；γ取1.4。需要注意的是，只有當飛行器滿足細長體假設［29］時式（22）才成立。

2 計算程序驗證

基于增廣Burgers 方程及逆向增廣Burgers方程，可以編寫計算程序對聲爆過壓信號在近場與遠場之間的傳播進行研究。

正向增廣Burgers 方程計算程序參考筆者所在實驗室自研聲爆程序［22］，圖 2 展示了正向增廣Burgers 方程計算程序的主要框架。過壓信號由近場傳播向遠場，首先使用聲管計算程序輸出對應的聲管路徑和大氣參數，并將其提供給求解器。求解器使用算子分裂法求解，輸入過壓信號依次通過非線性效應項求解器、經典耗散項求解器、分子弛豫項求解器、幾何擴散項求解器和大氣分層效應項求解器，然后推進至下一高度，繼續(xù)以上述順序進行迭代，直至高度推進至目標高度。

圖2 正向程序結構Fig.2 Forward program structure

圖3 展示了逆向增廣Burgers 方程計算程序的主要框架。與正向求解程序類似，在聲管計算程序輸出相應計算參數后，使輸入過壓信號依次通過逆向非線性效應項求解器、逆向經典耗散項求解器、逆向分子弛豫項求解器、幾何擴散項求解器和大氣分層效應項求解器，然后推進至下一高度，繼續(xù)以上述順序進行迭代，直至高度推進至近場高度。求解獲得近場過壓信號后，通過Abel 逆變換計算得到對應的逆向等效面積分布，完成計算任務。

圖3 逆向程序結構Fig.3 Reverse program structure

為驗證計算程序的有效性，選取第2 次AIAA 聲爆預測研討會（SBPW-2）提供的官方算例Axibody 來進行驗證工作。Axibody 是一個幾何回轉體，其外形大致如圖 4 所示［30］。

圖4 Axibody 外形［30］Fig.4 Appearance of Axibody［30］

2.1 正向計算程序驗證

選取飛行高度H=15 849.6 m，飛行馬赫數Ma=1.6，近場與機體距離R=128.93 m。Axibody 在周向角0°時的標準近場過壓信號如圖 5所示。

圖5 Axibody 近場過壓分布Fig.5 Near-field over-pressure distribution of Axibody

由正向增廣Burgers 計算程序計算所得的地面波形和NASA 的sBOOM 軟件的計算結果對比圖如圖 6 所示，驗證了編寫的正向增廣Burgers計算程序計算準確度符合要求。

圖6 地面信號對比Fig.6 Ground signal comparison

2.2 逆向計算程序驗證

由于經典耗散效應和分子弛豫效應的影響，近場過壓信號中的部分分量在向遠場傳遞的過程中會不可避免地出現耗散和消亡的現象［31］；此外，在逆向求解數值方法中，逆向過程會趨向于抹平波形，這使得使用逆向增廣Burgers 方程難以將遠場信號準確地反演至近場［24］。但倘若縮減逆向計算程序的反演距離，從一個距離近場較近的中場過壓信號開始反演，也可以獲得大致準確的近場過壓信號。

選取3 個不同的中場高度Hm=15 500 m、15 000 m、14 000 m 進行驗證。使用正向增廣Burgers 計算程序計算獲得對應的中場過壓信號如圖 7 所示，再使用逆向增廣Burgers 計算程序反演至近場，其計算結果與標準數據對比圖如圖 8 所示，反演過壓信號對應的等效面積分布與標準近場過壓信號的等效面積分布的對比圖如圖 9 所示?？梢杂^察到，雖然隨著反演距離的增長計算精度有所降低，但逆向增廣Burgers 方程計算程序的計算結果與標準近場數據大致相同，可以認為逆向Burgers 方程反演得到的信號在一定精度范圍內可以代表近場信號。

圖7 中場過壓信號Fig.7 Mid-field over-pressure signal

圖8 近場過壓分布對比Fig.8 Comparison of near-field over-pressure distributions

圖9 等效面積分布對比Fig.9 Comparison of equivalent-area distributions

2.3 逆向計算程序網格收斂性

在進行逆向增廣Burgers 方程計算時，不僅需要對聲管傳播路徑進行離散，還要對過壓信號進行離散，即在空間和時間上各有一個離散維度。

采用控制變量法對網格收斂性進行研究。選取中場高度Hm=14 000 m，首先固定空間網格數為10 000，使用不同的網格數量離散時間維度，其對應的計算結果對比圖如圖 10 所示。然后固定時間網格數為10 000，使用不同的網格數量離散空間維度，其對應的計算結果對比圖如圖 11所示。

圖10 不同時間網格數量計算結果對比Fig.10 Comparison of results with different grid numbers in time dimension

圖11 不同空間網格數量計算結果對比Fig.11 Comparison of results with different grid numbers in spatial dimension

觀察計算結果，認為時間采樣率提升至25 kHz 時，計算結果趨向于收斂，而空間網格的密度對結果無顯著影響。比較兩組計算結果，可以發(fā)現時間網格密度對計算結果的影響較大，而空間網格密度對計算結果影響較小。這和正向Burgers 方程離散求解的特點相似［22］。

3 數值結果

在第2 節(jié)中，研究已經針對正向和逆向的增廣Burgers 方程計算程序進行了準確度驗證和網格收斂性研究，認為正向和逆向的計算程序均能較為準確地完成計算任務?；谏鲜鲇嬎愠绦?，本節(jié)對方程求解實際問題中各因素的影響進行了較為系統的研究，對在不同中場高度與不同周向角條件下使用逆向增廣Burgers 方程計算程序進行反演計算，并對計算結果進行分析。

3.1 反演距離

前文已經提到，由于近場過壓信號在向遠場傳遞的過程中，不可避免地受到經典耗散效應和分子弛豫效應的影響，使得逆向增廣Burgers 方程難以準確地將遠場聲壓波形反演為對應的近場過壓信號。而若縮短反演距離，則中場過壓信號中因經典耗散效應與分子弛豫效應而消亡的頻率分量將減少，因此從一個距離較近的中場利用逆向增廣Burgers 方程進行反演計算能夠取得更高的準確度。反演距離對于逆向增廣Burgers方程的反演具有重要影響。

研究采用LM1021、C608 兩個標準算例進行研究。其外形如圖 12［30］與圖 13［32］所示。LM1021為SBPW-1 的可選算例和SBPW-2 的標準算例之一，是洛克希德·馬丁公司設計的新一代超聲速客機；C608 是SBPW-3 的標準算例之一，是NASA的低聲爆試驗飛行器X59的初步改形方案。

圖12 LM1021 飛行器［30］Fig.12 LM1021 aircraft［30］

圖13 C608 飛行器［32］Fig.13 C608 aircraft［32］

對于LM1021 算例，選定的計算條件為飛行高度H=16 764 m，近場與機體距離R=222.59 m，飛行馬赫數Ma=1.6。圖 14 與圖 15展示了LM1021 算例在不同中場高度開始反演計算所得到的近場過壓信號對比與對應的等效面積分布。

圖14 不同中場高度反演結果對比（LM1021）Fig.14 Comparison of inversion results at different mid-field heights（LM1021）

圖15 不同中場高度反演等效面積對比（LM1021）Fig.15 Comparison of equivalent-areas at different mid-field heights（LM1021）

對于C608 算例，選定的計算條件為飛行高度H=16 215.36 m，近場與機體距離R=82.296 m，飛行馬赫數Ma=1.4。圖 16 與圖 17展示了C608 算例在不同中場高度開始反演計算所得到的近場過壓信號對比與對應的等效面積分布。

圖16 不同中場高度反演結果對比（C608）Fig.16 Comparison of inversion results at different mid-field heights（C608）

圖17 不同中場高度反演等效面積對比（C608）Fig.17 Comparison of equivalent-areas at different mid-field heights（C608）

觀察計算結果，可以觀察到隨著中場高度的增大，逆向增廣Burgers 方程計算程序的計算結果越來越逼近標準近場過壓信號，這與理論分析的結果相一致。同時可以注意到，雖然不同中場高度的選取對反演的近場過壓信號波形有較大影響，但相應的等效面積分布變化相對較小，這使得基于逆向增廣Burgers 方程的低聲爆反設計方法能夠提供較為準確的低聲爆等效面積分布。

此外，從圖 14 與圖 16 中可以觀察到，反演結果的誤差主要來自于無法對過壓信號中的高頻分量進行準確的還原，且當中場位置靠近機體時，反演準確度會迅速增大。圖 18 與圖 19 展示了由正向增廣Burgers 方程計算程序得出的LM1021 與C608 算例距離近場約1 000 m 內聲爆過壓信號的演變過程?？梢杂^察到，在聲爆傳播過程初期，非線性效應的扭轉作用和大氣分層效應的作用效果尚不明顯，過壓信號的演化主要受經典耗散和分子弛豫效應對高頻分量的耗散作用以及聲管幾何擴散對聲爆信號的線性衰減作用的影響。由于逆向增廣Burgers 方程計算程序和正向程序在處理幾何散布效應的方式相同，不會對準確度產生較大影響，因此可以認為在距離機體較近的情況下，中場位置對反演精度的影響主要是由經典耗散效應和分子弛豫效應引起的。

圖18 LM1021 機體附近過壓信號演化過程Fig.18 Evolution process of over-pressure signals near LM1021 body

圖19 C608 機體附近過壓信號演化過程Fig.19 Evolution process of over-pressure signals near C608 body

3.2 不同周向角下反演效果

周向角用以表示近場過壓信號位置相對于機體的位置，其示意圖如圖 20 所示。通常情況下，飛機在不同周向角下的近場過壓信號是不同的，且其向遠場傳播的聲管路徑也不同。但值得注意的是，若飛機的外形為軸對稱的回轉體，則其不同周向角下的近場過壓信號相同。研究選取LM1021 與Axibody 兩個標準算例進行研究。

圖20 周向角相對于機體位置Fig.20 Rolling angles relative to body position

圖21～圖 25 展示了LM1021 算例在不同周向角下逆向增廣Burgers 方程計算程序的反演結果。飛行高度、近場距離和飛行馬赫數與3.1 節(jié)中相同，選取的中場高度為Hm=16 400 m。

圖21 周向角0°下近場波形對比（LM1021）Fig.21 Comparison of near-field waveforms at rolling angle 0°（LM1021）

圖22 周向角10°下近場波形對比（LM1021）Fig.22 Comparison of near-field waveforms at rolling angle 10°（LM1021）

圖23 周向角20°下近場波形對比（LM1021）Fig.23 Comparison of near-field waveforms at rolling angle 20°（LM1021）

圖24 周向角30°下近場波形對比（LM1021）Fig.24 Comparison of near-field waveforms at rolling angle 30°（LM1021）

圖25 周向角50°下近場波形對比（LM1021）Fig.25 Comparison of near-field waveforms at rolling angle 50°（LM1021）

圖26 展示了Axibody 算例在不同周向角下逆向增廣Burgers 方程計算程序的反演結果。飛行高度、近場距離和飛行馬赫數與第2.1 節(jié)中相同，選取的中場高度為Hm=15 500 m。由于Axibody為回轉體外形，其不同周向角下的近場過壓信號相同，因此可以直接比較不同周向角方向反演過壓信號的不同。

圖26 不同周向角下近場波形對比（Axibody）Fig.26 Comparison of near-field waveforms at different rolling angles（Axibody）

由結果可以觀察到，逆向增廣Burgers 方程計算程序的反演準確度隨著周向角的減小而有少量的提升，在周向角為0°時反演準確度最高?？紤]到在低聲爆反設計過程中，往往更為關注低周向角下聲爆波形的傳播，因此認為周向角對反演準確度的影響并不會對反設計任務造成過大的影響。

4 基于中場目標波形的優(yōu)化設計方法可行性

在低聲爆反設計方法中，通常先選取一個低噪聲的地面遠場目標波形，然后尋找到與其對應的近場波形，再利用Abel 逆變換獲得對應的等效面積分布，并以此來指導飛行器的外形優(yōu)化，其流程圖如圖 27 所示。由于前文提到的原因，逆向增廣Burgers 方程難以將地面聲壓波形準確地反演至近場。因此，可以嘗試在中場設計目標波形，以減輕反演過程帶來的誤差，其流程圖如圖 28 所示。本節(jié)對基于增廣Burgers 方程的中場目標波形方法展開研究，驗證中場目標波形方法的可行性，并對聲爆信號中的高頻分量對優(yōu)化設計過程的影響進行了探討。

圖27 基于地面目標波形的設計方法Fig.27 Design method based on ground target waveform

圖28 基于中場目標波形的設計方法Fig.28 Design method based on mid-field target waveform

基于Plotkin 等［33］的研究成果，將低聲爆地面目標波形修改為一個類正弦波形，能夠有效地抑制聲波中的高頻分量的能量，從而降低人耳聽覺頻率范圍內的噪聲分量。研究嘗試將這一結論用于中場目標波形的設計，使用一個類正弦的中場波形作為目標波形。選取LM1021 作為算例，飛行高度、近場與機體距離和飛行馬赫數的選取與3.1 節(jié)相同，選取中場高度為Hm=14 500 m。利用正向增廣Burgers 方程計算程序獲取LM1021 在中場高度的聲壓波形，在此基礎上設計了一個類正弦的目標中場波形，如圖 29 所示。并將中場目標波形使用逆向增廣Burgers 方程反演至近場，生成近場目標波形。

圖29 中場目標波形Fig.29 Mid-field target waveform

由于物理約束的限制，通常很難通過對飛行器外形的優(yōu)化設計，將對應的近場過壓信號波形優(yōu)化至和近場目標波形一致。通常情況下，近場過壓信號在峰值附近會出現不可避免的高頻振蕩。為研究這一特性對優(yōu)化過程的影響，研究在目標近場信號的峰值附近添加了頻率約為100 Hz、200 Hz 和300 Hz 的振蕩信號，以及200 Hz下最大 振幅為1.8×10-3、2.7×10-3與3.6×10-3的振蕩信號。將合成后的新波形作為假想的優(yōu)化后近場波形，合成近場波形和目標近場波形如圖 30 與圖 31 所示。

圖30 不同頻率的合成近場波形Fig.30 Synthetic near-field waveforms of different frequencies

圖31 不同幅值的合成近場波形Fig.31 Synthetic near-field waveforms of different amplitudes

將合成近場波形和目標近場波形使用正向增廣Burgers 方程計算程序傳播至遠場，其計算結果如圖 32～圖 34 所示。觀察計算結果，可以發(fā)現中場目標波形在傳遞至遠場后，其形狀依然是一個類正弦波形，說明Plotkin 等［33］的研究結論也可用于中場目標波形的設計。而合成近場波形經過傳播，可以發(fā)現同頻率不同幅值的附加信號也并非按照相同的比例線性衰減，由于非線性效應對波形的扭轉作用，大振蕩幅度的高頻分量會改變波形的整體形狀。此外，也能觀察到波形信號內高頻分量的衰減幅度遠大于低頻基波分量，由于經典耗散效應和分子弛豫效應的存在，大氣對高頻聲波分量的吸收作用更為明顯，使得大氣具有類似低通濾波器的性質［34］，而程序計算結果也與這一理論分析相吻合。

圖32 遠場目標波形Fig.32 Far-field target waveform

圖33 不同頻率的合成遠場波形Fig.33 Synthetic far-field waveforms of different frequencies

圖34 不同幅值的合成遠場波形Fig.34 Synthetic far-field waveforms of different amplitudes

合成近場波形傳播至遠場后，其波形與目標遠場波形大致相同。這說明對優(yōu)化時近場波形中的高頻分量無需過度關注，只需約束近場波形中的低頻和中頻分量，并對高頻分量的幅值進行一定的限制，使近場波形與目標近場波形相似，便能較好地完成優(yōu)化任務。

5 結論

1）使用基于算子分裂法的逆向增廣Burgers 方程數值求解方法在一定反演距離范圍內能夠取得較為準確的結果，且對于時間方向的網格密度較為敏感。

2）不同中場高度的選取對于逆向增廣Burgers 方程的數值求解有較大影響，中場高度越高，反演距離越短，反演的準確度越高。中場高度的選取對等效面積分布計算的影響較近場過壓信號計算的影響更小。

3）不同的周向角對于逆向增廣Burgers 方程的數值求解具有一定影響，在周向角較小時反演精度較高。

4）采用中場目標波形代替地面目標波形進行低聲爆反設計工作，能減輕反演距離對反演計算過程的影響。

5）近場過壓信號中的高頻分量對地面噪聲具有一定的影響，應在優(yōu)化設計過程中對高頻分量的幅值和頻率進行一定的限制。