基于流場/聲爆耦合伴隨方程的超聲速公務機聲爆優化

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

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

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

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

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

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

1 優化系統中的學科分析模塊

1.1 并行化CFD求解器PMB3D

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

(1)

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

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

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

(2)

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

(3)

(4)

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

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

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

圖2 聲壓計算與試驗數據對比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 流場伴隨方程求解梯度

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

(5)

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

L=I+ΛTR

(6)

式中：L為聲爆目標函數。

對式(6)進行求導，有

(7)

(8)

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

(9)

(10)

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

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

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

2.2 聲爆伴隨方程求解梯度

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

(12)

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

(13)

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

(14)

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

圖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

數采用以下形式：

式中：pT為聲爆設計目標特征?？梢钥闯雎暠殡S方程梯度計算結果與差分結果較為一致，可以為耦合伴隨系統提供準確的地面聲爆目標函數對近場聲壓的梯度。

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

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

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

p0=T(W)

(15)

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

(16)

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

(17)

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

(18)

(19)

(20)

3 并行環境下近場聲壓提取與變分結果裝配

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

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

圖6 并行環境下的網格分布與“聲壓盒”Fig.6 Grid distribution and “sound pressure box” in parallel environment

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

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

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

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

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

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

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

7) 轉換為CFD坐標系。

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

9) 流場伴隨方程求解。

4 聲爆優化設計體系

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

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

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

5 典型超聲速公務機聲爆優化

基于巡航馬赫數為1.5的小型超聲速公務機的基本構型，開展超聲速巡航狀態下的聲爆優化設計研究。

其優化數學模型為

Constraints:

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

圖9給出了基于NURBS基函數的自由式變形參數化示意圖，共采用60個控制頂點實現機身、機翼氣動外形參數化建模；設計變量為圖9中x方向4～6個截面，以及機翼展向4個截面，近場“包圍盒”處于機身下方60 m處。圖10為初始外形的對稱面壓力云圖，可以看出按照馬赫角拉伸的網格拓撲能夠較為清晰地捕捉空間激波形態。

圖8 CFD網格分布Fig.8 Distribution of CFD grid

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

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

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

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

圖11 耦合伴隨系統收斂歷程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 目標函數優化收斂歷程Fig.14 Convergence history of objective function optimization

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

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

圖16 地面過壓優化過程Fig.16 Optimization process of ground pressure

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

6 結 論

基于并行化結構網格RANS求解器PMB3D以及伴隨方程求解器PADJ3D，結合增廣Burgers聲爆預測與伴隨方程，開展了流場/聲爆耦合伴隨方程構造、求解方法以及聲爆優化研究。

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

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

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