劉紅陽,宋超,羅驍,周鑄,呂廣亮
中國空氣動力研究與發展中心 計算空氣動力研究所,綿陽 621000
經濟性和環保性一直是民航客機永恒不變的追求之一,氣動阻力與之緊密相關。客機各部件的阻力一般分為壓差阻力和摩擦阻力,湍流狀態下的摩擦阻力遠高于層流狀態下的摩擦阻力,因此使客機表面保持大部分層流區域將會大大減小客機的總阻力。早期已有研究表明,如果在機翼和尾翼保持一定的層流區,總阻力將降低30%[1]。然而機翼后掠角較大,且對升力和力矩有一定要求,自然層流機翼的設計較難實現。
隨著渦扇發動機技術的發展,大涵道比發動機以其噪聲低、燃油效率高、使用壽命長等優勢被更多的應用于民航客機。一般尺寸的民用渦扇發動機短艙的阻力占客機總阻力的5%左右,隨著發動機涵道比的不斷增大,短艙的尺寸隨之增大,短艙阻力的占比也不斷提升。由于飛機設計中對短艙的升力和力矩等特性不作要求,且沒有橫流不穩定和附著線轉捩等問題,自然層流短艙的設計更容易實現,進而達到減小阻力的目的[2]。
國外學者對層流短艙的研究工作開展的比較早。從20 世紀80年代中期開始,美國NASA蘭利研究中心開展了自然層流短艙的風洞和飛行試驗研究[3]。1984年,Younghans 等[4]對設計的單獨自然層流短艙和翼吊自然層流短艙開展了風洞試驗驗證研究,獲得了較大范圍的層流區域。1990年,Radespiel 等[5]研究了層流短艙設計的可行性,使用數值流動方法實現了層流短艙設計,在巡航飛行狀態下自然層流長度占短艙長度的60%。1998年,德國宇航院DLR、羅爾斯·羅伊斯公司和Motorenund Turbinen-Union Friendrichshafen(MTU)合作開展了發動機短艙自然層 流 和 混 合 層 流 飛 行 試 驗[6],在VFW614/ATTAS 飛機上獲得了短艙長度60%的層流區。2011年,Lin[2]對自然層流短艙開展了氣動設計、風洞試驗和轉捩預測方法的驗證研究。2014年,Vermeersch 等[7]通過求解歐拉方程并分析邊界層穩定性,對自然層流短艙和混合層流短艙開展了研究,分析了除冰和加熱等對轉捩的影響。隨著自然層流短艙技術的逐步發展和工業制造水平的不斷提升,波音公司已成功將其應用于波音787 系列客機上。
國內學者對自然層流短艙的研究起步較晚。2014年,何小龍等[8]基于EFFD 方法實現了軸對稱自然層流單獨通氣短艙和帶動力短艙的優化設計,獲得了短艙長度40%以上的層流覆蓋區。2016—2017年,Zhong 和Li 發展了短艙截面造型設計方法,根據二維型線生成三維短艙并進行評估和優化[9-10]。2019年,杜璽等[11]開展了自然層流短艙的氣動設計、數值計算和風洞試驗驗證研究,獲得了短艙長度30%~55%的層流區,Wang和Sun 等[12]使用有限差分算法對跨聲速自然層流短艙進行了優化設計研究,層流區面積增加16.64%,總阻力系數減小11.6 counts,孟曉軒等[13]基于線性穩定性分析方法,將雙eN方法同RANS 方程求解器耦合,研究了來流馬赫數、雷諾數、湍流度以及攻角對短艙轉捩的影響。2020—2021年,胡驍等[14-15]通過設計軸對稱層流短艙基準面,并二次開發CATIA 工具,實現了三維非軸對稱自然層流短艙設計。
實際上短艙屬于非軸對稱構型,受機翼機身的氣流影響,具有復雜的流動狀態。傳統的設計方法以二維或二維假設的形式開展,沒有考慮短艙較強的三維流動效應,這種方式會損失設計過程中獲得的層流區面積;如果直接以非軸對稱短艙為對象進行優化設計,則面臨設計變量數目多、CFD 求解次數多、計算成本高的問題。
本文提出一種基于生成拓撲映射(Generative Topographic Mapping,GTM)[16-17]的壓力分布反設計方法,充分考慮復雜三維流動效應,直接對短艙三維外形進行設計,通過對短艙外形及壓力分布組成的數據集進行降維處理,有效減少所需樣本數目,且設計過程中無需反復計算流場,能夠大大提高設計效率。
結合自由曲面變形技術(Free Form Deformation, FFD)[18]、RANS 方程、基于SST(Shear Stress Transport)湍 流 模 型 的γ-轉 捩 模型[19]、基于徑向基函數插值的網格變形方法[20-21]、GTM 模型和遺傳算法,構建了考慮三維流動效應的自然層流(Natural Laminar Flow,NLF)短艙優化設計方法,利用該方法對單獨通氣NLF 短艙進行了優化設計研究。
生成拓撲映射是一種非線性隱變量模型,能夠實現高維數據到低維空間的映射,并保持原有的拓撲相對關系,目前廣泛應用于數據分析及可視化分析。
設 有D維 數 據 集T={t1,t2,…,tN},但 數 據本質是L維的(L<D)。映射函數y(x,W)將L維空間的K個隱變量點x={x1,x2,…,xK}映射到D維數據集T。采用中心位于y(x,W)且方差為β的高斯函數作為D維數據集T的分布:
式中:p(t|x,W,β)為條件概率密度函數,對x積分可得:
其中:p(x)為隱變量空間維度的先驗概率密度函數:
表示均勻分布在隱變量空間內的K個點上。
由此得到:
通過最大對數似然函數可求得參數W和β:
分布p(t|W,β)為中心固定的帶約束混合高斯模型。GTM 算法可歸納為:已知隱變量空間中x的分布和數據空間的數據集合T,使用最大似然法求解參數W和β。常采用期望最大化(Expectation Maximization,EM)算法訓練GTM 模型參數[22]。
圖1 定義了非軸對稱短艙的4 個基準面,以短艙的頂部截面為0°基準面,沿順時針方向依次為90°基準面、180°基準面和270°基準面。

