聶晗，宋文萍，韓忠華,*，陳堅強，段茂昌，萬兵兵

1. 西北工業(yè)大學 航空學院 環(huán)保型超音速客機研究中心/氣動與多學科優(yōu)化設計研究所，西安 710072

2. 西北工業(yè)大學 翼型、葉柵空氣動力學國家級重點實驗室，西安 710072

3. 中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000

對超聲速民機而言，減阻是繼低聲爆技術基本攻克之后下一個亟需突破的核心關鍵問題[1]，是決定其能否持續(xù)商業(yè)運營的“攔路虎”，而層流減阻技術則是超聲速民機氣動減阻中的最重要途徑之一[2]。層流減阻技術是一種通過精細化外形設計、流動控制等手段推遲飛行器表面流動轉捩的發(fā)生，從而維持大范圍穩(wěn)定層流的氣動減阻技術。由于層流邊界層的摩擦阻力遠小于湍流邊界層，采用層流減阻技術可顯著降低飛行器表面摩擦阻力，從而大幅改善全機氣動性能。例如，對于典型超聲速民機構型，在巡航狀態(tài)下黏性阻力占總阻力40%左右，通過外形設計在機翼上表面實現(xiàn)40%范圍的自然層流，可使全機升阻比提升約5%[3]。由于巨大的減阻潛力，層流技術被國際航空運輸協(xié)會[4]在2050年飛機技術路線圖中評選入下一代飛機最有發(fā)展前景的技術。

轉捩過程的物理機理復雜、影響因素眾多，是流體力學中尚未完全解決的前沿問題之一。為了滿足工程設計應用的需求，需要結合現(xiàn)有經驗和理論研究，發(fā)展魯棒、準確的轉捩預測方法。現(xiàn)有轉捩預測方法可分為以下幾類：① 基于穩(wěn)定性理論的方法。穩(wěn)定性理論是最早用來研究轉捩現(xiàn)象的理論，其具備物理基礎，并經過了長期發(fā)展，在工程問題中得到廣泛的應用。這類方法包括：基于LST[5-9]/Bi-global[10]/PSE[11-12]的eN方法、DNS方法[13-14]等。② 轉捩經驗關系。轉捩經驗關系是通過理論研究、試驗分析等歸納得到一系列經驗性的轉捩判斷準則，包括：包絡線擬合法[15]、數(shù)據庫方法[16]、代理模型方法[17]、AHD準則[18]、Gleyzes準則[19]和C1準則[20]等。③ 轉捩模式。轉捩模式是以輸運方程等形式將轉捩關系式和CFD求解進行高效耦合的方法，包括：γ-Reθ模型[21-22]、k-ω-γ模型[23-24]、C-γ-Reθ模型[25]、基于AHD/C1準則的擾動放大因子輸運模型等[26]。

針對三維后掠層流機翼設計，現(xiàn)有轉捩預測方法中，目前在工程上應用最廣泛的是基于線性穩(wěn)定性理論的eN方法[27]。德國宇航院(DLR)[28-31]、法國宇航院(ONERA)[32]、荷蘭宇航院(NLR)[33]等均發(fā)表了比較成熟的耦合eN方法和RANS求解器的轉捩自動判斷方法，并在三維機翼、翼身組合體等構型上進行了應用。國內以周恒[34]為代表的天津大學力學系早在20世紀60年代就開展了邊界層穩(wěn)定性和轉捩的基礎研究。唐登斌[35]開展了后掠機翼可壓縮邊界層穩(wěn)定性分析研究。張坤[36-37]、朱震[38]和左歲寒[39]等分別基于線性穩(wěn)定性理論和線性拋物化穩(wěn)定性方程開展了三維機翼/翼身組合體邊界層轉捩預測研究。徐國亮和符松[40]研究了可壓縮流動中橫流擾動的失穩(wěn)及其控制。黃章峰等[41-42]研究了后掠機翼邊界層的橫流不穩(wěn)定特性以及轉捩預測，并總結了穩(wěn)定性理論在復雜外形流、非平行流、局部突變流、強三維流等工程問題中的應用。徐家寬等[43-44]通過線性穩(wěn)定性理論得到轉捩關系式，并發(fā)展了相應的轉捩模式方法，用于預測亞/跨/超聲速邊界層中流向和橫流不穩(wěn)定性誘導的轉捩。楊體浩等[45]通過eN方法對自然層流機翼翼套飛行試驗中的機翼轉捩位置進行了預測。韓忠華等[46-51]基于耦合雙eN方法的RANS求解器，開展了自然層流后掠機翼優(yōu)化設計的工作。

