冷巖，錢戰(zhàn)森,*，楊龍

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

自“協(xié)和號”飛機首飛以來，聲爆問題很大程度上決定了超聲速客機能否進入商業(yè)運營及取得商業(yè)成功，因而成為研制超聲速民用飛機必須解決的關鍵難題[1]。為了實現(xiàn)新一代環(huán)保型超聲速客機的成功運營，地面聲爆特征必須降低到航線附近居民可以接受范圍內。

超聲速飛機在實際飛行過程中，機體對大氣造成的壓力擾動(即聲爆)，沿著聲射線在真實大氣環(huán)境中傳播并最終到達地面，因而真實大氣環(huán)境對聲爆傳播有很大影響。然而，真實大氣環(huán)境千變萬化，包含溫度、濕度、季風、大氣湍流等多種干擾效應。其中溫度、濕度、季風等大氣不確定屬性雖隨著高度和季節(jié)變化，但其變化頻率相對較低，短時間內可以認為是確定性的影響，可根據(jù)氣象觀測數(shù)據(jù)來分析其在聲爆傳播過程中對傳播路徑和幅值的影響；而大氣湍流屬于相對高頻的變化因素，即便在短時間內也呈現(xiàn)為不確定性的影響，故長期以來針對大氣湍流對聲爆傳播路徑和幅值影響的數(shù)值模擬研究相對較少。但是根據(jù)NASA于1964年開展的飛行試驗所得數(shù)據(jù)可知[2]，大氣湍流對于聲爆的特征有著十分顯著的影響。在大氣湍流作用下，地面聲爆特征發(fā)生了顯著變化，并且呈現(xiàn)出一定隨機性，既有可能發(fā)生增強，也有可能發(fā)生衰減，在飛行架次有限的情況下，難以給出定量的判斷。因此，在聲爆數(shù)值模擬預測中考慮大氣湍流效應具有重要意義。

目前國外對于聲爆數(shù)值模擬預測模型中考慮大氣湍流影響研究方面尚未有成熟的方法。一種直接的思路為采用類似氣體動力學湍流研究的模擬方法，主要包括直接數(shù)值模擬(DNS)和大渦模擬(LES)方法，即針對Navier-Stokes方程直接開展空間數(shù)值模擬，該方法理論上是可行的，但是對于聲爆傳播這樣的超大尺度空間上的空氣動力學問題，其計算量巨大，遠遠超出了目前計算能力。另一種思路即為采用簡化的氣體動力模型，如Takeno等[3-4]發(fā)展的基于KZK方程的模型和Piacsek[5]、Locey[6]、Luquet[7]等發(fā)展的基于NPE方程的模型，這兩種方法在數(shù)學上是同源的，都根據(jù)一定的假設將全Navier-Stokes方程降階為關于擾動壓力的標量方程，進一步借助窄角近似，將方程簡化為具有主傳播方向的弱非線性聲學傳播方程，該類模型仍在發(fā)展之中，目前尚不成熟，且其計算量仍然相對較大。Yamashita和Obayashi[8]近期提出了在聲射線法中直接添加隨機湍流影響的思路，該方法操作相對簡便易行，可方便地添加到現(xiàn)有射線法程序中，且計算量也相對小很多，Yamashita和Obayashi雖然僅考慮了均勻各向同性湍流的影響，但是其研究成果具有重要參考意義。

國內在近場聲爆CFD計算與基于非線性Burgers方程的遠場聲爆預測方面取得了較大進展。西北工業(yè)大學[9-11]、中國空氣動力研究與發(fā)展中心[12]、航空工業(yè)空氣動力研究院[13-15]、中國航空研究院[16]、北京航空航天大學[17]等多家單位都開展了相關研究，發(fā)展了一系列聲爆的近、遠場預測方法，并針對AIAA聲爆預測大會發(fā)布的標模進行了計算研究。但是，國內在大氣湍流對聲爆特性影響研究方面，尚未見公開發(fā)表的文章。

