低聲爆超聲速客機聲爆預測及不確定度量化分析

陳樹生，顧奕然，楊 華，黃江濤，高正紅，*

(1.西北工業大學 航空學院，西安 710072；2.中航工業 第一飛機設計研究院，西安 710089；3.中國空氣動力研究與發展中心，綿陽 621000)

0 引 言

隨著全球化的發展，利用超聲速客機可以極大縮短國際航線的飛行時間，積極發展新一代超聲速客機，有望極大促進我國同世界各國在經濟、文化等各領域的交流與合作，促進全球化進程不斷深入。

為了充分發揮超聲速客機的速度優勢，未來的超聲速客機必須具備陸地上空超聲速巡航能力。但超聲速客機投入運營仍存在種種阻礙因素，其中聲爆噪聲顯得尤為突出。因此低聲爆/低噪聲成為未來超聲速客機的重要設計指標之一。

20世紀50年代 ，聲爆現象首次得到關注，相關研究人員對其開展了初步研究。Whitham針對此現象，首次提出了基于線化理論的聲爆預測方法[1]。在1965～1970年間，NASA舉辦了三屆聲爆研討會[2-4]，借此推動深刻理解聲爆的產生、傳播、預測和最小化問題。

21世紀，聲爆預測和低聲爆超聲速客機的研究在國際上重新掀起熱潮。2014年，AIAA舉辦了第一屆聲爆預測研討會[5]，評估了當前近場聲爆信號的計算精度，為遠場聲爆預測打下基礎。2017年，AIAA舉辦了第二屆聲爆預測研討會[6]，該研討會同時將近場過壓值信號的計算和遠場地面波形的預測作為重點。2020年，AIAA舉辦了第三屆聲爆預測研討會[7-8]，該研討會主要將C608低聲爆超聲速飛行器作為標準算例，討論先進的近場CFD聲爆信號預測方法及遠場地面波形預測方法，及其在具有復雜推進和環境控制系統條件下的低聲爆外形計算中的適應性。

目前對聲爆的研究主要包括兩方面，聲爆預測以及低聲爆設計。其中聲爆預測是指對飛機超聲速飛行時產生的聲爆水平進行評估，包括對近場、中場以及遠場地面的聲爆水平進行計算[7-9]。而低聲爆設計是指通過一定的方式對飛行器超聲速飛行時的聲爆水平進行控制，使其處于能夠接受的范圍[10]。

聲爆預測方法大致分為三類：基于線性理論的聲爆快速預測方法[11]；基于近場CFD與遠場非線性傳播方程的聲爆高精度預測方法[12]以及全場CFD數值模擬方法[13]。其中，基于近場CFD與遠場非線性傳播方程的聲爆高精度預測方法是當前更為合適的聲爆預測方法，其優勢在于：一是相比于線性理論，該方法能夠計算復雜外形的聲爆特征，同時結合增廣Burgers方程還可以較準確地計算聲爆的上升時間；二是相對于全場CFD數值模擬，該方法大大降低了計算量，同時還合理地考慮了真實大氣環境中的分子弛豫、熱黏性吸收等效應。

近年來，聲爆的相關研究在國內也逐漸得到重視。朱自強和蘭世隆[14]指出，高精度流場模擬和基于非線性聲學的增廣Burgers方程等方法在聲爆預測中的優勢。張繹典等[15]建立了增廣Burgers方程的數值求解方法，并提出了一些方法來提高過壓值、上升時間等關鍵參數的計算精度。

超聲速客機在飛行過程中，由于飛行高度和大氣密度等經常發生波動，大氣來流條件也隨之不斷發生變化。為得到更可靠的低聲爆近/遠場信號預測結果，在考慮計算網格、物理模型及相對濕度等對聲爆近/遠場信號預測的影響的基礎上， 也需要考慮來流參數的不確定性和由此帶來的聲爆近/遠場信號預測的不確定度。目前，國內外對于考慮來流不確定性的低聲爆超聲速客機近/遠場信號預測的相關研究還較少。