圖1 非軸對稱NLF 短艙各個基準面的定義Fig.1 Reference planes of non-axisymmetric NLF nacelles
相隔180°的2 個基準面之間的層流區相互影響很小,因此可以先對其中一組對向的2 個基準面分別進行反設計,然后在此基礎上設計另一組對向的2 個基準面。2 組基準面的設計順序不影響設計結果。圖2 給出了NLF 短艙的設計思路:

圖2 非軸對稱NLF 短艙的設計思路Fig.2 Design idea of non-axisymmetric NLF nacelles
步驟1在短艙初始外形baseline1 的基礎上,分別對90°基準面和270°基準面進行壓力分布反設計,獲得相應的短艙外形design1.1 和design1.2。
步驟2應用FFD 技術,將design1.1 的90°基準面和design1.2 的270°基準面匹配到baseline1 外形上,獲得新的短艙初始外形baseline2。
步驟3在baseline2 的基礎上,分別對0°基準面和180°基準面進行壓力分布反設計,獲得相應的短艙外形design2.1 和design2.2。
步驟4應用FFD 方法,將design2.1 的0°基準面和design2.2 的180°基準面匹配到baseline2外形上,獲得最終的短艙反設計外形design3。
圖3 給出了NLF 短艙各個基準面的壓力分布反設計流程,具體如下:

圖3 NLF 短艙各個基準面的壓力分布反設計流程Fig.3 Pressure distribution inverse design process of NLF nacelle reference planes
步驟1針對短艙外形和基準面上的設計變量進行實驗設計,采用拉丁超立方方法選取樣本點S,并計算對應的壓力分布Cp。
步驟2根據高維數據集{S,Cp}創建并訓練GTM 模型,對數據集進行降維處理,獲得相應隱空間映射關系。
步驟3結合預設目標壓力分布,采用全局優化算法在隱空間尋優,GTM 模型輸出預測的最優設計變量與壓力分布。
步驟4對最優設計變量進行CFD 校驗,如果相應壓力分布和層流區長度滿足設計要求,則設計結束,否則需要重新設定目標壓力分布或更新樣本數據集,重復以上步驟。
基于GTM 模型的壓力分布反設計有如下優點:
1) 目標壓力分布的設定方式靈活。不必給出繞短艙基準面一周的完整壓力分布,按照設計需求僅給定基準面上一段或多段壓力分布即可進行設計。
2) 尋優效率高。GTM 模型將短艙構型及其壓力分布組成的高維數據集進行降維處理,全局優化算法在低維隱空間尋優,極大提高了優化效率。
3) 計算成本低。GTM 模型通過樣本訓練,將高維數據集映射到低維隱空間,這種高精度映射關系避免了CFD 求解器的反復調用,縮短了設計時間,降低了計算成本。
初始的非軸對稱NLF 短艙模型尺寸參數如圖4 所示,短艙外罩軸向最大長度Lmax=5.544 m,周向最大直徑Dmax=3.719 m,最大直徑的軸向位置L1=1.865 m,下垂角Φ1=2.49°,外罩0°和180°基準面的船尾角分別為θ1=16.11°和θ2=13.75°,短艙外罩出口直徑Doutlet=2.642 m,后緣厚度thicknessoutlet=0.004 m。對短艙的0°,90°,180°和270°基準面(如圖1 所示)進行壓力分布反設計,目的是推遲邊界層轉捩的發生,擴大自然層流區面積,進而減小短艙的摩擦阻力。

