基于廣義Burgers 方程的聲爆傳播特性大氣湍流影響

王迪，冷巖，楊龍，韓忠華，錢戰(zhàn)森，*

1.西北工業(yè)大學 航空學院 翼型、葉柵空氣動力學國家級重點實驗室，西安 710072

2.中國航空工業(yè)空氣動力研究院 高超聲速氣動力/熱技術重點實驗室，沈陽 110034

3.中國航空工業(yè)空氣動力研究院 高速高雷諾數(shù)氣動力航空科技重點實驗室，沈陽 110034

聲爆作為發(fā)展超聲速民機不可回避的關鍵瓶頸問題，是未來新一代綠色超聲速民機能否取得商業(yè)成功的重要羈絆［1-5］。當飛機進入超聲速飛行時，其周圍流場產生的激波與膨脹波系經(jīng)過大氣環(huán)境中的各種宏、微觀效應的影響，到達地面后通常會演化成一道N 型波。然而，在地表附近的大氣邊界層內往往存在有大氣湍流等非靜止擾動，嚴重干擾聲爆特征的演化過程，從而對聲爆預測帶來了一定的不確定性。對于未來中國超聲速民機飛行軌跡規(guī)劃，精準的聲爆預測方法是其中的關鍵技術［6-8］。因此，開展大氣湍流效應對聲爆特性的影響規(guī)律研究具有重要意義。

大氣邊界層通常是指離地球表面約0.3～3 km高度的底層大氣［9］，由于受地球表面粗糙度和熱對流效應影響較為強烈，該區(qū)域的大氣運動具有明顯的湍流性質［10］，在短期內呈現(xiàn)出高度非定常性，導致聲爆信號的幅值發(fā)生隨機性增強或衰弱現(xiàn)象。圖 1 展示了NASA 曾開展的飛行試驗測量結果［11］，可以發(fā)現(xiàn)大氣湍流對聲爆特征的影響主要呈現(xiàn)出兩個特點［12-14］：一是激波后可能存在“釘狀細致結構”，導致聲爆的過壓峰值較大，最終演化成P 型波（Peaked Waveform）；二是激波峰值可能被抹平，導致聲爆上升時間的隨機性增加，最終演化成R 型波（Rounded Waveform）。

圖1 不同大氣條件下測量得到的地面聲爆波形［11］Fig.1 Sonic boom waveforms measured under different atmospheric conditions on ground［11］

目前的遠場聲爆預測技術主要是以線性化和擬線性化為基礎［15-17］發(fā)展而來的幾何聲學方法，該方法以Fermat 理論來確定聲射線的軌跡，建立了聲壓與聲射線管橫截面積之間的關系。1968 年，Randall［18］考慮分 層大氣 的影響，提出了 大氣物性沿高度變化對聲爆特性影響的處理方法。1972 年，Thomas［19］以射線法為框架發(fā)展的波形參數(shù)法成為最為經(jīng)典的快速聲爆預測方法，但是該方法對壓力波傳播過程中的幅值變化仍采用線性化方法處理，導致波形預測結果仍不理想。Robinson［20］、Cleveland［21］、Crow［22］考慮到真實大氣環(huán)境中的熱黏性吸收和分子弛豫效應的影響，建立了基于廣義Burgers 方程的遠場聲爆傳播模型，能夠更全面地涵蓋大氣的宏觀和微觀效應，得到了國內外的廣泛應用。近年來，國內在聲爆預測方面也積累了一些成果。王剛［23］和冷巖［24］等采用波形參數(shù)法進行了典型標模的數(shù)值預測和影響因素分析。張繹典［25］和喬建領［26］等使用廣義Burgers方程法開展了聲爆遠場計算和應用研究。王迪等［27］開展了廣義Burgers方程的高階格式離散方法研究。總體來說，靜止大氣條件下的聲爆預測方法發(fā)展的相對比較成熟，但上述方法并未考慮到大氣湍流效應的影響，無法解釋聲爆波形受大氣湍流作用產生的幅值和相位的隨機性變化。