針對CFD不確定度量化分析問題，多項式混沌方法[16]得到了廣泛應用。本文擬采用基于非嵌入式多項式的混沌方法（non-intrusive polynomial chaos, NIPC）開展不確定度量化分析研究，即選擇一組正交多項式作為空間的無限基函數，將變量分解為確定和隨機兩部分，從而構造出多項式混沌展開式。該方法將CFD求解器當做黑盒子，基于確定性的解來近似多項式混沌展開式的系數，并在非嵌入式多項式混沌方法的基礎上，采用Sobol指數來衡量輸入變量的敏感性。

本文擬以AIAA第三屆聲爆預測會標模C608低聲爆超聲速飛行器[8,17]為研究對象，開展復雜外形聲爆近/遠場信號特征分析及不確定度量化分析，以期為超聲速客機低聲爆設計提供參考依據，并進一步用于超聲速客機的魯棒優化設計及可靠性分析。

1 數值方法

1.1 聲爆近場預測方法

在聲爆計算體系下，采用基于Naiver-Stokes方程的CFD方法來計算近場過壓值分布，從而為聲爆傳播方程提供輸入信號以計算遠場地面信號。

在笛卡爾坐標系下，同時忽略體積力（如重力、電磁力等）和外部熱源，守恒形式的三維非定常可壓縮Naiver-Stokes方程可表示為[18]：

式中：Q為 守恒變量；f、g、h分別為三個坐標方向的通量，可表示為：

式中：ρ為密度；u、v、w為笛卡爾坐標系下的速度分量； τxx、τyy、τzz、τxy、τxz、τyz為黏性應力項；e為單位質量氣體的總能量；p為壓力；qx、qy、qz分別為三個坐標方向的熱流量。

1.2 聲爆遠場預測方法

在聲爆計算體系下，采用基于非線性聲學理論的增廣Burgers方程方法將近場波形推進到遠場，從而求解地面聲爆信號。

20世紀60年代，Blackstock[19]最早使用Burgers方程來模擬波在有損耗介質中的傳播，提出了經典Burgers方程[20]，其形式如下：

在經典Burgers方程中加入非均勻介質、幾何擴散以及分子弛豫效應的影響，即可得到增廣Burgers方程[21]：

式中：S為聲管面積[22]；(Δc)v為分子弛豫效應造成的聲速變化量； τv為弛豫時間，下標v表示不同的大氣組分（如 O2、N2）的弛豫過程。式（7）等號右邊 5 項分別對應非線性效應、經典吸收、大氣分層、幾何擴散以及分子弛豫效應對聲爆的影響。

為了便于求解，對式（7）進行無量綱化處理：

式中：P=p′/p0，p0為參考氣壓；無量綱距離 σ =ˉ，參考距離為無量綱時間 τ =ω0t′，參考時間為1 / ω0，可由輸入信號采樣率決定；無量綱耗散參數 Γ =bω0/(2βp0)；無量綱弛豫時間 θv= ω0τv；無量綱弛豫系數

1.3 混沌多項式不確定性分析

在考慮來流不確定性的低聲爆超聲速客機近/遠場信號不確定度量化分析中，以地面波形過壓值Δp等聲爆重要參數作為隨機輸出變量，表示為：

式中：α*是CFD方法或增廣Burgers方程計算結果；αj(x)是計算結果的確定部分耦合系數；隨機部分ψj(ξ)是以n維隨機變量 ξ =(ξ1，···，ξn)為自變量的隨機函數，為 ξ的正交多項式。式（9）可將任一隨機變量分解為確定和隨機兩個獨立部分。

出于計算量的考慮，對 ψj(ξ)多項式采取r階截斷，設隨機參數的維數為n，則混沌多項式（PCE）項數可以表示為[23]：

