黃江濤，張繹典，高正紅，余婧，周鑄，余雷

1.中國空氣動力研究與發(fā)展中心，綿陽 621000 2.西北工業(yè)大學(xué) 航空學(xué)院，西安 710072

隨著氣動設(shè)計技術(shù)、新能源技術(shù)的發(fā)展和未來市場需要，在各國民航對超聲速聲爆問題嚴(yán)格限制的條件下，民航業(yè)界普遍認(rèn)為，發(fā)展小型超聲速公務(wù)機(jī)的技術(shù)條件以及市場時機(jī)已經(jīng)基本成熟。至少在未來的幾年內(nèi)，小型超聲速公務(wù)機(jī)的研制、試飛將會被提上日程，實際上美國、俄羅斯、法國以及日本等國家的航空公司均已經(jīng)推出一系列50座以下的超聲速公務(wù)機(jī)設(shè)計方案，并提出了三代超聲速民機(jī)的技術(shù)要求[1-2]，例如灣流Boom、Aerion、Spike等公司，并進(jìn)一步制造出了縮比原型機(jī)進(jìn)行拓展試驗[3-5]。

超聲速公務(wù)機(jī)面臨的最大挑戰(zhàn)之一就是民航對其超聲速飛行時聲爆水平的嚴(yán)格限制，聲爆強(qiáng)度水平的影響因素主要包含了質(zhì)量、飛行高度、飛行速度等。在總體方案選型以及布局優(yōu)化過程中，計算流體力學(xué)以及相應(yīng)的優(yōu)化設(shè)計手段起著至關(guān)重要的作用，大幅度降低了設(shè)計成本。先進(jìn)超聲速公務(wù)機(jī)在氣動性能上最明顯的特點是高巡航效率、低聲爆，需要在保證工程約束條件下，充分挖掘氣動外形的升阻比、聲爆設(shè)計潛力，是典型的多目標(biāo)精細(xì)化設(shè)計問題，這對氣動外形的綜合設(shè)計方法提出了苛刻要求，傳統(tǒng)的優(yōu)化將面臨計算量龐大、維度障礙等瓶頸問題。此時基于伴隨方程的梯度優(yōu)化是較為合理的選擇。

國外在氣動聲爆優(yōu)化方面的起步較早，主要研究工作包含了梯度、非梯度優(yōu)化，大部分研究工作基于伴隨方程的梯度優(yōu)化進(jìn)行，基于伴隨方程的優(yōu)化分為兩個方向：近場聲壓變分伴隨與流場/聲爆伴隨方程，最具有代表性是，Jameson等基于近場變分形式進(jìn)行氣動力/聲爆優(yōu)化[6]， Rallabhandi基于聲爆預(yù)測方程耦合變分進(jìn)行超聲速飛機(jī)聲爆優(yōu)化[7]。國內(nèi)在聲爆預(yù)測、優(yōu)化設(shè)計方面也開展了研究，取得了一定的進(jìn)展，但大多工作基于進(jìn)化算法以及波形參數(shù)方法等進(jìn)行[8-10]，基于伴隨系統(tǒng)的可微型聲爆信號優(yōu)化上的研究較少。基于伴隨方法的優(yōu)化設(shè)計盡管在全局性優(yōu)化問題上存在不足，但在高維設(shè)計變量精細(xì)化優(yōu)化問題上具有傳統(tǒng)方法不具備的天然優(yōu)勢，由于伴隨系統(tǒng)具有計算代價小、梯度計算量與各個學(xué)科設(shè)計變量個數(shù)均無關(guān)等優(yōu)點，因此，在氣動/聲爆綜合優(yōu)化領(lǐng)域具有不可替代的優(yōu)勢，是一個值得發(fā)展的研究方向。

在地面聲爆信號設(shè)計中，盡管近場變分實現(xiàn)方式比較簡單，卻無法直接設(shè)計地面聲爆信號的形態(tài)，不利于聲爆信號上升時間、過壓峰值等綜合特征的有效抑制。因此，本文基于中國空氣動力研究與發(fā)展中心自主研發(fā)的大型并行結(jié)構(gòu)化網(wǎng)格RANS(Reynolds-Averaged Navier-Stokes)求解器PMB3D、并行化伴隨方程求解器PADJ3D以及聲爆預(yù)測軟件，在地面過壓分布目標(biāo)函數(shù)變分條件下，構(gòu)建了耦合伴隨優(yōu)化系統(tǒng)，開展了求解過程中關(guān)鍵環(huán)節(jié)的變分與裝配方法研究，并將進(jìn)一步應(yīng)用于超聲速公務(wù)機(jī)聲爆優(yōu)化設(shè)計。