由于聲爆波和大氣湍流相關作用的過程本質上是湍流與激波等復雜結構在近無限大空間范圍的多尺度相互作用問題。為了簡化求解過程，研究人員試圖從聲爆波和大氣湍流的物理特征出發(fā)，通過建立合理高效的數(shù)學模型來描述大氣湍流導致的散射和衍射等復雜效應。其中，最具代表的是由Zabolotskaya和Khokhlov［28］及Kuznetsov［29］提出的KZK 方程（Khokhlov-Zabolotskaya-Kuznetsov Eqution），Piacsek［30］、Locey［31］和Luquet［32］發(fā)展的NPE 方程（Nonlinear Progressive Equation）以及Dagrau［33］等推導 的HOWARD 方 程（Heterogeneous One-Way Approximation for the Resolution of Diffraction ）等形式的多維弱非線性方程，它們能夠較為準確地模擬聲爆信號受湍流效應作用后產生的形態(tài)變化，并與飛行試驗測量結果吻合較好。但是，基于多維弱非線性方程的預測方法仍面臨物理機理不甚明確，模型發(fā)展不夠成熟等問題，且對于聲爆從機體產生到傳播至地面的這類大尺度空間問題，求解多維弱非線性方程在計算量上依舊是難以接受的。1993 年，Robinson［34］提出了通過在射線法中添加隨機湍流速度場的方法模擬大氣湍流對聲爆傳播的影響。該方法思路簡單易行，可以快速預測聲爆的傳播過程，且能夠在一定程度上反映聲爆波形幅值受湍流效應影響的變化規(guī)律。2009 年，Yamashita 等［35］通過該方法進一步探究了均勻各向同性大氣湍流對聲爆傳播的不確定性影響，并與XB-70 飛機的飛行實測結果進行比較，獲得了良好的一致性規(guī)律。2020 年，冷巖等［36］受Yamashita 工作的啟發(fā)，將此方法應用于航空工業(yè)空氣動力研究院（以下簡稱航空工業(yè)氣動院）研發(fā)的ARI-Boom 聲爆預測程序中，詳細分析了大氣湍流對地面聲爆過壓峰值、聲射線傳播路徑及地面到達點位置的影響。然而上述研究工作主要是基于波形參數(shù)法開展的，該方法由于沒有考慮到大氣分子對聲能量的吸收和弛豫作用，導致遠場聲爆波形過壓峰值相對偏大且無法計算聲爆的上升時間，尤其是在大氣湍流的擾動下無法對地面聲爆波形做出準確判斷。相比之下，基于廣義Burgers 方程的遠場聲爆預測方法［20-22，25-27］由于考慮到真實大氣環(huán)境中的熱黏性吸收和分子弛豫等物理效應，能夠更加準確地反映大氣對聲爆傳播特性的影響，目前還未見到基于廣義Burgers 方程開展大氣湍流影響的研究工作。

本文在前期波形參數(shù)法的基礎上，開展了基于廣義Burgers 方程的均勻各向同性大氣湍流效應對聲爆傳播特性的影響研究。通過對中國航空工業(yè)空氣動力研究院發(fā)展的ARI-Boom 聲爆預測程序［24，27，36-37］進行改進，采用基于廣義Burgers 方程的遠場聲爆預測方法，結合射線法框架，在其中添加隨機湍流速度場，建立了一套可考慮熱黏性吸收、分子弛豫等物理效應的大氣湍流聲爆影響快速預測方法，并分析了湍流強度、大氣邊界層厚度等大氣湍流參數(shù)對典型遠程超聲速民機的聲爆特性影響規(guī)律。

1 計算模型與方法

考慮大氣湍流效應影響的聲爆快速預測方法的整體研究思路可以分為3 個步驟：

步驟 1通過求解三維Navier-Stokes 方程獲得飛機繞流流場的近場空間壓力分布。使用中國航空工業(yè)空氣動力研究院開發(fā)的ARI-Overset數(shù)值模擬軟件［37-39］，采用AUSMDV 通量差分格式［40］求解Navier-Stokes方程得到聲爆的近場壓力分布，作為遠場計算的輸入信號。根據(jù)經(jīng)驗［37，39］，為了準確且高效地捕捉激波，在近場區(qū)域的計算中采用混合型計算網(wǎng)格，詳細的參數(shù)設置及驗證信息這里不再給出。

步驟 2在大氣邊界層內，通過von Karman能量譜［41］，采用離散Fourier 模態(tài)有限和生成均勻各向同性脈動速度場，模擬大氣湍流對聲爆傳播過程中的速度擾動。von Karman能量譜可定義為