鑒于此，受Yamashita和Obayashi[8]工作的啟發(fā)，本文將波形參數(shù)方法與各向同性隨機大氣湍流場相結合，建立了考慮大氣湍流效應的修正波形參數(shù)方法模型，通過對自主研發(fā)的基于Thomas波形參數(shù)法的ARI_Boom聲爆預測程序進行改進，開展了均勻各向同性大氣湍流對聲爆傳播特性的影響分析，重點分析了大氣湍流對地面聲爆過壓峰值、聲射線傳播路徑及地面到達點位置的影響。

圖1 數(shù)值模型及湍流影響區(qū)域示意圖

1 數(shù)值模型和方法

圖1為數(shù)值模型及湍流影響區(qū)域示意圖。整個模型包括以下3步：① 基于航空工業(yè)空氣動力研究院自主研發(fā)的ARI_Overset高精度數(shù)值模擬軟件[18-22]，采用CFD手段數(shù)值求解三維Navier-Stokes方程得到近場流場結構并提取模型下方設定位置處壓力分布;② 利用離散Fourier模態(tài)有限和生成表征均勻各向同性大氣湍流的隨機速度場[8,23-24]; ③ 以均勻各向同性大氣湍流隨機速度場和近場壓力分布為初始值，基于修正波形參數(shù)方法[25-26]傳播得到地面聲爆特征。

為了體現(xiàn)大氣湍流效應對地面聲爆特征的不確定性影響，本文利用上述方法共生成100個隨機湍流速度場，CFD計算完成后調用100次修正的射線法，得到可供統(tǒng)計分析的聲爆地面信號和傳播路徑。需要聲明的是，本文在研究過程中假設湍流場是凍結的，即在聲爆傳播過程中假定湍流分布不再發(fā)生變化。從時間尺度上分析，與自然界中大氣湍流結構演化時間相比，聲爆傳播所需時間為小量，因此這一假設是合理的[27-28]。下文將對圖1中各步驟所用方法進行詳細介紹。

1.1 近場CFD數(shù)值模擬

在飛行馬赫數(shù)Ma=1.6、飛行迎角α=2°、飛行高度H=14 km的飛行條件下，采用航空工業(yè)空氣動力研究院ARI_Overset軟件在三維空間求解Navier-Stokes方程得到簡化超聲速公務機模型的空間流場信息。該超聲速公務機模型如圖2所示，其特征長度L=37.6 m，半展長b/2=8.1 m。近場CFD模擬的重點是提取空間壓力分布，因此對于空間網(wǎng)格精度要求較高。根據(jù)前期研究經驗[19,22]，綜合考慮提高激波捕捉精度和計算效率的要求，采用圖3所示混合型計算網(wǎng)格，網(wǎng)格數(shù)量約為3 000萬。ARI_Overset數(shù)值模擬軟件的精準度前期已經過大量算例驗證[18-22]，為表述簡潔起見，這里就不再給出其參數(shù)設置和結果驗證細節(jié)。

圖2 簡化超聲速公務機模型

圖3 CFD模擬采用的混合網(wǎng)格(對稱面上)

1.2 生成均勻大氣湍流隨機速度場

基于Von Karman能量譜[29]，采用離散Fourier模態(tài)有限和形式生成均勻各向同性大氣湍流隨機速度場。Von Karman能量譜的定義為

(1)

式中：k=[k1,k2,k3]為波矢量；L0=0.032 m為長度尺度；η=0.005為Kolmogorov尺度。

能量譜確定后，利用N階離散Fourier模態(tài)有限和生成均勻各向同性大氣湍流隨機速度場，其表達式為

(2)

式中：x=[x,y,z]為湍流場中給定點的位置坐標。

圖4 波數(shù)向量示意圖[30-31]

