基于連續伴隨方法的高超聲速飛行器高精度氣動優化

高昌，李正洲，黃江濤，賀元元，吳穎川，樂嘉陵，桂豐

1.中國空氣動力研究與發展中心 高超聲速沖壓發動機技術重點實驗室，綿陽 621000 2.中國空氣動力研究與發展中心，綿陽 621000 3.中國航發四川燃氣渦輪研究院，成都 610500

高超聲速飛行器具有實現全球快速到達和低成本進入空間的潛力，目前各主要航天大國都正在開展對高超聲速飛行器的深入研究。高超聲速飛行環境下，飛行器受到氣動力、熱載荷嚴酷，高超聲速飛行器機體和推進系統的設計面臨推阻平衡、升重平衡等難題，機體/推進一體化設計方法是解決這些難題的有效方法[1-2]。目前的一體化設計方法主要是基于無黏理論，高超聲速來流條件下黏性對于飛行器氣動性能有重要影響，實際應用中還要權衡幾何修型導致的性能惡化[3]。氣動設計作為概念設計的起始點，性能優良的氣動外形有利于整個總體設計方案的封閉性和可實現性。在考慮黏性的情況下，針對高超聲速飛行器開展高精度氣動優化可以在一體化設計的基礎上進一步提高飛行器的氣動性能。

高精度氣動優化需要借助精細化CFD(Computational Fluid Dynamics)工具開展氣動性能預測，優化過程通常采用進化策略或梯度策略[4]。進化策略基于遺傳算法、粒子群算法等優化算法對整個設計空間進行全局尋優，對于多設計變量優化問題，龐大的CFD計算開銷是面臨的主要困難，代理模型在一定程度上可以降低CFD計算開銷，但仍存在精度和泛化性問題[5-6]。梯度策略根據目標函數梯度信息對設計空間進行局部尋優，伴隨方法具有計算量與設計變量數目無關的優勢，是氣動優化中高效獲取目標函數梯度的方法，已成熟應用于亞、跨和低超聲速氣動優化問題[7-9]。但對于高超聲速流動，伴隨方程推導需要考慮黏性輸運系數變分的貢獻；同時高超聲速流動中存在的強激波間斷對伴隨方程對流項離散格式的精度和魯棒性提出了更高的要求。

針對伴隨方法在高超聲速優化問題中的應用，國內外研究人員進行了一系列探索。Copeland等在考慮高溫化學非平衡效應的情況下推導了連續伴隨方程，基于Steger-Warming格式對伴隨方程對流項進行離散，在高馬赫數來流下針對二維翼型對連續伴隨方法的計算精準度進行了驗證[10]；Kline等采用連續伴隨方法針對高超聲速進氣道開展了優化研究，但忽略黏性輸運系數變分對伴隨方程的貢獻，伴隨方程對流項離散采用JST(Jameson-Schmidt-Turkel)中心差分格式[11]。宋紅超等基于離散伴隨方法對高超聲速單邊膨脹噴管開展了優化研究[12]；高昌等基于連續伴隨方法對二維高超聲速進氣道開展了優化研究[13]。目前對于伴隨方法在高超聲速優化問題中的應用還處于初步探索階段，已有的研究工作主要針對二維或簡單流場，對伴隨方法在三維復雜高超聲速流動中的應用缺乏深入研究和驗證。

為進一步探索伴隨方法在高超聲速優化問題中的應用，建立高超聲速飛行器高精度氣動優化方法，本文基于Navier-Stokes(NS)方程推導了連續伴隨方程和氣動力目標函數對應的邊界條件和靈敏度公式，考慮層流黏性輸運系數對伴隨方程的貢獻，探究基于二階熵修正Roe格式的伴隨方程對流項離散格式，構建適用于高超聲速流動的連續伴隨求解器；基于連續伴隨求解器構建針對高超聲速飛行器的高精度氣動優化框架；通過二維翼型和Sanger飛行器機翼優化算例驗證連續伴隨求解器的計算精度，以Sanger飛行器升阻比為優化目標開展基于連續伴隨方法的高超聲速高精度氣動優化方法的應用研究。