式中：k=(k1，k2，k3)為波矢量；L0=0.032 m 為長度尺度；η=0.005 為Kolmogorov 尺度。

能量譜確定后，采用N階離散Fourier 模態(tài)有限和生成的隨機速度擾動可表示為

式中：x=(x，y，z)為湍流場中給定點的位置坐標為振幅，kn為每一個波數(shù)向量kn對應的模量；Ψn為相位的隨機數(shù)；σn為單位法向量。由前期經(jīng)驗［36］發(fā)現(xiàn)，離散Fourier 模態(tài)階數(shù)N取200 時較為合適，其他變量的取值可參考文獻［35-36］。由于大氣邊界層的厚度通常介于0.3～3 km 的范圍內，模擬大氣邊界層的區(qū)域在x（飛行方向）、y（飛行側向）、（z高度方向）3 個方向的尺度分別為6.0 km×0.24 km×3.0 km，并確保湍流場網(wǎng)格的位置能夠覆蓋聲爆信號的傳播路徑。每個方向的網(wǎng)格點等距分布，間距為40 m，總網(wǎng)格量為67 500。

步驟 3以近場聲爆壓力特征作為初始條件，通過沿聲射線求解廣義Burgers 方程，模擬聲爆信號從飛機近場穿過包含湍流邊界層在內的真實大氣環(huán)境直至傳播到地面的整個過程。整體研究思路如圖 2 所示。

圖2 整體研究思路示意圖Fig.2 Schematic diagram of overall research

無量綱化后的廣義Burgers 方程可表示為［19］

式中：P為無量綱聲壓；σ和τ分別為無量綱聲射線長度和波形延遲時間；ρ0和c0分別為環(huán)境大氣密度和聲速；S為聲線管面積；Γ為熱黏性吸收系數(shù)；θυ和Cυ分別為無量綱分子弛豫時間和系數(shù)。

采用算子分裂法將式（3）分解成5 項獨立方程依次進行求解［27］。其中，前3 項方程在波形方向上分別使用五階WENO 格式、四階Pade 格式和Crank-Nicolson 格式，沿射線傳播方向使用四步四階Runge-Kutta 格式；后2 項方程具有解析解。聲爆的傳播路徑及射線管面積可由幾何聲學的射線法［42］技術獲得，其表達式為

式中：R(i)為射線路徑；N(i)為波陣面單位法向量。增量ΔR(i)和ΔN(i)可通過式（6）和式（7）計算得到：

