基于增廣Burgers方程的超聲速客機遠場聲爆預測研究

崔青，白俊強，宋源，余培汛

（1.航空工業(yè)空氣動力研究院氣動研究與試驗一部，沈陽110000）

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

0 引言

超聲速客機能夠有效地縮短長距離飛行所花費的時間，提高出行效率，其優(yōu)點在越洋出行中更加明顯。因此各大國和國際機構對于超聲速客機的研究一直沒有中斷，而如何降低超聲速客機的聲爆水平和超聲速客機的航線規(guī)劃也一直是限制超聲速客機投入大規(guī)模應用的瓶頸［1］。如何精確預測聲爆強度是其中的關鍵技術，全場直接計算法因其過于耗費計算資源，目前沒有工程實用價值［2］。近場CFD計算與遠場聲學傳播計算相結合的方法是目前國際上常采用的高精度聲爆預測方法。其中遠場聲學傳播計算的方法主要分為波形參數(shù)法［3］和Burgers方程法［4］。

波形參數(shù)法所求得的遠場過壓波形中，激波是間斷的，且不包含激波上升時間，進行聲強級計算前需依據(jù)實驗數(shù)據(jù)進行人工激波增厚；波形參數(shù)法并沒有考慮大氣濕度的影響，模型具有理論缺陷［5］。但Burgers方程在求解過程中考慮了大氣濕度帶來的影響以及還有其他效應，求解遠場過壓的同時能求解激波上升時間。

國外對Burgers方程的研究較早，由于Burgers方程模型比較簡單，只能預測有限振幅平面波在熱黏性耗散介質中的傳播。這種簡單模型只能模擬傳播過程中的非線性效應和熱黏性吸收效應，被稱作經典Burgers方程［6］，之后有研究者提出了廣義Burgers方程［7］。和經典Burgers方程相比，廣義Burgers方程加入了傳播過程中幾何擴散效應和大氣分層效應的模擬。廣義Burgers方程的進一步發(fā)展就是增廣Burgers方程，增廣Burgers方程在廣義Burgers方程的基礎上加入了對分子遲豫效應的模擬［8］。至此，Burgers方程對聲爆傳播過程的模擬基本完善。現(xiàn)在常說的增廣Burgers方程是一個復雜的偏微分方程，右端有五項，包含之前提到的五種效應，對其進行整體的離散和數(shù)值求解有一定難度。直到1995年，R.O.Cleveland［9］提出了算子分裂法，該方法最早被Y.S.Lee和M.F.Hamil?ton用于求解KZK方程［10］。算子分裂法的 核心便是在小推進步長上依次單獨離散求解右邊的每一項，每一項求解的結果作為下一項求解的輸入，經驗證這種算法的計算結果收斂于整體離散求解的結果［11］。大量基于增廣Burgers方程的聲爆預測程序均采用了Cleveland的算子分裂法［5，12］，但對于整體求解的探索并沒有停止。2015年，M.Yama?moto等［13］提出了一種使用內部變量的變換，經過這種變換，右側項可以聯(lián)合求解，和算子分裂法相比，這種算法更加穩(wěn)定且計算速度更快；S.Rallab?handi［14］主 要 將 精 力 集 中 在Burgers方 程 和 外 形 優(yōu)化相結合，同時探究了傳播過程中因地球三維效應帶來的影響［15］，研究結果顯示三維效應對于無風大氣中的傳播沒有太大影響，但當大氣風速存在時，三維效應對于傳播路徑會有較大影響，甚至會得到與無風大氣相反的結論。

近年來，國內也有不少研究者對增廣Burgers方程開展了較深入的研究。張繹典等［4］構建了基于增廣Burgers方程的遠場求解方法，利用第二屆國際聲爆研討會（SBPW-2，簡稱：大會）給出的標準算例對其精度進行了驗證，同時探究了近場過壓信號的時間采樣頻率和遠場信號波形失真程度的關系，為了保證波形不發(fā)生過分的失真，信號的采樣頻率不應低于85 kHz，對大氣濕度與遠場過壓峰值的關系進行了初步研究；喬建領等［16-17］構建了基于增廣Burgers方程的近遠場結合的聲爆預測方法，利用大會給出的標模算例對近場、遠場的求解精度均進行了驗證，最后基于該預測方法，探究了風速對聲爆幅值和地面影響區(qū)域的影響，但是在相同推進步長情況下，激波插入精度有待提高。