對于馬赫數(shù)1.3～2.0超聲速民機的層流機翼設計，基于NTS-NCF雙N因子積分策略的eN方法還存在如下問題[31]：雙N因子方法主要考慮邊界層二維Tollmien-Schlichting(TS)波和橫流(Crossflow, CF)駐波誘導的轉捩，而超聲速機翼邊界層還需要考慮三維TS斜波和CF行波不穩(wěn)定性。高雷諾數(shù)、大后掠等特點，使得超聲速民機層流機翼和傳統(tǒng)低速、跨聲速層流機翼在邊界層穩(wěn)定性上存在很大區(qū)別。一方面，隨著馬赫數(shù)的增加，二維TS波的擾動放大率大幅降低，三維TS斜波成為更不穩(wěn)定模態(tài)[7-8]，其中最不穩(wěn)定三維TS斜波波角約為45°～60°[52]；另一方面，隨著馬赫數(shù)和后掠角的增加，CF不穩(wěn)定性也有所增強，且除了頻率為0 Hz的橫流駐波不穩(wěn)定性，一般還要考慮橫流行波不穩(wěn)定性[53]。牛海波等[54]在馬赫數(shù)為6的三角翼靜音風洞試驗中發(fā)現(xiàn)了機翼表面存在橫流行波主導的轉捩。萬兵兵等[55]對飛行試驗測得的低頻信號進行分析，首次證實了在真實飛行條件下高超聲速飛行器邊界層存在橫流行波。此外，日本宇宙航空研究開發(fā)機構(JAXA)[53]和美國國家航空航天局(NASA)[56]等將eN方法應用于超聲速大后掠層流機翼設計，并通過風洞試驗、飛行試驗等開展了轉捩判據標定，也提出了在轉捩預測中要考慮TS斜波和CF行波。

針對超聲速民機機翼層流減阻設計中考慮TS斜波和CF行波的轉捩自動判斷需求，本文發(fā)展了耦合RANS求解器的eN轉捩預測方法，主要工作包括3方面：① 改進了擾動放大因子積分策略，考慮了超聲速機翼邊界層TS斜波和CF行波誘導的轉捩；② 開展了NASA超聲速試驗機翼的邊界層穩(wěn)定性分析，驗證了本文方法的正確性；③ 開展了馬赫數(shù)2.0、雷諾數(shù)1.39×107、后掠角60°的無限翼展層流機翼初步設計，探索了超聲速大后掠機翼邊界層橫流不穩(wěn)定性的抑制方法，驗證了本文方法對于超聲速民機層流機翼設計的適用性。

1 耦合RANS求解器的eN轉捩預測方法

針對超聲速機翼邊界層轉捩預測需要考慮TS斜波和CF行波的問題，本文發(fā)展了耦合RANS求解器的eN轉捩預測方法。

1.1 eN轉捩預測方法

eN方法進行轉捩預測的思想是，在邊界層中引入一個微小擾動，并通過線性穩(wěn)定性分析，得到擾動的空間演化情況；當擾動幅值增長到初始幅值的eN倍時，認為轉捩發(fā)生。其中，線性穩(wěn)定性理論假設的三維小擾動形式為

(1)

對于三維擾動應用線性穩(wěn)定性理論，為了求解特征關系式，需要補充擾動展向放大率的計算或選取方式。Mack提出[57]，可以假設擾動沿群速度方向(簡化為無黏流線方向)增長，引入流線方向和弦向的夾角θ0，得到展向擾動放大率為βi=αitanθ0，展向距離為dz=dxtanθ0。對于無限翼展后掠機翼，該條件可進一步簡化為βi=0，即假設擾動在展向的增長可以忽略[58]。Arnal等[59]通過數(shù)值試驗發(fā)現(xiàn)βi是否為0不會給擾動放大因子N的計算帶來較大差異。因此，為了簡化穩(wěn)定性方程的求解，本文對于任意的機翼構型均假設βi=0成立。

考慮到超聲速機翼邊界層中可能導致轉捩的TS斜波和CF行波，本文發(fā)展了固定波角和固定頻率方法來尋找不穩(wěn)定TS和CF波擾動，進而計算中性曲線，再通過固定展向波數(shù)/頻率方法或包絡線方法積分求解擾動放大因子N。下面將詳細介紹TS波和CF波擾動的區(qū)分、中性曲線計算、擾動放大因子積分和轉捩判據Ntr等。

1.1.1 TS和CF波擾動的區(qū)分和中性曲線計算

中性曲線是擾動放大率為0的位置，即不穩(wěn)定擾動幅值增長的起點。為了計算中性曲線，首先應區(qū)分超聲速機翼邊界層中可能存在的TS和CF波擾動，包括二維TS波(2D TS wave)、TS斜波(Oblique TS wave, OTS)、CF駐波(Stationary CF wave, SCF)和CF行波(Traveling CF wave, TCF)等，在流線坐標系下給出其頻率/展向波數(shù)分布示意圖(見圖1)。在低速和跨聲速層流機翼設計中主要考慮的是二維TS波和CF駐波，可以看出，兩者的區(qū)分方式是：二維TS波沿無黏流線方向，波角為0°，而頻率f不為0 Hz；CF駐波的頻率為0 Hz，而波角不為0°。因此對于二維TS波，可通過固定波角φw=0°計算中性曲線，并得到不穩(wěn)定TS波擾動的頻率范圍；對于CF駐波，則通過固定頻率f=0 Hz計算中性曲線，從而得到不穩(wěn)定CF波擾動的展向波數(shù)范圍。