能量譜為表達不同頻率湍流速度漲落對湍動能貢獻的函數(shù)。圖5給出了離散Fourier模態(tài)階數(shù)N取不同值時，能量譜的分布曲線。可以看到，隨著N值的增加能量譜曲線峰值附近趨于光滑，并且當N值超過200之后，能量譜曲線基本無變化，故本文后續(xù)工作中N的取值均為200。

大氣邊界層厚度隨氣象條件、地形、地面粗糙度而變化，高度范圍一般介于300～3 000 m之間[32]，考慮典型情況，本文模擬大氣邊界層的區(qū)域在x(飛行方向)、y(飛行側向)、z(高度方向)3個方向的尺度分別為8.0 km×0.4 km×3.0 km。每個方向上的網(wǎng)格點分布為等距形式，間距取為30 m，總網(wǎng)格點數(shù)大約為356 000。圖6給出了100次湍流隨機速度場的均方根速度值Vrms，可以看出上述方法生成的湍流隨機速度場的均方根速度值基本都在設置值Vrms=2.5 m/s附近波動，分布較為理想。

圖5 不同N值下能量譜曲線對比

圖6 100次隨機湍流速度場的均方根速度值

1.3 聲爆遠場傳播模型

聲爆遠場傳播模型采用基于幾何聲學和等熵波理論的修正波形參數(shù)方法[25-26]，可以得到波振幅、射線路徑和持續(xù)時間等信息。在聲爆遠場傳播過程中，射線路徑為

R(i+1)=R(i)+ΔR(i)

(3)

N(i+1)=N(i)+ΔN(i)

(4)

式中：R(i)為射線路徑；N(i)為波陣面單位法向量。

增量ΔR(i)和 ΔN(i)為

ΔR(i)=[a0(i)N(i)+V0(i)]Δt

(5)

(6)

(7)

式中：Δt為聲爆在空間傳播過程中第i步到第i+1步的時間間隔；a0為當?shù)芈曀佟２豢紤]大氣湍流效應時，V0(i)=(V0x,V0y,0)僅代表所在高度平均風速，一般只有兩個水平分量；考慮大氣湍流速度之后，V0(i)=(V0x+u′,V0y+v′,w′)，其中u′、v′、w′為大氣湍流速度分量。

圖7給出了湍流隨機速度場及聲射線路徑示意圖。從圖中可以看出，當聲信號(紅色實線)在隨機湍流速度場中傳播時，沿著聲射線每一時間步射線所在位置可能與生成的湍流隨機速度不在同一網(wǎng)格點，本文采用三線性插值方法將隨機湍流場速度插值到射線路徑點，如圖8所示。

圖7 湍流隨機速度場及聲射線路徑示意圖

圖8 三線性插值示意圖

2 結果分析

聲爆遠場傳播模型需要提取近場空間壓力分布作為初始值，因此近場空間壓力特征計算精度直接決定地面聲爆特征的預測精度。圖9給出了對稱面上近場壓力云圖及h/L=1.0,1.5,2.0處壓力分布曲線，其中h/L為空間壓力提取位置(h)與模型特征長度(L)比值，dp為過壓值，p∞為來流靜壓值。圖10給出了不考慮大氣湍流條件下，以不同h/L處壓力分布作為遠場傳播初始值計算所得地面聲爆特征。從圖中可以看出，不同h/L情況下，地面聲爆特征基本一致，進一步說明本文近場計算網(wǎng)格分辨率和空間壓力特征精度達到聲爆數(shù)值模擬要求。因此，本文選取h/L=1.5處空間壓力分布作為遠場傳播模型的初始值。

圖9 近場壓力云圖及不同h/L處壓力分布曲線

圖10 以不同h/L處壓力特征作為遠場傳播初值所得地面聲爆特征(1 psf=47.848 6 Pa)

2.1 大氣湍流對地面聲爆過壓影響

圖11 有/無湍流條件下地面聲爆特征