根據隨機變量 ξ的分布，即輸入參數的分布類型不同，對隨機變量 ψj(ξ)選取不同形式的正交多項式。當輸入參數服從正態分布時，選取Hermite正交多項式；當輸入參數服從均勻分布時，選取Legendre正交多項式[23]。

目前，主要有兩種形式的基于非嵌入式混沌多項式（NIPC）的不確定度分析方法：隨機響應面法（stochastic response surface method, SRSM）和基于Galerkin的投影法。本文采用隨機響應面法來求解正交多項式系數，開展不確定度分析。

通常采用過采樣策略。文獻[24]中推薦選用PCE系數兩倍的過采樣方法，即N=2Nt=2(R+1)。根據Schaefer等[25]的比較結果，采用精度和收斂性均表現較好的拉丁超立方（Latin hypercube design,LHD）抽樣方法，對模型參數進行抽樣選取。

定義

利用線性最小二次回歸，使J(α)最小，求解得到PCE系數αj的值，同時可以由展開式系數直接計算輸出響應的均值 μ和方差D等概率統計量，如下所示：

使用敏感度指標來表征不同隨機輸入對模型響應輸出的影響程度。本文采用可以直接從多項式混沌展開式中得到的Sobol敏感度指標，開展基于方差分解的全局敏感度分析。

將多項式混沌展開式重新分解，得到總方差為：

部分方差為：

則Sobol指數作為敏感性指數定義為部分方差與總方差的比值：

另外，針對輸入參數i的Sobol指數（STi）則定義為包含變量i的所有部分Sobol指數之和：

Sobol指數的詳細推導過程和具體形式可以參閱文獻[26]。

2 模型介紹及程序驗證

2.1 C608飛行器外形及來流條件

采用美國AIAA第三屆聲爆預測研討會（SBPW-3，2020）提供的C608飛行器標準算例開展研究[7]。C608是NASA低聲爆試驗飛行器X59的初步改形方案，其來自于洛克希德·馬丁公司對NASA N+2超聲速客機概念設計研究。三維視圖如圖1所示。

圖1 C608三維視圖Fig.1 Three-dimensional views of C608

C608飛行器參考長度L= 27.432 m，半模參考面積s= 37.16 m2，來流條件和邊界條件如表1所示。C608具有推進和環境控制系統（environmental control system, ECS）。

表1 C608來流和邊界條件Table 1 Freestream and boundary conditions of C608

2.2 網格收斂性研究

聲爆近場輸入信號計算采用自研CFD軟件[27-28]，湍流模型為SA模型，近場距離為3個C608飛行器長度（z= 82.29 m）。網格使用SBPW-3提供的非結構混合網格[7]，半模網格節點數約5 021萬，特征網格尺度h定義為：

當此特征網格尺度h= 0.58時，網格具有較好的收斂性。具體網格細節見表2。

表2 非結構混合網格Table 2 Unstructured hybrid grid

遠場地面波形計算利用自研聲爆程序[15]。選取上述CFD計算所得到的近場過壓值分布作為輸入信號，利用增廣Burgers方程計算遠場地面波形。由于增廣Burgers方程具有兩個維度（時間維度和空間維度），在進行網格收斂性研究時需要依次對這兩個維度展開分析。在均勻網格下，采用不同空間網格密度和時間網格密度，具體網格量見表3和表4。表中，NZ為空間網格點數，NT為時間網格點數。也可用時間采樣率表示時間網格密度。

表3 不同空間網格密度Table 3 Different spatial grid sizes

表4 不同時間網格密度Table 4 Different temporal grid sizes

圖2和圖3分別給了出遠場地面波形在空間和時間維度上的網格收斂性。由圖可見，地面波形對空間和時間網格密度都有一定敏感性，但需要加密到一定程度后，遠場地面波形信號才能逐漸收斂。對于該算例，當空間網格密度和時間網格密度分別達到20 000點和30 000點以后，網格達到收斂，此后地面波形基本不隨網格密度變化。

