基于增廣Burgers方程的超聲速客機(jī)遠(yuǎn)場(chǎng)聲爆預(yù)測(cè)研究

崔青，白俊強(qiáng)，宋源，余培汛

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

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

0 引言

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

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

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

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

本文探究基于增廣Burgers方程的聲爆預(yù)測(cè)方法的求解過(guò)程，提出一種有效的非線性效應(yīng)求解方法；在算子分裂法的基礎(chǔ)上構(gòu)建聲爆遠(yuǎn)場(chǎng)預(yù)測(cè)方法，基于該方法編寫(xiě)聲爆預(yù)測(cè)程序，通過(guò)該程序計(jì)算大會(huì)的兩個(gè)標(biāo)準(zhǔn)算例［18-19］以驗(yàn)證該程序的正確性；分析空間、時(shí)間網(wǎng)格密度及緩沖信號(hào)對(duì)預(yù)測(cè)結(jié)果的影響。

1 方程求解

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

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

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

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

在廣義Burgers方程中添加分子弛豫效應(yīng)項(xiàng)，便得到增廣Burgers方式，形式為

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

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

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

式（5）第一項(xiàng)為非線性效應(yīng)項(xiàng)，用于計(jì)算由于傳播中非線性效應(yīng)累積導(dǎo)致的波形畸變，該效應(yīng)是聲爆波形在傳播中變化的主要原因。

式（5）第二項(xiàng)為熱黏性吸收效應(yīng)項(xiàng)，用于計(jì)算由于分子間摩擦導(dǎo)致的信號(hào)強(qiáng)度衰減。

式（5）第三項(xiàng)分子弛豫項(xiàng)與其他幾項(xiàng)形式上有一點(diǎn)差別，因?yàn)椴煌瑲怏w分子弛豫效應(yīng)不同，大氣整體弛豫效應(yīng)與不同種類氣體占比有關(guān)，因此該項(xiàng)需要將所有與大氣整體弛豫效應(yīng)有關(guān)的氣體進(jìn)行計(jì)算。真實(shí)大氣中對(duì)弛豫效應(yīng)影響較大的氣體有氮?dú)狻⒀鯕狻⒍趸己退魵狻Ｆ渲卸趸己退魵庥捎谠诖髿庵泻窟h(yuǎn)小于氮?dú)馀c氧氣含量，在弛豫效應(yīng)計(jì)算中可被略去，但需要注意大氣濕度對(duì)氧氣與氮?dú)獬谠バ?yīng)的影響沒(méi)有被忽略。因此第三項(xiàng)實(shí)際上是兩項(xiàng)求和的形式，出于簡(jiǎn)化離散求解的目的，求解弛豫效應(yīng)時(shí)依舊采取算子分裂法，即將氮?dú)馀c氧氣的弛豫效應(yīng)分開(kāi)計(jì)算，上一項(xiàng)的結(jié)果作為下一項(xiàng)的輸入，經(jīng)算例驗(yàn)證該算法能正確模擬大氣整體弛豫效應(yīng)。

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

式（5）最后一項(xiàng)為大氣分層效應(yīng)項(xiàng)，用于計(jì)算由于大氣參數(shù)如密度、溫度等隨高度變化導(dǎo)致的信號(hào)強(qiáng)度衰減。

其中除去非線性效應(yīng)外，文獻(xiàn)［9］中已有詳細(xì)的離散和求解方法，本文不再贅述。

1.1 非線性求解

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

為了便于求解，將波形轉(zhuǎn)換到扭曲時(shí)間坐標(biāo)下。

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

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

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

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

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

圖3 實(shí)際插入激波情況Fig.3 Insertion in reality

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

圖4 推進(jìn)步長(zhǎng)過(guò)大導(dǎo)致失效Fig.4 Long step causing failure

1.2 聲線追蹤及聲線管面積

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

式中：R為聲線位置矢量；W為當(dāng)?shù)仫L(fēng)速向量；N為波陣面單位法向量；I為單位矩陣；?為克羅克內(nèi)積。

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

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

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

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

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

1.3 計(jì)算流程