1 優化方法

1.1 連續伴隨方程

伴隨方法的基本思想是將流動控制方程作為優化問題的約束條件，對于氣動力優化問題，目標函數可以表示為統一形式：

(1)

式中：J為伴隨方程目標函數；Sw為壁面邊界；j為伴隨方程目標函數積分項；p為壓力；σ為應力；d為氣動力目標函數方向向量；I為單位矩陣；n為壁面法向量；dS為積分面元。

對于不同的氣動力目標函數：

(2)

(3)

式中：R為流動控制方程殘差；U為流動控制方程守恒變量；Fc為對流項；Fv為黏性項；μ為黏性輸運系數。通過Lagrange乘子Ψ將NS方程引入目標函數：

(4)

式中：D為計算域；dV為積分體積微元。對式(4)進行變分運算：

(5)

式中：δ為變分運算符；δSw為壁面變分。目標函數的變分項根據曲面變分公式可以簡化為

(6)

對于式(5)中NS方程的變分項，由于黏性項中包含流場變量的梯度，為了簡化推導過程，將NS方程視為流場變量和流場變量梯度的函數：

(7)

對式(7)進行變分運算，根據求導交換法則：

(8)

式中：各Jacobian矩陣可表示為

(9)

式中：Ac為無黏通量Jacobian矩陣；Av為黏性通量Jacobian矩陣；Aμ為黏性輸運系數變分項對伴隨方程的貢獻；D為黏性項關于守恒變量梯度的Jacobian矩陣。對于高超聲速流動，黏性輸運系數對近壁面流動起主導作用，對壁面壓力和摩阻分布有顯著影響，進而影響整個飛行器的氣動性能。湍流輸運系數變分項需要對湍流模型方程進行變分推導，先只考慮層流輸運系數變分項的貢獻，其中層流黏性系數μlam和層流熱傳導系數Klam的變分項可以通過Sutherland公式和層流Prandtl數Prlam計算：

(10)

式中：下標lam表示層流參數；T為溫度；μref和Tref為Sutherland公式中的參考黏性系數和參考溫度；Sref為Sutherland公式中的常數；cp為定壓比熱。將式(8)代入式(5)的第3項，采用Gauss積分定理和壁面無滑移條件簡化可以得到：

(11)

式中：u為速度矢量；φ為與動量對應的伴隨變量；S∞為遠場邊界；K為熱傳導系數；Σφ和η為中間變量；ψρE為與能量對應的伴隨變量；μ為黏性系數;?n表示法向梯度，式(11)中最后一項體積分項即為與NS方程對應的伴隨方程，設ψρ為與密度ρ對應的伴隨變量，Lagrange乘子Ψ=[ψρ,φ,ψρE]是與流場守恒變量相對應伴隨變量。其中：

(12)

式中：h為焓值。推導式(11)的目的是建立目標函數與壁面變分的導數關系，進一步推導需要對NS方程的壁面邊界條件進行變分運算，建立目標函數變分δJ′和壁面變分δSw的關系。

1.2 遠場邊界條件

對于氣動力目標函數，在遠場邊界可以忽略黏性的影響，式(11)第2項表示在遠場邊界?！奚系姆e分。高超聲速計算通常采用小遠場邊界，伴隨方程邊界條件需要結合流動控制方程邊界條件進行分析。伴隨變量在遠場邊界需滿足

ΨT(Ac·n)δU=0

(13)

將式(13)展開，寫為初始變量的形可得

(14)

式中：un為邊界法向速度分量；γ為比熱比。

對于超聲速入口邊界，邊界流場變量密度、速度和壓力均由來流計算條件指定，此時

(15)

已滿足式(13)，不需要對伴隨變量施加邊界條件。

對于超聲速出口邊界，邊界流場變量均由流場內點決定，此時

(16)

伴隨變量需滿足

(17)

在超聲速出口邊界Jacobian矩陣均為非零特征值，此時伴隨方程邊界條件為

(18)

對于亞聲速出口邊界，需要給定出口背壓，此時：