本節(jié)計算選取了100個湍流隨機速度場樣本。圖11給出了有/無湍流條件下地面聲爆特征對比結果。從圖中可以看出，雖然地面壓力特征均以h/L=1.5處壓力分布為遠場傳播模型初始值，但在湍流條件下地面聲爆特征發(fā)生了較大變化。圖11(a)表明，在無湍流條件下，地面聲爆過壓正峰值和負峰值分別為Δpmax(+)=86.93 Pa、Δpmax(-)=-70.75 Pa；在湍流條件下，地面聲爆過壓正峰值和負峰值均受到一定范圍的影響。

在本算例中，正峰值變化范圍為Δpmax(+)=51.67～260.73 Pa，負峰值變化范圍為Δpmax(-)=-40.10～-212.41 Pa。從圖中還可以看到，與無湍流狀態(tài)相比，均勻各向同性大氣湍流對地面過壓正峰值和負峰值的影響是正相關的，即對于給定的湍流速度場，如果均勻各項同性大氣湍流使得正峰值增強，也會使得負峰值增強，如果使得正峰值減弱，則會使得負峰值減弱。這一規(guī)律與大量飛行試驗觀測結果是一致的[2,33-35]。

圖12 地面聲爆過壓峰值概率柱狀圖

圖12給出了有/無湍流條件下地面聲爆過壓峰值變化的概率分布，圖12(a)和圖12(b)分別對應過壓正、負峰值。圖中橫坐標Δpmax-tur(+/-)對應100個湍流隨機速度場下的過壓峰值，Δpmax-no-tur(+/-)對應無湍流條件下的過壓峰值。從圖中看出，與無湍流的狀態(tài)相比，均勻各向同性大氣湍流使地面聲爆特征增強的概率約為55%，使地面聲爆特征減弱的概率約為45%，故總體上來看大氣湍流效應更傾向于增強地面聲爆特征。

在大氣湍流對地面聲爆過壓信號的影響方面，比值Δpmax-tur/Δpmax-no-tur在0.8～1.2范圍內的概率為73%。也就是說，在大部分情形下(73%的概率條件)，大氣湍流效應對地面聲爆特征的影響是相對緩和的，但是仍有27%的概率使得地面過壓峰值出現(xiàn)陡增或者陡降。特別是陡增條件下，比值可能達到3.1，這樣的放大倍數(shù)將使得聲爆的影響大大加強。

圖13給出了射線管面積隨高度變化曲線，紅色曲線為無湍流時的射線管面積隨高度的變化曲線?？傮w來看，射線管面積隨高度降低而逐漸增大，這表明聲爆信號總體上是隨著傳播距離增大而減小的，這主要是由于幾何擴散作用所致。圖13(a)和圖13(b)中的黑色曲線分別對應地面聲爆過壓峰值增加和減小時的射線管面積隨高度的變化，可以看到，在無湍流的情形下，射線管面積隨高度變化曲線基本是光滑的；在湍流條件下，射線管面積隨高度變化曲線會出現(xiàn)較大的波動，這主要是由湍流的衍射效應導致的；根據(jù)幾何聲學射線法基本原理，聲爆過壓值與射線管面積變化直接相關，當湍流效應使得射線管面積減小時地面聲爆過壓峰值呈增大趨勢(如圖13(a)所示)，當湍流效應使得射線管面積增加時地面聲爆過壓峰值呈減小趨勢(如圖13(b)所示)。

圖13 射線管面積隨高度變化曲線

2.2 大氣湍流對射線傳播路徑和地面到達點影響

圖14 不同條件下聲射線傳播路徑