圖2 不同空間網格密度收斂性Fig.2 Spatial grid convergence

圖3 不同時間網格密度收斂性Fig.3 Temporal grid convergence

綜上所示，利用增廣Burgers方程計算遠場地面波形時，需適當加密網格。基于上述分析，本文選用網格NZ= 20 000和NT= 50 000開展研究。

3 基準狀態分析

3.1 基準狀態近/遠場信號分析

幾何模型為C608飛行器，三維外形見圖1，基準狀態飛行條件如表1所示。CFD計算采用SA湍流模型，空間格式為二階精度，時間推進方法為GMRES 方法[29]。圖4（a，b）給出了基準狀態下馬赫數和壓力云圖，圖4（c）為飛行器表面壓力分布圖。機頭產生激波和膨脹波對，在未受干擾前，細長的機身使流動平緩壓縮。機翼前緣、艙蓋后面的渦流發生器和ECS進氣道產生了一系列快速變化的激波和膨脹波。入口溢出在機翼上形成一個高壓區域。機翼后緣、水平尾翼、T型尾翼、發動機羽流等產生極其復雜的流動干擾，同時引發另一系列快速變化的激波和膨脹波。

圖4 飛行器基準狀態云圖Fig.4 Aircraft in the reference state

圖5和圖6分別對比了當前計算的近/遠場波形與SBPW-3公布數據[8]的差異，其中紅色和黑色線為SBPW-3數據，綠色實線為自研程序計算數據。可以發現，預測的近/遠場波形與SBPW-3數據基本重合，說明自研CFD程序和聲爆程序計算準確度符合要求，達到了較好的精度。整體近場過壓值分布變化頻率較高，存在多個波峰和波谷。由于經典吸收、分子弛豫效應等，使得大氣具有類似低通濾波器的性質，聲波能量及高頻脈動在大氣傳播過程中被消耗成內能，因而遠場地面波形較為平滑。

圖5 近場過壓值分布與SBPW-3數據對比Fig.5 Comparison of near-field over-pressure values with SBPW-3 data

圖6 遠場地面波形與SBPW-3數據對比Fig.6 Comparison of far-field ground waveform with SBPW-3 data

3.2 物理模型影響

本文選取層流模型、SA模型、SA-QCR模型[30]和SST模型來研究物理模型對近場過壓值分布及遠場地面波形的影響。

首先，給出層流模型下計算所得的馬赫數和壓力云圖（見圖7）。與SA湍流模型（見圖4）相比，層流模型馬赫數和壓力云圖在上下部分和飛行器前部與SA模型計算結果類似，在飛行器尾流部分兩者相差較大，其主要原因在于尾部區域噴流與主流之間有較強的剪切、湍流摻混，該區域用層流模型是不合適的。

圖7 層流模型下馬赫數和壓力分布Fig.7 Mach number and pressureobtained by a laminar model

接著，圖8和圖9對比了不同物理模型下近場過壓值分布和遠場地面波形的差異。從圖中可以看出：層流模型與湍流模型計算結果差別較明顯（尤其尾部區域），同時層流模型對地面波形過壓值預測偏低；SA湍流模型與SA-QCR湍流模型計算結果幾乎重合；SST湍流模型與SA湍流模型，前部波形幾乎一致，只有在尾部區域有較明顯的差異。綜合來看，不同湍流模型對近/遠場波形的影響相對較小，層湍模型對近/遠場波形的影響相對較大。

圖8 不同物理模型下近場過壓值分布Fig.8 Distributions of near-field over-pressure values obtained by different physical models

圖9 不同物理模型下遠場地面波形Fig.9 Far-field ground waveforms obtained bydifferent physical models

3.3 相對濕度影響

未來超聲速客機飛行將跨越海洋和大陸，兩者之間的相對濕度相差較大，而分子弛豫效應的強度與空氣濕度極為相關，因此本節研究不同濕度對地面波形的影響。