本文探究基于增廣Burgers方程的聲爆預測方法的求解過程，提出一種有效的非線性效應求解方法；在算子分裂法的基礎上構建聲爆遠場預測方法，基于該方法編寫聲爆預測程序，通過該程序計算大會的兩個標準算例［18-19］以驗證該程序的正確性；分析空間、時間網格密度及緩沖信號對預測結果的影響。

1 方程求解

增廣Burgers方程來源于經典Burgers方程，其有量綱形式如式（1）所示。式中：x為從出射點開始的聲線傳播距離；p′為聲爆信號有量綱過壓；t′為延遲時間，是為了得到行波解而建立的新時間坐標系下的時間；ρ0，c0分別為巡航高度處于受擾動大氣密度與聲速；β，b分別為非線性效應系數(shù)和熱黏性吸收效應系數(shù)。

式中：σ為無量綱聲線傳播距離；-x為激波形成距離為無量綱壓力；p0為參考壓強，一般取巡航高度處的靜壓；τ為無量綱延遲時間；ω0為輸入波形的角頻率，通過將輸入的波形看作一個周期求得；Γ為無量綱熱黏性系數(shù)，當Γ變大時，非線性效應逐漸占主導地位，當Γ變小時，熱黏性吸收效應逐漸占主導地位。

經典Burgers方程只包含非線性效應項與熱黏性效應項，添加幾何擴散效應項和大氣分層效應項后便得到廣義Burgers方程，其無量綱形式如式（3）所示。

式中：S為聲線管截面積。

在廣義Burgers方程中添加分子弛豫效應項，便得到增廣Burgers方式，形式為

式中：Cv為第v個弛豫過程的無量綱色散系數(shù)；θv為第v個弛豫過程的無量綱弛豫時間。

式（4）中右側項依次為非線性效應項，熱黏性吸收項，分子弛豫項，幾何擴散項效應和大氣分層效應項。根據(jù)算子分裂法，式（4）被分解為下列形式：

按照算子分裂法，求解時由上至下依次求解式（5）每一項，將式（5）中上一項的結果作為下一項的輸入，每將上述五項遍歷一次等同于在真實空間中推進一步。推進時沿聲線傳播方向推進直至到達指定高度即可。

式（5）第一項為非線性效應項，用于計算由于傳播中非線性效應累積導致的波形畸變，該效應是聲爆波形在傳播中變化的主要原因。

式（5）第二項為熱黏性吸收效應項，用于計算由于分子間摩擦導致的信號強度衰減。

式（5）第三項分子弛豫項與其他幾項形式上有一點差別，因為不同氣體分子弛豫效應不同，大氣整體弛豫效應與不同種類氣體占比有關，因此該項需要將所有與大氣整體弛豫效應有關的氣體進行計算。真實大氣中對弛豫效應影響較大的氣體有氮氣、氧氣、二氧化碳和水蒸氣。其中二氧化碳和水蒸氣由于在大氣中含量遠小于氮氣與氧氣含量，在弛豫效應計算中可被略去，但需要注意大氣濕度對氧氣與氮氣弛豫效應的影響沒有被忽略。因此第三項實際上是兩項求和的形式，出于簡化離散求解的目的，求解弛豫效應時依舊采取算子分裂法，即將氮氣與氧氣的弛豫效應分開計算，上一項的結果作為下一項的輸入，經算例驗證該算法能正確模擬大氣整體弛豫效應。

式（5）第四項為幾何擴散效應項，用于計算聲爆傳播過程中由于波振面擴大導致的衰減。受大氣分層的影響，聲爆信號的傳播路徑一般不是直線，因此衡量幾何擴散效應與信號的傳播路徑有關。常用方法是選取信號發(fā)生處相鄰的四條聲線，這四條聲線圍成一個管道，用與其中一條聲線垂直的平面與聲管所截得的面積作為該處聲線管面積。每推進一步會得到每一步的聲線管面積，用每一步的聲線管面積變化即可衡量傳播中的幾何擴散效應。

式（5）最后一項為大氣分層效應項，用于計算由于大氣參數(shù)如密度、溫度等隨高度變化導致的信號強度衰減。

其中除去非線性效應外，文獻［9］中已有詳細的離散和求解方法，本文不再贅述。

1.1 非線性求解