1 優(yōu)化系統(tǒng)中的學(xué)科分析模塊

1.1 并行化CFD求解器PMB3D

PMB3D是中國空氣動力研究與發(fā)展中心自主研發(fā)的大型并行CFD代碼，可求解任意曲線坐標(biāo)系下的Navier-Stokes方程實現(xiàn)流場精細(xì)化數(shù)值模擬：

(1)

PMB3D求解器具備S-A (Spalart-Allmaras)一方程、剪切應(yīng)力輸運(SST)兩方程湍流模型以及Langtry-Menter轉(zhuǎn)捩預(yù)測模型，具備JST (Jameson-Schmidt-Turkel)、Roe、Vanleer等空間離散方法，具備MUSCL(Monotonic Upwind centered Scheme for Conservation Laws)迎風(fēng)插值與WENO(Weighted Essentially Non-Oscillatory)格式，可實現(xiàn)LU-SGS隱式時間推進(jìn)，支持FAS(Full Approximation Scheme)形式多重網(wǎng)格技術(shù)、多塊對接、拼接以及重疊網(wǎng)格技術(shù)，以及具備基于MPI(Message-Passing-Interface)通信協(xié)議的大規(guī)模并行計算能力，廣泛應(yīng)用于常規(guī)氣動力計算、多體分離、空中加油安全性分析、進(jìn)排氣系統(tǒng)、旋轉(zhuǎn)部件氣動特性計算以及火箭發(fā)射、級間分離等領(lǐng)域[11]。

1.2 基于Burgers方程的地面聲爆預(yù)測

在進(jìn)行地面聲爆預(yù)測時，傳統(tǒng)的直接利用CFD計算的方法帶來網(wǎng)格耗散、網(wǎng)格需求量大、不適用于工程快速設(shè)計要求以及無法模擬大氣分層特性等問題，因此，選擇合理的聲爆預(yù)測方程是評估設(shè)計結(jié)果的重要環(huán)節(jié)。目前用于預(yù)測地面聲爆信號的方法主要包含波形參數(shù)法與Burgers方程，兩者在聲爆預(yù)測中具有良好的表現(xiàn)。但波形參數(shù)法[12-13]存在無法預(yù)測激波上升時間、預(yù)測信號存在間斷導(dǎo)致聲爆信號不可微等問題，無法進(jìn)行快速傅里葉變換(FFT)和感覺噪聲級分析，且在梯度優(yōu)化體系中應(yīng)用受限。因此，文中基于Burgers方程進(jìn)行聲爆預(yù)測[14-15]：

(2)

在時間步長足夠小的前提下，可以通過算子分裂法對式(3)右端各項所代表的非線性項、吸收影響、分子松弛、射線管擴(kuò)張以及大氣分層項進(jìn)行求解：

(3)

(4)

式(2)～式(4)中各參數(shù)的具體含義參見文獻(xiàn)[14]。

圖1～圖3為課題組優(yōu)化軟件的CFD模塊、Burgers方程聲爆預(yù)測模塊[16]對雙圓錐、1021標(biāo)模聲爆[17-18]的預(yù)測與風(fēng)洞試驗以及NASA計算結(jié)果的對比，圖中：δP、t、H/L分別為過壓、時間和聲爆預(yù)測位置坐標(biāo)?？梢钥闯觯糜趦?yōu)化的聲爆預(yù)測模塊精度較高，為聲爆伴隨方程及優(yōu)化提供了基礎(chǔ)平臺。

圖1 雙圓錐壓力云圖Fig.1 Pressure contour of cone

圖2 聲壓計算與試驗數(shù)據(jù)對比Fig.2 Comparison of sound pressure calculation and experiment data

圖3 NASA1021標(biāo)模聲爆預(yù)測對比Fig.3 Sonic boom prediction comparison of NASA1021 standard model

2 基于流場/聲爆耦合伴隨思想的梯度計算

2.1 流場伴隨方程求解梯度