簡單起見，考慮5個不同的相對濕度（10%、30%、50%、70%、90%）進行計算，其中10%相對濕度近似模擬沙漠等干旱地區濕度，90%相對濕度近似模擬沿海等濕潤地區濕度。圖10給出了在不同相對濕度的大氣中的遠場地面波形。可見，相對于高濕度，低濕度對波形幅值的影響較大。相對濕度越低的大氣，其遠場地面波形的過壓值峰值更低，這是由于干燥的空氣對于聲波的耗散更強。文獻[31]表明，對于聲波在大氣中的傳播過程，因分子弛豫效應而產生的衰減，可以用一個衰減因子來衡量衰減的快慢，而聲爆波形頻率分量大多處于約化聲波頻率較低處，約化衰減因子會隨相對濕度的減小而增大，因此干燥的空氣對聲爆波形的耗散作用更強。此外需要說明的是，除本小節外，本文其余計算均假定相對濕度為70%。

圖10 不同相對濕度下遠場地面波形Fig.10 Far-field ground waveforms under different relative humidity

4 不確定度量化分析

超聲速客機在高空飛行過程中，來流條件會不斷發生波動。因此需要考慮來流參數存在不確定性時對聲爆預測結果的影響，并甄別來流參數中影響聲爆預測的關鍵性因素，為超聲速客機工程實際應用提供有價值的參考。本文選用C608低聲爆超聲速飛行器為研究對象，采用SA湍流模型和增廣Burgers方程開展聲爆近/遠場信號不確定度量化分析。

4.1 樣本點選取

超聲速客機巡航飛行環境一般位于大氣平流層，平流層中來流溫度基本不變。本文選取三個來流參數變量（來流馬赫數、來流攻角和單位雷諾數）作為隨機不確定性輸入參數，同時將地面波形過壓值Δp作為輸出響應，開展聲爆近/遠場波形的不確定傳播和參數敏感性分析。三個變量的取值范圍參考SBPW-3官網[7]和相關文獻[8，17]，并假設均服從均勻分布，不確定范圍設置如表5所示。其中：來流馬赫數基準值1.4，變化范圍為±3%，對應參數變化范圍為[1.358, 1.442]；來流攻角基準值0°，參數變化范圍為[-0.2°, 0.2°]；單位雷諾數 4 321 890 m-1，變化范圍±10%，對應參數變化范圍為[3 889 701 m-1, 4 754 079 m-1]。

表5 輸入變量參數變化范圍Table 5 Ranges of input variables

點配置非嵌入式多項式混沌方法（NIPC）采用二階截斷，為了更高精度地計算混沌多項式的系數，采用兩倍過采樣樣本。因此，本文利用拉丁超立方抽樣方法抽取20個樣本點。如圖11所示，本文選取的樣本點分布均勻。在每一個選取的樣本點上開展確定性的近場CFD計算和遠場增廣Burgers方程計算。后續將根據樣本點上的計算結果開展聲爆近/遠場信號不確定性量化分析。

圖11 拉丁超立方抽樣方法抽取樣本點Fig.11 Latin hypercube design method for sampling points

4.2 近場波形不確定度量化分析

為了便于開展不確定度量化分析，將所有樣本近場波形第一個波峰起點取為一致。圖12給出所有樣本點近場過壓值分布。從中可以發現，在近場過壓值的波峰、波谷等處存在較大的不確定度。

圖12 所有樣本點近場過壓值分布Fig.12 Near-field over-pressure values at all sampling points

在95%置信區間內，聲爆信號的不確定區間為(μ-1.96σ,μ+1.96σ)，區間大小為 2 ×1.96σ，其中 μ為平均值，σ為標準差。圖13為NIPC方法計算的近場過壓值均值與基準狀態近場過壓值。從圖中可以看到，基準狀態下的近場過壓值和NIPC近場過壓值均值基本一致，在波峰、波谷處有一定差異。