(19)

伴隨變量需滿足

(20)

求解得

(21)

式中：c為聲速；ψρE可由流場內點插值決定。

1.3 壁面邊界條件

等溫壁面上流動控制方程的邊界條件為

(22)

式中：Tw為給定的壁面溫度。線化處理后的變分形式為

(23)

將式(23)代入式(11)，為消除包含流場變量的積分項，得到連續伴隨方程在等溫壁面的邊界條件：

(24)

式(11)可以轉化為

(25)

式(25)中積分項表示氣動力目標函數關于局部壁面變分的靈敏度：

giso=μn·Σφ·?nu+K?nψρE?nT

(26)

式中：下標iso表示等溫壁面。

對于絕熱壁面條件，流動控制方程的邊界條件為

(27)

線化處理后的變分形式為

(28)

(29)

氣動力目標函數在絕熱壁面關于局部壁面變分的靈敏度為

(30)

式中：下標adi表示絕熱壁面。

1.4 數值方法

采用基于非結構網格的有限體積法對流動控制方程和連續伴隨方程進行求解。對于流動控制方程，對流項離散格式決定了流場波系結構的求解精度，對流項離散格式在強激波間斷等位置需要具備強魯棒性和低耗散性，在保證計算收斂的同時提高激波捕獲精度。為此采用MUSCL(Monotonic Upstream-centered Scheme for Conservation Laws)格式[14]對界面變量進行重構，網格界面上對流項采用穩定性好、計算精度高的Roe格式[15-16]進行離散：

(31)

式中：下標ij表示網格控制體j對i的貢獻，下標R、L表示重構后界面兩側流場變量對應的值；下標Roe表示界面兩側流場變量Roe平均后對應的值；PRoe為Roe平均后特征向量構成的矩陣；ΛRoe為Roe平均后對應的特征值。為在強激波和聲速點位置滿足熵條件，對Roe平均特征采用Harten-Yee型熵修正[17]：

(32)

式中：ε*為常數，通常取0～0.4；cRoe為Roe平均的聲速。

(33)

1.5 優化框架

單次優化循環的主要計算耗時集中在流動控制方程和連續伴隨方程的求解，連續伴隨方程求解計算量與流動控制方程相當，通過一次流動控制方程和連續伴隨方程的求解可以得到目標函數關于所有設計變量的梯度，因此伴隨方法具有計算量與設計變量數目無關的優勢，適用于多設計變量優化問題。

圖1給出了基于伴隨方法的梯度優化流程，首先對優化對象建立幾何參數化模型，采用FFD(Free Form Deformation)方法[21-22]對優化對象壁面進行參數控制。FFD方法可以較少的設計變量實現對變形物體的光滑連續控制；對于三維復雜外形可以方便地建立參數化模型；具備較強的局部變形控制能力[23]。通過網格變形方法將壁面變形傳遞到整個計算網格。優化對象氣動性能通過求解RANS(Reynolds-Averaged Navier-Stokes)方程進行評估，伴隨求解器讀取流場變量對連續伴隨方程進行求解得到目標函數壁面靈敏度，根據FFD方法將壁面靈敏度映射為目標函數關于設計變量的梯度。梯度優化算法采用SQP(Sequential Quadratic Programming)算法[24]，SQP算法根據目標函數梯度信息開展優化循環。

圖1 優化流程

2 流場求解器驗證

為驗證流場求解器精度，對Sanger飛行器在高超聲速來流條件下的氣動性能計算[25]，并與AHL3D(Airbreathing Hypersonic Lab 3D)計算結果進行對比。AHL3D是中國空氣動力研究與發展中心開發的大規模三維數值模擬軟件，經過了大量研究和工程應用驗證[26]。AHL3D計算采用結構網格如圖2所示，結構網格分為49個塊，第1層網格高度為10-6，網格量約為1 080萬，對流項離散采用AUSM格式，湍流模型采用k-ω剪切應力傳輸(Shear Stress Transter,SST)模型。求解器采用非結構混合網格，壁面邊界層采用三棱柱網格，第1層網格高度為10-5，流場中間采用四面體網格，對流項離散為Roe格式，湍流模型采用帶可壓縮修正的SA模型，計算條件為來流馬赫數6.0、高度35 km，壁面條件為1 000 K等溫壁面。圖3給出了不同攻角下升阻力系數對比，本文求解器得到的升阻力系數相對AHL3D偏小，阻力系數相差在5.0%以內，升力系數相差在3.0%以內。

