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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 簡化超聲速公務(wù)機(jī)模型

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

1.2 生成均勻大氣湍流隨機(jī)速度場

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

(1)

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

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

(2)

式中：x=[x,y,z]為湍流場中給定點(diǎn)的位置坐標(biāo)。

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

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

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

圖5 不同N值下能量譜曲線對(duì)比

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

1.3 聲爆遠(yuǎn)場傳播模型

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

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步的時(shí)間間隔；a0為當(dāng)?shù)芈曀?。不考慮大氣湍流效應(yīng)時(shí)，V0(i)=(V0x,V0y,0)僅代表所在高度平均風(fēng)速，一般只有兩個(gè)水平分量；考慮大氣湍流速度之后，V0(i)=(V0x+u′,V0y+v′,w′)，其中u′、v′、w′為大氣湍流速度分量。

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

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

圖8 三線性插值示意圖

2 結(jié)果分析

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

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

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

2.1 大氣湍流對(duì)地面聲爆過壓影響

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

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

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

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

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

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

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

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

2.2 大氣湍流對(duì)射線傳播路徑和地面到達(dá)點(diǎn)影響

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

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

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

圖15 有/無湍流條件下聲爆地面到達(dá)點(diǎn)位置分布

圖16 聲爆地面到達(dá)點(diǎn)分布

3 結(jié) 論

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

1) 均勻各向同性大氣湍流對(duì)地面過壓正峰值和負(fù)峰值的影響是正相關(guān)的，即使得正峰值增強(qiáng)的湍流條件也會(huì)使得負(fù)峰值增強(qiáng)，使得正峰值減弱的湍流條件也會(huì)使得負(fù)峰值減弱。

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

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

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

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