簡單回顧一下伴隨方程的推導(dǎo)，對于目標(biāo)函數(shù)最小化問題：

(5)

在流場收斂解的條件下，殘差約束R(W,X,D)=0,W、X、D分別代表流場守恒變量、網(wǎng)格變量、設(shè)計變量。引入拉格朗日算子Λ可以構(gòu)造以下目標(biāo)函數(shù)：

L=I+ΛTR

(6)

式中：L為聲爆目標(biāo)函數(shù)。

對式(6)進(jìn)行求導(dǎo)，有

(7)

(8)

式(8)就是流場伴隨方程，通過隱式迭代方法求解Λ之后，可以通過式(9)和式(10)進(jìn)行目標(biāo)函數(shù)梯度信息求解。

(9)

(10)

最小化問題中的目標(biāo)函數(shù)既可以是氣動力，也可以是聲壓、總壓恢復(fù)系數(shù)、流量、壓比等參數(shù)，由以上推導(dǎo)過程可以看出，對于伴隨方程來講，不同的目標(biāo)函數(shù)只需要改變伴隨方程右端的變分項。其中目標(biāo)函數(shù)I可以是升力、阻力、力矩、流量、壓比等參數(shù)，R為殘差約束，Λ為拉格朗日算子。

本文采用二階精度的中心格式、人工黏性以及SST兩方程湍流模型，采用手工推導(dǎo)方式進(jìn)行黏性離散伴隨方程構(gòu)造，最終表達(dá)式為

Rc(λ)j,k,l-RD(λ)j,k,l-Rv(λ)j,k,l=0

(11)

式中：Vj,k,l、Rc(λ)j,k,l、RD(λ)j,k,l、Rv(λ)j,k,l分別代表網(wǎng)格體積、伴隨方程的對流項、人工黏性項以及物理黏性項，對式(11)的迭代求解，文中采用了LU-SGS方法隱式時間推進(jìn)方法，其中，離散伴隨方程的邊界條件采用矩陣形式處理方式，黏性項采用薄層近似，離散伴隨方程求解時，并行機(jī)制依然采用單元數(shù)衡量的負(fù)載平衡、對等式計算以及MPI消息傳遞模式。本文求解器采用了多塊對接網(wǎng)格技術(shù)，MPI傳遞的信息是各個進(jìn)程中分割面上的兩層虛網(wǎng)格上的伴隨變量信息，詳細(xì)推導(dǎo)可以參考文獻(xiàn)[19]。

2.2 聲爆伴隨方程求解梯度

對于聲爆設(shè)計，可以對地面聲爆強(qiáng)度進(jìn)行變分，也可以對近場聲壓進(jìn)行變分。近場變分[6]實現(xiàn)方式更為簡單，但無法直接設(shè)計地面聲爆信號特征，因此，本文采用了地面變分形式，該方法需要推導(dǎo)聲爆預(yù)測伴隨方程，求解地面聲爆目標(biāo)函數(shù)對近場聲壓的梯度。聲爆伴隨方程的具體形式如下，詳細(xì)推導(dǎo)可以參考文獻(xiàn)[7]：

(12)

式中：λ、β、γ為中間伴隨變量；Ib為地面聲爆目標(biāo)函數(shù)；A、B與A2、B2分別對應(yīng)氮氣、氧氣分子弛豫矩陣；A3、B3為吸收過程矩陣。與聲爆預(yù)測方程不同，式(12)的求解過程是聲傳播的一個反向過程，利用最終的伴隨變量可以很方便地獲取地面聲爆目標(biāo)函數(shù)對近場聲壓的梯度：

(13)

式中：pin為均勻坐標(biāo)系下的過壓分布。需要注意的是，式(13)求解的是地面聲爆目標(biāo)函數(shù)對均勻坐標(biāo)系下的近場輸入聲壓的梯度，向網(wǎng)格單元裝配需要將該梯度轉(zhuǎn)化為CFD網(wǎng)格非均勻坐標(biāo)系下，依據(jù)網(wǎng)格非均勻坐標(biāo)系與聲爆均勻坐標(biāo)系的轉(zhuǎn)換關(guān)系，可以方便推導(dǎo)出對角稀疏化的坐標(biāo)轉(zhuǎn)換雅克比矩陣χ[7]，依據(jù)分段線性插值表達(dá)式可以實現(xiàn)均勻坐標(biāo)系與非均勻坐標(biāo)系的導(dǎo)數(shù)轉(zhuǎn)換：