圖2 Sanger飛行器結構網格

圖3 Sanger飛行器升阻力系數對比

為對流場求解器網格無關性進行研究，通過調整非結構網格的表面網格尺寸和空間網格增長率，得到了如圖4所示的3套網格，網格量分別為267萬、494萬和1 108萬，來流攻角為4°。表1給出了3套網格升阻力系數計算結果，其中相對誤差為相對網格3計算結果的誤差，升力系數最大相對誤差為1.9%，阻力系數最大相對誤差為0.8%。

圖4 Sanger飛行器非結構網格

表1 Sanger飛行器氣動力系數網格靈敏度

3 伴隨求解器驗證

采用NACA64A-204翼型對伴隨方法求解結果進行驗證，計算網格采用結構網格的形式(圖5)，采用FFD方法對翼型進行參數化控制(圖6)，設計變量為FFD控制點沿法向的位移量，來流條件和計算設置與第2節中相同，攻角為4°。

圖5 NACA64A-204翼型計算網格

圖6 NACA64A-204翼型FFD模型

圖7 NACA64A-204翼型伴隨方程收斂歷程

圖8 NACA64A-204翼型ψρE云圖

圖9 NACA64A-204翼型伴隨方法與有限差分法梯度對比

4 Sanger飛行器機翼優化

采用伴隨優化框架對Sanger飛行器進行驗證和優化研究，Sanger飛行器為德國二級入軌空天運輸系統的一級飛行器，采用后掠三角翼兼顧中低速氣動性能。機翼為影響整個飛行器氣動性能的主要部件，采用FFD方法對飛行器機翼進行參數化控制，上下表面共設置64個設計變量(圖10)，設計變量取值范圍以1～8設計變量所在剖面翼型厚度為基準，為避免機翼前緣過度變薄，控制機翼前緣變形的設計變量的取值范圍為基準的20%，其余設計變量的取值范圍設定為基準的50%。計算網格采用第2節中的網格2(圖4(b))。

圖10 Sanger飛行器機翼FFD控制模型

為對連續伴隨方法進行應用驗證，選擇Sanger飛行器巡航狀態進行優化，巡航來流狀態為馬赫數6.0、高度35 km、攻角4°，計算設置與第2節中相同，目標函數選擇升阻比。圖11給出了壁面靈敏度積分值和伴隨變量殘差的收斂歷程，圖中ψρw為與守恒變量ρw對應的伴隨變量，計算過程伴隨變量殘差一致收斂，經過2萬次迭代求解，翼型壁面靈敏度積分值基本穩定，伴隨變量殘差降低了2個數量級以上。圖12給出了伴隨變量ψρE在對稱面的分布云圖，可見采用的數值離散格式較好地捕捉到了伴隨流場的波系結構。圖13 給出了升阻比壁面靈敏度分布，紅色區域表示壁面內凹時飛行器升阻比增大，藍色區域反之。圖14給出了伴隨方法計算得到的升阻比關于設計變量梯度與有限差分結果的對比，由于設計變量數目較多，只選擇同一控制面上的設計變量9～16(機翼下表面)、41～48(機翼上表面)進行有限差分驗證，差分步長選擇與NACA64A-204算例相同，可見不同步長得到的有限差分結果一致性較好，同時前向差分和后向差分結果符合較好。從對比結果看，相對二維情況下此時伴隨方法與有限差分法結果在數值上的差異更明顯，這是由于三維情況下流場激波結構更加復雜，同時機翼下表面氣流受到的壓縮更強，這種差異在設計變量9～16上更加明顯；但兩種方法得到的結果變化規律和正負分布一致，能夠為梯度優化算法提供正確的優化方向。