圖4 初始短艙模型尺寸參數Fig.4 Size of nacelle baseline model
短艙模型的計算網格如圖5 所示,網格節點數為633 萬,邊界層內第1 層網格高度為1×10-5m,邊界層內網格增長率為1.2。設計狀態為高度H=11 km,自由來流馬赫數Ma=0.85,雷諾數Re=3.1×107,來流初始湍流度Tu=0.1%,黏性比Rt=5.0,攻角α=5.0°,側滑角β=0°。

圖5 非軸對稱層流短艙計算網格Fig.5 Computation grid of non-axisymmetric nacelles
采用自由曲面變形方法對短艙進行參數化,短艙各個基準面上的設計變量分布如圖6 中實心圓點所示,設計變量數目為8,設定降維后隱空間維度L=2,實驗設計選取的樣本數目為30。通過求解RANS 方程和γ-轉捩模型,獲得短艙壓力分布,實現邊界層轉捩的預測。采用遺傳算法在隱空間尋優,種群數目取為200,最大迭代步數取為50 步,目標函數為短艙基準面的實際壓力系數與目標壓力系數分布之間的差異:

圖6 短艙各個基準面的設計變量分布Fig.6 Distribution of design variables on nacelle reference planes
式中:M為壓力系數分布的點數;Cip為第i個點的壓力系數值;CipT為第i個點的目標壓力系數值。
順壓區是保持層流的前提,順壓區長度和順壓區內壓力梯度是影響層流發展的2 個關鍵因素,這種影響關系還與實際的流動狀態緊密相關。一般情況下,順壓區越長,層流區越長,順壓區內壓力梯度越大,越有利于推遲邊界層轉捩的發生,延長自然層流區長度,這2 個因素相互耦合,共同作用。
在短艙初始外形baseline1 基礎上,對短艙90°和270°基準面進行壓力分布反設計。圖7 給出了短艙90°基準面反設計結果,其中c為基準面當地弦長。圖7(a)中黑色圓點表示基準面弦長10%~30%區域內給定的目標壓力分布,目的是增大順壓區內的壓力梯度,紅色實線代表短艙設計外形design1.1 經CFD 校驗后提取的90°基準面壓力系數分布,其趨勢向目標壓力分布靠近,由于給定的目標壓力分布不一定具有真實的物理意義,即不存在對應的物理解,因此也不要求設計所得的壓力分布與目標壓力分布完全吻合。圖7(b)對比了短艙90°基準面設計前后的摩阻系數分布和構型,此處以摩阻系數作為判斷邊界層轉捩的依據,可以看出,90°基準面的層流區長度從當地弦長的23%延長到26.8%。圖8 是設計過程中真實壓力與目標壓力分布的差異Ep收斂示意圖,藍色三角形表示每一代所有樣本的Ep,紅色實心圓點線圖表示每一代樣本的最小Ep。可以看出,在第15 代時設計基本收斂,Ep<0.01。

圖7 短艙90°基準面反設計結果Fig.7 Design results of nacelle reference plane 90°

圖8 實際壓力與目標壓力分布的差異收斂示意圖Fig.8 Convergence of error between real and target pressure distribution
如果進行直接調用CFD 求解器的全局優化設計,其計算時間與迭代步數、樣本數目成正比,即使所有樣本并行計算,設計周期仍以天為單位。本文采用的反設計方法優化過程不調用CFD 求解器,因此尋優效率很高,在樣本計算完成后,只需要分鐘量級的時間就可獲得設計外形,一輪優化在個人工作電腦(CPU 主頻3.4 GHz)只需要50 s。
圖9 給出了短艙270°基準面反設計結果。圖9(a)中黑色圓點表示基準面弦長10%~30%區域內給定的目標壓力分布,目的是增大順壓區內的壓力梯度,紅色實線代表短艙設計外形design1.2 經CFD 校驗后提取的270°基準面壓力系數分布,其與目標壓力吻合較好;圖9(b)對比了短艙270°基準面設計前后的摩阻系數分布和構型,此處以摩阻系數作為判斷邊界層轉捩的依據,可以看出,270°基準面的層流區長度從當地弦長的25%延長到29.3%。