在馬赫數(shù)1.3～2.0的超聲速機翼邊界層中，還需要考慮三維TS斜波和CF行波。從圖1可以看出，在流線坐標系下，TS斜波和CF行波的展向波數(shù)和頻率均不為0，僅從展向波數(shù)和頻率上無法作出顯式區(qū)分。一種區(qū)分TS斜波和CF行波的方式是：觀察擾動放大率隨波角的變化(圖2)，CF行波只在波角為正或負的單一方向上增長(αi<0)，且最不穩(wěn)定CF行波波角相對較大，一般大于70°；TS斜波在波角為正和負的2個方向上均有所增長，且最不穩(wěn)定TS斜波波角相對較小，一般為45°～60°。該方法適用于由TS斜波或CF行波主導的典型邊界層流動，而當兩種擾動均存在時，可進一步結合波角等信息來判斷放大因子包絡線上的不穩(wěn)定擾動類型。上述方法是一種隱式區(qū)分最不穩(wěn)定擾動波類型的方法，其對于一般問題的普適性還需后續(xù)驗證。

圖1 超聲速邊界層TS和CF波擾動頻率/波數(shù)分布示意圖

圖2 超聲速邊界層中TS斜波和CF行波擾動放大率隨波角的變化

在此基礎上，對于超聲速機翼，可分別選取固定波角φw和固定頻率f兩種方法來確定不穩(wěn)定擾動。首先，計算不同波角上的擾動放大率分布情況，并判斷擾動類型；然后，為了確定不穩(wěn)定擾動的頻率范圍，采用固定波角φw方法來計算擾動中性曲線，在其上尋找不同頻率的TS斜波和CF波擾動；最后，為了找出各頻率下不穩(wěn)定擾動的波數(shù)范圍，采用固定頻率f的方法來計算擾動中性曲線，在中性曲線上找出各頻率下不穩(wěn)定TS斜波和CF波擾動的展向波數(shù)，從而確定不穩(wěn)定擾動的頻率/波數(shù)組合。

1.1.2 擾動放大因子N的計算

對于可能導致轉捩的不同波數(shù)/頻率的擾動，從中性曲線出發(fā)，可采用兩種方法積分計算擾動放大因子。一種是固定展向波數(shù)/頻率方法，該方法最早由Mack提出[58]，是通過對擾動波引入無旋假設，并采用無限翼展機翼條件簡化得到的。該方法的好處是可以根據擾動波角、頻率等區(qū)分不同類型的擾動，如TS斜波、CF行波、CF駐波等，但由于要針對許多不同頻率/波數(shù)組合的擾動進行分析，計算量相對較大。類似的方法還有固定波長/頻率、固定波角/頻率方法等[7]。另一種方法是包絡線方法，假定擾動總是以各方向上擾動放大率的最大值進行累加，即在求解放大因子時通過改變展向波數(shù)或波角得到擾動放大率的最大值。該方法計算量相對較小，但不能判斷是哪種擾動導致轉捩發(fā)生。采用固定展向波數(shù)/頻率和包絡線方法計算擾動放大因子的公式分別為

(2)

(3)

在馬赫數(shù)1.3～2.0超聲速民機的機翼邊界層中，由于馬赫數(shù)的增加，二維TS波受到顯著抑制，應主要考慮三維TS斜波；而邊界層轉捩由CF駐波還是由CF行波占主導，則取決于來流湍流度、壁面粗糙度等條件[40]。例如，在來流湍流度大或機翼前緣粗糙時，轉捩一般由CF行波主導；在來流湍流度小、機翼表面光滑的情況下，轉捩一般由CF駐波主導。這一結論在超聲速流動下還需進一步確認。靜音風洞和飛行試驗表明，在低湍流度、高馬赫數(shù)環(huán)境下，飛行器表面也存在橫流行波主導轉捩現(xiàn)象[54-55]。而在超聲速層流機翼設計中，考慮到CF行波的擾動放大率一般高于CF駐波，按照CF行波進行設計更加保險[60]。

1.1.3 轉捩臨界N因子的選取

對于計算得到的不同頻率/波數(shù)擾動的放大因子，結合相應的轉捩判據，就可以預測得到機翼邊界層轉捩位置。eN方法是一種半經驗方法，通常需要借助風洞試驗和飛行試驗等手段對轉捩判據進行標定。

針對不同類型的TS和CF波擾動，給出相應的轉捩判據。根據與DLR-F4機翼風洞試驗的結果對比[38]，二維TS波和CF駐波的臨界擾動放大因子可分別選用Ntr_2DTS=10.5,Ntr_SCF=7.5。根據NASA的靜音風洞試驗結果[61]和JAXA的飛行試驗結果[3]標定得到，在采用eN-包絡線方法時，可選用Ntr_ENV=14[61]或Ntr_ENV=12.5[3]作為超聲速機翼轉捩預測的臨界N因子；根據NEXST-1超聲速民機構型的飛行試驗轉捩位置與穩(wěn)定性分析結果的對比可得[53]，在采用固定展向波數(shù)/頻率方法時，可選用Ntr_OTS=9.5和Ntr_TCF=9.8分別作為TS斜波和CF行波臨界N因子。應注意的是，超聲速邊界層轉捩受到來流擾動條件、來流雷諾數(shù)、壁面粗糙度等許多因素的影響，而ONERA和JAXA[60]等僅對來流雷諾數(shù)、來流擾動的影響進行了定性分析。為了保證eN方法轉捩判據的泛化能力，還需要針對各種因素的影響開展深入研究。

1.2 耦合流場求解的轉捩自動判斷流程