圖11 Sanger飛行器伴隨方程收斂歷程

圖12 Sanger飛行器對稱面ψρE云圖

圖13 升阻比壁面靈敏度

圖14 Sanger飛行器目標函數梯度對比

圖15給出了飛行器升阻比和3個站位翼型最大厚度的變化歷程，優化循環共迭代29步，3個站位的位置分別為0.57(站位1)、0.72(站位2)和0.86(站位3)展寬位置，優化過程中站位1處翼型最大厚度變化較小，站位2和站位3處翼型最大厚度出現明顯下降。飛行器升阻比在第11～20迭代步出現震蕩收斂，這是由于優化歷程趨近收斂時，伴隨梯度精度的不足在一定程度上對優化歷程產生影響，此時SQP算法通過減小試探步長在一定程度上克服了這種影響，從而導致飛行器升阻比出現上下震蕩。表2給出了優化前后飛行器升阻比、升阻力系數和3個站位翼型相對厚度的變化，優化后飛行器升阻比增加5.0%，升力系數增加1.8%，阻力系數減小2.9%，3個站位翼型厚度分別減小1.8%、23.6%和25.4%。

表2 優化結果

圖15 優化歷程

圖16給出了優化后飛行器機翼構型變化，優化后機翼上表面局部出現了小幅隆起，機翼前緣變薄，下表面出現明顯的內凹，機翼構型變化與圖13 中壁面靈敏度的分布一致，驗證了FFD方法具備較強的局部變形控制能力。圖17給出了3個站位剖面翼型的變化，機翼下表面內凹使得剖面翼型前緣厚度減小、彎度增加。

圖16 優化后Sanger飛行器機翼構型

圖17 3個站位機翼構型變化

圖18和圖19分別給出了優化前后3個站位壓力系數沿流場截面和翼型表面的變化，優化后機翼前緣迎風面積減小，頭部激波強度減弱，減小了飛行器阻力；氣流經過頭部激波后出現明顯膨脹，在內凹拐點處經過二次激波壓縮，大幅增加了機翼下表面中后部壓力，顯著提高了飛行器的升力；優化后站位2上表面壓力由于機翼隆起出現了小幅度下降，但機翼上表面壓力變化較小，下表面壓力變化是飛行器升阻比提升的主要因素。圖20給出了優化前后機翼下表面的流場變化，從三維流動分析來看，優化后機翼前緣壓力明顯減小，緊貼機翼下表面的二次激波封閉了高壓氣流，減少了機翼前緣高壓氣流的泄流，機翼下表面高壓區范圍增大，達到了減阻增升的效果。

圖19 3個站位壓力系數分布變化

圖20 機翼附近流場變化

5 結 論

本文構建了基于連續伴隨方法的高超聲速飛行器高精度氣動優化方法，通過對二維翼型和Sanger飛行器機翼優化驗證和應用得到以下結論。

1)基于熵修正的Roe格式構造了伴隨方程對流項離散格式，在二維和三維伴隨流場求解過程中伴隨變量殘差一致收斂，較好地捕捉到了與初始流場共軛的伴隨流場波系結構，表明采用的數值格式具有較好的魯棒性和低耗散性。

2)基于構造的數值格式，建立了適用于高超聲速流動的連續伴隨求解器，通過二維翼型算例和Sanger飛行器機翼優化算例的驗證分析可知，在高超聲速來流條件下構造的連續伴隨求解器能夠較好地給出氣動力目標函數的梯度信息，但計算精度一定程度上仍受到強激波的影響。

3)FFD方法具備較強的局部變形控制能力，與基于連續伴隨方法的梯度優化相結合，實現了針對高超聲速飛行器的高精度優化。

4)Sanger飛行器機翼經過優化后，機翼前緣內凹，彎度增加，前緣迎風面積較小，產生的二次激波壓縮使機翼下表面壓力顯著增加，飛行器升阻比增加5.0%，達到了減阻增升的效果。

5)經過本文的驗證和應用研究可知，連續伴隨方法是實現高超聲速飛行器高精度氣動優化的可行方法。