計(jì)算聲爆的程序框架如圖6所示，其中輸入分為三部分：過(guò)壓信號(hào)輸入，風(fēng)速輸入和大氣參數(shù)輸入。如果不輸入風(fēng)速和大氣參數(shù)，計(jì)算將采用內(nèi)置的無(wú)風(fēng)標(biāo)準(zhǔn)大氣［23］。輸入過(guò)壓信號(hào)必須是無(wú)量綱化的過(guò)壓信號(hào)。

圖6 程序計(jì)算流程圖Fig.6 framework of the program

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

采用指定有量綱推進(jìn)距離再預(yù)估—修正的方法確定每次推進(jìn)步長(zhǎng)。當(dāng)剩余高度不足指定距離時(shí)，將剩余高度無(wú)量綱化作為最后的推進(jìn)步長(zhǎng)代入程序計(jì)算。

2 標(biāo)模驗(yàn)證

在展示標(biāo)模驗(yàn)證結(jié)果之前，需要明確信號(hào)提取位置與指向的定義和風(fēng)速與機(jī)頭指向的定義，再展示標(biāo)模驗(yàn)證的結(jié)果。

2.1 提取位置與指向

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

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

飛機(jī)飛行的方向與正北方的夾角α如圖8所示，該角度用于將聲速表中沿正北、正東向的風(fēng)速轉(zhuǎn)化到以飛機(jī)飛行方向?yàn)閤正向的坐標(biāo)系下。更新聲線坐標(biāo)等計(jì)算均在該坐標(biāo)下進(jìn)行，該坐標(biāo)為右手坐標(biāo)系，當(dāng)飛機(jī)向正東方飛行時(shí)，x軸正向?yàn)闁|，y軸正向?yàn)楸保瑉軸正向向上。

2.2 LM1021標(biāo)模驗(yàn)證

LM1021是SPBW-2中遠(yuǎn)場(chǎng)計(jì)算的標(biāo)模算例之一，原型為洛克希德馬丁公司的一款超聲速客機(jī)，實(shí)驗(yàn)?zāi)Ｐ屯庑稳鐖D9所示。

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

第二屆國(guó)際聲爆研討會(huì)提供了風(fēng)洞實(shí)驗(yàn)數(shù)據(jù)［26］和遠(yuǎn)場(chǎng)數(shù)據(jù)［27］，提取位置H L=3.129 9，分別提取出射角為0°，30°的信號(hào)（如圖10所示）。

圖10 LM1021近場(chǎng)輸入信號(hào)Fig.10 Farfield signature of LM1021

對(duì)應(yīng)的巡航高度為16 764 m，巡航馬赫數(shù)為1.6Ma，地面反射因數(shù)取為1.9，側(cè)滑角為0°，大會(huì)提供了多個(gè)大氣情況和風(fēng)速表供選擇，本次選擇標(biāo)準(zhǔn)無(wú)風(fēng)大氣進(jìn)行驗(yàn)證，0°和30°遠(yuǎn)場(chǎng)波形對(duì)比如圖11～圖12所示。

圖11 0°遠(yuǎn)場(chǎng)波形對(duì)比Fig.11 Comparison of farfield waveform of 0°

圖12 30°遠(yuǎn)場(chǎng)波形對(duì)比Fig.12 Comparison of farfield waveform of 30°

圖11 ～圖12中，cal表示本程序計(jì)算結(jié)果，sboom表示大會(huì)程序提供的標(biāo)準(zhǔn)結(jié)果，可以看出：計(jì)算結(jié)果與大會(huì)提供標(biāo)模值吻合良好。

2.3 C25D標(biāo)模驗(yàn)證

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

圖13 C25D外形Fig.13 Shape of C25D

圖14 C25D近場(chǎng)輸入信號(hào)Fig.14 Extracted signature of C25D

圖15 0°遠(yuǎn)場(chǎng)波形對(duì)比（C25D）Fig.15 Comparison of farfield waveform of 0°（C25D）

圖16 30°遠(yuǎn)場(chǎng)波形對(duì)比（C25D）Fig.16 Comparison of farfield waveform of 30°（C25D）

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

3 緩沖信號(hào)的設(shè)置

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

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

圖17 LM1021添加緩沖信號(hào)與否對(duì)比Fig.17 Comparison of LM1021 before and after adding signature

圖18 C25D添加緩沖信號(hào)與否對(duì)比Fig.18 Comparison of C25D before and after adding signature