本文方法將eN方法和RANS求解器相結合，其總體思路是通過從Navier-Stokes方程的層流邊界層解中提取邊界層信息，找出各位置的不穩(wěn)定擾動，計算相應的擾動放大因子，并結合一定的轉捩判據預測轉捩位置，最后將預測的轉捩位置返回流場求解器，實現(xiàn)閉合循環(huán)。流程框架如圖3所示。下面具體介紹該方法中的各個步驟。

圖3 面向超聲速民機層流機翼設計的轉捩預測方法框架

步驟1給定預估的轉捩位置，用于激活湍流模型中的生成項，從而保證預估轉捩位置的上游為層流邊界層，以便開展穩(wěn)定性分析。

步驟2在給定轉捩位置下，進行流場求解，直至收斂。

步驟3提取機翼層流邊界層信息，用于穩(wěn)定性分析。一種提取方式是首先從流場解中得到機翼表面壓力分布、邊界層外邊界速度等流動參數(shù)，然后求解可壓縮邊界層方程，從而獲取光滑的邊界層速度、溫度剖面及其一階、二階導數(shù)。另一種方式是直接從流場解中提取層流邊界層信息。根據NLR的研究結果[33]，采用后一種方式的話，至少需要在邊界層內法向布置約60～70個網格點才能保證準確提取邊界層信息。

步驟4采用固定波角φw和固定頻率f方法，求解線性穩(wěn)定性方程，尋找不穩(wěn)定TS和CF波擾動，計算得到不同頻率/波數(shù)擾動的中性曲線，作為擾動幅值增長積分計算的起點。

步驟5從中性曲線上選取不穩(wěn)定擾動，采用固定展向波數(shù)/頻率方法或包絡線方法，求解線性穩(wěn)定性方程，得到不同有量綱波數(shù)/頻率的擾動放大率，積分得到相應的擾動放大因子N。

步驟6根據計算得到的不同波數(shù)/頻率擾動的放大因子N，并結合相應的轉捩判據，預測新的轉捩位置。

步驟7將預測得到的轉捩位置返回流場求解過程，并重復進行步驟2～步驟7，直至預測的轉捩位置收斂。

本文的轉捩自動判斷方法框架，和DLR[28]、NLR[33]等發(fā)表方法的主要區(qū)別在于，對其中eN方法的擾動放大因子積分策略進行了改進，考慮了超聲速機翼邊界層TS斜波和CF行波誘導轉捩的預測。

2 超聲速機翼邊界層轉捩預測驗證算例

為了驗證所發(fā)展的轉捩預測方法的正確性，本文選取了NASA的Owens等[62]研究的超聲速后掠機翼模型，在不同來流雷諾數(shù)下進行了流場求解和邊界層穩(wěn)定性分析，并將計算得到的CF波擾動N因子和文獻中的計算結果進行了比較。

2.1 NASA試驗機翼模型和計算狀態(tài)

所選用于驗證轉捩預測方法的機翼模型源自NASA的Swept-Wing Laminar Flow項目[62]，目的是研究超聲速機翼邊界層橫流不穩(wěn)定性及控制。NASA近年內依次開展了CFD分析、風洞試驗和飛行試驗研究。風洞試驗在NASA Langley研究中心的20 in(1 in=2.54 cm)超聲速風洞中進行，飛行試驗在NASA Armstrong飛行研究中心的F-15B載機平臺上進行。

該機翼的平面形狀和剖面翼型如圖4[62]所示，機翼前緣后掠角為65°，翼尖切角35°。根弦長croot=566.8 mm，展長b=266.7 mm，垂直前緣的弦線長度cn=239.5 mm,xn為垂直前緣距離。機翼剖面翼型為雙弧形翼型，沿垂直機翼前緣方向布置。本文所用翼型由文獻[62]結果擬合得到。

圖4 NASA試驗模型[62]的平面形狀和剖面翼型

NASA開展風洞試驗和飛行試驗的來流馬赫數(shù)為Ma=2.0，來流單位雷諾數(shù)為Re=4.6×106～13.2×106m-1，攻角為α=0°。風洞試驗中觀測到機翼表面的轉捩前沿呈光滑狀，因此認為是橫流行波誘導的轉捩；飛行試驗中觀測到的轉捩前沿則是鋸齒狀，因此認為是橫流駐波誘導的轉捩。

2.2 NASA試驗機翼流場求解

為了獲得用于穩(wěn)定性分析的基本流，在Re=6.2×106m-1和Re=11.5×106m-1兩個試驗狀態(tài)下，對試驗機翼的繞流流場和邊界層流動進行分析。流場求解采用課題組自主開發(fā)的結構多塊網格RANS方程求解程序——“PMNS3D”[63]。計算總網格量390萬，機翼表面沿流向布置400個網格單元、展向布置60個網格單元、法向布置100個網格單元。第1層網格高度5×10-6m，法向增長率1.1。上下表面初始轉捩固定在90%弦長位置，進行流場求解，然后根據從流場中提取的機翼表面壓力分布、邊界層外邊界速度等信息，求解三維可壓縮邊界層方程，獲取邊界層速度剖面、溫度剖面及其一階、二階導數(shù)。