與普通聲波的傳播不同，聲爆的過壓值并不遠小于大氣靜壓值。聲爆在傳播過程中受到累積的非線性效應作用，波形會發(fā)生扭曲變形［20］。描述這種情況的解析解為Poisson解，具體形式如式（6）所示。

為了便于求解，將波形轉換到扭曲時間坐標下。

將式（7）代入式（6），得到：

即只要將波形轉換到扭曲時間坐標下，就能得到推進一步后的波形。為了后續(xù)計算方便，可將波形線性插值到均勻時間坐標。當推進步長過大或累積推進一定距離后，會出現(xiàn)多值現(xiàn)象，如圖1所示。

圖1 波形的多值現(xiàn)象Fig.1 Multivalue form of the wave shape

在真實傳播過程中并不會產生多值現(xiàn)象，非線性效應導致的波形扭曲，堆積會產生激波，進而消除多值現(xiàn)象。計算中需要插入激波來消除多值現(xiàn)象，激波插入位置可以使用面積相等原則［20］求得，如圖2所示，激波插入的位置使得激波所在直線和波形圍成的兩部分面積相等。實際插入激波情況和圖2類似，如圖3所示。

圖2 激波插入位置示意Fig.2 Position of the inserted shock wave

圖3 實際插入激波情況Fig.3 Insertion in reality

由面積相等原則可以得到S1=S2，雖然波形多值部分一般不會是折線，但用折線和三角形面積進行近似能滿足精度需求。推進步長不能過大，否則多值部分波形會更復雜以致不能使用折線和三角形面積進行近似，如圖4所示。如果繼續(xù)按三角形面積計算，所得過壓值會很快發(fā)散。

圖4 推進步長過大導致失效Fig.4 Long step causing failure

1.2 聲線追蹤及聲線管面積

根據(jù)文獻［9］的理論，在均勻無風介質中聲爆沿直線傳播，聲速沿高度線性分布時，聲爆沿一圓弧傳播。真實大氣存在分層效應和風速梯度，且聲速一般不會沿高度線性分布，此時聲線將會是一條任意曲線。如果不存在垂直方向風速，根據(jù)幾何聲學理論［21-22］聲線追蹤可按式（9）～式（10）計算。

式中：R為聲線位置矢量；W為當?shù)仫L速向量；N為波陣面單位法向量；I為單位矩陣；?為克羅克內積。

式中：Δt為相鄰兩推進步之間傳播所用的時間，Δt=Δz c0Nz；Δz為相鄰兩推進步之間聲線傳播距離；Nz為在這兩部間聲線波陣面單位法向量。

N需要在每次推進過程中更新，N更新需要Δt。但沒有N就不能求得Δt，如果推進步長足夠小，可以將上一步的N近似作這一步的N。本文采取一種預估—校正的方法進行聲線追蹤。首先設定一Δz，用上一步N和設定的Δz求得Δt；然后用該Δt計算出真實推進距離Δz，再進行高度更新；重復此步驟直到聲爆傳播到希望高度。使用該方法可以避免波陣面單位法向量N的近似帶來的誤差，同時可以適當放大推進步長。

聲線管是由四條相鄰的聲線所圍成的“管”，這 四 條 聲 線 分 別 從(ta，φ)，(ta，φ+Δφ)，(ta+Δta，φ)，(ta+Δta，φ+Δφ)這四個位置出射，位置關系如圖5所示。

圖5 聲線及聲線管面積示意Fig.5 Definition of ray and ray tube area

分別計算出這四條相鄰聲線在同一推進步時的位置坐標后，可以根據(jù)式（13）求得聲線管面積。

1.3 計算流程

計算聲爆的程序框架如圖6所示，其中輸入分為三部分：過壓信號輸入，風速輸入和大氣參數(shù)輸入。如果不輸入風速和大氣參數(shù)，計算將采用內置的無風標準大氣［23］。輸入過壓信號必須是無量綱化的過壓信號。

圖6 程序計算流程圖Fig.6 framework of the program

將過壓信號插值到均勻、無量綱的時間網格上，根據(jù)當前高度計算出對應大氣參數(shù)，再按照上述射線追蹤與面積計算方法更新聲線坐標和聲線管面積。接下來依次單獨計算非線性效應，熱黏性吸收，分子弛豫效應，幾何擴散效應和大氣分層效應。如果推進到指定高度（比如地面），程序停止，輸出結果；如果沒有推進到指定高度，重復上述過程直到指定高度。