在本文的計算中，僅考慮水平飛行條件，x軸為飛行方向(以x=0 km為飛機機頭位置)，y軸為側向(以y=0 km作為飛行軌跡正下方)，z軸為高度方向。圖14給出了聲射線的傳播路徑對比圖，其中黑實線表示無湍流條件下的射線路徑，紅虛線表示本文所采用的100個湍流隨機速度場中y>0 km方向偏移量最大時所對應的射線路徑，藍虛線表示采用的100個湍流隨機速度場中y<0 km方向偏移量最大時所對應的射線路徑。從圖中可以看出，這3條曲線的差異很小，即3條射線傳播路徑幾乎是相同的。這是因為本文所生成的100個湍流隨機速度場的湍流速度均方根值均在Vrms=2.5 m/s附近(如圖6所示)，該值遠遠小于真實大氣環(huán)境中當?shù)芈曀?約為a0=300 m/s)，故而聲爆射線傳播路徑仍由當?shù)芈曀偎鲗А?/p>

圖15給出了有/無湍流條件下聲射線的地面到達點位置，可以看出，在無湍流時聲爆地面到達點位于飛行軌跡正下方(y=0 km、x=13.01 km處)，這與預期結果一致；在湍流條件下，聲爆地面到達點隨機分布于無湍流條件下到達點的周圍，在本文計算中其地面覆蓋范圍為12.98 km≤x≤13.03 km,-7.11 m≤y≤6.84 m。如果將聲射線地面到達點按照過壓峰值增加或降低(與無湍流相比)兩種情況分類展示，則結果如圖16所示。圖16(a)給出了過壓峰值增加情形的聲射線地面到達點分布，圖16(b)給出了過壓峰值減小情形的聲射線地面到達點分布?？梢钥吹?，地面聲爆過壓峰值變化與聲射線地面到達點分布之間無明顯關聯(lián)，即與無湍流條件相比，地面聲爆過壓峰值的變化并不受聲射線地面到達點位置變化的影響。

圖15 有/無湍流條件下聲爆地面到達點位置分布

圖16 聲爆地面到達點分布

3 結 論

本文基于離散Fourier模態(tài)有限和生成的隨機大氣湍流場，采用修正波形參數(shù)方法，以簡化超聲速公務機模型為例，開展了均勻各向同性大氣湍流對地面聲爆特征影響的分析研究。應用航空工業(yè)空氣動力研究院自主研發(fā)的數(shù)值模擬軟件ARI_Overset在三維空間求解Navier-Stokes方程，得到作為聲爆遠場傳播初始值的近場壓力分布；基于Von Karman能量譜，生成了均勻各向同性大氣湍流場；基于修正波形參數(shù)方法，模擬了聲爆信號在隨機湍流速度場中的傳播過程。主要結論如下：

1) 均勻各向同性大氣湍流對地面過壓正峰值和負峰值的影響是正相關的，即使得正峰值增強的湍流條件也會使得負峰值增強，使得正峰值減弱的湍流條件也會使得負峰值減弱。

2) 即便是均勻各向同性大氣湍流，對于地面聲爆特征也有重要影響。與無湍流的狀態(tài)相比，均勻各向同性大氣湍流使得地面聲爆特征增強的概率約為55%，使得地面聲爆特征減弱的概率約為45%，故總體上來看大氣湍流效應更傾向于增強地面聲爆特征。

3) 大部分情形下(73%的概率條件)大氣湍流效應對地面聲爆特征的影響是相對緩和的(比值Δpmax-tur/Δpmax-no-tur在0.8～1.2范圍內)，但是仍有相當一部分情形(27%的概率條件)使得地面過壓峰值出現(xiàn)陡增或者陡降，特別是陡增時比值Δpmax-tur/Δpmax-no-tur可能達到3.1。

4) 在湍流條件下，射線管面積變化趨勢與地面聲爆過壓峰值有直接關聯(lián)，當聲射線臨近地面時，如果射線管面積減小，則地面聲爆過壓峰值增加，如果射線管面積增大，則地面聲爆過壓峰值減小。

5) 均勻各向同性大氣湍流對于聲爆傳播路徑(從飛行高度至地面)影響相對較小，但是仍會導致地面信號接收點(聲射線到達點)的位置發(fā)生一定范圍的不確定變化，且這種變化與聲爆地面過壓峰值增強或減小沒有明顯關聯(lián)。