圖13 近場過壓值均值Fig.13 Mean value of near-field over-pressure

4.3 遠場波形不確定度量化分析

本文將所有樣本遠場地面波形起點取為一致，以便開展不確定度量化分析。圖14給出所有樣本點遠場地面波形分布。可以發現，整體地面波形在第一波峰、中間轉折及波谷處存在較大的不確定度。

圖14 所有樣本遠場地面波形分布Fig.14 Far-field ground waveform distributions of all samples

圖15 地面波形均值和標準差Fig.15 Mean and standard deviation of ground waveforms

圖16分別給出了不同來流參數的敏感性（Sobol）指數分布。通過前述三個來流參數不同時刻對應的敏感性指數分析，可以觀察到來流參數對聲爆傳播的影響程度。由圖16可知，來流單位雷諾數對地面波形過壓值的影響顯著低于來流馬赫數和來流攻角的影響，后兩者對地面波形不確定度貢獻較大。從流動機理上來看，聲爆近場及遠場過壓分布取決于脫體膨脹波和壓縮波，其中來流馬赫數和迎角對脫體空間波系結構特征影響較大，而雷諾數主要作用于物面附近，對脫體空間波系結構特征影響相對較小。

圖16 地面波形敏感性指數Fig.16 Sobol index of ground waveform

不確定度量化和敏感性分析表明，在低聲爆超聲速飛行器聲爆傳播預測數值模擬中，需要重點關注來流馬赫數和來流攻角的變化。不確定性輸入變量會引起遠場地面波形、波峰與波谷處的過壓值產生明顯變化。來流馬赫數變化（±3%）和來流攻角變化（[-0.2°, 0.2°]）對地面波形 Sobol指數貢獻最大的分別是t= 0.074 14和t= 0.085 84時刻，對應Sobol指數分別為0.984 4和0.986 8。

5 結 論

本文對AIAA第三屆聲爆預測研討會提供的C608低聲爆超聲速飛行器開展了聲爆傳播預測及不確定度量化分析，開展了遠場地面信號空間維度和時間維度上的網格收斂性研究，并對比基準狀態下不同湍流模型以及相對濕度對聲爆預測的影響。此外，選取來流馬赫數、來流攻角及來流單位雷諾數等三個不確定性輸入變量，采用非嵌入式多項式混沌方法對聲爆近/遠場信號進行了不確定度量化分析和敏感性分析，得到的結論如下：

1） C608飛行器標準算例聲爆近/遠場波形與SBPW-3公布數據相吻合，表明自研CFD軟件及聲爆程序計算準確度符合要求，具有較好預測精度。

2） 層流模型相較于湍流模型，在飛行器尾流部分預測較差，同時對地面波形過壓值預測偏低，不同湍流模型之間對聲爆近/遠場波形的影響較小。當相對濕度大于30%時，相對濕度對聲爆近/遠場波形的影響基本可以忽略不計。

3） 采用基于非嵌入式多項式混沌方法的不確定度量化分析方法，得到了近場波形均值、遠場波形均值和標準差等。其中近/遠場波形過壓值均值與基準狀態均值基本一致，僅在波峰和波谷處有一定差異。

4） 采用基于Sobol指數的敏感性方法，對地面波形過壓值 Δp進行了敏感性分析。其中來流單位雷諾數對地面波形過壓值的影響顯著低于來流馬赫數和來流攻角的影響。此外，聲爆地面感覺噪聲級是對不同頻率噪聲強度的量化，與遠場過壓分布之間呈非線性關系，今后應開展基于感覺噪聲級的不確定度分析。

5） 來流條件的波動影響聲爆近/遠場信號，其中氣流變化頻率或飛行器運動頻率對聲爆信號傳播產生復雜影響，如超聲速加速飛行發生聲爆聚集現象[32]，有待開展深入研究。