采用指定有量綱推進距離再預估—修正的方法確定每次推進步長。當剩余高度不足指定距離時，將剩余高度無量綱化作為最后的推進步長代入程序計算。

2 標模驗證

在展示標模驗證結果之前，需要明確信號提取位置與指向的定義和風速與機頭指向的定義，再展示標模驗證的結果。

2.1 提取位置與指向

提取位置及指向角示意圖如圖7所示，即從飛機飛行正前方向后看去的視角。圖中虛線圓弧代表空間中距離飛機一定距離的圓柱面，近場信號在該面上提取［24］。提取距離一般用H L來衡量，H為提取位置距離飛機的距離，L為飛機長度。H L選取既要保證近場信號收斂，又不能過大，避免浪費時間。

圖7 提取位置及指向角示意圖Fig.7 Definition of extract position and pointing angle

飛機飛行的方向與正北方的夾角α如圖8所示，該角度用于將聲速表中沿正北、正東向的風速轉化到以飛機飛行方向為x正向的坐標系下。更新聲線坐標等計算均在該坐標下進行，該坐標為右手坐標系，當飛機向正東方飛行時，x軸正向為東，y軸正向為北，z軸正向向上。

2.2 LM1021標模驗證

LM1021是SPBW-2中遠場計算的標模算例之一，原型為洛克希德馬丁公司的一款超聲速客機，實驗模型外形如圖9所示。

圖9 LM1021外形［25］Fig.9 Shape of LM1021［25］

第二屆國際聲爆研討會提供了風洞實驗數(shù)據(jù)［26］和遠場數(shù)據(jù)［27］，提取位置H L=3.129 9，分別提取出射角為0°，30°的信號（如圖10所示）。

圖10 LM1021近場輸入信號Fig.10 Farfield signature of LM1021

對應的巡航高度為16 764 m，巡航馬赫數(shù)為1.6Ma，地面反射因數(shù)取為1.9，側滑角為0°，大會提供了多個大氣情況和風速表供選擇，本次選擇標準無風大氣進行驗證，0°和30°遠場波形對比如圖11～圖12所示。

圖11 0°遠場波形對比Fig.11 Comparison of farfield waveform of 0°

圖12 30°遠場波形對比Fig.12 Comparison of farfield waveform of 30°

圖11 ～圖12中，cal表示本程序計算結果，sboom表示大會程序提供的標準結果，可以看出：計算結果與大會提供標模值吻合良好。

2.3 C25D標模驗證

C25D是SPBW-2中的近場計算算例［28］（如圖13所示），但各參與者均進行了遠場傳播計算，大會也給出了重新插值后的標準過壓提取值。本文選用Aftosmis［29］的計算結果作為標準進行算例驗證，計算結果如圖14～圖16所示，其中，巡航高度為15 760 m，巡航馬赫數(shù)為1.6Ma，地面反射因數(shù)為1.9，側滑角為0°。

圖13 C25D外形Fig.13 Shape of C25D

圖14 C25D近場輸入信號Fig.14 Extracted signature of C25D

圖15 0°遠場波形對比（C25D）Fig.15 Comparison of farfield waveform of 0°（C25D）

圖16 30°遠場波形對比（C25D）Fig.16 Comparison of farfield waveform of 30°（C25D）

從圖14～圖16可以看出：計算結果與大會提供結果吻合良好，程序計算精度滿足要求。

3 緩沖信號的設置

在計算LM1021的遠場過壓信號時，如果采用大會提供的原始輸入信號，會發(fā)生激波上升時間不足和過壓幅值不足的情況。張繹典等［4］指出，如果在原始信號前加入一段過壓值為0的緩沖信號，情況會得到解決。

本文對該現(xiàn)象進行分析，在LM1021，C25D的原始信號前分別加入了長度為0.1 s，0.04 s的緩沖信號（長度約為原始信號長度1/5）作為對比，如圖17～圖18所示，圖中extended表示添加緩沖信號后的結果，original表示未添加緩沖信號的結果。需要說明的是，兩者原始信號在首尾已有一定長度的緩沖信號，即在添加緩沖信號前便是延展信號。

圖17 LM1021添加緩沖信號與否對比Fig.17 Comparison of LM1021 before and after adding signature