從圖17可以看出：在加入緩沖時(shí)間后LM1021遠(yuǎn)場(chǎng)信號(hào)大有改善，激波上升時(shí)間和過(guò)壓幅值達(dá)到與大會(huì)提供結(jié)果一致的水平。從圖18可以看出：緩沖信號(hào)的加入與否并不影響激波上升時(shí)間和過(guò)壓幅值。

在模擬聲爆傳播過(guò)程中，需要人工插入激波，這一過(guò)程會(huì)造成波形在時(shí)間方向上延長(zhǎng)，真實(shí)傳播過(guò)程中波形也會(huì)在時(shí)間方向上延長(zhǎng)。如果最終長(zhǎng)度超過(guò)原始波形并且沒(méi)有添加緩沖信號(hào)就會(huì)造成波形在時(shí)間軸左端“堆積”，導(dǎo)致結(jié)果不正確。

與C25D相比，LM1021尺寸更大，產(chǎn)生聲爆過(guò)壓值也更大，同樣無(wú)量綱推進(jìn)距離下波形扭曲更嚴(yán)重，波形在時(shí)間軸上延長(zhǎng)情況也更嚴(yán)重，因此LM1021需要加入緩沖信號(hào)而C25D則不需要。

4 時(shí)間網(wǎng)格與空間網(wǎng)格密度

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

圖19 不同時(shí)間網(wǎng)格密度結(jié)果對(duì)比Fig.19 Comparison of farfield signatures of different time grid densities

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

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

對(duì)于不同的信號(hào)，尤其是不同時(shí)間長(zhǎng)度的信號(hào)，在選取采樣點(diǎn)個(gè)數(shù)時(shí)應(yīng)關(guān)注兩相鄰采樣點(diǎn)間的無(wú)量綱時(shí)間間隔和采樣頻率。由于采樣頻率是一個(gè)消除信號(hào)輸入頻率和時(shí)間影響的無(wú)量綱量，不同長(zhǎng)度的信號(hào)輸入頻率不相同，但輸入信號(hào)長(zhǎng)度與頻率相關(guān)，采用采樣頻率來(lái)描述更具有普適性。而無(wú)量綱時(shí)間間隔會(huì)直接參與計(jì)算，因此該響應(yīng)與計(jì)算的方法與過(guò)程有較密切的關(guān)系。

由于在模擬過(guò)程中采用不同算法也會(huì)導(dǎo)致結(jié)果對(duì)頻率的響應(yīng)不同，如果采用本文的程序，C25D信號(hào)在9×104個(gè)采樣點(diǎn)時(shí)過(guò)壓峰值已經(jīng)開(kāi)始下降，對(duì)應(yīng)的采樣頻率約為430 kHz，而文獻(xiàn)［17］的計(jì)算結(jié)果則認(rèn)為采樣頻率不應(yīng)低于400 kHz。

空間網(wǎng)格的情況與時(shí)間網(wǎng)格類似，如圖20所示。

圖20 不同空間網(wǎng)格密度結(jié)果對(duì)比Fig.20 Comparison of farfield signatures of different spe?cial grid densities

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

5 結(jié)論

（1）本文提出了一種非線性效應(yīng)求解算法，并基于算子分裂法構(gòu)建了基于增廣Burgers方程的聲爆遠(yuǎn)場(chǎng)預(yù)測(cè)程序，采用SBPW-2中的標(biāo)模算例驗(yàn)證了其有效性，其計(jì)算精度滿足設(shè)計(jì)要求。

（2）為了達(dá)到所需計(jì)算精度，時(shí)間網(wǎng)格和空間網(wǎng)格密度不需取到滿足網(wǎng)格收斂性的密度，能達(dá)到波形固定和過(guò)壓峰值增長(zhǎng)不超過(guò)預(yù)定值即可。

（3）緩沖信號(hào)對(duì)于小型飛機(jī)聲爆計(jì)算來(lái)說(shuō)不是必要的，但對(duì)于大型飛機(jī)來(lái)說(shuō)是必要的。

（4）無(wú)量綱時(shí)間間隔會(huì)直接參與計(jì)算，具體算法會(huì)影響到波形對(duì)時(shí)間網(wǎng)格和空間網(wǎng)格的敏感程度。

（5）非線性算法受空間網(wǎng)格密度影響較大，為了提高計(jì)算速度應(yīng)采用對(duì)空間網(wǎng)格密度不敏感的算法。