以Re=11.5×106m-1狀態(tài)為例，給出試驗機翼表面壓力云圖、邊界層外邊界勢流流線如圖5 所示?？梢姡瑒萘髁骶€首先沿機翼前緣的方向，然后迅速恢復到接近自由來流方向，且和機翼表面壓力梯度方向存在一定的夾角。沿垂直機翼前緣的剖面(見圖5)，分別給出機翼剖面壓力分布、當?shù)睾舐咏欠植迹鐖D6所示，其中c為機翼弦長。此處，當?shù)睾舐咏铅誷w定義為邊界層外邊界勢流速度方向和當?shù)貕毫μ荻确较虻膴A角。從圖6 可以看出，試驗機翼表面存在較強的順壓梯度，而順壓梯度會增強橫流不穩(wěn)定性；機翼表面當?shù)睾舐咏菑臋C翼前緣的90°迅速減小到65°左右,接著緩慢減小，逐漸接近橫流最不穩(wěn)定的45°當?shù)睾舐咏?。給出試驗機翼不同流向位置的邊界層橫流速度形(如圖7所示，橫流速度w為邊界層內弦向和展向速度投影到垂直外邊界勢流方向的速度分量，ue為邊界層外邊界流向速度，x1、x2和x3為不同的流向位置)，可見機翼上表面橫流速度方向均為正(從翼根到翼梢)，且其最大值沿向下游的方向而增大；不同流向位置的邊界層流向速度形如圖8所示(圖中u為流向速度)，可以看出邊界層沿向下游的方向厚度呈增加的趨勢。下面將對NASA試驗機翼的邊界層擾動特性進行具體分析。

圖5 NASA試驗機翼表面壓力系數(shù)Cp分布和邊界層外邊界勢流流線(Ma=2.0,Re=11.5×106 m-1,α=0°)

圖6 NASA試驗機翼剖面壓力系數(shù)和當?shù)睾舐咏欠植?Ma=2.0, Re=11.5×106 m-1, α=0°)

圖7 NASA試驗機翼不同流向位置的邊界層橫流速度形(Ma=2.0, Re=11.5×106 m-1, α=0°, z/b=50%)

圖8 NASA試驗機翼不同流向位置的邊界層流向速度形(Ma=2.0, Re=11.5×106 m-1, α=0°, z/b=50%)

2.3 NASA試驗機翼邊界層穩(wěn)定性分析

作為對轉捩預測方法的驗證，分別在Re=6.2×106m-1和Re=11.5×106m-1兩個試驗狀態(tài)下對NASA試驗機翼進行邊界層穩(wěn)定性分析。首先，給出不同頻率下NASA機翼邊界層有量綱擾動放大率隨著波角的變化情況(見圖9)，可見擾動僅在波角為正的方向增長。結合機翼表面壓力分布、橫流速度型分布等信息分析可得，NASA試驗機翼邊界層轉捩由CF波占主導。然后，比較了不同頻率、不同雷諾數(shù)下的CF波擾動放大因子，并和文獻的計算結果進行了對比。

圖9 NASA試驗機翼邊界層內不同頻率和波角下的擾動放大率(Ma=2.0, Re=11.5×106 m-1, α=0°, z/b=50%, xn=23.5 mm)

2.3.1 不同頻率/波數(shù)下的CF波擾動放大因子

以Re=11.5×106m-1、風洞試驗狀態(tài)為例，試驗測得的轉捩位置在垂直前緣距離xn=23.5 mm 附近。在0～50 kHz的頻率范圍內，通過固定頻率方法找出中性曲線，并采用固定展向波數(shù)/頻率方法計算擾動放大率。圖10為不同頻率下的CF波擾動放大因子包絡線(圖中Nenvelope為各流向位置下，固定展向波數(shù)的擾動放大因子曲線族的最大值)，可見隨著頻率首先從0 Hz增加到33 kHz，各流向位置的CF波擾動放大因子均有所增長，最不穩(wěn)定橫流行波頻率約為33 kHz；而隨著頻率進一步增加到50 kHz，在xn=35 mm附近下游的擾動放大因子開始下降，而在xn=35 mm 附近上游的擾動放大因子略有增長。

圖10 NASA試驗機翼邊界層內不同頻率的CF波N因子包絡線計算結果(Ma=2.0, Re=11.5×106 m-1, α=0°, z/b=50%)

圖11分別為頻率0 Hz的CF駐波和頻率33 kHz 的CF行波的擾動放大因子，可見最不穩(wěn)定的CF行波放大因子高于CF駐波，且最不穩(wěn)定CF行波的展向波數(shù)約在2 424～2 727 m-1之間，低于最不穩(wěn)定CF駐波的3 540～4 215 m-1之間。

圖11 NASA試驗機翼邊界層內不同頻率/波數(shù)的CF波N因子計算結果(Ma=2.0, Re=11.5×106 m-1, α=0°, z/b=50%)

2.3.2 不同雷諾數(shù)下的CF波擾動放大因子

對Re=6.2×106m-1和Re=11.5×106m-1兩個狀態(tài)下的CF駐波與CF行波擾動放大因子進行比較，如圖12所示?？梢钥闯?，保持頻率為0 Hz，最不穩(wěn)定CF駐波的擾動放大因子NSCF的最大值均隨著雷諾數(shù)的增加而增加；給出不同頻率的CF波N因子包絡線，得到最不穩(wěn)定CF行波的放大因子包絡線Nenvelope的最大值也隨著雷諾數(shù)的增加而增加。