圖18 C25D添加緩沖信號與否對比Fig.18 Comparison of C25D before and after adding signature

從圖17可以看出：在加入緩沖時間后LM1021遠場信號大有改善，激波上升時間和過壓幅值達到與大會提供結果一致的水平。從圖18可以看出：緩沖信號的加入與否并不影響激波上升時間和過壓幅值。

在模擬聲爆傳播過程中，需要人工插入激波，這一過程會造成波形在時間方向上延長，真實傳播過程中波形也會在時間方向上延長。如果最終長度超過原始波形并且沒有添加緩沖信號就會造成波形在時間軸左端“堆積”，導致結果不正確。

與C25D相比，LM1021尺寸更大，產生聲爆過壓值也更大，同樣無量綱推進距離下波形扭曲更嚴重，波形在時間軸上延長情況也更嚴重，因此LM1021需要加入緩沖信號而C25D則不需要。

4 時間網格與空間網格密度

由于程序中采用指定采樣點個數(shù)的方式確定時間網格密度，因此在探究時間網格密度時為了方便考慮，直接改變采樣點個數(shù)，再轉化為采樣頻率的方式。本文分析從小到大8個采樣點對遠場信號波形的影響，結果如圖19所示。

圖19 不同時間網格密度結果對比Fig.19 Comparison of farfield signatures of different time grid densities

從圖19可以看出：在5×104和8×104個采樣點的情況下，波形負值部分的形狀仍在變化，但正值部分尤其是與上升時間相關的波峰形狀基本固定；之后波形形狀與位置基本保持不變，只會發(fā)生小范圍浮動；過壓峰值隨著采樣點個數(shù)上升而上升，上升逐漸減緩直到2.5×105個采樣點后開始下降。

如果按照上述分析，在計算時應采用2.5×105個時間采樣點對輸入信號進行離散。但這并沒有實際意義，由于計算的時間代價難以接受。對于8×104和105個采樣點，過壓峰值增量只相當于峰值的不到0.04%。而且從8×104個采樣點后波形形狀基本固定，尤其是與響度密切相關的正值部分形狀在6×104個采樣點后便基本固定。因此在實際計算中應根據(jù)對精度的要求選取合適的采樣點個數(shù)。

對于不同的信號，尤其是不同時間長度的信號，在選取采樣點個數(shù)時應關注兩相鄰采樣點間的無量綱時間間隔和采樣頻率。由于采樣頻率是一個消除信號輸入頻率和時間影響的無量綱量，不同長度的信號輸入頻率不相同，但輸入信號長度與頻率相關，采用采樣頻率來描述更具有普適性。而無量綱時間間隔會直接參與計算，因此該響應與計算的方法與過程有較密切的關系。

由于在模擬過程中采用不同算法也會導致結果對頻率的響應不同，如果采用本文的程序，C25D信號在9×104個采樣點時過壓峰值已經開始下降，對應的采樣頻率約為430 kHz，而文獻［17］的計算結果則認為采樣頻率不應低于400 kHz。

空間網格的情況與時間網格類似，如圖20所示。

圖20 不同空間網格密度結果對比Fig.20 Comparison of farfield signatures of different spe?cial grid densities

從圖20可以看出：隨著網格逐漸加密，波形形狀先開始固定，之后開始小幅度擾動。因此空間網格并不是越密越好，只要滿足精度需求即可。與時間網格的選取相同，空間網格的密度受到計算方法的影響。如果計算方法在推進步長較大時發(fā)散便只能采用較小的推進步長，時間網格被迫加密。

5 結論

（1）本文提出了一種非線性效應求解算法，并基于算子分裂法構建了基于增廣Burgers方程的聲爆遠場預測程序，采用SBPW-2中的標模算例驗證了其有效性，其計算精度滿足設計要求。

（2）為了達到所需計算精度，時間網格和空間網格密度不需取到滿足網格收斂性的密度，能達到波形固定和過壓峰值增長不超過預定值即可。

（3）緩沖信號對于小型飛機聲爆計算來說不是必要的，但對于大型飛機來說是必要的。

（4）無量綱時間間隔會直接參與計算，具體算法會影響到波形對時間網格和空間網格的敏感程度。

（5）非線性算法受空間網格密度影響較大，為了提高計算速度應采用對空間網格密度不敏感的算法。