式中：V0(i)表示聲爆信號所在高度的平均風速。在無湍流條件下，V0(i)=(V0x，V0y，0)僅需考慮兩個水平方向的風速分量；在大氣湍流的作用下，V0(i)=(V0x+u'，V0y+v'，w')，其 中u'、v'和w'表示大氣湍流的速度分量，即步驟2 中利用von Karman 能量譜生成的隨機速度擾動。

圖3 給出了聲爆信號沿聲射線穿過湍流速度場的示意圖，本文生成的隨機湍流場主要作用于聲射線傳播路徑和聲線管面積，繼而改變聲爆特征的演化過程。由于離散的聲射線位置與生成的湍流速度場不一定剛好在同一網(wǎng)格點上，需要通過三線性插值的方法將湍流速度分量插值到射線路徑點。另外，聲爆的傳播時間相對于湍流速度的響應時間尺度來說非常小，因此可以假設生成的湍流場是凍結的。

圖3 聲爆信號穿過湍流場示意圖Fig.3 Schematic diagram of sonic boom passing through turbulent field

2 計算結果與分析

本文使用的計算模型是航空工業(yè)氣動院提出的一類遠程超聲速民機［43］，其采用大長細比、大后掠角箭形翼和鴨式氣動布局方案，可有效降低超聲速飛行條件下的激波阻力和聲爆強度，幾何外形如圖 4所示。該模型的特征長度為67 m，飛行高度為15 km，巡航馬赫數(shù)為1.6，飛行攻角為4°，提取離機體機身長度1.5 倍距離處（h/L=1.5）的空間壓力分布曲線作為遠場的輸入波形，圖 5 給出了該算例的近場壓力云圖和提取的空間壓力分布曲線［43］。

圖4 一類典型遠程超聲速民機構型［43］Fig.4 Illustration of supersonic civil aircraft configuration［43］

圖5 近場壓力云圖及空間壓力分布曲線［43］Fig.5 Near-field pressure contour and pressure distribution［43］

圖6 給出了使用廣義Burgers 方程法與波形參數(shù)法在無湍流條件下計算得到的地面聲爆波形對比結果，取地面反射因子為1.9。可以發(fā)現(xiàn)，基于廣義Burgers 方程的遠場聲爆預測方法由于考慮到大氣分子對于聲能量的熱黏性吸收和分子弛豫作用，得到了更符合物理特性的聲爆信號，即波形中的激波不再是理想化的強間斷，而是有限厚度的強壓縮波，波形的最大過壓峰值也相對有所下降。另外，與波形參數(shù)法相比，廣義Burgers 方程模型能夠較準確地計算聲爆的上升時間，在本算例中約為1.80 ms。上升時間是計算聲爆感覺噪聲級的重要參數(shù)，在相同的過壓峰值下，上升時間越短，聲爆的感覺噪聲級越高，人們感受到的聲音越尖銳。

圖6 無湍流條件下廣義Burgers 方程法與波形參數(shù)法計算結果對比Fig.6 Comparison of calculation results of augmented Burgers equation and waveform parameter method without turbulence

圖7 為湍流條件下，分別采用廣義Burgers 方程法和波形參數(shù)法獲得的地面聲爆波形。每種方法分別選取100 組計算結果作為統(tǒng)計樣本，這里假設湍流邊界層厚度為1.5 km，生成的湍流隨機速度場的均方根速度值Vrms約為2.5 m/s。其中，黑色曲線表示不受湍流影響的地面聲爆特征；綠色和紅色曲線分別表示聲爆過壓峰值受湍流效應影響后出現(xiàn)增強或減弱的現(xiàn)象。根據(jù)統(tǒng)計結果表明，采用廣義Burgers方程法得到的地面聲爆過壓峰值范圍介于19.44～41.85 Pa 之間；而波形參數(shù)法的過壓峰值范圍介于12.89～228.03 Pa 之間。顯然，基于廣義Burgers方程的預測方法由于能夠表征熱黏性吸收、分子弛豫等大氣物理效應，其過壓峰值的波動范圍處于合理的區(qū)間內，與NASA 的飛行試驗測量結果規(guī)律相一致［11］。由于波形參數(shù)法對聲爆傳播過程中的幅值變化采用了線性化假設，通過該方法獲得的聲爆過壓峰值不夠準確，尤其是經(jīng)過大氣湍流的干擾后，其變化趨勢與實際的飛行試驗測量結果有所偏差，一般情況下均有明顯的過度放大。因此，本文建立的方法能夠更加真實地反映大氣湍流對聲爆傳播特性的影響。

圖7 湍流條件下廣義Burgers 方程法與波形參數(shù)法計算結果對比Fig.7 Comparison of calculation results of augmented Burgers equation and waveform parameter method with turbulence

另外，相比于前期的典型超聲速公務機［36］，本文采用的遠程超聲速民機聲爆波形更加復雜，該方法依然能給出復雜波系的大氣湍流影響規(guī)律。雖然廣義Burgers 方程仍是一維標量方程，但是其能夠在一定程度上更加合理地模擬大氣湍流對聲爆特征幅值產生的隨機性影響。同時，與KZK［28-29］、NPE［30-32］和HOWARD［33］等多維方程模型相比，一維廣義Burgers 方程具有計算量小的優(yōu)勢，能夠快速判斷大氣邊界層內的湍流效應對聲爆特性的影響規(guī)律，在未來超聲速民機飛行軌跡規(guī)劃的早期階段，可方便快速地預估地面聲爆強度受湍流影響的閾值。

2.1 湍流強度對聲爆特性的影響

為了探究大氣湍流強度對地面聲爆特性的影響規(guī)律，假設大氣邊界層厚度為1.5 km，通過設置隨機湍流速度場的均方根速度值Vrms分別為1.0、1.8、2.5 m/s，模擬低、中、高3 種湍流強度的影響。圖 8 為3 種湍流條件下的地面聲爆特征對比圖，每組條件分別統(tǒng)計了100 組計算樣本。可以發(fā)現(xiàn)，隨著湍流強度的提升，聲爆的過壓峰值范圍及最大過壓峰值均顯著增大，而最小過壓峰值有變得更低的趨勢，說明地面聲爆特征發(fā)生的隨機性變化會更加明顯，具體的統(tǒng)計結果如表 1所示。

圖8 不同湍流條件下的地面聲爆特征Fig.8 Ground sonic boom signatures for different turbulence intensities

表1 不同湍流強度下地面聲爆特征統(tǒng)計Table 1 Distribution statistics of ground sonic boom signatures for different turbulence intensities

圖9 給出了不同湍流強度下地面聲爆過壓峰值的概率分布柱狀圖。對于低湍流強度條件，大部分統(tǒng)計結果與無湍流條件下聲爆過壓峰值的比值（Δpmax-tur/Δpmax-no-tur）主要集中在0.9～1.1 之間，說明大氣湍流對聲爆傳播過程中的擾動較弱，因此地面聲爆特征的幅值變化幅度較小。當湍流強度持續(xù)增大后，聲爆的過壓峰值不再穩(wěn)定在較小的區(qū)間，根據(jù)統(tǒng)計，Δpmax-tur/Δpmax-no-tur介于0.9～1.1 之間的概率分別為93%、66%和38%，聲爆過壓峰值的變化范圍有進一步擴大的趨勢，可見高湍流強度會使聲爆強度增強的概率顯著提升。當湍流擾動速度達到2.5 m/s時，湍流條件下的最大聲爆過壓峰值達到無湍流條件下的1.54倍。

圖9 不同湍流強度下地面聲爆過壓峰值的概率分布柱狀圖Fig.9 Ground peak overpressure probability distribution histogram for different turbulence intensities

圖10 給出了所有條件下聲爆到達點的位置分布。可以觀察到，湍流條件下的聲爆到達點位置（黑色方塊）隨機分布在無湍流條件下聲爆到達點的周圍（紅色圓形，x=14.10 km，y=0 km處），且隨著湍流強度的增大，聲爆到達點圍成的范圍也呈現(xiàn)出擴大的趨勢。其中，涵蓋95%到達點圍成的圓的半徑分別為1.65 m、3.05 m 和3.55 m，說明高湍流強度不僅可能會增強對聲爆特性產生的隨機性變化，還有可能使聲爆到達點的位置更加分散，繼而影響聲爆毯的面積。表 2列出了更詳細的統(tǒng)計結果，聲爆到達點沿飛行方向（x方向）的波動范圍要明顯大于沿飛行側向（y方向）的波動范圍，該規(guī)律與波形參數(shù)法［35-36］的預測結果相一致。但總體上達到點分布區(qū)域仍在比較小的范圍，說明大氣湍流對聲爆到達點影響整體較小。

圖10 不同湍流條件下聲爆到達點位置分布Fig.10 Distribution of sonic boom arrival points on ground for different turbulence intensities

表2 不同湍流強度下聲爆到達點位置分布統(tǒng)計Table 2 Distribution statistics of sonic boom arrival points for different turbulence intensities

2.2 湍流邊界層厚度對聲爆特性的影響

氣象觀測結果表明，地球表面的湍流邊界層厚度通常介于0.3～3 km 的范圍，會伴隨地形和氣象條件有所變化。為了探究湍流邊界層厚度對聲爆特性影響的一般性規(guī)律，本節(jié)假設隨機湍流場的均方根速度值Vrms為2.5 m/s，設置邊界層厚度z從0.5 km 依次增加到3 km，每個狀態(tài)條件下同樣統(tǒng)計了100 組數(shù)據(jù)樣本，計算得到的地面聲爆特征如圖 11 所示。很明顯，聲爆特征的隨機性變化范圍隨著湍流邊界層厚度的增加有擴大的趨勢。另外，根據(jù)表 3 的統(tǒng)計結果顯示，聲爆特征增強或減弱的概率與湍流邊界層厚度沒有明顯的規(guī)律可循，而聲爆過壓峰值的標準差持續(xù)增大，說明聲爆幅值的波動程度將更加劇烈。

表3 不同邊界層厚度的地面聲爆特征統(tǒng)計Table 3 Distribution statistics of ground sonic boom signatures for different turbulence boundary layer heights

圖11 不同邊界層厚度的地面聲爆特征Fig.11 Ground sonic boom signatures for different turbulence boundary layer heights

圖12 給出了不同湍流邊界層厚度條件下，地面聲爆過壓峰值與無湍流條件的對比結果。對于4 種不同厚度的邊界層，Δpmax-tur/Δpmax-no-tur介于0.9～1.1 之間的概率分別為81%、61%、44% 和38%，該變化規(guī)律與圖 11 的觀察結果相一致，進一步說明隨著邊界層增厚，大氣湍流效應對聲爆特性產生的隨機性影響更加明顯，其最大過壓峰值將逐漸升高。當z=3 km 時，聲爆到達地面的最大過壓峰值達到了無湍流條件下的2.04 倍，這將對人們的正常生活產生嚴重的影響，因此在飛行軌跡規(guī)劃時應盡量避免湍流邊界層較厚的區(qū)域。

圖12 不同湍流邊界層厚度下地面聲爆過壓峰值的概率分布柱狀圖Fig.12 Ground peak overpressure probability distribution histogram for different turbulence boundary layer heights

圖13 給出了4 種湍流邊界層厚度對應的聲爆地面到達點的位置分布。根據(jù)表4，當邊界層厚度為0.5 km 時，有95%的聲射線落點都在以無湍流條件下聲爆到達點的位置為圓心、半徑為2.30 m 的圓內。當湍流邊界層厚度依次增加到1，2，3 km 時，其涵蓋95%聲爆到達點圍成的圓的半徑依次增加至3.61，5.09，6.70 m，說明隨著湍流邊界層厚度的增加，聲爆到達地面的位置同樣呈現(xiàn)出向周圍分散的趨勢，但達到點仍分散在比較小的范圍。

表4 不同湍流邊界層厚度下的聲爆到達點統(tǒng)計Table 4 Distribution statistics of sonic boom arrival points on ground for different turbulence boundary layer heights

圖13 不同湍流邊界層厚度下的聲爆到達點位置分布Fig.13 Distribution of sonic boom arrival points on ground for different turbulence boundary layer heights

3 結論

基于von Karman 能量譜，采用離散Fourier模態(tài)有限和生成均勻各向同性大區(qū)湍流隨機速度場，結合廣義Burgers 方程遠場聲爆預測方法，建立了可考慮熱黏性吸收、分子弛豫等物理效應的大氣湍流聲爆影響快速預測方法，并開展了典型遠程超聲速民機聲爆傳播特性的大氣湍流效應影響研究。得到的主要結論如下：

1）相比于前期基于波形參數(shù)法框架的聲爆預測方法，建立的方法能夠更合理地表征熱黏性吸收和分子弛豫等大氣物理效應，可以更真實地反映大氣湍流對聲爆傳播特性的影響規(guī)律，該方法的計算結果也更符合真實飛行試驗測量結果的變化趨勢。

2）使用的算例相比于前期的典型超聲速公務機生成的聲爆波形更加復雜，但該預測方法仍能給出復雜波系的大氣湍流影響規(guī)律；同時，該方法可以對地面聲爆信號進行快速預測，在未來超聲速民機的飛行軌跡規(guī)劃中具有獨到優(yōu)勢。

3）湍流強度和大氣邊界層厚度是影響聲爆特性的重要因素，隨著湍流強度和邊界層厚度的增加，聲爆過壓峰值的波動范圍有明顯增大的趨勢。在極端條件下，聲爆信號受大氣湍流影響后的最大過壓峰值可能達到無湍流條件下的2 倍以上。因此，在規(guī)劃飛行路線時，應該仔細評估當?shù)赝牧鲝姸群痛髿馔牧鬟吔鐚雍穸鹊挠绊憽?/p>

4）受大氣湍流效應的影響，聲爆到達地面的位置呈現(xiàn)出隨機性分布，且隨著湍流強度和邊界層高度的增加更加分散；同時，聲爆到達點沿飛行方向的波動范圍要明顯大于沿飛行側向的波動范圍，該規(guī)律與波形參數(shù)法的預測結果相一致；但總體上達到點分布區(qū)域仍在比較小的范圍，說明大氣湍流對聲爆到達點影響整體較小。