圖9 短艙270°基準面反設計結果Fig.9 Design results of nacelle reference plane 270°
利用FFD 參數化方法,將短艙外形design1.1 對應的90°基準面構型和design1.2 對應的270°基準面構型匹配到初始外形baseline1 上,獲得新的短艙初始外形baseline2,并利用CFD對其進行校驗。圖10 為計算得到的壓力系數分布和摩阻系數分布,并與baseline1 外形的相應數據進行對比,由于baseline1 外形是左右對稱的,因此其90°和270°基準面的壓力分布基本重合,摩阻分布略有差異,但層流區長度基本一致。baseline2 外形在90°和270°基準面上均保持了與design1.1 和design1.2 一致的較長層流區,此時baseline2 左右不再對稱,且270°基準面的順壓區更長,因此層流區更長。

圖10 baseline2 的CFD 校 驗 結 果Fig.10 Validation of baseline2 with CFD solver
在新的短艙初始外形baseline2 基礎上,對短艙0°和180°基準面進行壓力分布反設計。圖11 給出了短艙0°基準面反設計結果,與其他3 個基準面相比,由于0°基準面處在背風區,其外表面順壓區很短,邊界層轉捩發生的比較靠前,這里給定基準面弦長5%~30%區域內的目標壓力分布,如圖11(a)中黑色圓點所示,目的是減小該區域內的壓力梯度,增長順壓區,紅色實線代表短艙設計外形design2.1 經CFD 校驗后提取的0°基準面壓力系數分布,其趨勢向目標壓力分布靠近;圖11(b)對比了短艙0°基準面設計前后的摩阻系數分布和構型,此處以摩阻系數作為判斷邊界層轉捩的依據,可以看出,0°基準面的層流區長度從當地弦長的15.8%延長到20.3%。

圖11 短艙0°基準面反設計結果Fig.11 Design results of nacelle reference plane 0°
圖12給出了短艙180°基準面反設計結果,圖12(a)中黑色圓點表示基準面弦長10%~30%區域內給定的目標壓力分布,目的是增大順壓區內的壓力梯度,紅色實線代表短艙設計外形design2.2 經CFD 校驗后提取的180°基準面壓力系數分布,其趨勢向目標壓力分布靠近;圖12(b)對比了短艙0°基準面設計前后的摩阻系數分布和構型,此處以摩阻系數作為判斷邊界層轉捩的依據,可以看出,180°基準面的層流區長度從當地弦長的28.3%延長到40.5%。

圖12 短艙180°基準面反設計結果Fig.12 Design results of nacelle reference plane 180°
利用FFD 參數化方法,將短艙外形design2.1 對應的0°基準面構型和design2.2 對應的180°基準面構型匹配到外形baseline2 上,獲得最終的非軸對稱短艙設計外形design3,并利用CFD 對其進行校驗。圖13 對比了不同視圖下短艙設計前后的摩阻系數云圖,可以看出,與初始外形baseline1 相比,設計外形design3 各個基準面上的層流區長度均有不同程度的延長,其中180°基準面的延長量最大,0°,90°和270°基準面的延長量相當,這與各個基準面所處的當地流動狀態有關,還與其樣本及目標壓力分布的差異性相關,不同的樣本集和目標壓力分布驅動優化算法朝著不同的方向進化,因此取得的優化效果也不同。此外,由于0°基準面處在背風區域,流動特性決定其無法獲得較長的順壓區來保持層流。總的來說,本文所提方法設計的非軸對稱短艙能夠增大外表面自然層流區域,與設計過程中的各階段外形相比,不會損失最終成型后的短艙層流區面積。

圖13 不同視圖下短艙設計前后摩阻系數Cf云圖對比Fig.13 Comparison of Cf contour between design3 and baseline1 with different views
本文提出了一種基于GTM 模型的壓力分布反設計方法,利用該方法對非軸對稱單獨通氣短艙進行了自然層流優化設計,得到以下結論:
1) GTM 模型能夠對氣動外形及壓力分布組成的高維數據集進行有效的降維處理,僅需較少的樣本點數即可實現精確的壓力分布反設計,且無需迭代計算流場,所提設計方法具有高效便捷的優點。
2) 傳統的非軸對稱NLF 短艙設計在二維或二維假設下開展,并在三維情形下校驗設計結果,由于未充分考慮復雜三維流場的影響,使得最終短艙的層流區域有所損失。提出的設計方法計算成本較低,直接在三維外形上進行設計,過程中充分考慮三維流動效應,設計所得非軸對稱短艙能夠保持理想的層流區域。
3) 所提設計方法具有較高的工程應用價值,在下一步研究中,擬深入探究壓力分布對邊界層轉捩的影響規律,將該方法應用于帶動力單獨NLF短艙和考慮機翼機身干擾的NLF 短艙的設計中,進一步驗證該優化設計平臺的魯棒性和實用性。