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

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

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

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

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

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

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

1.1 并行化CFD求解器PMB3D

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

(1)

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

1.2 基于Burgers方程的地面聲爆預測

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

(2)

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

(3)

(4)

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

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

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

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

圖3 NASA1021標模聲爆預測對比Fig.3 Sonic boom prediction comparison of NASA1021 standard model

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

2.1 流場伴隨方程求解梯度

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

(5)

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

L=I+ΛTR

(6)

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

對式(6)進行求導，有

(7)

(8)

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

(9)

(10)

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

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

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方法隱式時間推進方法，其中，離散伴隨方程的邊界條件采用矩陣形式處理方式，黏性項采用薄層近似，離散伴隨方程求解時，并行機制依然采用單元數(shù)衡量的負載平衡、對等式計算以及MPI消息傳遞模式。本文求解器采用了多塊對接網(wǎng)格技術，MPI傳遞的信息是各個進程中分割面上的兩層虛網(wǎng)格上的伴隨變量信息，詳細推導可以參考文獻[19]。

2.2 聲爆伴隨方程求解梯度

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

(12)

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

(13)

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

(14)

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

圖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為聲爆設計目標特征。可以看出聲爆伴隨方程梯度計算結(jié)果與差分結(jié)果較為一致，可以為耦合伴隨系統(tǒng)提供準確的地面聲爆目標函數(shù)對近場聲壓的梯度。

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

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

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

p0=T(W)

(15)

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

(16)

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

(17)

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

(18)

(19)

(20)

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

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

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

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

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

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

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

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

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

5) 聲爆預測迭代推進計算。

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

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

8) 按編號、坐標文件將變分結(jié)果按對應的進程編號輸出，向各個進程裝配。

9) 流場伴隨方程求解。

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

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

本文發(fā)展的低聲爆超聲速公務機優(yōu)化技術屬于AMDEsign中的耦合伴隨模塊，其中CFD方法采用中心格式、SST兩方程湍流模型、多重網(wǎng)格加速收斂以及大規(guī)模并行計算；梯度求解采用對應的并行化伴隨方程與梯度求解器；參數(shù)化方法采用基于NURBS基函數(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 典型超聲速公務機聲爆優(yōu)化

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

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

Constraints:

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

圖9給出了基于NURBS基函數(shù)的自由式變形參數(shù)化示意圖，共采用60個控制頂點實現(xiàn)機身、機翼氣動外形參數(shù)化建模；設計變量為圖9中x方向4～6個截面，以及機翼展向4個截面，近場“包圍盒”處于機身下方60 m處。圖10為初始外形的對稱面壓力云圖，可以看出按照馬赫角拉伸的網(wǎng)格拓撲能夠較為清晰地捕捉空間激波形態(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é)合最終的導數(shù)計算式(9)，可以從伴隨變量云圖分布上定性看出流動敏感性區(qū)域。同樣，結(jié)合耦合伴隨變量云圖分布以及伴隨方程耦合變分項的裝配位置可以看出，隨著偽時間推進迭代，由空間近場向物面“傳播”，因此，耦合伴隨方程呈現(xiàn)先收斂，進而階段振蕩，最終完全收斂的特征。

圖14給出了地面聲爆目標函數(shù)優(yōu)化過程，經(jīng)過9代優(yōu)化，聲爆抑制效果逐漸收斂，圖15給出了地面聲爆目標函數(shù)對近場聲壓的梯度優(yōu)化前后對比，圖16給出了優(yōu)化前后地面過壓對比?？梢钥闯?，伴隨優(yōu)化效果較為明顯，第2道激波的峰值明顯降低。由于機頭的控制點沒有作為設計變量，第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 目標函數(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)化前后感覺噪聲級的頻譜特性，橫坐標為1/3倍頻程頻段，縱坐標為響度級。由于頻譜特性峰值與激波上升時間密切相關，而第1道激波上升時間與機頭形狀關系密切，如前面所述，文中沒有將機頭作為設計變量，圖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聲爆預測與伴隨方程，開展了流場/聲爆耦合伴隨方程構(gòu)造、求解方法以及聲爆優(yōu)化研究。

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

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

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