圖12 NASA試驗機翼邊界層內不同雷諾數(shù)的CF駐波和CF行波N因子計算結果(Ma=2.0, α=0°, z/b=50%)

2.3.3 計算結果和文獻結果的比較

將Re=6.2×106m-1和Re=11.5×106m-1兩個狀態(tài)下，計算得到的CF駐波與CF行波擾動放大因子，與NASA文獻[62]給出的計算結果進行比較，如圖13所示，圖中xtr_TCF和xtr_SCF分別為不同試驗中測量的橫流行波和橫流駐波誘導轉捩的位置，Ntr_CAL和Ntr_REF則是相應的計算與文獻給出的N值。可見，在不同來流雷諾數(shù)下，本文計算的CF行波的放大因子和文獻結果吻合良好，CF駐波的放大因子也基本吻合，驗證了本文發(fā)展的轉捩預測方法的正確性。其中一些微小差異可能是由采用的邊界層信息提取方法或擾動放大因子積分、穩(wěn)定性方程求解方法的差異所導致，對轉捩預測的影響不大。此外，還將飛行試驗狀態(tài)Ma=1.98、Re=10.9×106m-1、α=0°下預測的轉捩線和試驗測量的轉捩線進行了比較(見圖14)，飛行試驗中轉捩由橫流駐波誘導發(fā)生?？梢钥吹剑斶x取橫流駐波臨界轉捩放大因子Ntr_SCF=7.5時，預測的轉捩位置和試驗轉捩位置吻合良好，驗證了本文方法對于超聲速機翼轉捩預測的適用性。

圖13 NASA試驗機翼邊界層內CF駐波和CF行波N因子計算結果與文獻[62]結果的比較(Ma=2.0, α=0°, z/b=50%)

圖14 預測的NASA試驗機翼轉捩位置和飛行試驗測量轉捩位置[62]的比較(Ma=1.98, Re=10.9×106 m-1, α=0°)

3 超聲速層流機翼設計算例

為了進一步驗證本文發(fā)展的轉捩預測方法對于超聲速民機層流機翼設計的適用性，將其用于后掠角φsw=60°、來流馬赫數(shù)Ma=2.0、弦長雷諾數(shù)Rec=1.39×107的無限翼展層流機翼初步設計。首先，基于可壓縮三維邊界層相似解，探討了有利于抑制橫流不穩(wěn)定的理想壓力分布；然后，通過氣動反設計方法獲得了滿足理想壓力分布特性的機翼外形；最后，基于固定展向波數(shù)/頻率方法，給出了設計機翼的邊界層TS波和CF波擾動穩(wěn)定性分析結果與預測的轉捩位置。

3.1 理想壓力分布的探討

在高雷諾數(shù)、大后掠超聲速民機的層流機翼設計中，為了避免擾動在機翼前緣快速增長并誘發(fā)轉捩，關鍵在于抑制邊界層橫流不穩(wěn)定性。對此，通過對不同當?shù)睾舐咏铅誷w、壓力梯度因子βh和邊界層外邊界馬赫數(shù)Mae的三維可壓縮邊界層相似解(即Falkner-Skan-Cooke解)進行比較分析，探討了適合于抑制橫流不穩(wěn)定性的設計目標。由于橫流不穩(wěn)定性源自邊界層橫流速度剖面上的拐點，是一種無黏不穩(wěn)定性，因此在分析不同的邊界層剖面時，認為橫流速度最大值wmax越大，即橫流速度形越飽滿，則橫流不穩(wěn)定性越強。

針對上面提到的3種因素，分別開展了3組Falkner-Skan-Cooke方程解分析，其流動參數(shù)分別為：①Mae=2.0,βh=0.3,φsw=0～90°；②Mae=2.0,βh=-0.06～0.5,φsw=60°；③Mae=0.5～4.0,βh=0.3,φsw=60°。第1組解如圖15所示，其中η為無量綱法向坐標?？梢钥闯霎?shù)睾舐咏菑?°增加到45°時，橫流速度最大值wmax隨著后掠角的增大而增大；而隨著當?shù)睾舐咏堑倪M一步增大到90°，wmax反而逐漸減??；說明當?shù)睾舐咏窃浇咏?5°，橫流不穩(wěn)定性越強。第2組解如圖16所示，可見，當壓力梯度為0時，不存在橫流速度分量，橫流不穩(wěn)定性最弱；而不管壓力梯度因子>0還是<0，即邊界層為順壓或逆壓情況下，橫流不穩(wěn)定性都隨著壓力梯度的絕對值增大而增強。第3組解如圖17所示，可以看出，wmax隨馬赫數(shù)增加而增加，橫流不穩(wěn)定性也隨之增強。

圖15 不同當?shù)睾舐咏窍翭alkner-Skan-Cooke解的橫流速度形(Mae=2.0，βh=0.3，φsw=0°～90°)

圖16 不同壓力梯度因子下Falkner-Skan-Cooke解的橫流速度形(Mae=2.0，βh=-0.06～0.5，φsw=60°)

圖17 不同馬赫數(shù)下Falkner-Skan-Cooke解的橫流速度形(Mae=0.5～4.0，βh=0.3，φsw=60°)