(14)

圖4和圖5給出了基于聲爆伴隨方程中間伴隨變量的分布，以及地面聲爆目標(biāo)函數(shù)對近場非均勻坐標(biāo)系下聲壓的梯度驗證，地面聲爆目標(biāo)函

圖4 不同高度(H)聲爆伴隨變量Fig.4 Adjoint variables of sonic boom at different altitudes (H)

圖5 聲爆伴隨梯度與差分對比Fig.5 Comparison of gradients of sonic boom adjoint and finite difference method

數(shù)采用以下形式：

式中：pT為聲爆設(shè)計目標(biāo)特征。可以看出聲爆伴隨方程梯度計算結(jié)果與差分結(jié)果較為一致，可以為耦合伴隨系統(tǒng)提供準(zhǔn)確的地面聲爆目標(biāo)函數(shù)對近場聲壓的梯度。

2.3 流場/聲爆耦合伴隨方程

Rallabhandi在文獻(xiàn)[7]中采用非結(jié)構(gòu)網(wǎng)格求解器FUN3D進(jìn)行聲爆優(yōu)化，由于非結(jié)構(gòu)的不規(guī)則性以及采用了自適應(yīng)網(wǎng)格產(chǎn)生的拉伸影響，需要將網(wǎng)格格點單元的流場變量向同高度坐標(biāo)變換，因此，近場聲壓函數(shù)關(guān)系式是關(guān)于流場變量與網(wǎng)格坐標(biāo)的函數(shù)，雅克比矩陣χ就包含了對網(wǎng)格坐標(biāo)X的變分。

為簡化耦合伴隨系統(tǒng)的變分推導(dǎo)過程，降低變分難度，文中依據(jù)結(jié)構(gòu)網(wǎng)格拓?fù)涞目煽匦?，進(jìn)行以下操作規(guī)定：① 在近場過壓提取站位附近將網(wǎng)格單元分布劃分為規(guī)整格式，即高度、寬度方向均為直線，這樣近場過壓分布就不需要向同高度轉(zhuǎn)換；② 非均勻坐標(biāo)下沿X方向各個站位初始過壓的提取均從本單元選取。由上述規(guī)則，近場過壓的提取基本消除對X的依賴，雅克比矩陣χ不再包含對網(wǎng)格坐標(biāo)X的變分，且僅與自身單元守恒變量W相關(guān),即

p0=T(W)

(15)

式(15)大幅度簡化了近場聲壓雅克比轉(zhuǎn)換矩陣的變分難度。與文獻(xiàn)[7]不同，本文進(jìn)行變分的約束沒有網(wǎng)格伴隨方程的殘差項，而只有流場殘差R=0與對聲壓轉(zhuǎn)換關(guān)系(p0-T)=0，基于上述原則，下面給出耦合伴隨的推導(dǎo)過程，將聲爆目標(biāo)函數(shù)引入流場以及聲爆拉格朗日算子λf、λb：

(16)

對式(16)進(jìn)行變分展開：

(17)

綜合文中網(wǎng)格劃分以及近場聲壓提取原則，可以看出式(17)右端第1、第6項為零，變分表達(dá)式為

(18)

(19)

(20)

3 并行環(huán)境下近場聲壓提取與變分結(jié)果裝配

耦合伴隨方程的第1步是為聲爆傳播方程提供近場聲壓輸入，在該問題上面臨的主要關(guān)鍵技術(shù)是并行環(huán)境下提取多塊網(wǎng)格運算中的進(jìn)程號、網(wǎng)格塊編號以及單元編號。

為方便近場聲壓的提取，定義長方體“盒子”，該盒子由長方體兩個角點定義，用于方便選定近場網(wǎng)格單元，避免繁瑣的人工操作，如圖6所示。

圖6 并行環(huán)境下的網(wǎng)格分布與“聲壓盒”Fig.6 Grid distribution and “sound pressure box” in parallel environment

進(jìn)程號、網(wǎng)格塊編號和單元編號提取、聲爆預(yù)測以及變分結(jié)果裝配過程如下：