層流機翼反設計中，通過調整機翼外形來獲得所需要的表面壓力分布。根據上述因素影響分析可知，壓力梯度絕對值越小，即機翼表面壓力分布越平坦，橫流不穩(wěn)定性越弱。由此，可以提出一種“讓流動在機翼前緣迅速加速(約0.5%弦長內)，繼而維持微弱壓力梯度”的理想壓力分布，以此來抑制橫流不穩(wěn)定性的增長。

3.2 兩輪次壓力分布反設計

針對來流馬赫數(shù)Ma=2.0的無限翼展后掠機翼，為了獲得在機翼上表面滿足理想壓力分布的氣動外形，基于課題組自主開發(fā)的代理優(yōu)化程序——“SurroOpt”[64]，開展了機翼氣動反設計。設計機翼弦長為1 m，后掠角60°。設計狀態(tài)為：Ma=2.0,Rec=1.39×107,α=2°。反設計機翼幾何外形和上表面目標壓力分布如圖18所示。根據收斂性分析結果，選取計算網格量為520×100×16，沿翼面流向布置400個網格單元，后緣割縫處布置60個網格單元，法向布置100個網格單元，展向布置16個網格單元，并在兩側采用循環(huán)邊界條件。為了實現(xiàn)圖18中的目標壓力分布，需要讓流動在機翼前緣很小一段范圍(約0.5%弦長)內迅速加速，即反設計對機翼前緣的幾何精度具有很高的要求。而直接采用類/形函數(shù)變換(CST)參數(shù)化方法對機翼前緣的擬合精度有限。因此，在NACA0003翼型的基準機翼基礎上，本文采用了兩輪次設計的方法：第1輪采用直接CST參數(shù)化方法，獲得基本滿足目標壓力分布的中間外形；第2輪在該外形基礎上，采用自由變形(FFD)參數(shù)化方法獲得較好滿足目標壓力分布的設計外形。

圖18 反設計機翼平面形狀和上表面目標壓力分布(Ma=2.0, Rec=1.39×107, α=2°)

3.2.1 基于CST參數(shù)化方法的第1輪反設計

在代理優(yōu)化框架[65-67]下，以當前壓力分布和目標壓力分布的差量為設計目標，以機翼上表面外形為設計對象，開展了基于直接CST方法的氣動反設計。僅對機翼上表面外形進行參數(shù)化，機翼下表面維持NACA0003翼型不變。參數(shù)化方法采用8階直接CST，并將機翼前緣的CST類函數(shù)也作為1個設計變量，得到一共10個設計變量。設計空間通過對基準的NACA0003翼型法向坐標上下擾動50%得到。優(yōu)化數(shù)學模型為

(4)

式中：n為流向位置標號；n0為流向位置總數(shù)；Cp,n和Cp_target,n分別為設計機翼表面壓力系數(shù)和目標壓力系數(shù)；si為設計變量；si_min和si_max分別為設計變量的上下界，即設計空間。

優(yōu)化收斂曲線如圖19所示。首先通過拉丁超立方抽樣方法，在設計空間內選取50個初始樣本點，建立代理模型；然后通過改善期望(EI)和最小代理模型預測值(MSP)兩種加點準則，在固定設計空間內，不斷尋找新的樣本點，更新代理模型；最后根據樣本分布情況，自適應調整設計空間，繼續(xù)加點和更新代理模型，直至收斂。得到的第1輪設計機翼表面壓力分布如圖20所示，可以看到機翼上表面基本滿足設計目標，但壓力分布仍存在微小抖動。

圖19 第1輪基于CST參數(shù)化方法的機翼反設計收斂曲線和樣本分布

圖20 第1輪設計機翼表面壓力分布和基準、目標壓力分布的比較(Ma=2.0, Rec=1.39×107, α=2°)

3.2.2 基于FFD參數(shù)化方法的第2輪反設計

為了抹平第1輪設計機翼表面壓力分布的微小擾動，在其基礎上又開展了基于FFD方法的氣動反設計。在翼型表面壓力分布波動的位置布置FFD控制點，通過改變FFD控制點坐標獲得新的機翼上表面翼型，一共20個設計變量。設計空間由FFD控制點坐標上下各擾動10%得到。機翼下表面同樣維持NACA0003翼型不變。

優(yōu)化收斂過程如圖21所示，首先通過拉丁超立方抽樣方法，在設計空間內選取50個初始樣本點，建立代理模型；然后通過EI+MSP兩種加點準則，在固定設計空間內，不斷尋找新的樣本點，并更新代理模型，直至收斂。得到的第2輪設計機翼表面壓力分布如圖22所示，可以看到設計機翼能夠良好滿足給定的上表面目標壓力分布。

圖21 第2輪基于FFD參數(shù)化方法的機翼反設計收斂曲線和樣本分布

圖22 第2輪設計機翼表面壓力分布和基準、目標壓力分布的比較(Ma=2.0, Rec=1.39×107, α=2°)

3.3 設計機翼邊界層穩(wěn)定性分析和轉捩預測

在設計狀態(tài)下，分別給出配置NACA0003翼型的基準機翼和第2輪設計機翼的上表面邊界層內不同波數(shù)/頻率擾動的空間演化情況以及預測的轉捩位置。轉捩判據采用從JAXA文獻[53]中得到的Ntr_OTS=9.5和Ntr_TCF=9.8。

對于基準機翼，從計算得到的不同頻率、波角的有量綱擾動放大率分布(圖23)可以看出，雖然擾動在波角為正和負兩個方向都存在增長，但波角為正方向的擾動放大率明顯高于波角為負的方向，且正波角方向的最不穩(wěn)定擾動波角大于70°，符合CF行波特征，而負波角方向的最不穩(wěn)定擾動波角在40°～50°，符合TS斜波特征，說明在波角為正方向最不穩(wěn)定CF行波不穩(wěn)定性要強于TS斜波模態(tài)；因此認為，在波角為正的方向計算的不同波數(shù)下最不穩(wěn)定擾動放大率圍成的包絡線(圖24)屬于CF行波。根據前面給出的橫流行波轉捩臨界因子Ntr_TCF=9.8，得到轉捩位置在弦長x/c=0.27，即機翼上表面層流范圍為27%。

圖23 基準機翼邊界層不同頻率和波角下的擾動放大率(Ma=2.0, Rec=1.39×107, α=2°, x/c=0.2)

圖24 基準機翼邊界層內不同頻率下的CF行波擾動N因子包絡線計算結果(Ma=2.0, Rec=1.39×107, α=2°)

對于設計機翼，計算得到不同頻率下設計機翼邊界層內的有量綱擾動放大率隨著波角的變化情況(見圖25)，可見各頻率下不穩(wěn)定擾動在正負波角方向上均存在增長，可得設計機翼邊界層轉捩由TS斜波主導。然后，采用固定展向波數(shù)/頻率方法，計算得到設計機翼邊界層內不同頻率的TS斜波和CF行波擾動N因子包絡線(見圖26)，可見CF波擾動在機翼前緣附近得到了很好的抑制，而TS斜波擾動N因子從5%弦長位置附近開始增長，在邊界層終點達到最大值，但仍未超過轉捩臨界放大因子Ntr_OTS=9.5。最后，給出計算得到的不同頻率下TS和CF波擾動放大因子如圖27所示,可見，波角為負的最不穩(wěn)定TS波擾動放大因子均大于波角為正的TS波；隨著頻率從9 kHz增加到12 kHz，在x/c=0～0.75范圍內邊界層最不穩(wěn)定擾動放大因子Nmax增大，在x/c>0.75范圍內邊界層最不穩(wěn)定擾動放大因子Nmax有所減小。

圖25 設計機翼邊界層不同頻率和波角下的擾動放大率(Ma=2.0, Rec=1.39×107, α=2°, x/c=0.2)

圖26 設計機翼邊界層內不同頻率下的擾動N因子包絡線計算結果(Ma=2.0, Rec=1.39×107, α=2°)

圖27 設計機翼邊界層內不同頻率/展向波數(shù)下TS和CF波擾動N因子計算結果(Ma=2.0, Rec=1.39×107, α=2°)

綜上可得，設計機翼相比基準機翼能夠更好地抑制CF行波不穩(wěn)定性，在設計狀態(tài)Ma=2.0、Rec=1.39×107、α=2°下機翼上表面可以接近維持全層流，具有良好的層流特性。

4 結 論

針對超聲速民機機翼層流減阻設計對轉捩自動判斷的需求，本文發(fā)展了能考慮TS斜波和CF行波的耦合RANS求解器的eN轉捩預測方法。主要結論如下：

1) 改進了eN方法的擾動放大因子積分策略，考慮了超聲速邊界層中可能誘發(fā)轉捩的TS斜波、CF駐波、CF行波等不穩(wěn)定性。采用固定波角和固定頻率方法來尋找TS和CF波擾動中性曲線，通過固定展向波數(shù)/頻率方法或包絡線方法積分計算擾動放大因子，最后結合相應的轉捩判據預測轉捩位置。

2) 為驗證所發(fā)展的轉捩預測方法的正確性，選取NASA的Owens等發(fā)布的超聲速后掠機翼模型，在不同來流雷諾數(shù)下進行了流場求解和邊界層穩(wěn)定性分析，結果表明計算的CF駐波和CF行波擾動放大因子N和文獻結果基本吻合，驗證了本文方法的正確性。

3) 將發(fā)展的轉捩預測方法用于來流馬赫數(shù)Ma=2.0、雷諾數(shù)Rec=1.39×107和后掠角φsw=60°的無限翼展層流機翼設計，提出了一種使機翼表面流動在前緣附近迅速加速，然后維持弱壓力梯度的理想壓力分布。經評估，設計機翼邊界層內CF波在機翼前緣附近得到了很好的抑制，且TS斜波和CF行波擾動放大因子均未達到轉捩臨界放大因子Ntr，表明設計機翼上表面可以接近維持全層流，驗證了該方法對于超聲速民機層流機翼設計的適用性。

致 謝

感謝西北工業(yè)大學許建華、宋科、喬建領、張力文、丁玉臨、許朕銘等的意見和幫助，感謝中國空氣動力研究與發(fā)展中心涂國華、李曉虎、楊強、陳曦、董思衛(wèi)、向星皓等的大力支持和幫助。