1) 各個進(jìn)程的網(wǎng)格判斷是否有格心坐標(biāo)處于盒子內(nèi)，若有，記錄該進(jìn)程的編號。

2) 記錄該網(wǎng)格塊在當(dāng)前進(jìn)程中的編號。

3) 記錄格心在當(dāng)前網(wǎng)格中的編號及X坐標(biāo)。

4) 主進(jìn)程將各個進(jìn)程的編號、坐標(biāo)文件收集寫出，并將過壓按X順序進(jìn)行排列輸出近場文件。

5) 聲爆預(yù)測迭代推進(jìn)計算。

6) 聲爆伴隨方程反向迭代推進(jìn)求解。

7) 轉(zhuǎn)換為CFD坐標(biāo)系。

8) 按編號、坐標(biāo)文件將變分結(jié)果按對應(yīng)的進(jìn)程編號輸出，向各個進(jìn)程裝配。

9) 流場伴隨方程求解。

4 聲爆優(yōu)化設(shè)計體系

本文的研究工作基于課題組自主研發(fā)的優(yōu)化設(shè)計軟件AMDEsign進(jìn)行，AMDEsign是面向航空航天飛行器氣動外形多學(xué)科優(yōu)化研發(fā)的分布式優(yōu)化軟件，集成了氣動、氣動/結(jié)構(gòu)、氣動/隱身多目標(biāo)多學(xué)科設(shè)計體系，具有進(jìn)化算法/代理模型、耦合伴隨等模塊，適用于不同設(shè)計問題的各種優(yōu)化模型、參數(shù)化方法、網(wǎng)格變形技術(shù)以及學(xué)科分析。

本文發(fā)展的低聲爆超聲速公務(wù)機(jī)優(yōu)化技術(shù)屬于AMDEsign中的耦合伴隨模塊，其中CFD方法采用中心格式、SST兩方程湍流模型、多重網(wǎng)格加速收斂以及大規(guī)模并行計算；梯度求解采用對應(yīng)的并行化伴隨方程與梯度求解器；參數(shù)化方法采用基于NURBS基函數(shù)的自由變形技術(shù)(FFD)[20]；變形網(wǎng)格采用并行化RBF_TFI[21]，聲爆優(yōu)化基本流程如圖7所示。

圖7 流場/聲爆耦合伴隨優(yōu)化流程Fig.7 Flowchart of flow field/sonic boom coupled adjoint optimization

5 典型超聲速公務(wù)機(jī)聲爆優(yōu)化

基于巡航馬赫數(shù)為1.5的小型超聲速公務(wù)機(jī)的基本構(gòu)型，開展超聲速巡航狀態(tài)下的聲爆優(yōu)化設(shè)計研究。

其優(yōu)化數(shù)學(xué)模型為

Constraints:

主要部件包含機(jī)翼、機(jī)身以及立尾。網(wǎng)格劃分為239塊，半模網(wǎng)格規(guī)模為900萬量級，為比較準(zhǔn)確地捕捉空間激波形態(tài)，空間網(wǎng)格拓?fù)浒凑振R赫角進(jìn)行x方向拉伸且對下表面進(jìn)行加密，如圖8 所示。采用64核進(jìn)行并行計算。

圖9給出了基于NURBS基函數(shù)的自由式變形參數(shù)化示意圖，共采用60個控制頂點實現(xiàn)機(jī)身、機(jī)翼氣動外形參數(shù)化建模；設(shè)計變量為圖9中x方向4～6個截面，以及機(jī)翼展向4個截面，近場“包圍盒”處于機(jī)身下方60 m處。圖10為初始外形的對稱面壓力云圖，可以看出按照馬赫角拉伸的網(wǎng)格拓?fù)淠軌蜉^為清晰地捕捉空間激波形態(tài)。

圖8 CFD網(wǎng)格分布Fig.8 Distribution of CFD grid

圖9 自由式變形參數(shù)化Fig.9 Parametric lattice of free form deformation

圖10 初始外形巡航狀態(tài)壓力云圖Fig.10 Pressure contour of initial configuration in cruise state

圖11給出了耦合伴隨系統(tǒng)的收斂歷程，圖12和圖13分別給出了y=0與x=8站位截面耦合伴隨方程的第1伴隨變量云圖，云圖呈反向傳播，實際上，結(jié)合最終的導(dǎo)數(shù)計算式(9)，可以從伴隨變量云圖分布上定性看出流動敏感性區(qū)域。同樣，結(jié)合耦合伴隨變量云圖分布以及伴隨方程耦合變分項的裝配位置可以看出，隨著偽時間推進(jìn)迭代，由空間近場向物面“傳播”，因此，耦合伴隨方程呈現(xiàn)先收斂，進(jìn)而階段振蕩，最終完全收斂的特征。

圖14給出了地面聲爆目標(biāo)函數(shù)優(yōu)化過程，經(jīng)過9代優(yōu)化，聲爆抑制效果逐漸收斂，圖15給出了地面聲爆目標(biāo)函數(shù)對近場聲壓的梯度優(yōu)化前后對比，圖16給出了優(yōu)化前后地面過壓對比?？梢钥闯?，伴隨優(yōu)化效果較為明顯，第2道激波的峰值明顯降低。由于機(jī)頭的控制點沒有作為設(shè)計變量，第1道激波壓力峰值抑制以及激波上升時間控制效果不太明顯。

圖11 耦合伴隨系統(tǒng)收斂歷程Fig.11 Convergence history of coupled adjoint system

圖12 耦合伴隨方程第1伴隨變量云圖(y=0)Fig.12 First adjoint variable contour of coupled adjoint equations (y=0)

圖13 耦合伴隨方程第1伴隨變量云圖(x=8)Fig.13 First adjoint variable contour of coupled adjoint equations (x=8)

圖14 目標(biāo)函數(shù)優(yōu)化收斂歷程Fig.14 Convergence history of objective function optimization

圖15 地面信號對近場聲壓梯度的優(yōu)化前后對比Fig.15 Comparison of gradient of ground signal with near field sound pressure between initial and optimized configuration

圖17給出了優(yōu)化前后感覺噪聲級的頻譜特性，橫坐標(biāo)為1/3倍頻程頻段，縱坐標(biāo)為響度級。由于頻譜特性峰值與激波上升時間密切相關(guān)，而第1道激波上升時間與機(jī)頭形狀關(guān)系密切，如前面所述，文中沒有將機(jī)頭作為設(shè)計變量，圖16第1 道激波上升時間變化不大，因此，峰值變化不大，但在整個頻段內(nèi)，感覺噪聲級幅值得到有效抑制，從而起到降爆作用，驗證了文中耦合伴隨優(yōu)化方法的有效性。

圖16 地面過壓優(yōu)化過程Fig.16 Optimization process of ground pressure

圖17 優(yōu)化前后感覺噪聲級頻譜特性Fig.17 Characteristics of spectrum on perceived noise level of initial and optimized configuration

6 結(jié) 論

基于并行化結(jié)構(gòu)網(wǎng)格RANS求解器PMB3D以及伴隨方程求解器PADJ3D，結(jié)合增廣Burgers聲爆預(yù)測與伴隨方程，開展了流場/聲爆耦合伴隨方程構(gòu)造、求解方法以及聲爆優(yōu)化研究。

1) 文中提出的“包圍盒”、單元定位以及排序方法，大幅度提高了并行環(huán)境下CFD數(shù)值模擬結(jié)果的近場過壓分布提取效率，避免了繁瑣的人工操作。提出的并行環(huán)境下目標(biāo)函數(shù)變分結(jié)果向耦合伴隨方程對應(yīng)的網(wǎng)格單元裝配方式簡單有效。

3) 文中構(gòu)建的聲爆模塊計算精度較高，伴隨模塊計算的地面聲爆目標(biāo)函數(shù)對近場聲壓的梯度與有限差分結(jié)果較為一致；優(yōu)化測試算例表明，耦合伴隨系統(tǒng)具有較高的優(yōu)化設(shè)計效率，地面過壓優(yōu)化效果比較明顯。

本文基于結(jié)構(gòu)網(wǎng)格求解器，開展了流場/聲爆耦合伴隨優(yōu)化方法研究，驗證了文中提出的插值、變分原則的有效性，為氣動力、聲爆綜合一體化設(shè)計提供了研究基礎(chǔ)。在此基礎(chǔ)之上，將進(jìn)一步基于高效的耦合伴隨系統(tǒng)開展氣動力/聲爆綜合設(shè)計以及聲爆抑制原則等研究。